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

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