/[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 2620 - (show annotations)
Thu Aug 20 06:24:00 2009 UTC (9 years, 8 months ago) by gross
File MIME type: text/x-python
File size: 5656 byte(s)
some small additions to pycad to make life a bit easier.
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, reduced=False):
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 @param reduced: if True topography reduced order is used for reconstruction of internal velocity and topography
53 """
54 if eps<0:
55 raise ValueError("Smooting parameter eps must be non-negative.")
56 self.__domain = domain
57 self.__reduced=reduced
58 self.__DIM=domain.getDim()
59 z=domain.getX()[self.__DIM-1]
60
61 self.__PDE_W = LinearPDE(domain)
62 self.__PDE_W.setSymmetryOn()
63 A=kronecker(domain)*eps*0
64 A[self.__DIM-1,self.__DIM-1]=(sup(FunctionOnBoundary(self.__domain).getSize())/log(2.))**2
65 self.__PDE_W.setValue(D=1, A=A, q=whereZero(sup(z)-z))
66
67 self.__PDE_H = LinearPDE(domain)
68 self.__PDE_H.setSymmetryOn()
69 if reduced: self.__PDE_H.setReducedOrderOn()
70 A=kronecker(domain)*0
71 A[self.__DIM-1,self.__DIM-1]=0.1
72 self.__PDE_H.setValue(D=1.0)
73 # self.__PDE_H.getSolverOptions().setSolverMethod(SolverOptions.LUMPING)
74
75 self.setVelocity()
76 self.setTopography()
77 def getSolverOptionsForSmooting(self):
78 """
79 returns the solver options for the smoothing/extrapolation
80 """
81 return self.__PDE_W.getSolverOptions()
82
83 def getSolverOptionsForUpdate(self):
84 """
85 returns the solver options for the topograthy update
86 """
87 return self.__PDE_H.getSolverOptions()
88 def getDomain(self):
89 """
90 Returns the domain.
91 """
92 return self.__domain
93
94 def setVelocity(self,v=None):
95 """
96 set a new velocity. v is define on the entire domain but only the surface values are used.
97
98 @param v: velocity field. If None zero is used.
99 @type v: vector
100 """
101 self.__dt=None
102 self.__v=Vector(0.,Solution(self.getDomain()))
103 if not v == None:
104 for d in range(self.__DIM):
105 self.__PDE_W.setValue(r=v[d])
106 self.__v[d]=self.__PDE_W.getSolution()
107 def getVelocity(self):
108 """
109 returns the smoothed/extrapolated velocity
110 @rtype: vector L{Data}
111 """
112 return self.__v
113
114 def setTopography(self,H=None):
115 """
116 set the topography to H where H defines the vertical displacement. H is defined for the entire domain.
117
118 @param H: the topography. If None zero is used.
119 @type H: scalar
120 """
121 if self.__reduced:
122 fs=ReducedSolution(self.getDomain())
123 else:
124 fs=Solution(self.getDomain())
125
126 if H==None:
127 self.__H=Scalar(0.0, fs)
128 else:
129 self.__H=interpolate(H, fs)
130
131 def getTopography(self):
132 """
133 returns the current topography.
134 @rtype: scalar L{Data}
135 """
136 return self.__H
137
138 def getSafeTimeStepSize(self):
139 """
140 Returns the time step value.
141
142 @rtype: C{float}
143 """
144 if self.__dt == None:
145 h=self.getDomain().getSize()
146 self.__dt=0.5*inf(h/length(interpolate(self.getVelocity(),h.getFunctionSpace())))
147 return self.__dt
148 def update(self,dt=None, allow_substeps=True):
149 """
150 Sets a new W and updates the H function.
151
152 @param dt: time step forward. If None the save time step size is used.
153 @type dt: positve C{float} which is less or equal than the safe time step size.
154
155 """
156 if dt == None:
157 dt = self.getSafeTimeStepSize()
158 if dt<=0:
159 raise ValueError("Time step size must be positive.")
160 dt_safe=self.getSafeTimeStepSize()
161 n=max(int(math.ceil(dt/dt_safe)+0.5),1)
162 if n>1 and not allow_substeps:
163 raise SubSteppingException,"Substepping required."
164 dt/=n
165
166 H=self.getTopography()
167 w=self.getVelocity()
168 w_tilda=1.*w
169 w_tilda[self.__DIM-1]=0
170 w_z=w[self.__DIM-1]
171
172 t=0
173 for i in range(n):
174 L=integrate(w_z*dt+H)/vol(self.__PDE_H.getDomain())
175 self.__PDE_H.setValue(X=(inner(w_tilda,grad(H))*dt/2+H)*w_tilda*dt, Y=w_z*dt+H-L)
176 H=self.__PDE_H.getSolution()
177 print "DDD : ava = ",integrate(H)
178 t+=dt
179 self.setTopography(H)
180
181 return self.getTopography()

  ViewVC Help
Powered by ViewVC 1.1.26