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

  ViewVC Help
Powered by ViewVC 1.1.26