/[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 2563 - (show annotations)
Tue Jul 28 03:50:45 2009 UTC (10 years, 1 month ago) by gross
File MIME type: text/x-python
File size: 5063 byte(s)
reengineered version of the mountain class so it becomes easier to use for the convection code
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):
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 """
45 if eps<=0:
46 raise ValueError("Smmoting parameter eps must be positive.")
47 self.__domain = domain
48 self.__eps = eps
49 self.__DIM=domain.getDim()
50 z=domain.getX()[self.__DIM-1]
51
52 self.__PDE_W = LinearPDE(domain)
53 self.__PDE_W.setSymmetryOn()
54 A=kronecker(domain)
55 A[self.__DIM-1,self.__DIM-1]=1/self.__eps
56 self.__PDE_W.setValue(A=A, q=whereZero(sup(z)-z)+whereZero(inf(z)-z))
57
58 self.__PDE_H = LinearPDE(domain)
59 self.__PDE_H.setSymmetryOn()
60 self.__PDE_H.setValue(D=1.0)
61 self.__PDE_H.getSolverOptions().setSolverMethod(SolverOptions.LUMPING)
62 self.__PDE_H.setValue(q=whereZero(inf(z)-z))
63
64 self.setVelocity()
65 self.setTopography()
66 def getSolverOptionsForSmooting(self):
67 """
68 returns the solver options for the smoothing/extrapolation
69 """
70 return self.__PDE_W.getSolverOptions()
71
72 def getSolverOptionsForUpdate(self):
73 """
74 returns the solver options for the topograthy update
75 """
76 return self.__PDE_H.getSolverOptions()
77 def getDomain(self):
78 """
79 Returns the domain.
80 """
81 return self.__domain
82
83 def setVelocity(self,v=None):
84 """
85 set a new velocity. v is define on the entire domain but only the surface values are used.
86
87 @param v: velocity field. If None zero is used.
88 @type v: vector
89 """
90 self.__dt=None
91 self.__v=Vector(0.,Solution(self.getDomain()))
92 if not v == None:
93 z=self.getDomain().getX()[self.__DIM-1]
94 z_min=inf(z)
95 z_max=sup(z)
96 f=(z-z_min)/(z_max-z_min)
97 for d in range(self.__DIM):
98 self.__PDE_W.setValue(r=v[d]*f)
99 self.__v[d]=self.__PDE_W.getSolution()
100 def getVelocity(self):
101 """
102 returns the smoothed/extrapolated velocity
103 @rtype: vector L{Data}
104 """
105 return self.__v
106
107 def setTopography(self,H=None):
108 """
109 set the topography to H where H defines the vertical displacement. H is defined for the entire domain.
110
111 @param H: the topography. If None zero is used.
112 @type H: scalar
113 """
114
115 if H==None:
116 self.__H=Scalar(0.0, Solution(self.getDomain()))
117 else:
118 self.__H=interpolate(H, Solution(self.getDomain()))
119
120 def getTopography(self):
121 """
122 returns the current topography.
123 @rtype: scalar L{Data}
124 """
125 return self.__H
126
127 def getSafeTimeStepSize(self):
128 """
129 Returns the time step value.
130
131 @rtype: C{float}
132 """
133 if self.__dt == None:
134 self.__dt=0.5*inf(self.getDomain().getSize()/length(self.getVelocity()))
135 return self.__dt
136 def update(self,dt=None):
137 """
138 Sets a new W and updates the H function.
139
140 @param dt: time step forward. If None the save time step size is used.
141 @type dt: positve C{float} which is less or equal than the safe time step size.
142
143 """
144 if dt == None:
145 dt = self.getSafeTimeStepSize()
146 if dt<=0:
147 raise ValueError("Time step size must be positive.")
148 if dt>self.getSafeTimeStepSize():
149 raise ValueError("Time step must be less than safe time step size = %e."%self.getSafeTimeStepSize())
150
151 H=self.getTopography()
152 w=self.getVelocity()
153 w_tilda=1.*w
154 w_tilda[self.__DIM-1]=0
155 w_z=w[self.__DIM-1]
156
157 self.__PDE_H.setValue(X=((div(w_tilda*H)+w_z)*dt/2+H)*w_tilda*dt, Y=w_z*dt+H, r=H)
158 self.setTopography(self.__PDE_H.getSolution())
159
160 return self.getTopography()
161

  ViewVC Help
Powered by ViewVC 1.1.26