/[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 2627 - (show annotations)
Mon Aug 24 01:57:34 2009 UTC (10 years, 6 months ago) by gross
File MIME type: text/x-python
File size: 5728 byte(s)
some bugs in pycad fixed.

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]=(0.3*(sup(z)-inf(z))/log(2.))**2
65 # A[self.__DIM-1,self.__DIM-1]=(sup(FunctionOnBoundary(self.__domain).getSize())/log(2.))**2
66 self.__PDE_W.setValue(D=1, A=A, q=whereZero(sup(z)-z))
67
68 self.__PDE_H = LinearPDE(domain)
69 self.__PDE_H.setSymmetryOn()
70 if reduced: self.__PDE_H.setReducedOrderOn()
71 # A=kronecker(domain)*0
72 # A[self.__DIM-1,self.__DIM-1]=0.1
73 self.__PDE_H.setValue(D=1.0)
74 # self.__PDE_H.getSolverOptions().setSolverMethod(SolverOptions.LUMPING)
75
76 self.setVelocity()
77 self.setTopography()
78 def getSolverOptionsForSmooting(self):
79 """
80 returns the solver options for the smoothing/extrapolation
81 """
82 return self.__PDE_W.getSolverOptions()
83
84 def getSolverOptionsForUpdate(self):
85 """
86 returns the solver options for the topograthy update
87 """
88 return self.__PDE_H.getSolverOptions()
89 def getDomain(self):
90 """
91 Returns the domain.
92 """
93 return self.__domain
94
95 def setVelocity(self,v=None):
96 """
97 set a new velocity. v is define on the entire domain but only the surface values are used.
98
99 :param v: velocity field. If None zero is used.
100 :type v: vector
101 """
102 self.__dt=None
103 self.__v=Vector(0.,Solution(self.getDomain()))
104 if not v == None:
105 for d in range(self.__DIM):
106 self.__PDE_W.setValue(r=v[d])
107 self.__v[d]=self.__PDE_W.getSolution()
108 def getVelocity(self):
109 """
110 returns the smoothed/extrapolated velocity
111 :rtype: vector `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 `Data`
136 """
137 return self.__H
138
139 def getSafeTimeStepSize(self):
140 """
141 Returns the time step value.
142
143 :rtype: ``float``
144 """
145 if self.__dt == None:
146 h=self.getDomain().getSize()
147 self.__dt=0.5*inf(h/length(interpolate(self.getVelocity(),h.getFunctionSpace())))
148 return self.__dt
149 def update(self,dt=None, allow_substeps=True):
150 """
151 Sets a new W and updates the H function.
152
153 :param dt: time step forward. If None the save time step size is used.
154 :type dt: positve ``float`` which is less or equal than the safe time step size.
155
156 """
157 if dt == None:
158 dt = self.getSafeTimeStepSize()
159 if dt<=0:
160 raise ValueError("Time step size must be positive.")
161 dt_safe=self.getSafeTimeStepSize()
162 n=max(int(math.ceil(dt/dt_safe)+0.5),1)
163 if n>1 and not allow_substeps:
164 raise SubSteppingException,"Substepping required."
165 dt/=n
166
167 H=self.getTopography()
168 w=self.getVelocity()
169 w_tilda=1.*w
170 w_tilda[self.__DIM-1]=0
171 w_z=w[self.__DIM-1]
172
173 t=0
174 for i in range(n):
175 L=integrate(w_z*dt+H)/vol(self.__PDE_H.getDomain())
176 self.__PDE_H.setValue(X=(inner(w_tilda,grad(H))*dt/2+H)*w_tilda*dt, Y=w_z*dt+H-L)
177 H=self.__PDE_H.getSolution()
178 print "DDD : ava = ",integrate(H)
179 t+=dt
180 self.setTopography(H)
181
182 return self.getTopography()

  ViewVC Help
Powered by ViewVC 1.1.26