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 |
|