/[escript]/branches/restext/escript/py_src/rheologies.py
ViewVC logotype

Contents of /branches/restext/escript/py_src/rheologies.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2616 - (show annotations)
Wed Aug 19 06:33:08 2009 UTC (9 years, 10 months ago) by jfenwick
File MIME type: text/x-python
File size: 29164 byte(s)
Changed last M{} markup.
Had some issues working out how to deal with expressions which ended with *
ended up doing X* ->  *(X * )*

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 *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 *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: ``int``
59 :param verbose: if ``True`` some informations are printed.
60 :type verbose: ``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: ``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: ``int``
88 :return: ``True`` is the id is valid
89 :rtype: ``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 ``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 ``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 ``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 """
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 ``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 """
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 ``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 """
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: ``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 *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 ``float`` if present
240 :param iter_max: maximum number of iteration steps.
241 :type iter_max: ``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: `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 v: a vector value/field
321 :param p: initial pressure
322 :type p: a scalar value/field
323 :param t: initial time
324 :type t: ``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.setDeviatoricStress(stress)
337 self.setVelocity(v)
338 self.setPressure(p)
339 self.setDeviatoricStrain()
340 self.setTime(t)
341 #=============================================================
342 self.setExternals(F=Data(), f=Data(), fixed_v_mask=Data(), v_boundary=Data())
343
344 def getDomain(self):
345 """
346 returns the domain.
347
348 :return: the domain
349 :rtype: `Domain`
350 """
351 return self.__domain
352
353 def getTime(self):
354 """
355 Returns current time.
356
357 :return: current time
358 :rtype: ``float``
359 """
360 return self.__t
361
362 def setExternals(self, F=None, f=None, fixed_v_mask=None, v_boundary=None):
363 """
364 sets external forces and velocity constraints
365
366 :param F: external force
367 :type F: vector value/field
368 :param f: surface force
369 :type f: vector value/field on boundary
370 :param fixed_v_mask: location of constraints maked by positive values
371 :type fixed_v_mask: vector value/field
372 :param v_boundary: value of velocity at location of constraints
373 :type v_boundary: vector value/field
374 :note: Only changing parameters need to be specified.
375 """
376 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
381 def getForce(self):
382 """
383 Returns the external force
384
385 :return: external force
386 :rtype: `Data`
387 """
388 return self.__F
389
390 def getSurfaceForce(self):
391 """
392 Returns the surface force
393
394 :return: surface force
395 :rtype: `Data`
396 """
397 return self.__f
398
399 def getVelocityConstraint(self):
400 """
401 Returns the constraint for the velocity as a pair of the
402 mask of the location of the constraint and the values.
403
404 :return: the locations of fixed velocity and value of velocities at these locations
405 :rtype: ``tuple`` of `Data` s
406 """
407 return self.__fixed_v_mask, self.__v_boundary
408
409 def checkVerbose(self):
410 """
411 Returns True if verbose is switched on
412
413 :return: value of verbosity flag
414 :rtype: ``bool``
415 """
416 return self.__verbose
417
418 def setTime(self,t=0.):
419 """
420 Updates current time.
421
422 :param t: new time mark
423 :type t: ``float``
424 """
425 self.__t=t
426 #=======================================================================================
427 def getStress(self):
428 """
429 Returns current stress.
430
431 :return: current stress
432 :rtype: `Data` of rank 2
433 """
434 s=self.getDeviatoricStress()
435 p=self.getPressure()
436 k=util.kronecker(self.getDomain())
437 return s-p*(k/trace(k))
438
439 def getDeviatoricStress(self):
440 """
441 Returns current deviatoric stress.
442
443 :return: current deviatoric stress
444 :rtype: `Data` of rank 2
445 """
446 return self.__stress
447
448 def setDeviatoricStress(self, stress):
449 """
450 Sets the current deviatoric stress
451
452 :param stress: new deviatoric stress
453 :type stress: `Data` of rank 2
454 """
455 dom=self.getDomain()
456 s=util.interpolate(stress,Function(dom))
457 self.__stress=util.deviatoric(s)
458
459 def getPressure(self):
460 """
461 Returns current pressure.
462
463 :return: current stress
464 :rtype: scalar `Data`
465 """
466 return self.__p
467
468 def setPressure(self, p):
469 """
470 Sets current pressure.
471 :param p: new deviatoric stress
472 :type p: scalar `Data`
473 """
474 self.__p=util.interpolate(p,ReducedSolution(self.getDomain()))
475
476 def getVelocity(self):
477 """
478 Returns current velocity.
479
480 :return: current velocity
481 :rtype: vector `Data`
482 """
483 return self.__v
484
485 def setVelocity(self, v):
486 """
487 Sets current velocity.
488
489 :param v: new current velocity
490 :type v: vector `Data`
491 """
492 self.__v=util.interpolate(v,Solution(self.getDomain()))
493
494 def setDeviatoricStrain(self, D=None):
495 """
496 set deviatoric strain
497
498 :param D: new deviatoric strain. If D is not present the current velocity is used.
499 :type D: `Data` of rank 2
500 """
501 if D==None: D=util.deviatoric(util.symmetric(util.grad(2.*self.getVelocity())))
502 self.__D=util.deviatoric(util.interpolate(D,Function(self.getDomain())))
503
504 def getDeviatoricStrain(self):
505 """
506 Returns deviatoric strain of current velocity.
507
508 :return: deviatoric strain
509 :rtype: `Data` of rank 2
510 """
511 return self.__D
512
513 def getTau(self):
514 """
515 Returns current second invariant of deviatoric stress
516
517 :return: second invariant of deviatoric stress
518 :rtype: scalar `Data`
519 """
520 s=self.getDeviatoricStress()
521 return util.sqrt(0.5)*util.length(s)
522
523 def getGammaDot(self):
524 """
525 Returns current second invariant of deviatoric strain
526
527 :return: second invariant of deviatoric strain
528 :rtype: scalar `Data`
529 """
530 s=self.getDeviatoricStrain()
531 return util.sqrt(2)*util.length(s)
532
533 def setTolerance(self,tol=1.e-4):
534 """
535 Sets the tolerance used to terminate the iteration on a time step.
536 See the implementation of the rheology for details.
537
538 :param tol: relative tolerance to terminate iteration on time step.
539 :type tol: positive ``float``
540 """
541 if tol<=0.:
542 raise ValueError,"tolerance must be non-negative."
543 self.__tol=tol
544
545 def getTolerance(self):
546 """
547 Returns the set tolerance for terminate the iteration on a time step.
548
549 :rtype: positive ``float``
550 """
551 return self.__tol
552
553 #=======================================================================
554 def setFlowTolerance(self, tol=1.e-4):
555 """
556 Sets the relative tolerance for the flow solver
557
558 :param tol: desired relative tolerance for the flow solver
559 :type tol: positive ``float``
560 :note: Typically this method is overwritten by a subclass.
561 """
562 pass
563 def getFlowTolerance(self):
564 """
565 Returns the relative tolerance for the flow solver
566
567 :return: tolerance of the flow solver
568 :rtype: ``float``
569 :note: Typically this method is overwritten by a subclass.
570 """
571 pass
572 #====================================================================================================================================
573
574 class IncompressibleIsotropicFlowCartesian(PowerLaw,Rheology):
575 """
576 This class implements the rheology of an isotropic Kelvin material.
577
578 Typical usage::
579
580 sp = IncompressibleIsotropicFlow(domain, stress=0, v=0)
581 sp.setTolerance()
582 sp.initialize(...)
583 v,p = sp.solve()
584
585 :note: This model has been used in the self-consistent plate-mantle model
586 proposed in U{Hans-Bernd Muhlhaus<emailto:h.muhlhaus@uq.edu.au>}
587 and U{Klaus Regenauer-Lieb<mailto:klaus.regenauer-lieb@csiro.au>}:
588 I{Towards a self-consistent plate mantle model that includes elasticity: simple benchmarks and application to basic modes of convection},
589 see U{doi: 10.1111/j.1365-246X.2005.02742.x<http://www3.interscience.wiley.com/journal/118661486/abstract>}
590 """
591 def __init__(self, domain, stress=0, v=0, p=0, t=0, numMaterials=1, verbose=True, adaptSubTolerance=True):
592 """
593 Initializes the model.
594
595 :param domain: problem domain
596 :type domain: `Domain`
597 :param stress: initial (deviatoric) stress
598 :type stress: a tensor value/field of order 2
599 :param v: initial velocity field
600 :type v: a vector value/field
601 :param p: initial pressure
602 :type p: a scalar value/field
603 :param t: initial time
604 :type t: ``float``
605 :param numMaterials: number of materials
606 :type numMaterials: ``int``
607 :param verbose: if ``True`` some informations are printed.
608 :type verbose: ``bool``
609 :param adaptSubTolerance: If True the tolerance for subproblem is set automatically.
610 :type adaptSubTolerance: ``bool``
611 """
612 PowerLaw. __init__(self, numMaterials,verbose)
613 Rheology. __init__(self, domain, stress, v, p, t, verbose)
614 self.__solver=StokesProblemCartesian(self.getDomain(),verbose=verbose,adaptSubTolerance=adaptSubTolerance)
615 self.__eta_eff=None
616 self.setTolerance()
617 self.setFlowTolerance()
618
619 def update(self, dt, iter_max=100, inner_iter_max=20, verbose=False, usePCG=True):
620 """
621 Updates stress, velocity and pressure for time increment dt.
622
623 :param dt: time increment
624 :param inner_iter_max: maximum number of iteration steps in the
625 incompressible solver
626 :param verbose: prints some infos in the incompressible solver
627 """
628 if self.checkVerbose(): print "IncompressibleIsotropicFlowCartesian: start iteration for t = %s."%(self.getTime()+dt,)
629 v_last=self.getVelocity()
630 s_last=self.getDeviatoricStress()
631 F=self.getForce()
632 f=self.getSurfaceForce()
633 mask_v,v_b=self.getVelocityConstraint()
634 mu=self.getElasticShearModulus()
635 #=========================================================================
636 #
637 # we use velocity and pressure from the last time step as initial guess:
638 #
639 v=v_last
640 p=self.getPressure()
641 #
642 # calculate eta_eff if we don't have one or elasticity is present.
643 #
644 if self.__eta_eff == None or mu!=None:
645 D=self.__getDeviatoricStrain(v)
646 if mu==None:
647 gamma=util.sqrt(2.)*util.length(D)
648 else:
649 gamma=util.sqrt(2.)*util.length(D+s_last/(2*dt*mu))
650 if self.__eta_eff == None:
651 eta0=None
652 else:
653 eta0=self.__eta_eff
654 eta_eff=self.getEtaEff(gamma, pressure=p,dt=dt, eta0=eta0, iter_max=iter_max)
655 if self.checkVerbose(): print "IncompressibleIsotropicFlowCartesian: eta_eff has been initialied."
656 else:
657 eta_eff = self.__eta_eff
658 iter=0
659 converged=False
660 while not converged:
661 #
662 # intialize the solver
663 #
664 if mu==None:
665 stress0=Data()
666 else:
667 stress0=-(eta_eff/(dt*mu))*s_last
668 self.__solver.initialize(f=F,fixed_u_mask=mask_v,eta=eta_eff,surface_stress=f,stress=stress0)
669 #
670 # get a new velcocity and pressure:
671 #
672 if mask_v.isEmpty() or v_b.isEmpty():
673 v0=v
674 else:
675 v0=v_b*mask_v+v*(1.-mask_v)
676 v,p=self.__solver.solve(v0,p,verbose=self.checkVerbose(),max_iter=inner_iter_max,usePCG=usePCG)
677 #
678 # update eta_eff:
679 #
680 D=self.__getDeviatoricStrain(v)
681 if mu==None:
682 gamma=util.sqrt(2.)*util.length(D)
683 else:
684 gamma=util.sqrt(2.)*util.length(D+s_last/(2*dt*mu))
685 eta_eff_old ,eta_eff=eta_eff, self.getEtaEff(gamma, pressure=p,dt=dt, eta0=eta_eff, iter_max=iter_max)
686 if self.checkVerbose(): print "IncompressibleIsotropicFlowCartesian: eta_eff has been updated."
687 #
688 # check the change on eta_eff:
689 #
690 diff=util.Lsup(eta_eff_old-eta_eff)
691 n=util.Lsup(eta_eff)
692 if self.checkVerbose(): print "IncompressibleIsotropicFlowCartesian: step %s: max. rel. change in eta_eff is %s."%(iter,diff/n)
693 converged = diff <= self.getTolerance()* n
694 iter+=1
695 if iter >= iter_max:
696 raise MaxIterReached,"maximum number of iteration steps on time step %e reached."%(self.getTime()+dt)
697 #
698 # finally we can update the return values:
699 #
700 self.setPressure(p)
701 self.setVelocity(v)
702 self.setDeviatoricStrain(D)
703 if mu==None:
704 stress=(2*eta_eff)*D
705 else:
706 stress=(2.*eta_eff)*(D+s_last/(2*dt*mu))
707 self.setDeviatoricStress(stress)
708 self.__eta_eff = eta_eff
709 self.setTime(self.getTime()+dt)
710 if self.checkVerbose(): print "IncompressibleIsotropicFlowCartesian: iteration on time step %s completed after %s steps."%(self.getTime(),iter)
711 return self.getVelocity(), self.getPressure()
712
713 def getCurrentEtaEff(self):
714 """
715 returns the effective viscosity
716 """
717 return self.__eta_eff
718
719 def __getDeviatoricStrain(self, v):
720 """
721 Returns deviatoric strain of velocity v:
722 """
723 return util.deviatoric(util.symmetric(util.grad(v)))
724
725 def setFlowTolerance(self, tol=1.e-4):
726 """
727 Sets the relative tolerance for the flow solver. See `StokesProblemCartesian.setTolerance` for details.
728
729 :param tol: desired relative tolerance for the flow solver
730 :type tol: positive ``float``
731 """
732 self.__solver.setTolerance(tol)
733 def getFlowTolerance(self):
734 """
735 Returns the relative tolerance for the flow solver
736
737 :return: tolerance of the flow solver
738 :rtype: ``float``
739 """
740 return self.__solver.getTolerance()
741
742 def getSolverOptionsVelocity(self):
743 """
744 returns the solver options used solve the equation for velocity in the
745 incompressible solver, see `StokesProblemCartesian.getSolverOptionsVelocity` for details.
746
747 :rtype: `SolverOptions`
748 """
749 return self.__solver.getSolverOptionsVelocity()
750 def setSolverOptionsVelocity(self, options=None):
751 """
752 set the solver options for solving the equation for velocity in the
753 incompressible solver, see `StokesProblemCartesian.setSolverOptionsVelocity` for details.
754
755 :param options: new solver options
756 :type options: `SolverOptions`
757 """
758 self.__solver.setSolverOptionsVelocity(options)
759
760 def getSolverOptionsPressure(self):
761 """
762 returns the solver options used solve the equation for pressure in the
763 incompressible solver, see `StokesProblemCartesian.getSolverOptionsPressure` for details.
764 :rtype: `SolverOptions`
765 """
766 return self.__solver.getSolverOptionsPressure()
767 def setSolverOptionsPressure(self, options=None):
768 """
769 set the solver options for solving the equation for pressure in the
770 incompressible solver, see `StokesProblemCartesian.setSolverOptionsPressure` for details.
771 :param options: new solver options
772 :type options: `SolverOptions`
773 """
774 self.__solver.setSolverOptionsPressure(options)
775
776 def setSolverOptionsDiv(self, options=None):
777 """
778 set the solver options for solving the equation to project the divergence of
779 the velocity onto the function space of pressure in the
780 incompressible solver, see `StokesProblemCartesian.setSolverOptionsDiv` for details.
781
782 :param options: new solver options
783 :type options: `SolverOptions`
784 """
785 self.__solver.setSolverOptionsDiv(options)
786 def getSolverOptionsDiv(self):
787 """
788 returns the solver options for solving the equation to project the divergence of
789 the velocity onto the function space of presure in the
790 incompressible solver, see `StokesProblemCartesian.getSolverOptionsDiv` for details..
791
792 :rtype: `SolverOptions`
793 """
794 return self.__solver.getSolverOptionsDiv()
795

  ViewVC Help
Powered by ViewVC 1.1.26