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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 123 by jgs, Fri Jul 8 04:08:13 2005 UTC revision 126 by jgs, Fri Jul 22 03:53:08 2005 UTC
# Line 5  from esys.linearPDEs import * Line 5  from esys.linearPDEs import *
5  from esys import finley  from esys import finley
6    
7  g=9.81 * 1e-3 # km/sec  g=9.81 * 1e-3 # km/sec
8  C=2           # sqrt(km)/sec  C=2*1e25           # sqrt(km)/sec
9  g_c=1e-5      # Hz  g_c=0*1e-5      # Hz
10  H_level=1  H_level=1
11  #=====================================  #=====================================
12  L=1000. #km  L=1000. #km
13  mydomain=finley.Rectangle(l0=L,l1=L,n0=100,n1=100)  mydomain=finley.Rectangle(l0=L,l1=L,n0=50,n1=50)
14  #mydomain=finley.ReadMesh("test.msh")  #mydomain=finley.ReadMesh("test.msh")
15  e_size=mydomain.getSize()  e_size=mydomain.getSize()
16  #====================================  #====================================
17  #  profile:  #  profile:
18  #===================================  #===================================
19  x=mydomain.getX()  x=mydomain.getX()
20  H=H_level # *(1.-exp(-2.*x[0]/L))  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)  h_pde=LinearPDE(mydomain)
26  h_pde.setValue(D=1.)  h_pde.setValue(D=1.)
27  h_pde.setSymmetryOn()  h_pde.setSymmetryOn()
28    h_pde.setLumpingOn()
29  hv_pde=LinearPDE(mydomain)  hv_pde=LinearPDE(mydomain)
30  hv_pde.setValue(D=kronecker(mydomain))  hv_pde.setValue(D=kronecker(mydomain))
31  hv_pde.setSymmetryOn()  hv_pde.setSymmetryOn()
32    hv_pde.setLumpingOn()
33  #====================================  #====================================
34    
35  def getDh(h_pde,h,v,hv,H):  def getDh(h_pde,h,v,hv,H):
36     h_pde.setValue(Y=-trace(grad(hv)))     h_pde.setValue(X=hv)
37     return h_pde.getSolution()     return h_pde.getSolution()
38    
39  def getDhv(hv_pde,h,v,hv,H):  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())))     # 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     # 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     hv_pde.setValue(Y=-F-Q)     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()     return hv_pde.getSolution()
46    
47    
48  x=mydomain.getX()  x=mydomain.getX()
49  h=H-0.015*exp(-0.1*length(x-[L/2.,L/2.]))  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))  v=Vector(0.,ContinuousFunction(mydomain))
55  hv=h*v  hv=h*v
56    
57  t=0  t=0
58  dt=0.01  dt=0.01
59  t_count=0  t_count=0
60  while t_count<200:  while t_count<20000:
61    t_count+=1    t_count+=1
62    t+=dt    t+=dt
63    print "@ ",t_count,"-th time step t =",t," (step size = ",dt,")"      print "@ ",t_count,"-th time step t =",t," (step size = ",dt,")"  
# Line 58  while t_count<200: Line 69  while t_count<200:
69    hv=hv+dt*getDhv(hv_pde,h_half,v_half,hv_half,H)    hv=hv+dt*getDhv(hv_pde,h_half,v_half,hv_half,H)
70    v=hv/(h+1.e-15)    v=hv/(h+1.e-15)
71    
72    if t_count%10==0:    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        print "save ",t_count
76        
77    dt=min(inf(e_size/(sqrt(g*h)+length(v)+1.e-15)),10.)    dt=min(inf(e_size/(sqrt(g*h)+length(v)+1.e-15)),1.)
78    print "@v =",sup(v),inf(v)    print "@v =",sup(v),inf(v)
79    print "@h =",sup(h),inf(h)    print "@h =",sup(h),inf(h)
80    print "@c =",Lsup(sqrt(g*h))    print "@c =",Lsup(sqrt(g*h))

Legend:
Removed from v.123  
changed lines
  Added in v.126

  ViewVC Help
Powered by ViewVC 1.1.26