/[escript]/trunk/escript/py_src/rheologies.py
ViewVC logotype

Annotation of /trunk/escript/py_src/rheologies.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3360 - (hide annotations)
Thu Nov 18 00:20:21 2010 UTC (9 years, 11 months ago) by jfenwick
File MIME type: text/x-python
File size: 27203 byte(s)
Fix some epydoc errors
Fixed issue 564
removed minimize.py since it hasn't worked in who knows how long
1 ksteube 1809 ########################################################
2 gross 1639 #
3 jfenwick 2881 # Copyright (c) 2003-2010 by University of Queensland
4 ksteube 1809 # Earth Systems Science Computational Center (ESSCC)
5     # http://www.uq.edu.au/esscc
6 gross 1639 #
7 ksteube 1809 # Primary Business: Queensland, Australia
8     # Licensed under the Open Software License version 3.0
9     # http://www.opensource.org/licenses/osl-3.0.php
10 gross 1639 #
11 ksteube 1809 ########################################################
12 gross 1639
13 jfenwick 2881 __copyright__="""Copyright (c) 2003-2010 by University of Queensland
14 ksteube 1809 Earth Systems Science Computational Center (ESSCC)
15     http://www.uq.edu.au/esscc
16     Primary Business: Queensland, Australia"""
17     __license__="""Licensed under the Open Software License version 3.0
18     http://www.opensource.org/licenses/osl-3.0.php"""
19 jfenwick 2344 __url__="https://launchpad.net/escript-finley"
20 ksteube 1809
21 gross 1639 """
22     Some models for flow
23    
24 jfenwick 2625 :var __author__: name of author
25     :var __copyright__: copyrights
26     :var __license__: licence agreement
27     :var __url__: url entry point on documentation
28     :var __version__: version
29     :var __date__: date of the version
30 gross 1639 """
31    
32     __author__="Lutz Gross, l.gross@uq.edu.au"
33    
34     from escript import *
35     import util
36 gross 2415 from flows import StokesProblemCartesian
37     from pdetools import MaxIterReached
38 gross 1639
39 gross 2300 class PowerLaw(object):
40     """
41     this implements the power law for a composition of a set of materials where the viscosity eta of each material is given by a
42     power law relationship of the form
43    
44 jfenwick 2625 *eta=eta_N*(tau/tau_t)**(1./power-1.)*
45 gross 2300
46     where tau is equivalent stress and eta_N, tau_t and power are given constant. Moreover an elastic component can be considered.
47     Moreover tau meets the Drucker-Prager type yield condition
48    
49 jfenwick 2625 *tau <= tau_Y + friction * pressure*
50 gross 2300
51     where gamma_dot is the equivalent.
52     """
53     def __init__(self, numMaterials=1,verbose=False):
54     """
55     initializes a power law
56    
57 jfenwick 2625 :param numMaterials: number of materials
58     :type numMaterials: ``int``
59     :param verbose: if ``True`` some informations are printed.
60     :type verbose: ``bool``
61 gross 2300 """
62     if numMaterials<1:
63     raise ValueError,"at least one material must be defined."
64     self.__numMaterials=numMaterials
65     self.__eta_N=[None for i in xrange(self.__numMaterials)]
66     self.__tau_t=[1. for i in xrange(self.__numMaterials)]
67     self.__power=[1. for i in xrange(self.__numMaterials)]
68     self.__tau_Y=None
69 gross 2301 self.__friction=None
70 gross 2300 self.__mu=None
71     self.__verbose=verbose
72     self.setEtaTolerance()
73     #===========================================================================
74     def getNumMaterials(self):
75     """
76     returns the numebr of materials
77    
78 jfenwick 2625 :return: number of materials
79     :rtype: ``int``
80 gross 2300 """
81     return self.__numMaterials
82     def validMaterialId(self,id=0):
83     """
84     checks if a given material id is valid
85    
86 jfenwick 2625 :param id: a material id
87     :type id: ``int``
88     :return: ``True`` is the id is valid
89     :rtype: ``bool``
90 gross 2300 """
91     return 0<=id and id<self.getNumMaterials()
92 gross 2440 def setEtaTolerance(self,rtol=1.e-4):
93 gross 2300 """
94     sets the relative tolerance for the effectice viscosity.
95    
96 jfenwick 2625 :param rtol: relative tolerance
97     :type rtol: positive ``float``
98 gross 2300 """
99     if rtol<=0:
100     raise ValueError,"rtol needs to positive."
101     self.__rtol=rtol
102     def getEtaTolerance(self):
103     """
104     returns the relative tolerance for the effectice viscosity.
105    
106 jfenwick 2625 :return: relative tolerance
107     :rtype: positive ``float``
108 gross 2300 """
109     return self.__rtol
110     #===========================================================================
111 gross 2301 def setDruckerPragerLaw(self,tau_Y=None,friction=None):
112 gross 2300 """
113     Sets the parameters for the Drucker-Prager model.
114    
115 jfenwick 2625 :param tau_Y: yield stress
116     :param friction: friction coefficient
117 gross 2300 """
118     self.__tau_Y=tau_Y
119     self.__friction=friction
120     def getFriction(self):
121     """
122     returns the friction coefficient
123    
124 jfenwick 2625 :return: friction coefficient
125 gross 2300 """
126     return self.__friction
127     def getTauY(self):
128     """
129     returns the yield stress
130    
131 jfenwick 2625 :return: the yield stress
132 gross 2300 """
133     return self.__tau_Y
134     #===========================================================================
135     def getElasticShearModulus(self):
136     """
137     returns the elastic shear modulus.
138    
139 jfenwick 2625 :return: elastic shear modulus
140 gross 2300 """
141     return self.__mu
142     def setElasticShearModulus(self,mu=None):
143     """
144     Sets the elastic shear modulus.
145    
146 jfenwick 2625 :param mu: elastic shear modulus
147 gross 2300 """
148     self.__mu=mu
149     #===========================================================================
150     def getPower(self, id=None):
151     """
152     returns the power in the power law
153    
154 jfenwick 2625 :param id: if present, the power for material ``id`` is returned.
155     :type id: ``int``
156     :return: the list of the powers for all matrials is returned. If ``id`` is present only the power for material ``id`` is returned.
157 gross 2300 """
158     if id == None:
159     return self.__power
160     else:
161     if self.validMaterialId(id):
162     return self.__power[id]
163     else:
164     raise ValueError,"Illegal material id %s."%id
165     def getEtaN(self, id=None):
166     """
167     returns the viscosity
168    
169 jfenwick 2625 :param id: if present, the viscosity for material ``id`` is returned.
170     :type id: ``int``
171     :return: the list of the viscosities for all matrials is returned. If ``id`` is present only the viscosity for material ``id`` is returned.
172 gross 2300 """
173     if id == None:
174     return self.__eta_N
175     else:
176     if self.validMaterialId(id):
177     return self.__eta_N[id]
178     else:
179     raise ValueError,"Illegal material id %s."%id
180     def getTauT(self, id=None):
181     """
182     returns the transition stress
183    
184 jfenwick 2625 :param id: if present, the transition stress for material ``id`` is returned.
185     :type id: ``int``
186     :return: the list of the transition stresses for all matrials is returned. If ``id`` is present only the transition stress for material ``id`` is returned.
187 gross 2300 """
188     if id == None:
189     return self.__tau_t
190     else:
191     if self.validMaterialId(id):
192     return self.__tau_t[id]
193     else:
194     raise ValueError,"Illegal material id %s."%id
195    
196     def setPowerLaw(self,eta_N, id=0, tau_t=1, power=1):
197     """
198     Sets the power-law parameters for material id
199    
200 jfenwick 2625 :param id: material id
201     :type id: ``int``
202     :param eta_N: viscosity for tau=tau_t
203     :param tau_t: transition stress
204     :param power: power law coefficient
205 gross 2300 """
206     if self.validMaterialId(id):
207     self.__eta_N[id]=eta_N
208     self.__power[id]=power
209     self.__tau_t[id]=tau_t
210     else:
211     raise ValueError,"Illegal material id %s."%id
212    
213     def setPowerLaws(self,eta_N, tau_t, power):
214     """
215     Sets the parameters of the power-law for all materials.
216    
217 jfenwick 2625 :param eta_N: list of viscosities for tau=tau_t
218     :param tau_t: list of transition stresses
219     :param power: list of power law coefficient
220 gross 2300 """
221     if len(eta_N)!=self.__numMaterials or len(tau_t)!=self.__numMaterials or len(power)!=self.__numMaterials:
222     raise ValueError,"%s materials are expected."%self.__numMaterials
223     for i in xrange(self.__numMaterials):
224     self.setPowerLaw(id=i, eta_N=eta_N[i],tau_t=tau_t[i],power=power[i])
225    
226     #===========================================================================
227 gross 2438 def getEtaEff(self,gamma_dot, eta0=None, pressure=None,dt=None, iter_max=30):
228 gross 2300 """
229     returns the effective viscosity eta_eff such that
230    
231 jfenwick 2625 *tau=eta_eff * gamma_dot*
232 gross 2300
233     by solving a non-linear problem for tau.
234    
235 jfenwick 2625 :param gamma_dot: equivalent strain gamma_dot
236     :param eta0: initial guess for the effective viscosity (e.g from a previous time step). If not present, an initial guess is calculated.
237     :param pressure: pressure used to calculate yield condition
238     :param dt: time step size. only needed if elastic component is considered.
239     :type dt: positive ``float`` if present
240     :param iter_max: maximum number of iteration steps.
241     :type iter_max: ``int``
242     :return: effective viscosity.
243 gross 2300 """
244 gross 2843 if pressure == None:
245     p2 = None
246     else:
247     p2=(abs(pressure)+pressure)/2.
248 gross 2301 SMALL=1./(util.DBLE_MAX/100.)
249 gross 2300 numMaterial=self.getNumMaterials()
250 gross 2438 s=[p-1. for p in self.getPower() ]
251 gross 2300 eta_N=self.getEtaN()
252     tau_t=self.getTauT()
253     mu=self.getElasticShearModulus()
254     fric=self.getFriction()
255     tau_Y=self.getTauY()
256     if eta0==None:
257     theta=0.
258     for i in xrange(numMaterial):
259     inv_eta_i=0**s[i]/eta_N[i]
260     theta=theta+inv_eta_i
261     if util.inf(theta)<=0:
262     raise ValueError,"unable to set positive initial guess for eta_eff. Most likely no power law with power 1 set."
263     eta_eff=1./theta
264     else:
265     if util.inf(eta0)<=0:
266     raise ValueError,"initial guess for eta_eff is not positive."
267     eta_eff=eta0
268    
269 gross 2793 if mu !=None:
270     if dt == None: raise ValueError,"Time stepsize dt must be given."
271     if dt<=0: raise ValueError,"Time step size must be positive."
272 gross 2301 if tau_Y==None and fric==None:
273 gross 2300 eta_max=None
274     else:
275 gross 2843 if fric == None or p2==None:
276 gross 2301 eta_max=tau_Y/(gamma_dot+SMALL*util.whereZero(gamma_dot))
277     else:
278     if tau_Y==None: tau_Y==0
279     if util.inf(fric)<=0:
280     raise ValueError,"if friction present it needs to be positive."
281 gross 2793 eta_max=fric*util.clip(tau_Y/fric+p2,minval=0)/(gamma_dot+SMALL*util.whereZero(gamma_dot))
282 gross 2300 rtol=self.getEtaTolerance()
283     iter =0
284     converged=False
285     tau=eta_eff*gamma_dot
286 gross 2415 if self.__verbose: print "PowerLaw: Start calculation of eta_eff (tolerance = %s)\nPowerLaw: initial max eta_eff = %s, tau = %s."%(rtol,util.Lsup(eta_eff),util.Lsup(tau))
287 gross 2300 while not converged:
288     if iter>max(iter_max,1):
289 gross 2301 raise RuntimeError,"tolerance not reached after %s steps."%max(iter_max,1)
290 gross 2300 #===========================================
291     theta=0. # =1/eta
292     omega=0. # = tau*theta'= eta'*tau/eta**2
293     if mu !=None: theta=1./(dt*mu)
294     for i in xrange(numMaterial):
295     inv_eta_i=(tau/tau_t[i])**s[i]/eta_N[i]
296     theta=theta+inv_eta_i
297     omega=omega+s[i]*inv_eta_i
298     #===========================================
299     eta_eff, eta_eff_old=util.clip(eta_eff*(theta+omega)/(eta_eff*theta**2+omega),maxval=eta_max), eta_eff
300     tau=eta_eff*gamma_dot
301     d=util.Lsup(eta_eff-eta_eff_old)
302     l=util.Lsup(eta_eff)
303     iter+=1
304 gross 2415 if self.__verbose: print "PowerLaw: step %s: correction = %s, max eta_eff = %s, max tau= %s"%(iter, d, l,util.Lsup(tau))
305 gross 2300 converged= d<= rtol* l
306 gross 2415 if self.__verbose: print "PowerLaw: Start calculation of eta_eff finalized after %s steps."%iter
307 gross 2300 return eta_eff
308    
309     #====================================================================================================================================
310 gross 2415 class Rheology(object):
311 gross 1639 """
312 gross 2415 General framework to implement a rheology
313 gross 1639 """
314 gross 2415 def __init__(self, domain, stress=None, v=None, p=None, t=0, verbose=True):
315 gross 1639 """
316 gross 2415 Initializes the rheology
317 gross 1639
318 jfenwick 2625 :param domain: problem domain
319     :type domain: `Domain`
320     :param stress: initial (deviatoric) stress
321     :type stress: a tensor value/field of order 2
322     :param v: initial velocity field
323     :type v: a vector value/field
324     :param p: initial pressure
325     :type p: a scalar value/field
326     :param t: initial time
327     :type t: ``float``
328 gross 1639 """
329     self.__domain=domain
330     self.__t=t
331 gross 2415 self.__verbose=verbose
332 gross 1639 #=======================
333 gross 2415 #
334 gross 1639 # state variables:
335     #
336 gross 2415 if stress == None: stress=Tensor(0.,Function(self.__domain))
337     if v == None: v=Vector(0.,Solution(self.__domain))
338     if p == None: p=Vector(0.,ReducedSolution(self.__domain))
339 gross 2620 self.setStatus(t, v, p, stress)
340     self.setExternals(F=Data(), f=Data(), fixed_v_mask=Data(), v_boundary=Data(), restoration_factor=0)
341 gross 2415
342     def getDomain(self):
343     """
344     returns the domain.
345 caltinay 2169
346 jfenwick 2625 :return: the domain
347     :rtype: `Domain`
348 gross 2415 """
349     return self.__domain
350 gross 1639
351 gross 2415 def getTime(self):
352     """
353     Returns current time.
354 caltinay 2169
355 jfenwick 2625 :return: current time
356     :rtype: ``float``
357 gross 1639 """
358 gross 2415 return self.__t
359    
360 gross 2620 def setExternals(self, F=None, f=None, fixed_v_mask=None, v_boundary=None, restoration_factor=None):
361 caltinay 2169 """
362 gross 2415 sets external forces and velocity constraints
363 gross 2100
364 jfenwick 2625 :param F: external force
365     :type F: vector value/field
366     :param f: surface force
367     :type f: vector value/field on boundary
368     :param fixed_v_mask: location of constraints maked by positive values
369     :type fixed_v_mask: vector value/field
370     :param v_boundary: value of velocity at location of constraints
371     :type v_boundary: vector value/field
372     :param restoration_factor: factor for normal restoration force
373     :type restoration_factor: scalar values/field
374     :note: Only changing parameters need to be specified.
375 gross 1639 """
376 gross 2415 if F != None: self.__F=F
377     if f != None: self.__f=f
378     if fixed_v_mask != None: self.__fixed_v_mask=fixed_v_mask
379     if v_boundary != None: self.__v_boundary=v_boundary
380 gross 2620 if restoration_factor!=None: self.__restoration_factor=restoration_factor
381 gross 2415
382     def getForce(self):
383     """
384     Returns the external force
385 caltinay 2169
386 jfenwick 2625 :return: external force
387     :rtype: `Data`
388 gross 1639 """
389 gross 2415 return self.__F
390 gross 2100
391 gross 2415 def getSurfaceForce(self):
392 gross 1639 """
393 gross 2415 Returns the surface force
394    
395 jfenwick 2625 :return: surface force
396     :rtype: `Data`
397 gross 1639 """
398 gross 2415 return self.__f
399 gross 2100
400 gross 2415 def getVelocityConstraint(self):
401     """
402     Returns the constraint for the velocity as a pair of the
403     mask of the location of the constraint and the values.
404    
405 jfenwick 2625 :return: the locations of fixed velocity and value of velocities at these locations
406     :rtype: ``tuple`` of `Data` s
407 gross 2415 """
408     return self.__fixed_v_mask, self.__v_boundary
409 gross 2793
410 gross 2620 def getRestorationFactor(self):
411     """
412     Returns the restoring force factor
413 gross 2415
414 gross 2620 @return: restoring force factor
415     @rtype: C{float} or L{Data}
416     """
417     return self.__restoration_factor
418    
419    
420 gross 2415 def checkVerbose(self):
421     """
422     Returns True if verbose is switched on
423    
424 jfenwick 2625 :return: value of verbosity flag
425     :rtype: ``bool``
426 gross 2415 """
427     return self.__verbose
428    
429     def setTime(self,t=0.):
430     """
431     Updates current time.
432    
433 jfenwick 2625 :param t: new time mark
434     :type t: ``float``
435 gross 2415 """
436     self.__t=t
437     #=======================================================================================
438     def getStress(self):
439     """
440     Returns current stress.
441    
442 jfenwick 2625 :return: current stress
443     :rtype: `Data` of rank 2
444 gross 2415 """
445     s=self.getDeviatoricStress()
446     p=self.getPressure()
447     k=util.kronecker(self.getDomain())
448     return s-p*(k/trace(k))
449    
450     def getDeviatoricStress(self):
451     """
452     Returns current deviatoric stress.
453    
454 jfenwick 2625 :return: current deviatoric stress
455     :rtype: `Data` of rank 2
456 gross 2415 """
457     return self.__stress
458    
459     def setDeviatoricStress(self, stress):
460     """
461     Sets the current deviatoric stress
462    
463 jfenwick 2625 :param stress: new deviatoric stress
464     :type stress: `Data` of rank 2
465 gross 2415 """
466     dom=self.getDomain()
467     s=util.interpolate(stress,Function(dom))
468     self.__stress=util.deviatoric(s)
469    
470     def getPressure(self):
471     """
472     Returns current pressure.
473    
474 jfenwick 2625 :return: current stress
475     :rtype: scalar `Data`
476 gross 2415 """
477     return self.__p
478    
479     def setPressure(self, p):
480     """
481     Sets current pressure.
482 jfenwick 2625 :param p: new deviatoric stress
483     :type p: scalar `Data`
484 gross 2415 """
485     self.__p=util.interpolate(p,ReducedSolution(self.getDomain()))
486    
487     def getVelocity(self):
488     """
489     Returns current velocity.
490    
491 jfenwick 2625 :return: current velocity
492     :rtype: vector `Data`
493 gross 2415 """
494     return self.__v
495    
496     def setVelocity(self, v):
497     """
498     Sets current velocity.
499    
500 jfenwick 2625 :param v: new current velocity
501     :type v: vector `Data`
502 gross 2415 """
503     self.__v=util.interpolate(v,Solution(self.getDomain()))
504 gross 2620 def setStatus(self,t, v, p, stress):
505     """
506     Resets the current status given by pressure p and velocity v.
507    
508     @param t: new time mark
509     @type t: C{float}
510     @param v: new current velocity
511     @type v: vector L{Data}
512     @param p: new deviatoric stress
513     @type p: scalar L{Data}
514     @param stress: new deviatoric stress
515     @type stress: L{Data} of rank 2
516     """
517     self.setDeviatoricStress(stress)
518     self.setVelocity(v)
519     self.setPressure(p)
520     self.setDeviatoricStrain()
521 gross 2850 self.setGammaDot()
522 gross 2620 self.setTime(t)
523 gross 2415
524     def setDeviatoricStrain(self, D=None):
525     """
526     set deviatoric strain
527    
528 gross 2850 :param D: new deviatoric strain. If ``D`` is not present the current velocity is used.
529 jfenwick 2625 :type D: `Data` of rank 2
530 gross 2415 """
531 gross 2850 if D==None:
532     self.__D=self.getDeviatoricStrain(self.getVelocity())
533     else:
534     self.__D=util.deviatoric(util.interpolate(D,Function(self.getDomain())))
535 gross 2415
536 gross 2850 def getDeviatoricStrain(self, v=None):
537 gross 2415 """
538 gross 2850 Returns deviatoric strain of current velocity or if ``v`` is present the
539     deviatoric strain of velocity ``v``:
540 gross 2415
541 gross 2850 :param v: a velocity field
542     :type v: `Data` of rank 1
543     :return: deviatoric strain of the current velocity field or if ``v`` is present the deviatoric strain of velocity ``v``
544 jfenwick 2625 :rtype: `Data` of rank 2
545 gross 2415 """
546 gross 2850 if v==None:
547     return self.__D
548     else:
549     return util.deviatoric(util.symmetric(util.grad(v)))
550 gross 2415
551     def getTau(self):
552     """
553     Returns current second invariant of deviatoric stress
554    
555 jfenwick 2625 :return: second invariant of deviatoric stress
556     :rtype: scalar `Data`
557 gross 2415 """
558     s=self.getDeviatoricStress()
559     return util.sqrt(0.5)*util.length(s)
560    
561 gross 2850 def setGammaDot(self, gammadot=None):
562 gross 2415 """
563 gross 2850 set the second invariant of deviatoric strain rate. If ``gammadot`` is not present zero is used.
564 gross 2415
565 gross 2850 :param gammadot: second invariant of deviatoric strain rate.
566     :type gammadot: `Data` of rank 1
567     """
568     if gammadot == None:
569     self.__gammadot = Scalar(0.,Function(self.getDomain()))
570     else:
571     self.__gammadot=gammadot
572    
573     def getGammaDot(self, D=None):
574     """
575     Returns current second invariant of deviatoric strain rate or if ``D`` is present the second invariant of ``D``.
576    
577     :param D: deviatoric strain rate tensor
578     :type D: `Data` of rank 0
579 jfenwick 2625 :return: second invariant of deviatoric strain
580     :rtype: scalar `Data`
581 gross 2415 """
582 gross 2850 if D==None:
583     return self.__gammadot
584     else:
585     return util.sqrt(2.)*util.length(D)
586    
587 gross 2415
588     #====================================================================================================================================
589 gross 1639
590 gross 2793 class IncompressibleIsotropicFlowCartesian(PowerLaw,Rheology, StokesProblemCartesian):
591     """
592     This class implements the rheology of an isotropic Kelvin material.
593 caltinay 2169
594 gross 2793 Typical usage::
595 caltinay 2169
596 gross 2793 sp = IncompressibleIsotropicFlowCartesian(domain, stress=0, v=0)
597 gross 2415 sp.initialize(...)
598     v,p = sp.solve()
599 gross 2100
600 gross 2793 :note: This model has been used in the self-consistent plate-mantle model
601 gross 2415 proposed in U{Hans-Bernd Muhlhaus<emailto:h.muhlhaus@uq.edu.au>}
602     and U{Klaus Regenauer-Lieb<mailto:klaus.regenauer-lieb@csiro.au>}:
603     I{Towards a self-consistent plate mantle model that includes elasticity: simple benchmarks and application to basic modes of convection},
604 gross 2438 see U{doi: 10.1111/j.1365-246X.2005.02742.x<http://www3.interscience.wiley.com/journal/118661486/abstract>}
605 gross 2793 """
606     def __init__(self, domain, stress=0, v=0, p=0, t=0, numMaterials=1, verbose=True):
607 gross 2415 """
608     Initializes the model.
609    
610 jfenwick 2625 :param domain: problem domain
611     :type domain: `Domain`
612     :param stress: initial (deviatoric) stress
613     :type stress: a tensor value/field of order 2
614     :param v: initial velocity field
615     :type v: a vector value/field
616     :param p: initial pressure
617     :type p: a scalar value/field
618     :param t: initial time
619     :type t: ``float``
620     :param numMaterials: number of materials
621     :type numMaterials: ``int``
622     :param verbose: if ``True`` some informations are printed.
623     :type verbose: ``bool``
624 gross 2415 """
625 gross 2793 PowerLaw. __init__(self, numMaterials,verbose=verbose)
626     Rheology. __init__(self, domain, stress, v, p, t,verbose=verbose)
627     StokesProblemCartesian.__init__(self,domain,verbose=verbose)
628 gross 2415 self.__eta_eff=None
629    
630 gross 2793 def getCurrentEtaEff(self):
631 gross 1639 """
632 gross 2793 returns the effective viscosity used in the last iteration step of the last time step.
633     """
634     return self.__eta_eff
635    
636    
637     def updateStokesEquation(self, v, p):
638     """
639     updates the underlying Stokes equation to consider dependencies from ``v`` and ``p``
640     """
641     dt=self.__dt
642     mu=self.getElasticShearModulus()
643     F=self.getForce()
644     f=self.getSurfaceForce()
645     mask_v,v_b=self.getVelocityConstraint()
646     s_last=self.getDeviatoricStress()
647     #
648     # calculate eta_eff if we don't have one or elasticity is present.
649     #
650     if mu==None:
651 gross 2850 gamma=self.getGammaDot(self.getDeviatoricStrain(v))
652 gross 2793 else:
653 gross 2850 gamma=self.getGammaDot(self.getDeviatoricStrain(v)+s_last/(2*dt*mu))
654 gross 2793
655     self.__eta_eff_save=self.getEtaEff(gamma, pressure=p,dt=dt, eta0=self.__eta_eff_save, iter_max=self.__eta_iter_max)
656    
657     if self.checkVerbose(): print "IncompressibleIsotropicFlowCartesian: eta_eff has been updated."
658    
659     if mu==None:
660     stress0=Data()
661     else:
662     stress0=-(self.__eta_eff_save/(dt*mu))*s_last
663    
664     self.setStokesEquation(eta=self.__eta_eff_save,stress=stress0)
665    
666    
667     def initialize(self, F=None, f=None, fixed_v_mask=None, v_boundary=None, restoration_factor=None):
668     """
669     sets external forces and velocity constraints
670    
671     :param F: external force
672     :type F: vector value/field
673     :param f: surface force
674     :type f: vector value/field on boundary
675     :param fixed_v_mask: location of constraints maked by positive values
676     :type fixed_v_mask: vector value/field
677     :param v_boundary: value of velocity at location of constraints
678     :type v_boundary: vector value/field
679     :param restoration_factor: factor for normal restoration force
680     :type restoration_factor: scalar values/field
681     :note: Only changing parameters need to be specified.
682     """
683     self.setExternals(F, f, fixed_v_mask, v_boundary, restoration_factor)
684    
685     def update(self, dt, iter_max=10, eta_iter_max=20, verbose=False, usePCG=True, max_correction_steps=50):
686     """
687 caltinay 2169 Updates stress, velocity and pressure for time increment dt.
688 gross 1639
689 jfenwick 2625 :param dt: time increment
690 gross 2793 :param iter_max: maximum number of iteration steps in the incompressible solver
691     :param eta_iter_max: maximum number of iteration steps in the incompressible solver
692 jfenwick 2625 :param verbose: prints some infos in the incompressible solver
693 gross 1639 """
694 gross 2793 mu=self.getElasticShearModulus()
695     if mu != None:
696     if not dt > 0.:
697     raise ValueError,"dt must be positive."
698     else:
699     dt=max(0,dt)
700     self.__dt=dt
701     self.__eta_iter_max=max(eta_iter_max,1)
702     v_last=self.getVelocity()
703 gross 2415 s_last=self.getDeviatoricStress()
704     mask_v,v_b=self.getVelocityConstraint()
705 gross 2793 p_last=self.getPressure()
706     self.__eta_eff_save=self.getCurrentEtaEff()
707    
708     self.setStokesEquation(f=self.getForce(),fixed_u_mask=mask_v,surface_stress=self.getSurfaceForce(), restoration_factor=self.getRestorationFactor())
709    
710     if self.checkVerbose(): print "IncompressibleIsotropicFlowCartesian: start iteration for t = %s."%(self.getTime()+dt,)
711     #
712     # get a new velcocity and pressure:
713 gross 2415 #
714 gross 2850 if mask_v.isEmpty():
715 gross 2793 v0=v_last
716 gross 2415 else:
717 gross 2850 if v_b.isEmpty():
718     v0=v_last*(1.-mask_v)
719     else:
720     v0=v_b*mask_v+v_last*(1.-mask_v)
721 gross 2793
722     v,p=self._solve(v0,p_last,verbose=self.checkVerbose(),max_iter=iter_max,usePCG=usePCG, max_correction_steps=max_correction_steps)
723 gross 2415 #
724     # finally we can update the return values:
725     #
726     self.setPressure(p)
727     self.setVelocity(v)
728 gross 2850 self.setDeviatoricStrain(self.getDeviatoricStrain(v))
729     if mu==None:
730     D=self.getDeviatoricStrain(v)
731     else:
732     D=self.getDeviatoricStrain(v)+s_last/(2*dt*mu)
733     gamma=self.getGammaDot(D)
734     self.setGammaDot(gamma)
735 gross 2793 self.__eta_eff = self.getEtaEff(self.getGammaDot(), pressure=p,dt=dt, eta0=self.__eta_eff_save, iter_max=self.__eta_iter_max)
736 gross 2850 self.setDeviatoricStress(2.*self.__eta_eff*D)
737 gross 2415 self.setTime(self.getTime()+dt)
738 gross 2793 if self.checkVerbose(): print "IncompressibleIsotropicFlowCartesian: iteration on time step %s completed."%(self.getTime(),)
739 gross 2415 return self.getVelocity(), self.getPressure()
740 gross 1639
741 gross 2440

  ViewVC Help
Powered by ViewVC 1.1.26