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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1000 - (show annotations)
Wed Feb 28 00:34:42 2007 UTC (14 years, 7 months ago) by gross
File MIME type: text/x-python
File size: 9487 byte(s)
coulomb stress failure added
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, ParameterSet
22 from esys.escript import *
23 from esys.escript.linearPDEs import LinearPDE
24 from mines import parse
25 import numarray
26
27 class MiningHistory(Model):
28 """
29 manages the history of mines. It mandels the time steping according to the available
30 data and a dictionary of mass changes per year for all the mines.
31
32 @ivar history: mine site history file.
33 @type history: C{DataSource}
34 @ivar mass_changes: current mass change per year.
35 @type mass_changes: C{DataSource}
36
37 """
38 def __init__(self,**kwargs):
39 """
40 """
41 super(MiningHistory,self).__init__(**kwargs)
42 self.declareParameter(t=0.,
43 history=None,
44 mass_changes=None)
45
46 def doInitialization(self):
47 """
48 initialize time integration
49 """
50 self.__minesite=parse(open(self.history.getLocalFileName(),'r'))
51 self.mass_changes=self.__minesite.getMassChanges(self.t)
52
53 def doStepPreprocessing(self, dt):
54 self.mass_changes=self.__minesite.getMassChanges(self.t)
55
56 def getSafeTimeStepSize(self, dt):
57 return self.__minesite.getNextTimeMarker(self.t)-self.t
58
59
60 class DensityChange(Model):
61 """
62 translates local mass extraction into density changes.
63 "local" refers a taged region.
64
65 @ivar domain: mining region
66 @type domain: L{Domian}
67 @ivar tag_map: a tagmap
68 @type tag_map: L{TagMap}
69 @ivar mass_rate: rate of change mass in a tagged region
70 @type mass_rate: C{dict}
71 @ivar density_rate: density in each region
72 @type density_rate: C{Data}
73 """
74 def __init__(self,**kwargs):
75 """
76 """
77 super(DensityChange,self).__init__(**kwargs)
78 self.declareParameter(domain=None,
79 tag_map=None,
80 mass_rate={},
81 density_rate=None)
82
83 def doInitialization(self):
84 """
85 initialize time integration
86 """
87 self.__volumes={}
88 for i in self.tag_map.getNames():
89 c=Scalar(0.,Function(self.domain))
90 self.tag_map.insert(c,**{i : 1.})
91 c=integrate(c)
92 if c>0: self.__volumes[i] = c
93
94 def doStepPreprocessing(self, dt):
95 self.density_rate=Scalar(0.,Function(self.domain))
96 d={}
97 for i in self.__volumes:
98 if self.mass_rate.has_key(i): d[i]=-self.mass_rate[i]/self.__volumes[i]
99 self.tag_map.insert(self.density_rate,**d)
100
101 class LinearElasticStressChange(Model):
102 """
103 calculates the stress according to an initial garvity field and a consecutive
104 chenge of density based on linar elastic material model. It is assumed that
105 the lame coefficients don't change over time, the direction of gravity is
106 pointing into the x_{dim} direction.
107
108 @note: add stress changes due to tectonic loading and slip
109
110 @ivar domain: mining region
111 @type domain: L{Domian}
112 @ivar tag_map: a tagmap
113 @type tag_map: L{TagMap}
114 @ivar displacement: displacement field
115 @type displacement: L{Vector} or C{None}
116 @ivar stress: displacement field
117 @type stress: L{Vector} or C{None}
118 @ivar density: initial density distribution
119 @type density: C{float} or L{Scalar}
120 @ivar density_rate: density rate by tag (may be changed of time)
121 @type density_rate: C{dict}
122 @ivar lame_lambda: elasticity coefficient lambda (assumed to be constant over time)
123 @type lame_lambda: C{float} or L{Scalar}
124 @ivar lame_mu: elasticity coefficient mu (assumed to be constant over time)
125 @type lame_mu: C{float} or L{Scalar}
126 @ivar location_of_fixed_displacement: mask of locations and component with zero displacements
127 @type location_of_fixed_displacement: L{Vector} or C{None}
128 """
129 def __init__(self,**kwargs):
130 """
131 """
132 super(LinearElasticStressChange,self).__init__(**kwargs)
133 self.declareParameter(domain=None,
134 tag_map=None,
135 displacement=None, \
136 stress=None, \
137 density=1., \
138 lame_lambda=2., \
139 lame_mu=1., \
140 location_of_fixed_displacement=None, \
141 density_rate=None)
142
143 def doInitialization(self):
144 """
145 initialize time integration
146 """
147 if self.stress == None: self.stress=Tensor(0,Function(self.domain))
148 if self.displacement==None: self.displacement=Scalar(0,Solution(self.domain))
149 self.__pde=LinearPDE(self.domain)
150 self.__pde.setSymmetryOn()
151 A =Tensor4(0.,Function(self.domain))
152 for i in range(self.domain.getDim()):
153 for j in range(self.domain.getDim()):
154 A[i,j,j,i] += self.lame_mu
155 A[i,j,i,j] += self.lame_mu
156 A[i,i,j,j] += self.lame_lambda
157 self.__pde.setValue(A=A,q=self.location_of_fixed_displacement)
158
159 self.__stress=self.stress
160 self.__displacement=self.displacement
161 self.__first=True
162
163 def doInitialStep(self):
164 """
165 """
166 d=self.domain.getDim()
167 self.__pde.setValue(Y=-kronecker(Function(self.domain))[d-1]*9.81*self.density)
168 ddisp=self.__pde.getSolution()
169 deps=symmetric(grad(ddisp))
170 self.stress=self.stress+self.lame_mu*deps+self.lame_lambda*trace(deps)*kronecker(Function(self.domain))
171 self.displacement=self.displacement+ddisp
172 self.__first=False
173
174 def terminateInitialIteration(self):
175 return not self.__first
176
177 def terminateIteration(self):
178 return not self.__first
179
180 def doInitialPostprocessing(self):
181 self.__stress=self.stress
182 self.__displacement=self.displacement
183 self.trace("initial displacement range : %s: %s"%(inf(self.displacement[2]),sup(self.displacement[2])))
184 self.trace("initial stress range : %s: %s"%(inf(self.stress),sup(self.stress)))
185
186 def doStepPreprocessing(self,dt):
187 self.stress=self.__stress
188 self.displacement=self.__displacement
189 self.__first=True
190
191 def doStep(self,dt):
192 """
193 """
194 d=self.domain.getDim()
195 if not self.density_rate == None:
196 self.__pde.setValue(Y=-kronecker(Function(self.domain))[d-1]*9.81*self.density_rate)
197 ddisp=self.__pde.getSolution()
198 deps=symmetric(grad(ddisp))
199 self.stress=self.stress+dt*(self.lame_mu*deps+self.lame_lambda*trace(deps)*kronecker(Function(self.domain)))
200 self.displacement=self.displacement+dt*ddisp
201 self.__first=False
202
203 def doStepPostprocessing(self, dt):
204 self.__stress=self.stress
205 self.__displacement=self.displacement
206 self.trace("displacement range : %s: %s"%(inf(self.displacement[2]),sup(self.displacement[2])))
207 self.trace("stress range : %s: %s"%(inf(self.stress),sup(self.stress)))
208
209 class CoulombFailureStress(ParameterSet):
210 """
211 calculates the Coulomb failure stress an planes of a given orientation.
212
213 @ivar domain: mining region
214 @type domain: L{Domian}
215 @ivar tag_map: a tagmap
216 @type tag_map: L{TagMap}
217 @ivar displacement: displacement field
218 @type displacement: L{Vector} or C{None}
219 @ivar stress: displacement field
220 @type stress: L{Vector} or C{None}
221 @ivar density: initial density distribution
222 @type density: C{float} or L{Scalar}
223 @ivar density_rate: density rate by tag (may be changed of time)
224 @type density_rate: C{dict}
225 @ivar lame_lambda: elasticity coefficient lambda (assumed to be constant over time)
226 @type lame_lambda: C{float} or L{Scalar}
227 @ivar lame_mu: elasticity coefficient mu (assumed to be constant over time)
228 @type lame_mu: C{float} or L{Scalar}
229 @ivar location_of_fixed_displacement: mask of locations and component with zero displacements
230 @type location_of_fixed_displacement: L{Vector} or C{None}
231 """
232 def __init__(self,**kwargs):
233 """
234 """
235 super(CoulombFailureStress,self).__init__(**kwargs)
236 self.declareParameter(stress=numarray.zeros((3,3)),
237 friction_coefficient=0.,
238 normal=numarray.array([1.,0.,0.]))
239 def cfs(self):
240 """
241 returns Coulomb failure stress
242 """
243 sn=matrixmult(self.stress,self.normal)
244 nsn=inner(self.normal,sn)
245 nssn=inner(sn,sn)
246 return (sqrt(nssn-nsn**2)-self.friction_coefficient*nsn)/length(self.normal)
247

  ViewVC Help
Powered by ViewVC 1.1.26