# Contents of /trunk/downunder/py_src/forwardmodels/pressure.py

Revision 5707 - (show annotations)
Mon Jun 29 03:59:06 2015 UTC (3 years, 8 months ago) by sshaw
File MIME type: text/x-python
File size: 4037 byte(s)
```adding copyright headers to files without copyright info, moved header to top of file in some cases where it wasn't
```
 1 ############################################################################## 2 # 3 # Copyright (c) 2003-2015 by The University of Queensland 4 5 # 6 # Primary Business: Queensland, Australia 7 # Licensed under the Open Software License version 3.0 8 9 # 10 # Development until 2012 by Earth Systems Science Computational Center (ESSCC) 11 # Development 2012-2013 by School of Earth Sciences 12 # Development from 2014 by Centre for Geoscience Computing (GeoComp) 13 # 14 ############################################################################## 15 16 """Isostatic Pressure calculation""" 17 from __future__ import division, print_function 18 19 __copyright__="""Copyright (c) 2003-2015 by The University of Queensland 20 http://www.uq.edu.au 21 Primary Business: Queensland, Australia""" 22 __license__="""Licensed under the Open Software License version 3.0 23 24 __url__= 25 26 __all__ = ['IsostaticPressure'] 27 28 from esys.downunder.coordinates import makeTranformation 29 from esys.escript import unitsSI as U 30 from esys.escript import Scalar, Vector, Function 31 from esys.escript.linearPDEs import LinearSinglePDE 32 from esys.escript.util import * 33 from math import pi as PI 34 35 36 class IsostaticPressure(object): 37 """ 38 class to calculate isostatic pressure field correction due to gravity forces 39 """ 40 def __init__(self, domain, p0=0., level0=0, gravity0=-9.81*U.m*U.sec**(-2), 41 background_density=2670* U.kg*U.m**(-3), 42 gravity_constant=U.Gravitational_Constant, 43 coordinates=None, tol=1e-8): 44 """ 45 :param domain: domain of the model 46 :type domain: `Domain` 47 :param p0: pressure at level0 48 :type p0: scalar `Data` or ``float`` 49 :param background_density: defines background_density in kg/m^3 50 :type background_density: ``float`` 51 :param coordinates: defines coordinate system to be used 52 :type coordinates: ReferenceSystem` or `SpatialCoordinateTransformation` 53 :param tol: tolerance of underlying PDE 54 :type tol: positive ``float`` 55 :param level0: pressure for z>=`level0` is set to zero. 56 :type level0: ``float`` 57 :param gravity0: vertical background gravity at `level0` 58 :type gravity0: ``float`` 59 """ 60 DIM=domain.getDim() 61 self.__domain = domain 62 self.__trafo=makeTranformation(domain, coordinates) 63 self.__pde=LinearSinglePDE(domain) 64 self.__pde.getSolverOptions().setTolerance(tol) 65 self.__pde.setSymmetryOn() 66 67 z = domain.getX()[DIM-1] 68 self.__pde.setValue(q=whereNonNegative(z-level0), r=p0) 69 70 fw = self.__trafo.getScalingFactors()**2 * self.__trafo.getVolumeFactor() 71 A=self.__pde.createCoefficient("A") 72 for i in range(DIM): A[i,i]=fw[i] 73 self.__pde.setValue(A=A) 74 z = Function(domain).getX()[DIM-1] 75 self.__g_b= 4*PI*gravity_constant/self.__trafo.getScalingFactors()[DIM-1]*background_density*(level0-z) + gravity0 76 self.__rho_b=background_density 77 78 def getPressure(self, g = None, rho=None): 79 """ 80 return the pressure for gravity force anomaly `g` and 81 density anomaly `rho` 82 83 :param g: gravity anomaly data 84 :type g: ``Vector`` 85 :param rho: gravity anomaly data 86 :type rho: ``Scalar`` 87 :return: pressure distribution 88 :rtype: ``Scalar` 89 """ 90 if not g: g=Vector(0., Function(self.__domain)) 91 if not rho: rho=Scalar(0., Function(self.__domain)) 92 93 g2=(rho * self.__g_b)*[0,0,1] + self.__rho_b*g + rho*g 94 # Tests need to be updated before the following is uncommented: 95 #g2=((rho+self.__rho_b) * self.__g_b)*[0,0,1] + self.__rho_b*g + rho*g 96 d=self.__trafo.getScalingFactors() 97 V= self.__trafo.getVolumeFactor() 98 self.__pde.setValue(X = -g2*d*V) 99 #self.__pde.setValue(X = g2*d*V) 100 return self.__pde.getSolution() 101