/[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 126 - (show annotations)
Fri Jul 22 03:53:08 2005 UTC (13 years, 8 months ago) by jgs
File MIME type: text/x-python
File size: 2233 byte(s)
Merge of development branch back to main trunk on 2005-07-22

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*1e25 # sqrt(km)/sec
9 g_c=0*1e-5 # Hz
10 H_level=1
11 #=====================================
12 L=1000. #km
13 mydomain=finley.Rectangle(l0=L,l1=L,n0=50,n1=50)
14 #mydomain=finley.ReadMesh("test.msh")
15 e_size=mydomain.getSize()
16 #====================================
17 # profile:
18 #===================================
19 x=mydomain.getX()
20 m=(x[0]-30.).whereNegative()
21 H=(H_level+0.0015)*m+H_level*(1.-m)*exp(-0.05*x[0]**2)
22 H.saveVTK("H.xml")
23
24 #====================================
25 h_pde=LinearPDE(mydomain)
26 h_pde.setValue(D=1.)
27 h_pde.setSymmetryOn()
28 h_pde.setLumpingOn()
29 hv_pde=LinearPDE(mydomain)
30 hv_pde.setValue(D=kronecker(mydomain))
31 hv_pde.setSymmetryOn()
32 hv_pde.setLumpingOn()
33 #====================================
34
35 def getDh(h_pde,h,v,hv,H):
36 h_pde.setValue(X=hv)
37 return h_pde.getSolution()
38
39 def getDhv(hv_pde,h,v,hv,H):
40 # F=trace(grad(outer(hv,v)+0.5*g*(h**2-H**2)*kronecker(hv_pde.getDomain())))
41 # Q=h*g_c*matrixmult(numarray.array([[0,-1],[1,0]]),v)-g*(h-H)*grad(H)+g*length(v)/(C*h**2)*v
42 F=outer(hv,v)+0.5*g*(h**2-H**2)*kronecker(hv_pde.getDomain())
43 Q=-g*(h-H)*grad(H)
44 hv_pde.setValue(X=F,Y=-Q)
45 return hv_pde.getSolution()
46
47
48 x=mydomain.getX()
49 l=length(x-[L/2.,L/2.])
50 m=(l-30.).whereNegative()
51 lift=H_level+0.015*(m+(1.-m)*exp(-0.0005*l**2))
52
53 h=(lift-H).whereNegative()*H+(lift-H).whereNonNegative()*lift
54 v=Vector(0.,ContinuousFunction(mydomain))
55 hv=h*v
56
57 t=0
58 dt=0.01
59 t_count=0
60 while t_count<20000:
61 t_count+=1
62 t+=dt
63 print "@ ",t_count,"-th time step t =",t," (step size = ",dt,")"
64 # Taylor Galerkin method
65 h_half=h+dt/2*getDh(h_pde,h,v,hv,H)
66 hv_half=hv+dt/2*getDhv(hv_pde,h,v,hv,H)
67 v_half=hv_half/(h_half+1.e-15)
68 h=h+dt*getDh(h_pde,h_half,v_half,hv_half,H)
69 hv=hv+dt*getDhv(hv_pde,h_half,v_half,hv_half,H)
70 v=hv/(h+1.e-15)
71
72 if (t_count-1)%10==0:
73 h.saveVTK("h.%d.xml"%((t_count-1)/10))
74 # v.saveVTK("v.%d.xml"%((t_count-1)/10))
75 print "save ",t_count
76
77 dt=min(inf(e_size/(sqrt(g*h)+length(v)+1.e-15)),1.)
78 print "@v =",sup(v),inf(v)
79 print "@h =",sup(h),inf(h)
80 print "@c =",Lsup(sqrt(g*h))
81 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