/[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 2415 - (show annotations)
Wed May 13 02:48:39 2009 UTC (10 years, 5 months ago) by gross
File MIME type: text/x-python
File size: 28357 byte(s)
FileWriter added: this class takes care of writing data which are global in MPI to a file. It is recommended to use this class rather then the build in open as it takes care of the case of many processors.
1 ########################################################
2 #
3 # Copyright (c) 2003-2008 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-2008 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.-1./power)}
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=util.sqrt(util.EPSILON)):
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=10):
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=[1.-1./p 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.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: L{Domain}
350 """
351 return self.__domain
352
353 def getTime(self):
354 """
355 Returns current time.
356
357 @return: current time
358 @rtype: C{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: L{Data}
387 """
388 return self.__F
389
390 def getSurfaceForce(self):
391 """
392 Returns the surface force
393
394 @return: surface force
395 @rtype: L{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: C{tuple} of L{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: C{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: C{float}
424 """
425 self.__t=t
426 #=======================================================================================
427 def getStress(self):
428 """
429 Returns current stress.
430
431 @return: current stress
432 @rtype: L{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: L{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: L{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 L{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 L{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 L{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 L{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: L{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: L{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 L{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 L{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 C{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 C{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 C{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: C{float}
569 @note: Typically this method is overwritten by a subclass.
570 """
571 pass
572 def setFlowSubTolerance(self, tol=1.e-8):
573 """
574 Sets the relative tolerance for the subsolver of the flow solver
575
576 @param tol: desired relative tolerance for the subsolver
577 @type tol: positive C{float}
578 @note: Typically this method is overwritten by a subclass.
579 """
580 pass
581 def getFlowSubTolerance(self):
582 """
583 Returns the relative tolerance for the subsolver of the flow solver
584
585 @return: tolerance of the flow subsolver
586 @rtype: C{float}
587 @note: Typically this method is overwritten by a subclass.
588 """
589 pass
590
591
592 #====================================================================================================================================
593
594 class IncompressibleIsotropicFlowCartesian(PowerLaw,Rheology):
595 """
596 This class implements the rheology of an isotropic Kelvin material.
597
598 Typical usage::
599
600 sp = IncompressibleIsotropicFlow(domain, stress=0, v=0)
601 sp.setTolerance()
602 sp.initialize(...)
603 v,p = sp.solve()
604
605 @note: This model has been used in the self-consistent plate-mantle model
606 proposed in U{Hans-Bernd Muhlhaus<emailto:h.muhlhaus@uq.edu.au>}
607 and U{Klaus Regenauer-Lieb<mailto:klaus.regenauer-lieb@csiro.au>}:
608 I{Towards a self-consistent plate mantle model that includes elasticity: simple benchmarks and application to basic modes of convection},
609 see U{doi: 10.1111/j.1365-246X.2005.02742.x<http://www3.interscience.wiley.com/journal/118661486/abstract?CRETRY=1&SRETRY=0>}.
610
611 """
612 def __init__(self, domain, stress=0, v=0, p=0, t=0, numMaterials=1, verbose=True):
613 """
614 Initializes the model.
615
616 @param domain: problem domain
617 @type domain: L{Domain}
618 @param stress: initial (deviatoric) stress
619 @type stress: a tensor value/field of order 2
620 @param v: initial velocity field
621 @type stress: a vector value/field
622 @param p: initial pressure
623 @type p: a scalar value/field
624 @param t: initial time
625 @type t: C{float}
626 @param numMaterials: number of materials
627 @type numMaterials: C{int}
628 @param verbose: if C{True} some informations are printed.
629 @type verbose: C{bool}
630 """
631 PowerLaw. __init__(self, numMaterials,verbose)
632 Rheology. __init__(self, domain, stress, v, p, t, verbose)
633 self.__solver=StokesProblemCartesian(self.getDomain(),verbose=verbose)
634 self.__eta_eff=None
635 self.setTolerance()
636 self.setFlowTolerance()
637 self.setFlowSubTolerance()
638
639 def update(self, dt, iter_max=100, inner_iter_max=20, verbose=False, usePCG=True):
640 """
641 Updates stress, velocity and pressure for time increment dt.
642
643 @param dt: time increment
644 @param inner_iter_max: maximum number of iteration steps in the
645 incompressible solver
646 @param verbose: prints some infos in the incompressible solver
647 """
648 if self.checkVerbose(): print "IncompressibleIsotropicFlowCartesian: start iteration for t = %s."%(self.getTime()+dt,)
649 v_last=self.getVelocity()
650 s_last=self.getDeviatoricStress()
651 F=self.getForce()
652 f=self.getSurfaceForce()
653 mask_v,v_b=self.getVelocityConstraint()
654 mu=self.getElasticShearModulus()
655 #=========================================================================
656 #
657 # we use velocity and pressure from the last time step as initial guess:
658 #
659 v=v_last
660 p=self.getPressure()
661 #
662 # calculate eta_eff if we don't have one or elasticity is present.
663 #
664 if self.__eta_eff == None or mu!=None:
665 D=self.__getDeviatoricStrain(v)
666 if mu==None:
667 gamma=util.sqrt(2.)*util.length(D)
668 else:
669 gamma=util.sqrt(2.)*util.length(D+s_last/(2*dt*mu))
670 if self.__eta_eff == None:
671 eta0=None
672 else:
673 eta0=self.__eta_eff
674 eta_eff=self.getEtaEff(gamma, pressure=p,dt=dt, eta0=eta0, iter_max=iter_max)
675 if self.checkVerbose(): print "IncompressibleIsotropicFlowCartesian: eta_eff has been initialied."
676 else:
677 eta_eff = self.__eta_eff
678
679 iter=0
680 converged=False
681 while not converged:
682 #
683 # intialize the solver
684 #
685 if mu==None:
686 stress0=Data()
687 else:
688 stress0=-(eta_eff/(dt*mu))*s_last
689 self.__solver.initialize(f=F,fixed_u_mask=mask_v,eta=eta_eff,surface_stress=f,stress=stress0)
690 #
691 # get a new velcocity and pressure:
692 #
693 if mask_v.isEmpty() or v_b.isEmpty():
694 v0=v
695 else:
696 v0=v_b*mask_v+v*(1.-mask_v)
697 v,p=self.__solver.solve(v0,p,show_details=False,
698 verbose=self.checkVerbose(),max_iter=inner_iter_max,usePCG=usePCG)
699 util.saveVTK("v.xml",v=v)
700 util.saveVTK("p.xml",p=p)
701 #
702 # update eta_eff:
703 #
704 D=self.__getDeviatoricStrain(v)
705 print D[0,0]
706 print D[1,1]
707 print D[1,0]
708 if mu==None:
709 gamma=util.sqrt(2.)*util.length(D)
710 else:
711 gamma=util.sqrt(2.)*util.length(D+s_last/(2*dt*mu))
712 print gamma
713 1/0
714 eta_eff_old ,eta_eff=eta_eff, self.getEtaEff(gamma, pressure=p,dt=dt, eta0=eta_eff, iter_max=iter_max)
715 if self.checkVerbose(): print "IncompressibleIsotropicFlowCartesian: eta_eff has been updated."
716 #
717 # check the change on eta_eff:
718 #
719 diff=util.Lsup(eta_eff_old-eta_eff)
720 n=util.Lsup(eta_eff)
721 if self.checkVerbose(): print "IncompressibleIsotropicFlowCartesian: step %s: max. change in eta_eff is %s."%(iter,diff)
722 converged = diff <= self.getTolerance()* n
723 iter+=1
724 if iter >= iter_max:
725 raise MaxIterReached,"maximum number of iteration steps on time step %e reached."%(self.getTime()+dt)
726 #
727 # finally we can update the return values:
728 #
729 self.setPressure(p)
730 self.setVelocity(v)
731 self.setDeviatoricStrain(D)
732 if mu==None:
733 stress=(2*eta_eff)*D
734 else:
735 stress=(2.*eta_eff)*(D+s_last/(2*dt*mu))
736 self.setDeviatoricStress(stress)
737 self.__eta_eff = eta_eff
738 self.setTime(self.getTime()+dt)
739 if self.checkVerbose(): print "IncompressibleIsotropicFlowCartesian: iteration on time step %s completed after %s steps."%(self.getTime(),iter)
740 return self.getVelocity(), self.getPressure()
741
742 def __getDeviatoricStrain(self, v):
743 """
744 Returns deviatoric strain of velocity v:
745 """
746 return util.deviatoric(util.symmetric(util.grad(2.*v)))
747
748 def setFlowTolerance(self, tol=1.e-6):
749 """
750 Sets the relative tolerance for the flow solver. See L{StokesProblemCartesian.setTolerance} for details.
751
752 @param tol: desired relative tolerance for the flow solver
753 @type tol: positive C{float}
754 """
755 self.__solver.setTolerance(tol)
756 def getFlowTolerance(self):
757 """
758 Returns the relative tolerance for the flow solver
759
760 @return: tolerance of the flow solver
761 @rtype: C{float}
762 """
763 return self.__solver.getTolerance()
764 def setFlowSubTolerance(self, tol=1.e-12):
765 """
766 Sets the relative tolerance for the subsolver of the flow solver. See L{StokesProblemCartesian.setSubProblemTolerance} for details
767
768 @param tol: desired relative tolerance for the subsolver
769 @type tol: positive C{float}
770 """
771 self.__solver.setSubProblemTolerance(tol)
772 def getFlowSubTolerance(self):
773 """
774 Returns the relative tolerance for the subsolver of the flow solver
775
776 @return: tolerance of the flow subsolver
777 @rtype: C{float}
778 """
779 return self.__solver.getSubProblemTolerance()
780

  ViewVC Help
Powered by ViewVC 1.1.26