Contents of /trunk/escriptcore/py_src/mountains.py

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

 ViewVC Help Powered by ViewVC 1.1.26