py_src/crustal/setups.py

"""
set up of mining activities in modelframe
                                                                                                                                                                                                     
@var __author__: name of author
@var __licence__: licence agreement
@var __url__: url entry point on documentation
@var __version__: version
@var __date__: date of the version
"""

__copyright__="""  Copyright (c) 2006 by ACcESS MNRF
                    http://www.access.edu.au
                Primary Business: Queensland, Australia"""
__license__="""Licensed under the Open Software License version 3.0
             http://www.opensource.org/licenses/osl-3.0.php"""
__author__="Lutz Gross, l.gross@uq.edu.au"
__url__="http://www.iservo.edu.au/esys/escript"
__version__="$Revision$"
__date__="$Date$"

from esys.escript.modelframe import Model
from esys.escript import *
from esys.escript.linearPDEs import LinearPDE
from mines import parse

class MiningHistory(Model):
    """
    manages the history of mines. It mandels the time steping according to the available 
    data and a dictionary of mass changes per year for all the mines. 

    @ivar history: mine site history file.
    @type history: C{DataSource}
    @ivar mass_changes: current mass change per year.
    @type mass_changes: C{DataSource}

    """
    def __init__(self,**kwargs):
        """
        """
        super(MiningHistory,self).__init__(**kwargs)
        self.declareParameter(t=0.,
                              history=None,
                              mass_changes=None)

    def doInitialization(self):
        """
        initialize time integration
        """
        self.__minesite=parse(open(self.history.getLocalFileName(),'r'))
        self.mass_changes=self.__minesite.getMassChanges(self.t)

    def doStepPreprocessing(self, dt):
        self.mass_changes=self.__minesite.getMassChanges(self.t)

    def getSafeTimeStepSize(self, dt):
        return self.__minesite.getNextTimeMarker(self.t)-self.t


class DensityChange(Model):
    """
    translates local mass extraction into density changes.
    "local" refers a taged region.

    @ivar domain: mining region
    @type domain: L{Domian}
    @ivar tag_map: a tagmap 
    @type tag_map: L{TagMap}
    @ivar mass_rate: rate of change mass in a tagged region
    @type mass_rate: C{dict}
    @ivar density_rate: density in each region
    @type density_rate:  C{Data}
    """
    def __init__(self,**kwargs):
        """
        """
        super(DensityChange,self).__init__(**kwargs)
        self.declareParameter(domain=None,
                              tag_map=None,
                              mass_rate={},
                              density_rate=None)

    def doInitialization(self):
        """
        initialize time integration
        """
        self.__volumes={}
        for i in self.tag_map.getNames():
            c=Scalar(0.,Function(self.domain))
            self.tag_map.insert(c,**{i : 1.})
            c=integrate(c)
            if c>0: self.__volumes[i] = c

    def doStepPreprocessing(self, dt):
         self.density_rate=Scalar(0.,Function(self.domain))
         d={}
         for i in self.__volumes:
            if self.mass_rate.has_key(i): d[i]=-self.mass_rate[i]/self.__volumes[i]
         self.tag_map.insert(self.density_rate,**d)

class LinearElasticStressChange(Model):
    """
    calculates the stress according to an initial garvity field and a consecutive 
    chenge of density  based on linar elastic material model. It is assumed that 
    the lame coefficients don't change over time, the direction of gravity is
    pointing into the x_{dim} direction. 
  
    @note: add stress changes due to tectonic loading and slip

    @ivar domain: mining region
    @type domain: L{Domian}
    @ivar tag_map: a tagmap 
    @type tag_map: L{TagMap}
    @ivar displacement: displacement field
    @type displacement: L{Vector} or C{None}
    @ivar stress: displacement field 
    @type stress: L{Vector}  or C{None}
    @ivar density: initial density distribution
    @type density: C{float} or L{Scalar}
    @ivar density_rate: density rate by tag (may be changed of time)
    @type density_rate:  C{dict}
    @ivar lame_lambda: elasticity coefficient lambda (assumed to be constant over time)
    @type lame_lambda: C{float} or L{Scalar}
    @ivar lame_mu: elasticity coefficient mu (assumed to be constant over time)
    @type lame_mu: C{float} or L{Scalar}
    @ivar location_of_fixed_displacement: mask of locations and component with zero displacements
    @type location_of_fixed_displacement: L{Vector} or C{None}
    """
    def __init__(self,**kwargs):
        """
        """
        super(LinearElasticStressChange,self).__init__(**kwargs)
        self.declareParameter(domain=None,
                              tag_map=None,
                              displacement=None, \
                              stress=None, \
                              density=1., \
                              lame_lambda=2., \
                              lame_mu=1., \
                              location_of_fixed_displacement=None, \
                              density_rate=None)

    def doInitialization(self):
        """
        initialize time integration
        """
        if self.stress == None: self.stress=Tensor(0,Function(self.domain))
        if self.displacement==None: self.displacement=Scalar(0,Solution(self.domain))
        self.__pde=LinearPDE(self.domain)
        self.__pde.setSymmetryOn()
        A =Tensor4(0.,Function(self.domain))
        for i in range(self.domain.getDim()):
           for j in range(self.domain.getDim()):
              A[i,j,j,i] += self.lame_mu
              A[i,j,i,j] += self.lame_mu
              A[i,i,j,j] += self.lame_lambda
        self.__pde.setValue(A=A,q=self.location_of_fixed_displacement)

        self.__stress=self.stress
        self.__displacement=self.displacement
        self.__first=True
        
    def doInitialStep(self):
        """
        """
        d=self.domain.getDim()
        self.__pde.setValue(Y=-kronecker(Function(self.domain))[d-1]*9.81*self.density)
        ddisp=self.__pde.getSolution()
        deps=symmetric(grad(ddisp))
        self.stress=self.stress+self.lame_mu*deps+self.lame_lambda*trace(deps)*kronecker(Function(self.domain))
        self.displacement=self.displacement+ddisp
        self.__first=False

    def terminateInitialIteration(self):
        return not self.__first

    def terminateIteration(self):
        return not self.__first

    def doInitialPostprocessing(self):
        self.__stress=self.stress
        self.__displacement=self.displacement
        self.trace("initial displacement range : %s: %s"%(inf(self.displacement[2]),sup(self.displacement[2])))
        self.trace("initial stress range : %s: %s"%(inf(self.stress),sup(self.stress)))

    def doStepPreprocessing(self,dt):
        self.stress=self.__stress
        self.displacement=self.__displacement
        self.__first=True
         
    def doStep(self,dt):
        """
        """
        d=self.domain.getDim()
        if not self.density_rate == None:
           self.__pde.setValue(Y=-kronecker(Function(self.domain))[d-1]*9.81*self.density_rate)
        ddisp=self.__pde.getSolution()
        deps=symmetric(grad(ddisp))
        self.stress=self.stress+dt*(self.lame_mu*deps+self.lame_lambda*trace(deps)*kronecker(Function(self.domain)))
        self.displacement=self.displacement+dt*ddisp
        self.__first=False

    def doStepPostprocessing(self, dt):
        self.__stress=self.stress
        self.__displacement=self.displacement
        self.trace("displacement range : %s: %s"%(inf(self.displacement[2]),sup(self.displacement[2])))
        self.trace("stress range : %s: %s"%(inf(self.stress),sup(self.stress)))
1	"""
2	set up of mining activities in modelframe
3
4	@var __author__: name of author
5	@var __licence__: licence agreement
6	@var __url__: url entry point on documentation
7	@var __version__: version
8	@var __date__: date of the version
9	"""
10
11	__copyright__=""" Copyright (c) 2006 by ACcESS MNRF
12	http://www.access.edu.au
13	Primary Business: Queensland, Australia"""
14	__license__="""Licensed under the Open Software License version 3.0
15	http://www.opensource.org/licenses/osl-3.0.php"""
16	__author__="Lutz Gross, l.gross@uq.edu.au"
17	__url__="http://www.iservo.edu.au/esys/escript"
18	__version__="$Revision$"
19	__date__="$Date$"
20
21	from esys.escript.modelframe import Model
22	from esys.escript import *
23	from esys.escript.linearPDEs import LinearPDE
24	from mines import parse
25
26	class MiningHistory(Model):
27	"""
28	manages the history of mines. It mandels the time steping according to the available
29	data and a dictionary of mass changes per year for all the mines.
30
31	@ivar history: mine site history file.
32	@type history: C{DataSource}
33	@ivar mass_changes: current mass change per year.
34	@type mass_changes: C{DataSource}
35
36	"""
37	def __init__(self,**kwargs):
38	"""
39	"""
40	super(MiningHistory,self).__init__(**kwargs)
41	self.declareParameter(t=0.,
42	history=None,
43	mass_changes=None)
44
45	def doInitialization(self):
46	"""
47	initialize time integration
48	"""
49	self.__minesite=parse(open(self.history.getLocalFileName(),'r'))
50	self.mass_changes=self.__minesite.getMassChanges(self.t)
51
52	def doStepPreprocessing(self, dt):
53	self.mass_changes=self.__minesite.getMassChanges(self.t)
54
55	def getSafeTimeStepSize(self, dt):
56	return self.__minesite.getNextTimeMarker(self.t)-self.t
57
58
59	class DensityChange(Model):
60	"""
61	translates local mass extraction into density changes.
62	"local" refers a taged region.
63
64	@ivar domain: mining region
65	@type domain: L{Domian}
66	@ivar tag_map: a tagmap
67	@type tag_map: L{TagMap}
68	@ivar mass_rate: rate of change mass in a tagged region
69	@type mass_rate: C{dict}
70	@ivar density_rate: density in each region
71	@type density_rate: C{Data}
72	"""
73	def __init__(self,**kwargs):
74	"""
75	"""
76	super(DensityChange,self).__init__(**kwargs)
77	self.declareParameter(domain=None,
78	tag_map=None,
79	mass_rate={},
80	density_rate=None)
81
82	def doInitialization(self):
83	"""
84	initialize time integration
85	"""
86	self.__volumes={}
87	for i in self.tag_map.getNames():
88	c=Scalar(0.,Function(self.domain))
89	self.tag_map.insert(c,**{i : 1.})
90	c=integrate(c)
91	if c>0: self.__volumes[i] = c
92
93	def doStepPreprocessing(self, dt):
94	self.density_rate=Scalar(0.,Function(self.domain))
95	d={}
96	for i in self.__volumes:
97	if self.mass_rate.has_key(i): d[i]=-self.mass_rate[i]/self.__volumes[i]
98	self.tag_map.insert(self.density_rate,**d)
99
100	class LinearElasticStressChange(Model):
101	"""
102	calculates the stress according to an initial garvity field and a consecutive
103	chenge of density based on linar elastic material model. It is assumed that
104	the lame coefficients don't change over time, the direction of gravity is
105	pointing into the x_{dim} direction.
106
107	@note: add stress changes due to tectonic loading and slip
108
109	@ivar domain: mining region
110	@type domain: L{Domian}
111	@ivar tag_map: a tagmap
112	@type tag_map: L{TagMap}
113	@ivar displacement: displacement field
114	@type displacement: L{Vector} or C{None}
115	@ivar stress: displacement field
116	@type stress: L{Vector} or C{None}
117	@ivar density: initial density distribution
118	@type density: C{float} or L{Scalar}
119	@ivar density_rate: density rate by tag (may be changed of time)
120	@type density_rate: C{dict}
121	@ivar lame_lambda: elasticity coefficient lambda (assumed to be constant over time)
122	@type lame_lambda: C{float} or L{Scalar}
123	@ivar lame_mu: elasticity coefficient mu (assumed to be constant over time)
124	@type lame_mu: C{float} or L{Scalar}
125	@ivar location_of_fixed_displacement: mask of locations and component with zero displacements
126	@type location_of_fixed_displacement: L{Vector} or C{None}
127	"""
128	def __init__(self,**kwargs):
129	"""
130	"""
131	super(LinearElasticStressChange,self).__init__(**kwargs)
132	self.declareParameter(domain=None,
133	tag_map=None,
134	displacement=None, \
135	stress=None, \
136	density=1., \
137	lame_lambda=2., \
138	lame_mu=1., \
139	location_of_fixed_displacement=None, \
140	density_rate=None)
141
142	def doInitialization(self):
143	"""
144	initialize time integration
145	"""
146	if self.stress == None: self.stress=Tensor(0,Function(self.domain))
147	if self.displacement==None: self.displacement=Scalar(0,Solution(self.domain))
148	self.__pde=LinearPDE(self.domain)
149	self.__pde.setSymmetryOn()
150	A =Tensor4(0.,Function(self.domain))
151	for i in range(self.domain.getDim()):
152	for j in range(self.domain.getDim()):
153	A[i,j,j,i] += self.lame_mu
154	A[i,j,i,j] += self.lame_mu
155	A[i,i,j,j] += self.lame_lambda
156	self.__pde.setValue(A=A,q=self.location_of_fixed_displacement)
157
158	self.__stress=self.stress
159	self.__displacement=self.displacement
160	self.__first=True
161
162	def doInitialStep(self):
163	"""
164	"""
165	d=self.domain.getDim()
166	self.__pde.setValue(Y=-kronecker(Function(self.domain))[d-1]9.81self.density)
167	ddisp=self.__pde.getSolution()
168	deps=symmetric(grad(ddisp))
169	self.stress=self.stress+self.lame_mudeps+self.lame_lambdatrace(deps)*kronecker(Function(self.domain))
170	self.displacement=self.displacement+ddisp
171	self.__first=False
172
173	def terminateInitialIteration(self):
174	return not self.__first
175
176	def terminateIteration(self):
177	return not self.__first
178
179	def doInitialPostprocessing(self):
180	self.__stress=self.stress
181	self.__displacement=self.displacement
182	self.trace("initial displacement range : %s: %s"%(inf(self.displacement[2]),sup(self.displacement[2])))
183	self.trace("initial stress range : %s: %s"%(inf(self.stress),sup(self.stress)))
184
185	def doStepPreprocessing(self,dt):
186	self.stress=self.__stress
187	self.displacement=self.__displacement
188	self.__first=True
189
190	def doStep(self,dt):
191	"""
192	"""
193	d=self.domain.getDim()
194	if not self.density_rate == None:
195	self.__pde.setValue(Y=-kronecker(Function(self.domain))[d-1]9.81self.density_rate)
196	ddisp=self.__pde.getSolution()
197	deps=symmetric(grad(ddisp))
198	self.stress=self.stress+dt(self.lame_mudeps+self.lame_lambdatrace(deps)kronecker(Function(self.domain)))
199	self.displacement=self.displacement+dt*ddisp
200	self.__first=False
201
202	def doStepPostprocessing(self, dt):
203	self.__stress=self.stress
204	self.__displacement=self.displacement
205	self.trace("displacement range : %s: %s"%(inf(self.displacement[2]),sup(self.displacement[2])))
206	self.trace("stress range : %s: %s"%(inf(self.stress),sup(self.stress)))