/[escript]/trunk/modellib/py_src/crustal/setups.py
ViewVC logotype

Annotation of /trunk/modellib/py_src/crustal/setups.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1809 - (hide annotations)
Thu Sep 25 06:43:44 2008 UTC (10 years, 11 months ago) by ksteube
File MIME type: text/x-python
File size: 9344 byte(s)
Copyright updated in all python files

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

  ViewVC Help
Powered by ViewVC 1.1.26