# Contents of /trunk/escript/py_src/mountains.py

Revision 3432 - (show annotations)
Fri Jan 7 01:32:07 2011 UTC (10 years, 4 months ago) by jfenwick
File MIME type: text/x-python
File size: 6422 byte(s)
```Made import statements a bit more specific to clean up the epydoc
```
 1 2 ######################################################## 3 # 4 # Copyright (c) 2003-2010 by University of Queensland 5 # Earth Systems Science Computational Center (ESSCC) 6 7 # 8 # Primary Business: Queensland, Australia 9 # Licensed under the Open Software License version 3.0 10 11 # 12 ######################################################## 13 14 __copyright__="""Copyright (c) 2003-2010 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 20 __url__= 21 22 from esys.escript import inf, atan, Lsup, log, whereZero, Vector, Solution, Scalar, length, interpolate, vol, integrate 23 import esys.escript as escript 24 from esys.escript.linearPDEs import LinearPDE, SolverOptions 25 26 import math 27 import sys 28 29 class SubSteppingException(Exception): 30 """ 31 Thrown if the L{Mountains} class uses substepping. 32 """ 33 pass 34 35 36 class Mountains: 37 """ 38 The Mountains class is defined by the following equations: 39 40 (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. 41 42 (2) Integration of topography PDE using Taylor-Galerkin upwinding to stabilize the advection terms 43 H^(t+dt)=H^t+dt*w_3+w_hat*dt*[(div(w_hat*H^t)+w_3)+(dt/2)+H^t], 44 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 (?) 45 46 """ 47 def __init__(self,domain,eps=0.01): 48 """ 49 Sets up the level set method. 50 51 :param domain: the domain where the mountains is used 52 :param eps: the smoothing parameter for (1) 53 """ 54 order=escript.Solution(domain).getApproximationOrder() 55 if order>1: 56 reduced = True 57 if ReducedSolution(domain).getApproximationOrder()>1: raise ValueError,"Reduced order needs to be equal to 1." 58 else: 59 reduced = False 60 if eps<0: 61 raise ValueError("Smooting parameter eps must be non-negative.") 62 self.__domain = domain 63 self.__reduced=reduced 64 self.__DIM=domain.getDim() 65 z=domain.getX()[self.__DIM-1] 66 67 self.__PDE_W = LinearPDE(domain) 68 self.__PDE_W.setSymmetryOn() 69 A=escript.kronecker(domain)*eps*0 70 A[self.__DIM-1,self.__DIM-1]=(0.3*(escript.sup(z)-inf(z))/log(2.))**2 71 # A[self.__DIM-1,self.__DIM-1]=(sup(FunctionOnBoundary(self.__domain).getSize())/log(2.))**2 72 self.__PDE_W.setValue(D=1, A=A, q=whereZero(escript.sup(z)-z)+whereZero(inf(z)-z)) 73 74 self.__PDE_H = LinearPDE(domain) 75 self.__PDE_H.setSymmetryOn() 76 if reduced: self.__PDE_H.setReducedOrderOn() 77 # A=kronecker(domain)*0 78 # A[self.__DIM-1,self.__DIM-1]=0.1 79 self.__PDE_H.setValue(D=1.0, q=whereZero(inf(z)-z)) 80 # self.__PDE_H.getSolverOptions().setSolverMethod(SolverOptions.LUMPING) 81 82 self.setVelocity() 83 self.setTopography() 84 def getSolverOptionsForSmooting(self): 85 """ 86 returns the solver options for the smoothing/extrapolation 87 """ 88 return self.__PDE_W.getSolverOptions() 89 90 def getSolverOptionsForUpdate(self): 91 """ 92 returns the solver options for the topograthy update 93 """ 94 return self.__PDE_H.getSolverOptions() 95 def getDomain(self): 96 """ 97 Returns the domain. 98 """ 99 return self.__domain 100 101 def setVelocity(self,v=None): 102 """ 103 set a new velocity. v is define on the entire domain but only the surface values are used. 104 105 :param v: velocity field. If None zero is used. 106 :type v: vector 107 """ 108 self.__dt=None 109 self.__v=Vector(0.,Solution(self.getDomain())) 110 if not v == None: 111 xi=self.getDomain().getX()[self.getDomain().getDim()-1] 112 v=(xi-inf(xi))/(escript.sup(xi)-inf(xi))*v 113 for d in range(self.__DIM): 114 self.__PDE_W.setValue(r=v[d]) 115 self.__v[d]=self.__PDE_W.getSolution() 116 def getVelocity(self): 117 """ 118 returns the smoothed/extrapolated velocity 119 :rtype: vector `Data` 120 """ 121 return self.__v 122 123 def setTopography(self,H=None): 124 """ 125 set the topography to H where H defines the vertical displacement. H is defined for the entire domain. 126 127 :param H: the topography. If None zero is used. 128 :type H: scalar 129 """ 130 if self.__reduced: 131 fs=ReducedSolution(self.getDomain()) 132 else: 133 fs=Solution(self.getDomain()) 134 135 if H==None: 136 self.__H=Scalar(0.0, fs) 137 else: 138 self.__H=interpolate(H, fs) 139 140 def getTopography(self): 141 """ 142 returns the current topography. 143 :rtype: scalar `Data` 144 """ 145 return self.__H 146 147 def getSafeTimeStepSize(self): 148 """ 149 Returns the time step value. 150 151 :rtype: ``float`` 152 """ 153 if self.__dt == None: 154 h=self.getDomain().getSize() 155 self.__dt=0.5*inf(h/length(interpolate(self.getVelocity(),h.getFunctionSpace()))) 156 return self.__dt 157 def update(self,dt=None, allow_substeps=True): 158 """ 159 Sets a new W and updates the H function. 160 161 :param dt: time step forward. If None the save time step size is used. 162 :type dt: positve ``float`` which is less or equal than the safe time step size. 163 164 """ 165 if dt == None: 166 dt = self.getSafeTimeStepSize() 167 if dt<=0: 168 raise ValueError("Time step size must be positive.") 169 dt_safe=self.getSafeTimeStepSize() 170 n=max(int(math.ceil(dt/dt_safe)+0.5),1) 171 if n>1 and not allow_substeps: 172 raise SubSteppingException,"Substepping required." 173 dt/=n 174 175 H=self.getTopography() 176 w=self.getVelocity() 177 w_tilda=1.*w 178 w_tilda[self.__DIM-1]=0 179 w_z=w[self.__DIM-1] 180 V=vol(self.__PDE_H.getDomain()) 181 182 t=0 183 for i in range(n): 184 # L=integrate(w_z*dt+H)/vol(self.__PDE_H.getDomain()) 185 # self.__PDE_H.setValue(X=(inner(w_tilda,grad(H))*dt/2+H)*w_tilda*dt, Y=w_z*dt+H-L) 186 # H=self.__PDE_H.getSolution() 187 L=integrate(w_z*dt+H)/V 188 self.__PDE_H.setValue(X=w_tilda*H*(dt/2), Y=w_z*(dt/2)+H-L) 189 Hhalf=self.__PDE_H.getSolution() 190 self.__PDE_H.setValue(X=w_tilda*Hhalf*dt, Y=w_z*dt+H-L) 191 H=self.__PDE_H.getSolution() 192 print "DDD : ava = ",integrate(H) 193 t+=dt 194 self.setTopography(H) 195 196 return self.getTopography()