/[escript]/trunk/downunder/py_src/forwardmodels/pressure.py
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5707 - (show annotations)
Mon Jun 29 03:59:06 2015 UTC (3 years, 6 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 # http://www.uq.edu.au
5 #
6 # Primary Business: Queensland, Australia
7 # Licensed under the Open Software License version 3.0
8 # http://www.opensource.org/licenses/osl-3.0.php
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 http://www.opensource.org/licenses/osl-3.0.php"""
24 __url__="https://launchpad.net/escript-finley"
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

  ViewVC Help
Powered by ViewVC 1.1.26