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

  ViewVC Help
Powered by ViewVC 1.1.26