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

Annotation of /trunk/escript/py_src/mountains.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2564 - (hide annotations)
Tue Jul 28 04:46:38 2009 UTC (10 years, 1 month ago) by gross
File MIME type: text/x-python
File size: 5487 byte(s)
reduced PDE solutions added to topography
1 artak 2210
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 jfenwick 2344 __url__="https://launchpad.net/escript-finley"
21 artak 2210
22     from esys.escript import *
23 gross 2474 from esys.escript.linearPDEs import LinearPDE, SolverOptions
24 artak 2210 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 gross 2564 def __init__(self,domain,eps=0.01, reduced=False):
39 artak 2210 """
40     Sets up the level set method.
41    
42     @param domain: the domain where the mountains is used
43 gross 2563 @param eps: the smoothing parameter for (1)
44 gross 2564 @param reduced: if True topography reduced order is used for reconstruction of internal velocity and topography
45 artak 2210 """
46 gross 2563 if eps<=0:
47     raise ValueError("Smmoting parameter eps must be positive.")
48 artak 2210 self.__domain = domain
49     self.__eps = eps
50 gross 2564 self.__reduced=reduced
51 gross 2563 self.__DIM=domain.getDim()
52     z=domain.getX()[self.__DIM-1]
53    
54 artak 2210 self.__PDE_W = LinearPDE(domain)
55     self.__PDE_W.setSymmetryOn()
56 gross 2564 if reduced: self.__PDE_W.setReducedOrderOn()
57 gross 2563 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 artak 2210 self.__PDE_H = LinearPDE(domain)
62     self.__PDE_H.setSymmetryOn()
63 gross 2564 if reduced: self.__PDE_H.setReducedOrderOn()
64 artak 2210 self.__PDE_H.setValue(D=1.0)
65 gross 2474 self.__PDE_H.getSolverOptions().setSolverMethod(SolverOptions.LUMPING)
66 gross 2563 self.__PDE_H.setValue(q=whereZero(inf(z)-z))
67 artak 2210
68 gross 2563 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 artak 2210 """
83 gross 2563 Returns the domain.
84     """
85     return self.__domain
86 artak 2210
87 gross 2563 def setVelocity(self,v=None):
88 artak 2210 """
89 gross 2563 set a new velocity. v is define on the entire domain but only the surface values are used.
90 artak 2210
91 gross 2563 @param v: velocity field. If None zero is used.
92     @type v: vector
93     """
94     self.__dt=None
95 gross 2564 if self.__reduced:
96     fs=ReducedSolution(self.getDomain())
97     else:
98     fs=Solution(self.getDomain())
99     self.__v=Vector(0.,fs)
100 gross 2563 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 artak 2210
115 gross 2563 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 artak 2210
119 gross 2563 @param H: the topography. If None zero is used.
120     @type H: scalar
121     """
122 gross 2564 if self.__reduced:
123     fs=ReducedSolution(self.getDomain())
124     else:
125     fs=Solution(self.getDomain())
126 artak 2210
127 gross 2563 if H==None:
128 gross 2564 self.__H=Scalar(0.0, fs)
129 gross 2563 else:
130 gross 2564 self.__H=interpolate(H, fs)
131 gross 2563
132     def getTopography(self):
133     """
134     returns the current topography.
135     @rtype: scalar L{Data}
136     """
137     return self.__H
138 artak 2210
139 gross 2563 def getSafeTimeStepSize(self):
140 artak 2210 """
141 artak 2305 Returns the time step value.
142 gross 2563
143     @rtype: C{float}
144 artak 2210 """
145 gross 2563 if self.__dt == None:
146     self.__dt=0.5*inf(self.getDomain().getSize()/length(self.getVelocity()))
147 artak 2210 return self.__dt
148 gross 2563 def update(self,dt=None):
149     """
150     Sets a new W and updates the H function.
151 artak 2210
152 gross 2563 @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 artak 2210 """
156 gross 2563 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 artak 2210
169 gross 2563 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