/[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 998 - (hide annotations)
Mon Feb 26 07:09:37 2007 UTC (14 years, 9 months ago) by gross
File MIME type: text/x-python
File size: 7858 byte(s)
basic mine scenario
1 gross 987 """
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 gross 997 from esys.escript import *
23     from esys.escript.linearPDEs import LinearPDE
24 gross 987 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 gross 997 class DensityChange(Model):
60     """
61     translates local mass extraction into density changes.
62     "local" refers a taged region.
63 gross 987
64 gross 997 @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 gross 998 if self.mass_rate.has_key(i): d[i]=-self.mass_rate[i]/self.__volumes[i]
98 gross 997 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 gross 998 d=self.domain.getDim()
166 gross 997 self.__pde.setValue(Y=-kronecker(Function(self.domain))[d-1]*9.81*self.density)
167     ddisp=self.__pde.getSolution()
168     deps=symmetric(grad(ddisp))
169     self.stress=self.stress+self.lame_mu*deps+self.lame_lambda*trace(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 gross 998 def doInitialPostprocessing(self):
180 gross 997 self.__stress=self.stress
181     self.__displacement=self.displacement
182 gross 998 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 gross 997
185     def doStepPreprocessing(self,dt):
186     self.stress=self.__stress
187     self.displacement=self.__displacement
188     self.__first=True
189    
190 gross 998 def doStep(self,dt):
191 gross 997 """
192     """
193 gross 998 d=self.domain.getDim()
194 gross 997 if not self.density_rate == None:
195     self.__pde.setValue(Y=-kronecker(Function(self.domain))[d-1]*9.81*self.density_rate)
196     ddisp=self.__pde.getSolution()
197     deps=symmetric(grad(ddisp))
198     self.stress=self.stress+dt*(self.lame_mu*deps+self.lame_lambda*trace(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 gross 998 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)))

  ViewVC Help
Powered by ViewVC 1.1.26