/[escript]/trunk/doc/examples/usersguide/mount.py
ViewVC logotype

Diff of /trunk/doc/examples/usersguide/mount.py

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

revision 2562 by gross, Tue Jun 30 05:49:22 2009 UTC revision 2563 by gross, Tue Jul 28 03:50:45 2009 UTC
# Line 1  Line 1 
1    ########################################################
2    #
3    # Copyright (c) 2003-2009 by University of Queensland
4    # Earth Systems Science Computational Center (ESSCC)
5    # http://www.uq.edu.au/esscc
6    #
7    # Primary Business: Queensland, Australia
8    # Licensed under the Open Software License version 3.0
9    # http://www.opensource.org/licenses/osl-3.0.php
10    #
11    ########################################################
12    
13    __copyright__="""Copyright (c) 2003-2009 by University of Queensland
14    Earth Systems Science Computational Center (ESSCC)
15    http://www.uq.edu.au/esscc
16    Primary Business: Queensland, Australia"""
17    __license__="""Licensed under the Open Software License version 3.0
18    http://www.opensource.org/licenses/osl-3.0.php"""
19    __url__="https://launchpad.net/escript-finley"
20  from esys.escript import *  from esys.escript import *
21  from esys.escript.models import Mountains  from esys.escript.models import Mountains
22  from esys.finley import Brick,Rectangle  from esys.finley import Brick,Rectangle
# Line 7  NE=16 Line 26  NE=16
26  DIM=3  DIM=3
27  H=1.  H=1.
28  L=2*H  L=2*H
 TOL=1.e-4  
29  OMEGA=10  OMEGA=10
30  EPS=0.01  EPS=0.01
31  t=0  t=0
# Line 17  if DIM==2: Line 35  if DIM==2:
35    mydomain=Rectangle(int(ceil(L*NE/H)),NE,l0=L,l1=H,order=1, useFullElementOrder=True,optimize=True)    mydomain=Rectangle(int(ceil(L*NE/H)),NE,l0=L,l1=H,order=1, useFullElementOrder=True,optimize=True)
36  else:  else:
37    mydomain=Brick(int(ceil(L*NE/H)),int(ceil(L*NE/H)),NE,l0=L,l1=L,l2=H,order=1, useFullElementOrder=True,optimize=True)    mydomain=Brick(int(ceil(L*NE/H)),int(ceil(L*NE/H)),NE,l0=L,l1=L,l2=H,order=1, useFullElementOrder=True,optimize=True)
38    
39  x=mydomain.getX()  x=mydomain.getX()
40  v = Vector(0.0, Solution(mydomain))  v = Vector(0.0, Solution(mydomain))
41  if DIM==2:  if DIM==2:
# Line 38  else: Line 57  else:
57      v[2]=a2*cos(pi*n0*x[0])* cos(pi*n1*x[1])* sin(pi*n2*x[2])      v[2]=a2*cos(pi*n0*x[0])* cos(pi*n1*x[1])* sin(pi*n2*x[2])
58    
59    
60  H_t=Scalar(0.0, Solution(mydomain))  mts=Mountains(mydomain,eps=EPS)
 mts=Mountains(mydomain,v,eps=EPS,z=H)  
 dt=0.  
61  while t<T_END:  while t<T_END:
62      print "STEP ", t      print "STEP ", t
63      u=v*cos(OMEGA*t)      mts.setVelocity(v*cos(OMEGA*t))
64      u,Z=mts.update(u=u,H_t=H_t)      Z=mts.update()
65            
66      saveVTK("state.%d.vtu"%n,sol=Z)      saveVTK("state.%d.vtu"%n,sol=Z, v=mts.getVelocity())
67      print "Integral(Z)=",integrate(Z),Lsup(u[DIM-1])      print "Integral(Z)=",integrate(Z),Lsup(mts.getVelocity()[DIM-1])
68      n+=1      n+=1
69      H_t=Z      t+=mts.getSafeTimeStepSize()
     dt=mts.getDt()  
     t+=dt  
70    
Legend:
Removed from v.2562  
changed lines
  Added in v.2563
  ViewVC Help
Powered by ViewVC 1.1.26