/[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 3432 - (show annotations)
Fri Jan 7 01:32:07 2011 UTC (8 years, 7 months ago) by jfenwick
File MIME type: text/x-python
File size: 27308 byte(s)
Made import statements a bit more specific to clean up the epydoc
1 ########################################################
2 #
3 # Copyright (c) 2003-2010 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-2010 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 import escript
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: 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 if pressure == None:
245 p2 = None
246 else:
247 p2=(abs(pressure)+pressure)/2.
248 SMALL=1./(util.DBLE_MAX/100.)
249 numMaterial=self.getNumMaterials()
250 s=[p-1. for p in self.getPower() ]
251 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 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 if tau_Y==None and fric==None:
273 eta_max=None
274 else:
275 if fric == None or p2==None:
276 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 eta_max=fric*util.clip(tau_Y/fric+p2,minval=0)/(gamma_dot+SMALL*util.whereZero(gamma_dot))
282 rtol=self.getEtaTolerance()
283 iter =0
284 converged=False
285 tau=eta_eff*gamma_dot
286 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 while not converged:
288 if iter>max(iter_max,1):
289 raise RuntimeError,"tolerance not reached after %s steps."%max(iter_max,1)
290 #===========================================
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 if self.__verbose: print "PowerLaw: step %s: correction = %s, max eta_eff = %s, max tau= %s"%(iter, d, l,util.Lsup(tau))
305 converged= d<= rtol* l
306 if self.__verbose: print "PowerLaw: Start calculation of eta_eff finalized after %s steps."%iter
307 return eta_eff
308
309 #====================================================================================================================================
310 class Rheology(object):
311 """
312 General framework to implement a rheology
313 """
314 def __init__(self, domain, stress=None, v=None, p=None, t=0, verbose=True):
315 """
316 Initializes the rheology
317
318 :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 """
329 self.__domain=domain
330 self.__t=t
331 self.__verbose=verbose
332 #=======================
333 #
334 # state variables:
335 #
336 if stress == None: stress=Tensor(0.,escript.Function(self.__domain))
337 if v == None: v=Vector(0.,escript.Solution(self.__domain))
338 if p == None: p=Vector(0.,Reducedescript.Solution(self.__domain))
339 self.setStatus(t, v, p, stress)
340 self.setExternals(F=escript.Data(), f=escript.Data(), fixed_v_mask=escript.Data(), v_boundary=escript.Data(), restoration_factor=0)
341
342 def getDomain(self):
343 """
344 returns the domain.
345
346 :return: the domain
347 :rtype: `Domain`
348 """
349 return self.__domain
350
351 def getTime(self):
352 """
353 Returns current time.
354
355 :return: current time
356 :rtype: ``float``
357 """
358 return self.__t
359
360 def setExternals(self, F=None, f=None, fixed_v_mask=None, v_boundary=None, restoration_factor=None):
361 """
362 sets external forces and velocity constraints
363
364 :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 """
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 if restoration_factor!=None: self.__restoration_factor=restoration_factor
381
382 def getForce(self):
383 """
384 Returns the external force
385
386 :return: external force
387 :rtype: `Data`
388 """
389 return self.__F
390
391 def getSurfaceForce(self):
392 """
393 Returns the surface force
394
395 :return: surface force
396 :rtype: `Data`
397 """
398 return self.__f
399
400 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 :return: the locations of fixed velocity and value of velocities at these locations
406 :rtype: ``tuple`` of `Data` s
407 """
408 return self.__fixed_v_mask, self.__v_boundary
409
410 def getRestorationFactor(self):
411 """
412 Returns the restoring force factor
413
414 @return: restoring force factor
415 @rtype: C{float} or L{Data}
416 """
417 return self.__restoration_factor
418
419
420 def checkVerbose(self):
421 """
422 Returns True if verbose is switched on
423
424 :return: value of verbosity flag
425 :rtype: ``bool``
426 """
427 return self.__verbose
428
429 def setTime(self,t=0.):
430 """
431 Updates current time.
432
433 :param t: new time mark
434 :type t: ``float``
435 """
436 self.__t=t
437 #=======================================================================================
438 def getStress(self):
439 """
440 Returns current stress.
441
442 :return: current stress
443 :rtype: `Data` of rank 2
444 """
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 :return: current deviatoric stress
455 :rtype: `Data` of rank 2
456 """
457 return self.__stress
458
459 def setDeviatoricStress(self, stress):
460 """
461 Sets the current deviatoric stress
462
463 :param stress: new deviatoric stress
464 :type stress: `Data` of rank 2
465 """
466 dom=self.getDomain()
467 s=util.interpolate(stress,escript.Function(dom))
468 self.__stress=util.deviatoric(s)
469
470 def getPressure(self):
471 """
472 Returns current pressure.
473
474 :return: current stress
475 :rtype: scalar `Data`
476 """
477 return self.__p
478
479 def setPressure(self, p):
480 """
481 Sets current pressure.
482 :param p: new deviatoric stress
483 :type p: scalar `Data`
484 """
485 self.__p=util.interpolate(p,escript.ReducedSolution(self.getDomain()))
486
487 def getVelocity(self):
488 """
489 Returns current velocity.
490
491 :return: current velocity
492 :rtype: vector `Data`
493 """
494 return self.__v
495
496 def setVelocity(self, v):
497 """
498 Sets current velocity.
499
500 :param v: new current velocity
501 :type v: vector `Data`
502 """
503 self.__v=util.interpolate(v,escript.Solution(self.getDomain()))
504 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 self.setGammaDot()
522 self.setTime(t)
523
524 def setDeviatoricStrain(self, D=None):
525 """
526 set deviatoric strain
527
528 :param D: new deviatoric strain. If ``D`` is not present the current velocity is used.
529 :type D: `Data` of rank 2
530 """
531 if D==None:
532 self.__D=self.getDeviatoricStrain(self.getVelocity())
533 else:
534 self.__D=util.deviatoric(util.interpolate(D,escript.Function(self.getDomain())))
535
536 def getDeviatoricStrain(self, v=None):
537 """
538 Returns deviatoric strain of current velocity or if ``v`` is present the
539 deviatoric strain of velocity ``v``:
540
541 :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 :rtype: `Data` of rank 2
545 """
546 if v==None:
547 return self.__D
548 else:
549 return util.deviatoric(util.symmetric(util.grad(v)))
550
551 def getTau(self):
552 """
553 Returns current second invariant of deviatoric stress
554
555 :return: second invariant of deviatoric stress
556 :rtype: scalar `Data`
557 """
558 s=self.getDeviatoricStress()
559 return util.sqrt(0.5)*util.length(s)
560
561 def setGammaDot(self, gammadot=None):
562 """
563 set the second invariant of deviatoric strain rate. If ``gammadot`` is not present zero is used.
564
565 :param gammadot: second invariant of deviatoric strain rate.
566 :type gammadot: `Data` of rank 1
567 """
568 if gammadot == None:
569 self.__gammadot = escript.Scalar(0.,escript.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 :return: second invariant of deviatoric strain
580 :rtype: scalar `Data`
581 """
582 if D==None:
583 return self.__gammadot
584 else:
585 return util.sqrt(2.)*util.length(D)
586
587
588 #====================================================================================================================================
589
590 class IncompressibleIsotropicFlowCartesian(PowerLaw,Rheology, StokesProblemCartesian):
591 """
592 This class implements the rheology of an isotropic Kelvin material.
593
594 Typical usage::
595
596 sp = IncompressibleIsotropicFlowCartesian(domain, stress=0, v=0)
597 sp.initialize(...)
598 v,p = sp.solve()
599
600 :note: This model has been used in the self-consistent plate-mantle model
601 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 see U{doi: 10.1111/j.1365-246X.2005.02742.x<http://www3.interscience.wiley.com/journal/118661486/abstract>}
605 """
606 def __init__(self, domain, stress=0, v=0, p=0, t=0, numMaterials=1, verbose=True):
607 """
608 Initializes the model.
609
610 :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 """
625 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 self.__eta_eff=None
629
630 def getCurrentEtaEff(self):
631 """
632 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 gamma=self.getGammaDot(self.getDeviatoricStrain(v))
652 else:
653 gamma=self.getGammaDot(self.getDeviatoricStrain(v)+s_last/(2*dt*mu))
654
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=escript.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 Updates stress, velocity and pressure for time increment dt.
688
689 :param dt: time increment
690 :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 :param verbose: prints some infos in the incompressible solver
693 """
694 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 s_last=self.getDeviatoricStress()
704 mask_v,v_b=self.getVelocityConstraint()
705 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 #
714 if mask_v.isEmpty():
715 v0=v_last
716 else:
717 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
722 v,p=self._solve(v0,p_last,verbose=self.checkVerbose(),max_iter=iter_max,usePCG=usePCG, max_correction_steps=max_correction_steps)
723 #
724 # finally we can update the return values:
725 #
726 self.setPressure(p)
727 self.setVelocity(v)
728 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 self.__eta_eff = self.getEtaEff(self.getGammaDot(), pressure=p,dt=dt, eta0=self.__eta_eff_save, iter_max=self.__eta_iter_max)
736 self.setDeviatoricStress(2.*self.__eta_eff*D)
737 self.setTime(self.getTime()+dt)
738 if self.checkVerbose(): print "IncompressibleIsotropicFlowCartesian: iteration on time step %s completed."%(self.getTime(),)
739 return self.getVelocity(), self.getPressure()
740
741

  ViewVC Help
Powered by ViewVC 1.1.26