/[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 2564 - (show annotations)
Tue Jul 28 04:46:38 2009 UTC (10 years ago) by gross
File MIME type: text/x-python
File size: 5487 byte(s)
reduced PDE solutions added to topography
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 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 def __init__(self,domain,eps=0.01, reduced=False):
39 """
40 Sets up the level set method.
41
42 @param domain: the domain where the mountains is used
43 @param eps: the smoothing parameter for (1)
44 @param reduced: if True topography reduced order is used for reconstruction of internal velocity and topography
45 """
46 if eps<=0:
47 raise ValueError("Smmoting parameter eps must be positive.")
48 self.__domain = domain
49 self.__eps = eps
50 self.__reduced=reduced
51 self.__DIM=domain.getDim()
52 z=domain.getX()[self.__DIM-1]
53
54 self.__PDE_W = LinearPDE(domain)
55 self.__PDE_W.setSymmetryOn()
56 if reduced: self.__PDE_W.setReducedOrderOn()
57 A=kronecker(domain)
58 A[self.__DIM-1,self.__DIM-1]=1/self.__eps
59 self.__PDE_W.setValue(A=A, q=whereZero(sup(z)-z)+whereZero(inf(z)-z))
60
61 self.__PDE_H = LinearPDE(domain)
62 self.__PDE_H.setSymmetryOn()
63 if reduced: self.__PDE_H.setReducedOrderOn()
64 self.__PDE_H.setValue(D=1.0)
65 self.__PDE_H.getSolverOptions().setSolverMethod(SolverOptions.LUMPING)
66 self.__PDE_H.setValue(q=whereZero(inf(z)-z))
67
68 self.setVelocity()
69 self.setTopography()
70 def getSolverOptionsForSmooting(self):
71 """
72 returns the solver options for the smoothing/extrapolation
73 """
74 return self.__PDE_W.getSolverOptions()
75
76 def getSolverOptionsForUpdate(self):
77 """
78 returns the solver options for the topograthy update
79 """
80 return self.__PDE_H.getSolverOptions()
81 def getDomain(self):
82 """
83 Returns the domain.
84 """
85 return self.__domain
86
87 def setVelocity(self,v=None):
88 """
89 set a new velocity. v is define on the entire domain but only the surface values are used.
90
91 @param v: velocity field. If None zero is used.
92 @type v: vector
93 """
94 self.__dt=None
95 if self.__reduced:
96 fs=ReducedSolution(self.getDomain())
97 else:
98 fs=Solution(self.getDomain())
99 self.__v=Vector(0.,fs)
100 if not v == None:
101 z=self.getDomain().getX()[self.__DIM-1]
102 z_min=inf(z)
103 z_max=sup(z)
104 f=(z-z_min)/(z_max-z_min)
105 for d in range(self.__DIM):
106 self.__PDE_W.setValue(r=v[d]*f)
107 self.__v[d]=self.__PDE_W.getSolution()
108 def getVelocity(self):
109 """
110 returns the smoothed/extrapolated velocity
111 @rtype: vector L{Data}
112 """
113 return self.__v
114
115 def setTopography(self,H=None):
116 """
117 set the topography to H where H defines the vertical displacement. H is defined for the entire domain.
118
119 @param H: the topography. If None zero is used.
120 @type H: scalar
121 """
122 if self.__reduced:
123 fs=ReducedSolution(self.getDomain())
124 else:
125 fs=Solution(self.getDomain())
126
127 if H==None:
128 self.__H=Scalar(0.0, fs)
129 else:
130 self.__H=interpolate(H, fs)
131
132 def getTopography(self):
133 """
134 returns the current topography.
135 @rtype: scalar L{Data}
136 """
137 return self.__H
138
139 def getSafeTimeStepSize(self):
140 """
141 Returns the time step value.
142
143 @rtype: C{float}
144 """
145 if self.__dt == None:
146 self.__dt=0.5*inf(self.getDomain().getSize()/length(self.getVelocity()))
147 return self.__dt
148 def update(self,dt=None):
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 if dt>self.getSafeTimeStepSize():
161 raise ValueError("Time step must be less than safe time step size = %e."%self.getSafeTimeStepSize())
162
163 H=self.getTopography()
164 w=self.getVelocity()
165 w_tilda=1.*w
166 w_tilda[self.__DIM-1]=0
167 w_z=w[self.__DIM-1]
168
169 self.__PDE_H.setValue(X=((div(w_tilda*H)+w_z)*dt/2+H)*w_tilda*dt, Y=w_z*dt+H, r=H)
170 self.setTopography(self.__PDE_H.getSolution())
171
172 return self.getTopography()
173

  ViewVC Help
Powered by ViewVC 1.1.26