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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26