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

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 gross 2620 import math
26 artak 2210 import sys
27    
28 gross 2620 class SubSteppingException(Exception):
29     """
30     Thrown if the L{Mountains} class uses substepping.
31     """
32     pass
33    
34    
35 artak 2210 class Mountains:
36     """
37     The Mountains class is defined by the following equations:
38    
39 gross 2620 (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 artak 2210
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 gross 2564 def __init__(self,domain,eps=0.01, reduced=False):
47 artak 2210 """
48     Sets up the level set method.
49    
50 jfenwick 2625 :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 artak 2210 """
54 gross 2620 if eps<0:
55     raise ValueError("Smooting parameter eps must be non-negative.")
56 artak 2210 self.__domain = domain
57 gross 2564 self.__reduced=reduced
58 gross 2563 self.__DIM=domain.getDim()
59     z=domain.getX()[self.__DIM-1]
60    
61 artak 2210 self.__PDE_W = LinearPDE(domain)
62     self.__PDE_W.setSymmetryOn()
63 gross 2620 A=kronecker(domain)*eps*0
64 gross 2627 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 gross 2620 self.__PDE_W.setValue(D=1, A=A, q=whereZero(sup(z)-z))
67 gross 2563
68 artak 2210 self.__PDE_H = LinearPDE(domain)
69     self.__PDE_H.setSymmetryOn()
70 gross 2564 if reduced: self.__PDE_H.setReducedOrderOn()
71 gross 2627 # A=kronecker(domain)*0
72     # A[self.__DIM-1,self.__DIM-1]=0.1
73 artak 2210 self.__PDE_H.setValue(D=1.0)
74 gross 2620 # self.__PDE_H.getSolverOptions().setSolverMethod(SolverOptions.LUMPING)
75 artak 2210
76 gross 2563 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 artak 2210 """
91 gross 2563 Returns the domain.
92     """
93     return self.__domain
94 artak 2210
95 gross 2563 def setVelocity(self,v=None):
96 artak 2210 """
97 gross 2563 set a new velocity. v is define on the entire domain but only the surface values are used.
98 artak 2210
99 jfenwick 2625 :param v: velocity field. If None zero is used.
100     :type v: vector
101 gross 2563 """
102     self.__dt=None
103 gross 2620 self.__v=Vector(0.,Solution(self.getDomain()))
104 gross 2563 if not v == None:
105     for d in range(self.__DIM):
106 gross 2620 self.__PDE_W.setValue(r=v[d])
107 gross 2563 self.__v[d]=self.__PDE_W.getSolution()
108     def getVelocity(self):
109     """
110     returns the smoothed/extrapolated velocity
111 jfenwick 2625 :rtype: vector `Data`
112 gross 2563 """
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 jfenwick 2625 :param H: the topography. If None zero is used.
120     :type H: scalar
121 gross 2563 """
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 jfenwick 2625 :rtype: scalar `Data`
136 gross 2563 """
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 jfenwick 2625 :rtype: ``float``
144 artak 2210 """
145 gross 2563 if self.__dt == None:
146 gross 2620 h=self.getDomain().getSize()
147     self.__dt=0.5*inf(h/length(interpolate(self.getVelocity(),h.getFunctionSpace())))
148 artak 2210 return self.__dt
149 gross 2620 def update(self,dt=None, allow_substeps=True):
150 gross 2563 """
151     Sets a new W and updates the H function.
152 artak 2210
153 jfenwick 2625 :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 gross 2563
156 artak 2210 """
157 gross 2563 if dt == None:
158     dt = self.getSafeTimeStepSize()
159     if dt<=0:
160     raise ValueError("Time step size must be positive.")
161 gross 2620 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 gross 2563 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 artak 2210
173 gross 2620 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 gross 2563
182     return self.getTopography()

  ViewVC Help
Powered by ViewVC 1.1.26