/[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 2862 - (show annotations)
Thu Jan 21 04:45:39 2010 UTC (9 years, 7 months ago) by gross
File MIME type: text/x-python
File size: 6281 byte(s)
more review in the cookbook
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, SolverOptions
24 from esys import finley
25 import math
26 import sys
27
28 class SubSteppingException(Exception):
29 """
30 Thrown if the L{Mountains} class uses substepping.
31 """
32 pass
33
34
35 class Mountains:
36 """
37 The Mountains class is defined by the following equations:
38
39 (1) eps*w_i,aa+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.
40
41 (2) Integration of topography PDE using Taylor-Galerkin upwinding to stabilize the advection terms
42 H^(t+dt)=H^t+dt*w_3+w_hat*dt*[(div(w_hat*H^t)+w_3)+(dt/2)+H^t],
43 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 (?)
44
45 """
46 def __init__(self,domain,eps=0.01):
47 """
48 Sets up the level set method.
49
50 :param domain: the domain where the mountains is used
51 :param eps: the smoothing parameter for (1)
52 """
53 order=Solution(domain).getApproximationOrder()
54 if order>1:
55 reduced = True
56 if ReducedSolution(domain).getApproximationOrder()>1: raise ValueError,"Reduced order needs to be equal to 1."
57 else:
58 reduced = False
59 if eps<0:
60 raise ValueError("Smooting parameter eps must be non-negative.")
61 self.__domain = domain
62 self.__reduced=reduced
63 self.__DIM=domain.getDim()
64 z=domain.getX()[self.__DIM-1]
65
66 self.__PDE_W = LinearPDE(domain)
67 self.__PDE_W.setSymmetryOn()
68 A=kronecker(domain)*eps*0
69 A[self.__DIM-1,self.__DIM-1]=(0.3*(sup(z)-inf(z))/log(2.))**2
70 # A[self.__DIM-1,self.__DIM-1]=(sup(FunctionOnBoundary(self.__domain).getSize())/log(2.))**2
71 self.__PDE_W.setValue(D=1, A=A, q=whereZero(sup(z)-z)+whereZero(inf(z)-z))
72
73 self.__PDE_H = LinearPDE(domain)
74 self.__PDE_H.setSymmetryOn()
75 if reduced: self.__PDE_H.setReducedOrderOn()
76 # A=kronecker(domain)*0
77 # A[self.__DIM-1,self.__DIM-1]=0.1
78 self.__PDE_H.setValue(D=1.0, q=whereZero(inf(z)-z))
79 # self.__PDE_H.getSolverOptions().setSolverMethod(SolverOptions.LUMPING)
80
81 self.setVelocity()
82 self.setTopography()
83 def getSolverOptionsForSmooting(self):
84 """
85 returns the solver options for the smoothing/extrapolation
86 """
87 return self.__PDE_W.getSolverOptions()
88
89 def getSolverOptionsForUpdate(self):
90 """
91 returns the solver options for the topograthy update
92 """
93 return self.__PDE_H.getSolverOptions()
94 def getDomain(self):
95 """
96 Returns the domain.
97 """
98 return self.__domain
99
100 def setVelocity(self,v=None):
101 """
102 set a new velocity. v is define on the entire domain but only the surface values are used.
103
104 :param v: velocity field. If None zero is used.
105 :type v: vector
106 """
107 self.__dt=None
108 self.__v=Vector(0.,Solution(self.getDomain()))
109 if not v == None:
110 xi=self.getDomain().getX()[self.getDomain().getDim()-1]
111 v=(xi-inf(xi))/(sup(xi)-inf(xi))*v
112 for d in range(self.__DIM):
113 self.__PDE_W.setValue(r=v[d])
114 self.__v[d]=self.__PDE_W.getSolution()
115 def getVelocity(self):
116 """
117 returns the smoothed/extrapolated velocity
118 :rtype: vector `Data`
119 """
120 return self.__v
121
122 def setTopography(self,H=None):
123 """
124 set the topography to H where H defines the vertical displacement. H is defined for the entire domain.
125
126 :param H: the topography. If None zero is used.
127 :type H: scalar
128 """
129 if self.__reduced:
130 fs=ReducedSolution(self.getDomain())
131 else:
132 fs=Solution(self.getDomain())
133
134 if H==None:
135 self.__H=Scalar(0.0, fs)
136 else:
137 self.__H=interpolate(H, fs)
138
139 def getTopography(self):
140 """
141 returns the current topography.
142 :rtype: scalar `Data`
143 """
144 return self.__H
145
146 def getSafeTimeStepSize(self):
147 """
148 Returns the time step value.
149
150 :rtype: ``float``
151 """
152 if self.__dt == None:
153 h=self.getDomain().getSize()
154 self.__dt=0.5*inf(h/length(interpolate(self.getVelocity(),h.getFunctionSpace())))
155 return self.__dt
156 def update(self,dt=None, allow_substeps=True):
157 """
158 Sets a new W and updates the H function.
159
160 :param dt: time step forward. If None the save time step size is used.
161 :type dt: positve ``float`` which is less or equal than the safe time step size.
162
163 """
164 if dt == None:
165 dt = self.getSafeTimeStepSize()
166 if dt<=0:
167 raise ValueError("Time step size must be positive.")
168 dt_safe=self.getSafeTimeStepSize()
169 n=max(int(math.ceil(dt/dt_safe)+0.5),1)
170 if n>1 and not allow_substeps:
171 raise SubSteppingException,"Substepping required."
172 dt/=n
173
174 H=self.getTopography()
175 w=self.getVelocity()
176 w_tilda=1.*w
177 w_tilda[self.__DIM-1]=0
178 w_z=w[self.__DIM-1]
179 V=vol(self.__PDE_H.getDomain())
180
181 t=0
182 for i in range(n):
183 # L=integrate(w_z*dt+H)/vol(self.__PDE_H.getDomain())
184 # self.__PDE_H.setValue(X=(inner(w_tilda,grad(H))*dt/2+H)*w_tilda*dt, Y=w_z*dt+H-L)
185 # H=self.__PDE_H.getSolution()
186 L=integrate(w_z*dt+H)/V
187 self.__PDE_H.setValue(X=w_tilda*H*(dt/2), Y=w_z*(dt/2)+H-L)
188 Hhalf=self.__PDE_H.getSolution()
189 self.__PDE_H.setValue(X=w_tilda*Hhalf*dt, Y=w_z*dt+H-L)
190 H=self.__PDE_H.getSolution()
191 print "DDD : ava = ",integrate(H)
192 t+=dt
193 self.setTopography(H)
194
195 return self.getTopography()

  ViewVC Help
Powered by ViewVC 1.1.26