/[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 126 - (hide annotations)
Fri Jul 22 03:53:08 2005 UTC (15 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 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 jgs 126 C=2*1e25 # sqrt(km)/sec
9     g_c=0*1e-5 # Hz
10 jgs 123 H_level=1
11     #=====================================
12     L=1000. #km
13 jgs 126 mydomain=finley.Rectangle(l0=L,l1=L,n0=50,n1=50)
14 jgs 123 #mydomain=finley.ReadMesh("test.msh")
15     e_size=mydomain.getSize()
16     #====================================
17     # profile:
18     #===================================
19     x=mydomain.getX()
20 jgs 126 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 jgs 123 #====================================
25     h_pde=LinearPDE(mydomain)
26     h_pde.setValue(D=1.)
27     h_pde.setSymmetryOn()
28 jgs 126 h_pde.setLumpingOn()
29 jgs 123 hv_pde=LinearPDE(mydomain)
30     hv_pde.setValue(D=kronecker(mydomain))
31     hv_pde.setSymmetryOn()
32 jgs 126 hv_pde.setLumpingOn()
33 jgs 123 #====================================
34    
35     def getDh(h_pde,h,v,hv,H):
36 jgs 126 h_pde.setValue(X=hv)
37 jgs 123 return h_pde.getSolution()
38    
39     def getDhv(hv_pde,h,v,hv,H):
40 jgs 126 # 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 jgs 123 return hv_pde.getSolution()
46    
47    
48     x=mydomain.getX()
49 jgs 126 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 jgs 123 v=Vector(0.,ContinuousFunction(mydomain))
55     hv=h*v
56    
57     t=0
58     dt=0.01
59     t_count=0
60 jgs 126 while t_count<20000:
61 jgs 123 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 jgs 126 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 jgs 123 print "save ",t_count
76    
77 jgs 126 dt=min(inf(e_size/(sqrt(g*h)+length(v)+1.e-15)),1.)
78 jgs 123 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