/[escript]/trunk/esys2/finley/test/python/Tsunami.py
ViewVC logotype

Annotation of /trunk/esys2/finley/test/python/Tsunami.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 123 - (hide annotations)
Fri Jul 8 04:08:13 2005 UTC (15 years, 9 months ago) by jgs
File MIME type: text/x-python
File size: 1820 byte(s)
Merge of development branch back to main trunk on 2005-07-08

1 jgs 123 # $Id$
2    
3     from esys.escript import *
4     from esys.linearPDEs import *
5     from esys import finley
6    
7     g=9.81 * 1e-3 # km/sec
8     C=2 # sqrt(km)/sec
9     g_c=1e-5 # Hz
10     H_level=1
11     #=====================================
12     L=1000. #km
13     mydomain=finley.Rectangle(l0=L,l1=L,n0=100,n1=100)
14     #mydomain=finley.ReadMesh("test.msh")
15     e_size=mydomain.getSize()
16     #====================================
17     # profile:
18     #===================================
19     x=mydomain.getX()
20     H=H_level # *(1.-exp(-2.*x[0]/L))
21     #====================================
22     h_pde=LinearPDE(mydomain)
23     h_pde.setValue(D=1.)
24     h_pde.setSymmetryOn()
25     hv_pde=LinearPDE(mydomain)
26     hv_pde.setValue(D=kronecker(mydomain))
27     hv_pde.setSymmetryOn()
28     #====================================
29    
30     def getDh(h_pde,h,v,hv,H):
31     h_pde.setValue(Y=-trace(grad(hv)))
32     return h_pde.getSolution()
33    
34     def getDhv(hv_pde,h,v,hv,H):
35     F=trace(grad(outer(hv,v)+0.5*g*(h**2-H**2)*kronecker(hv_pde.getDomain())))
36     Q=h*g_c*matrixmult(numarray.array([[0,-1],[1,0]]),v)-g*(h-H)*grad(H)+g*length(v)/(C*h**2)*v
37     hv_pde.setValue(Y=-F-Q)
38     return hv_pde.getSolution()
39    
40    
41     x=mydomain.getX()
42     h=H-0.015*exp(-0.1*length(x-[L/2.,L/2.]))
43     v=Vector(0.,ContinuousFunction(mydomain))
44     hv=h*v
45    
46     t=0
47     dt=0.01
48     t_count=0
49     while t_count<200:
50     t_count+=1
51     t+=dt
52     print "@ ",t_count,"-th time step t =",t," (step size = ",dt,")"
53     # Taylor Galerkin method
54     h_half=h+dt/2*getDh(h_pde,h,v,hv,H)
55     hv_half=hv+dt/2*getDhv(hv_pde,h,v,hv,H)
56     v_half=hv_half/(h_half+1.e-15)
57     h=h+dt*getDh(h_pde,h_half,v_half,hv_half,H)
58     hv=hv+dt*getDhv(hv_pde,h_half,v_half,hv_half,H)
59     v=hv/(h+1.e-15)
60    
61     if t_count%10==0:
62     print "save ",t_count
63    
64     dt=min(inf(e_size/(sqrt(g*h)+length(v)+1.e-15)),10.)
65     print "@v =",sup(v),inf(v)
66     print "@h =",sup(h),inf(h)
67     print "@c =",Lsup(sqrt(g*h))
68     print "@dt =",inf(e_size/(sqrt(g*h)+length(v)+1.e-15))

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26