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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 123 - (show annotations)
Fri Jul 8 04:08:13 2005 UTC (15 years 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 # $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