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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26