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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2620 - (show annotations)
Thu Aug 20 06:24:00 2009 UTC (9 years, 8 months ago) by gross
File MIME type: text/x-python
File size: 30162 byte(s)
some small additions to pycad to make life a bit easier.
1 ########################################################
2 #
3 # Copyright (c) 2003-2009 by University of Queensland
4 # Earth Systems Science Computational Center (ESSCC)
5 # http://www.uq.edu.au/esscc
6 #
7 # 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 #
11 ########################################################
12
13 __copyright__="""Copyright (c) 2003-2009 by University of Queensland
14 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 __url__="https://launchpad.net/escript-finley"
20
21 """
22 Some models for flow
23
24 @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 """
31
32 __author__="Lutz Gross, l.gross@uq.edu.au"
33
34 from escript import *
35 import util
36 from flows import StokesProblemCartesian
37 from pdetools import MaxIterReached
38
39 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 M{eta=eta_N*(tau/tau_t)**(1./power-1.)}
45
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 M{tau <= tau_Y + friction * pressure}
50
51 where gamma_dot is the equivalent.
52 """
53 def __init__(self, numMaterials=1,verbose=False):
54 """
55 initializes a power law
56
57 @param numMaterials: number of materials
58 @type numMaterials: C{int}
59 @param verbose: if C{True} some informations are printed.
60 @type verbose: C{bool}
61 """
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 self.__friction=None
70 self.__mu=None
71 self.__verbose=verbose
72 self.setEtaTolerance()
73 #===========================================================================
74 def getNumMaterials(self):
75 """
76 returns the numebr of materials
77
78 @return: number of materials
79 @rtype: C{int}
80 """
81 return self.__numMaterials
82 def validMaterialId(self,id=0):
83 """
84 checks if a given material id is valid
85
86 @param id: a material id
87 @type id: C{int}
88 @return: C{True} is the id is valid
89 @rtype: C{bool}
90 """
91 return 0<=id and id<self.getNumMaterials()
92 def setEtaTolerance(self,rtol=1.e-4):
93 """
94 sets the relative tolerance for the effectice viscosity.
95
96 @param rtol: relative tolerance
97 @type rtol: positive C{float}
98 """
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 @return: relative tolerance
107 @rtype rtol: positive C{float}
108 """
109 return self.__rtol
110 #===========================================================================
111 def setDruckerPragerLaw(self,tau_Y=None,friction=None):
112 """
113 Sets the parameters for the Drucker-Prager model.
114
115 @param tau_Y: yield stress
116 @param friction: friction coefficient
117 """
118 self.__tau_Y=tau_Y
119 self.__friction=friction
120 def getFriction(self):
121 """
122 returns the friction coefficient
123
124 @return: friction coefficient
125 """
126 return self.__friction
127 def getTauY(self):
128 """
129 returns the yield stress
130
131 @return: the yield stress
132 """
133 return self.__tau_Y
134 #===========================================================================
135 def getElasticShearModulus(self):
136 """
137 returns the elastic shear modulus.
138
139 @return: elastic shear modulus
140 """
141 return self.__mu
142 def setElasticShearModulus(self,mu=None):
143 """
144 Sets the elastic shear modulus.
145
146 @param mu: elastic shear modulus
147 """
148 self.__mu=mu
149 #===========================================================================
150 def getPower(self, id=None):
151 """
152 returns the power in the power law
153
154 @param id: if present, the power for material C{id} is returned.
155 @type id: C{int}
156 @return: the list of the powers for all matrials is returned. If C{id} is present only the power for material C{id} is returned.
157 """
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 @param id: if present, the viscosity for material C{id} is returned.
170 @type id: C{int}
171 @return: the list of the viscosities for all matrials is returned. If C{id} is present only the viscosity for material C{id} is returned.
172 """
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 @param id: if present, the transition stress for material C{id} is returned.
185 @type id: C{int}
186 @return: the list of the transition stresses for all matrials is returned. If C{id} is present only the transition stress for material C{id} is returned.
187 """
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 @param id: material id
201 @type id: C{int}
202 @param eta_N: viscosity for tau=tau_t
203 @param tau_t: transition stress
204 @param power: power law coefficient
205 """
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 @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 """
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 def getEtaEff(self,gamma_dot, eta0=None, pressure=None,dt=None, iter_max=30):
228 """
229 returns the effective viscosity eta_eff such that
230
231 M{tau=eta_eff * gamma_dot}
232
233 by solving a non-linear problem for tau.
234
235 @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 C{float} if present
240 @param iter_max: maximum number of iteration steps.
241 @type iter_max: C{int}
242 @return: effective viscosity.
243 """
244 SMALL=1./(util.DBLE_MAX/100.)
245 numMaterial=self.getNumMaterials()
246 s=[p-1. for p in self.getPower() ]
247 eta_N=self.getEtaN()
248 tau_t=self.getTauT()
249 mu=self.getElasticShearModulus()
250 fric=self.getFriction()
251 tau_Y=self.getTauY()
252 if eta0==None:
253 theta=0.
254 for i in xrange(numMaterial):
255 inv_eta_i=0**s[i]/eta_N[i]
256 theta=theta+inv_eta_i
257 if util.inf(theta)<=0:
258 raise ValueError,"unable to set positive initial guess for eta_eff. Most likely no power law with power 1 set."
259 eta_eff=1./theta
260 else:
261 if util.inf(eta0)<=0:
262 raise ValueError,"initial guess for eta_eff is not positive."
263 eta_eff=eta0
264
265 if mu !=None and dt == None:
266 raise ValueError,"Time stepsize dt must be given."
267 if dt !=None:
268 if dt<=0: raise ValueError,"time step size must be positive."
269 if tau_Y==None and fric==None:
270 eta_max=None
271 else:
272 if fric == None:
273 eta_max=tau_Y/(gamma_dot+SMALL*util.whereZero(gamma_dot))
274 else:
275 if tau_Y==None: tau_Y==0
276 if util.inf(fric)<=0:
277 raise ValueError,"if friction present it needs to be positive."
278 eta_max=fric*util.clip(tau_Y/fric+pressure,minval=0)/(gamma_dot+SMALL*util.whereZero(gamma_dot))
279 rtol=self.getEtaTolerance()
280 iter =0
281 converged=False
282 tau=eta_eff*gamma_dot
283 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))
284 while not converged:
285 if iter>max(iter_max,1):
286 raise RuntimeError,"tolerance not reached after %s steps."%max(iter_max,1)
287 #===========================================
288 theta=0. # =1/eta
289 omega=0. # = tau*theta'= eta'*tau/eta**2
290 if mu !=None: theta=1./(dt*mu)
291 for i in xrange(numMaterial):
292 inv_eta_i=(tau/tau_t[i])**s[i]/eta_N[i]
293 theta=theta+inv_eta_i
294 omega=omega+s[i]*inv_eta_i
295 #===========================================
296 eta_eff, eta_eff_old=util.clip(eta_eff*(theta+omega)/(eta_eff*theta**2+omega),maxval=eta_max), eta_eff
297 tau=eta_eff*gamma_dot
298 d=util.Lsup(eta_eff-eta_eff_old)
299 l=util.Lsup(eta_eff)
300 iter+=1
301 if self.__verbose: print "PowerLaw: step %s: correction = %s, max eta_eff = %s, max tau= %s"%(iter, d, l,util.Lsup(tau))
302 converged= d<= rtol* l
303 if self.__verbose: print "PowerLaw: Start calculation of eta_eff finalized after %s steps."%iter
304 return eta_eff
305
306 #====================================================================================================================================
307 class Rheology(object):
308 """
309 General framework to implement a rheology
310 """
311 def __init__(self, domain, stress=None, v=None, p=None, t=0, verbose=True):
312 """
313 Initializes the rheology
314
315 @param domain: problem domain
316 @type domain: L{Domain}
317 @param stress: initial (deviatoric) stress
318 @type stress: a tensor value/field of order 2
319 @param v: initial velocity field
320 @type stress: a vector value/field
321 @param p: initial pressure
322 @type p: a scalar value/field
323 @param t: initial time
324 @type t: C{float}
325 """
326 self.__domain=domain
327 self.__t=t
328 self.__verbose=verbose
329 #=======================
330 #
331 # state variables:
332 #
333 if stress == None: stress=Tensor(0.,Function(self.__domain))
334 if v == None: v=Vector(0.,Solution(self.__domain))
335 if p == None: p=Vector(0.,ReducedSolution(self.__domain))
336 self.setStatus(t, v, p, stress)
337 self.setExternals(F=Data(), f=Data(), fixed_v_mask=Data(), v_boundary=Data(), restoration_factor=0)
338
339 def getDomain(self):
340 """
341 returns the domain.
342
343 @return: the domain
344 @rtype: L{Domain}
345 """
346 return self.__domain
347
348 def getTime(self):
349 """
350 Returns current time.
351
352 @return: current time
353 @rtype: C{float}
354 """
355 return self.__t
356
357 def setExternals(self, F=None, f=None, fixed_v_mask=None, v_boundary=None, restoration_factor=None):
358 """
359 sets external forces and velocity constraints
360
361 @param F: external force
362 @type F: vector value/field
363 @param f: surface force
364 @type f: vector value/field on boundary
365 @param fixed_v_mask: location of constraints maked by positive values
366 @type fixed_v_mask: vector value/field
367 @param v_boundary: value of velocity at location of constraints
368 @type v_boundary: vector value/field
369 @param restoration_factor: factor for normal restoration force
370 @type restoration_factor: scalar values/field
371 @note: Only changing parameters need to be specified.
372 """
373 if F != None: self.__F=F
374 if f != None: self.__f=f
375 if fixed_v_mask != None: self.__fixed_v_mask=fixed_v_mask
376 if v_boundary != None: self.__v_boundary=v_boundary
377 if restoration_factor!=None: self.__restoration_factor=restoration_factor
378
379 def getForce(self):
380 """
381 Returns the external force
382
383 @return: external force
384 @rtype: L{Data}
385 """
386 return self.__F
387
388 def getSurfaceForce(self):
389 """
390 Returns the surface force
391
392 @return: surface force
393 @rtype: L{Data}
394 """
395 return self.__f
396
397 def getVelocityConstraint(self):
398 """
399 Returns the constraint for the velocity as a pair of the
400 mask of the location of the constraint and the values.
401
402 @return: the locations of fixed velocity and value of velocities at these locations
403 @rtype: C{tuple} of L{Data}s
404 """
405 return self.__fixed_v_mask, self.__v_boundary
406 def getRestorationFactor(self):
407 """
408 Returns the restoring force factor
409
410 @return: restoring force factor
411 @rtype: C{float} or L{Data}
412 """
413 return self.__restoration_factor
414
415
416 def checkVerbose(self):
417 """
418 Returns True if verbose is switched on
419
420 @return: value of verbosity flag
421 @rtype: C{bool}
422 """
423 return self.__verbose
424
425 def setTime(self,t=0.):
426 """
427 Updates current time.
428
429 @param t: new time mark
430 @type t: C{float}
431 """
432 self.__t=t
433 #=======================================================================================
434 def getStress(self):
435 """
436 Returns current stress.
437
438 @return: current stress
439 @rtype: L{Data} of rank 2
440 """
441 s=self.getDeviatoricStress()
442 p=self.getPressure()
443 k=util.kronecker(self.getDomain())
444 return s-p*(k/trace(k))
445
446 def getDeviatoricStress(self):
447 """
448 Returns current deviatoric stress.
449
450 @return: current deviatoric stress
451 @rtype: L{Data} of rank 2
452 """
453 return self.__stress
454
455 def setDeviatoricStress(self, stress):
456 """
457 Sets the current deviatoric stress
458
459 @param stress: new deviatoric stress
460 @type stress: L{Data} of rank 2
461 """
462 dom=self.getDomain()
463 s=util.interpolate(stress,Function(dom))
464 self.__stress=util.deviatoric(s)
465
466 def getPressure(self):
467 """
468 Returns current pressure.
469
470 @return: current stress
471 @rtype: scalar L{Data}
472 """
473 return self.__p
474
475 def setPressure(self, p):
476 """
477 Sets current pressure.
478 @param p: new deviatoric stress
479 @type p: scalar L{Data}
480 """
481 self.__p=util.interpolate(p,ReducedSolution(self.getDomain()))
482
483 def getVelocity(self):
484 """
485 Returns current velocity.
486
487 @return: current velocity
488 @rtype: vector L{Data}
489 """
490 return self.__v
491
492 def setVelocity(self, v):
493 """
494 Sets current velocity.
495
496 @param v: new current velocity
497 @type v: vector L{Data}
498 """
499 self.__v=util.interpolate(v,Solution(self.getDomain()))
500 def setStatus(self,t, v, p, stress):
501 """
502 Resets the current status given by pressure p and velocity v.
503
504 @param t: new time mark
505 @type t: C{float}
506 @param v: new current velocity
507 @type v: vector L{Data}
508 @param p: new deviatoric stress
509 @type p: scalar L{Data}
510 @param stress: new deviatoric stress
511 @type stress: L{Data} of rank 2
512 """
513 self.setDeviatoricStress(stress)
514 self.setVelocity(v)
515 self.setPressure(p)
516 self.setDeviatoricStrain()
517 self.setTime(t)
518
519 def setDeviatoricStrain(self, D=None):
520 """
521 set deviatoric strain
522
523 @param D: new deviatoric strain. If D is not present the current velocity is used.
524 @type D: L{Data} of rank 2
525 """
526 if D==None: D=util.deviatoric(util.symmetric(util.grad(2.*self.getVelocity())))
527 self.__D=util.deviatoric(util.interpolate(D,Function(self.getDomain())))
528
529 def getDeviatoricStrain(self):
530 """
531 Returns deviatoric strain of current velocity.
532
533 @return: deviatoric strain
534 @rtype: L{Data} of rank 2
535 """
536 return self.__D
537
538 def getTau(self):
539 """
540 Returns current second invariant of deviatoric stress
541
542 @return: second invariant of deviatoric stress
543 @rtype: scalar L{Data}
544 """
545 s=self.getDeviatoricStress()
546 return util.sqrt(0.5)*util.length(s)
547
548 def getGammaDot(self):
549 """
550 Returns current second invariant of deviatoric strain
551
552 @return: second invariant of deviatoric strain
553 @rtype: scalar L{Data}
554 """
555 s=self.getDeviatoricStrain()
556 return util.sqrt(2)*util.length(s)
557
558 def setTolerance(self,tol=1.e-4):
559 """
560 Sets the tolerance used to terminate the iteration on a time step.
561 See the implementation of the rheology for details.
562
563 @param tol: relative tolerance to terminate iteration on time step.
564 @type tol: positive C{float}
565 """
566 if tol<=0.:
567 raise ValueError,"tolerance must be non-negative."
568 self.__tol=tol
569
570 def getTolerance(self):
571 """
572 Returns the set tolerance for terminate the iteration on a time step.
573
574 @rtype: positive C{float}
575 """
576 return self.__tol
577
578 #=======================================================================
579 def setFlowTolerance(self, tol=1.e-4):
580 """
581 Sets the relative tolerance for the flow solver
582
583 @param tol: desired relative tolerance for the flow solver
584 @type tol: positive C{float}
585 @note: Typically this method is overwritten by a subclass.
586 """
587 pass
588 def getFlowTolerance(self):
589 """
590 Returns the relative tolerance for the flow solver
591
592 @return: tolerance of the flow solver
593 @rtype: C{float}
594 @note: Typically this method is overwritten by a subclass.
595 """
596 pass
597 #====================================================================================================================================
598
599 class IncompressibleIsotropicFlowCartesian(PowerLaw,Rheology):
600 """
601 This class implements the rheology of an isotropic Kelvin material.
602
603 Typical usage::
604
605 sp = IncompressibleIsotropicFlow(domain, stress=0, v=0)
606 sp.setTolerance()
607 sp.initialize(...)
608 v,p = sp.solve()
609
610 @note: This model has been used in the self-consistent plate-mantle model
611 proposed in U{Hans-Bernd Muhlhaus<emailto:h.muhlhaus@uq.edu.au>}
612 and U{Klaus Regenauer-Lieb<mailto:klaus.regenauer-lieb@csiro.au>}:
613 I{Towards a self-consistent plate mantle model that includes elasticity: simple benchmarks and application to basic modes of convection},
614 see U{doi: 10.1111/j.1365-246X.2005.02742.x<http://www3.interscience.wiley.com/journal/118661486/abstract>}
615 """
616 def __init__(self, domain, stress=0, v=0, p=0, t=0, numMaterials=1, verbose=True, adaptSubTolerance=True):
617 """
618 Initializes the model.
619
620 @param domain: problem domain
621 @type domain: L{Domain}
622 @param stress: initial (deviatoric) stress
623 @type stress: a tensor value/field of order 2
624 @param v: initial velocity field
625 @type stress: a vector value/field
626 @param p: initial pressure
627 @type p: a scalar value/field
628 @param t: initial time
629 @type t: C{float}
630 @param numMaterials: number of materials
631 @type numMaterials: C{int}
632 @param verbose: if C{True} some informations are printed.
633 @type verbose: C{bool}
634 @param adaptSubTolerance: If True the tolerance for subproblem is set automatically.
635 @type adaptSubTolerance: C{bool}
636 """
637 PowerLaw. __init__(self, numMaterials,verbose)
638 Rheology. __init__(self, domain, stress, v, p, t, verbose)
639 self.__solver=StokesProblemCartesian(self.getDomain(),verbose=verbose,adaptSubTolerance=adaptSubTolerance)
640 self.__eta_eff=None
641 self.setTolerance()
642 self.setFlowTolerance()
643
644 def update(self, dt, iter_max=100, inner_iter_max=20, verbose=False, usePCG=True):
645 """
646 Updates stress, velocity and pressure for time increment dt.
647
648 @param dt: time increment
649 @param inner_iter_max: maximum number of iteration steps in the
650 incompressible solver
651 @param verbose: prints some infos in the incompressible solver
652 """
653 if self.checkVerbose(): print "IncompressibleIsotropicFlowCartesian: start iteration for t = %s."%(self.getTime()+dt,)
654 v_last=self.getVelocity()
655 s_last=self.getDeviatoricStress()
656 F=self.getForce()
657 f=self.getSurfaceForce()
658 rf=self.getRestorationFactor()
659 mask_v,v_b=self.getVelocityConstraint()
660 mu=self.getElasticShearModulus()
661 #=========================================================================
662 #
663 # we use velocity and pressure from the last time step as initial guess:
664 #
665 v=v_last
666 p=self.getPressure()
667 #
668 # calculate eta_eff if we don't have one or elasticity is present.
669 #
670 if self.__eta_eff == None or mu!=None:
671 D=self.__getDeviatoricStrain(v)
672 if mu==None:
673 gamma=util.sqrt(2.)*util.length(D)
674 else:
675 gamma=util.sqrt(2.)*util.length(D+s_last/(2*dt*mu))
676 if self.__eta_eff == None:
677 eta0=None
678 else:
679 eta0=self.__eta_eff
680 eta_eff=self.getEtaEff(gamma, pressure=p,dt=dt, eta0=eta0, iter_max=iter_max)
681 if self.checkVerbose(): print "IncompressibleIsotropicFlowCartesian: eta_eff has been initialied."
682 else:
683 eta_eff = self.__eta_eff
684 iter=0
685 converged=False
686 while not converged:
687 #
688 # intialize the solver
689 #
690 if mu==None:
691 stress0=Data()
692 else:
693 stress0=-(eta_eff/(dt*mu))*s_last
694 self.__solver.initialize(f=F,fixed_u_mask=mask_v,eta=eta_eff,surface_stress=f,stress=stress0, restoration_factor=rf)
695 #
696 # get a new velcocity and pressure:
697 #
698 if mask_v.isEmpty() or v_b.isEmpty():
699 v0=v
700 else:
701 v0=v_b*mask_v+v*(1.-mask_v)
702 v,p=self.__solver.solve(v0,p,verbose=self.checkVerbose(),max_iter=inner_iter_max,usePCG=usePCG)
703 #
704 # update eta_eff:
705 #
706 D=self.__getDeviatoricStrain(v)
707 if mu==None:
708 gamma=util.sqrt(2.)*util.length(D)
709 else:
710 gamma=util.sqrt(2.)*util.length(D+s_last/(2*dt*mu))
711 eta_eff_old ,eta_eff=eta_eff, self.getEtaEff(gamma, pressure=p,dt=dt, eta0=eta_eff, iter_max=iter_max)
712 if self.checkVerbose(): print "IncompressibleIsotropicFlowCartesian: eta_eff has been updated."
713 #
714 # check the change on eta_eff:
715 #
716 diff=util.Lsup(eta_eff_old-eta_eff)
717 n=util.Lsup(eta_eff)
718 if self.checkVerbose(): print "IncompressibleIsotropicFlowCartesian: step %s: max. rel. change in eta_eff is %s."%(iter,diff/n)
719 converged = diff <= self.getTolerance()* n
720 iter+=1
721 if iter >= iter_max:
722 raise MaxIterReached,"maximum number of iteration steps on time step %e reached."%(self.getTime()+dt)
723 #
724 # finally we can update the return values:
725 #
726 self.setPressure(p)
727 self.setVelocity(v)
728 self.setDeviatoricStrain(D)
729 if mu==None:
730 stress=(2*eta_eff)*D
731 else:
732 stress=(2.*eta_eff)*(D+s_last/(2*dt*mu))
733 self.setDeviatoricStress(stress)
734 self.__eta_eff = eta_eff
735 self.setTime(self.getTime()+dt)
736 if self.checkVerbose(): print "IncompressibleIsotropicFlowCartesian: iteration on time step %s completed after %s steps."%(self.getTime(),iter)
737 return self.getVelocity(), self.getPressure()
738
739 def getCurrentEtaEff(self):
740 """
741 returns the effective viscosity
742 """
743 return self.__eta_eff
744
745 def __getDeviatoricStrain(self, v):
746 """
747 Returns deviatoric strain of velocity v:
748 """
749 return util.deviatoric(util.symmetric(util.grad(v)))
750
751 def setFlowTolerance(self, tol=1.e-4):
752 """
753 Sets the relative tolerance for the flow solver. See L{StokesProblemCartesian.setTolerance} for details.
754
755 @param tol: desired relative tolerance for the flow solver
756 @type tol: positive C{float}
757 """
758 self.__solver.setTolerance(tol)
759 def getFlowTolerance(self):
760 """
761 Returns the relative tolerance for the flow solver
762
763 @return: tolerance of the flow solver
764 @rtype: C{float}
765 """
766 return self.__solver.getTolerance()
767
768 def getSolverOptionsVelocity(self):
769 """
770 returns the solver options used solve the equation for velocity in the
771 incompressible solver, see L{StokesProblemCartesian.getSolverOptionsVelocity} for details.
772
773 @rtype: L{SolverOptions}
774 """
775 return self.__solver.getSolverOptionsVelocity()
776 def setSolverOptionsVelocity(self, options=None):
777 """
778 set the solver options for solving the equation for velocity in the
779 incompressible solver, see L{StokesProblemCartesian.setSolverOptionsVelocity} for details.
780
781 @param options: new solver options
782 @type options: L{SolverOptions}
783 """
784 self.__solver.setSolverOptionsVelocity(options)
785
786 def getSolverOptionsPressure(self):
787 """
788 returns the solver options used solve the equation for pressure in the
789 incompressible solver, see L{StokesProblemCartesian.getSolverOptionsPressure} for details.
790 @rtype: L{SolverOptions}
791 """
792 return self.__solver.getSolverOptionsPressure()
793 def setSolverOptionsPressure(self, options=None):
794 """
795 set the solver options for solving the equation for pressure in the
796 incompressible solver, see L{StokesProblemCartesian.setSolverOptionsPressure} for details.
797 @param options: new solver options
798 @type options: L{SolverOptions}
799 """
800 self.__solver.setSolverOptionsPressure(options)
801
802 def setSolverOptionsDiv(self, options=None):
803 """
804 set the solver options for solving the equation to project the divergence of
805 the velocity onto the function space of pressure in the
806 incompressible solver, see L{StokesProblemCartesian.setSolverOptionsDiv} for details.
807
808 @param options: new solver options
809 @type options: L{SolverOptions}
810 """
811 self.__solver.setSolverOptionsDiv(options)
812 def getSolverOptionsDiv(self):
813 """
814 returns the solver options for solving the equation to project the divergence of
815 the velocity onto the function space of presure in the
816 incompressible solver, see L{StokesProblemCartesian.getSolverOptionsDiv} for details..
817
818 @rtype: L{SolverOptions}
819 """
820 return self.__solver.getSolverOptionsDiv()
821

  ViewVC Help
Powered by ViewVC 1.1.26