/[escript]/trunk/escript/py_src/mountains.py
ViewVC logotype

Contents of /trunk/escript/py_src/mountains.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2344 - (show annotations)
Mon Mar 30 02:13:58 2009 UTC (10 years, 5 months ago) by jfenwick
File MIME type: text/x-python
File size: 3635 byte(s)
Change __url__ to launchpad site

1
2 ########################################################
3 #
4 # Copyright (c) 2003-2009 by University of Queensland
5 # Earth Systems Science Computational Center (ESSCC)
6 # http://www.uq.edu.au/esscc
7 #
8 # Primary Business: Queensland, Australia
9 # Licensed under the Open Software License version 3.0
10 # http://www.opensource.org/licenses/osl-3.0.php
11 #
12 ########################################################
13
14 __copyright__="""Copyright (c) 2003-2009 by University of Queensland
15 Earth Systems Science Computational Center (ESSCC)
16 http://www.uq.edu.au/esscc
17 Primary Business: Queensland, Australia"""
18 __license__="""Licensed under the Open Software License version 3.0
19 http://www.opensource.org/licenses/osl-3.0.php"""
20 __url__="https://launchpad.net/escript-finley"
21
22 from esys.escript import *
23 from esys.escript.linearPDEs import LinearPDE
24 from esys import finley
25 import sys
26
27 class Mountains:
28 """
29 The Mountains class is defined by the following equations:
30
31 (1) w_i,aa+(1/eps)*w_i,33=0 where 0<eps<<1 and a=1,2 and w_i is the extension of the surface velocity where w_i(x_3=1)=v_i.
32
33 (2) Integration of topography PDE using Taylor-Galerkin upwinding to stabilize the advection terms
34 H^(t+dt)=H^t+dt*w_3+w_hat*dt*[(div(w_hat*H^t)+w_3)+(dt/2)+H^t],
35 where w_hat=w*[1,1,0], dt<0.5*d/max(w_i), d is a characteristic element size; H(x_3=1)=lambda (?)
36
37
38 """
39 def __init__(self,domain,v,eps=0.01,z=1):
40 """
41 Sets up the level set method.
42
43 @param domain: the domain where the mountains is used
44 @param eps: the parameter for (1)
45 @param z: the height of the box
46 """
47 self.__domain = domain
48 self.__eps = eps
49 self.__PDE_W = LinearPDE(domain)
50 self.__PDE_W.setSymmetryOn()
51 self.__PDE_H = LinearPDE(domain)
52 self.__PDE_H.setSymmetryOn()
53 self.__PDE_H.setValue(D=1.0)
54 self.__PDE_H.setSolverMethod(solver=LinearPDE.LUMPING)
55 self.__h = inf(domain.getSize())
56 self.__v=v
57 self.__z=z
58 self.__x=domain.getX()
59 self.__DIM=domain.getDim()
60 self.__fixed_w_mask=whereZero(self.__z-self.__x[self.__DIM-1])+whereZero(self.__x[self.__DIM-1])
61 self.__fixed_H_mask=whereZero(self.__x[self.__DIM-1])
62 self.__A=kronecker(domain)
63 self.__A[self.__DIM-1,self.__DIM-1]=1/self.__eps
64 self.__PDE_W.setValue(A=self.__A, q=self.__fixed_w_mask)
65 self.__PDE_H.setValue(q=self.__fixed_H_mask)
66 self.__H_t=Scalar(0.0, Solution(domain))
67 self.__dt=0.
68
69 def update(self,u=None,H_t=None,dt=None,verbose=False):
70 """
71 Sets a new W and updates the H function.
72
73 @param dt: time step forward
74 """
75 if H_t==None:
76 H_t=self.__H_t
77
78 if u==None:
79 u=self.__v
80
81 for d in range(self.__DIM):
82 self.__PDE_W.setValue(r=u[d]*self.__x[self.__DIM-1]/self.__z)
83 u[d]=self.__PDE_W.getSolution(verbose=verbose)
84
85 if dt==None:
86 dt=0.5*self.__h/sup(u)
87 self.__dt=dt
88 else:
89 self.__dt=dt
90
91 w=u
92 w_tilda=1.*w
93 w_tilda[self.__DIM-1]=0
94
95 self.__PDE_H.setValue(X=((div(w_tilda*H_t)+w[self.__DIM-1])*dt/2+H_t)*w_tilda*dt, Y=w[self.__DIM-1]*dt+H_t)
96 self.__H_t=self.__PDE_H.getSolution()
97
98 self.__v=w
99
100 return w,self.__H_t
101
102 def getH(self):
103 """
104 Returns the mesh size.
105 """
106 return self.__h
107
108 def getDt(self):
109 """
110 Returns the time step value.
111 """
112 return self.__dt
113
114
115 def getDomain(self):
116 """
117 Returns the domain.
118 """
119 return self.__domain
120
121 def setTolerance(self,tolerance=1e-3):
122 self.__PDE_W.setTolerance(tolerance)
123 self.__PDE_H.setTolerance(tolerance)
124
125

  ViewVC Help
Powered by ViewVC 1.1.26