| 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 |
|
Parent Directory
| 
Revision Log