/[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 5024 - (show annotations)
Tue Jun 10 08:46:07 2014 UTC (5 years, 2 months ago) by jfenwick
File MIME type: text/x-python
File size: 27533 byte(s)
Fixing more == None

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

  ViewVC Help
Powered by ViewVC 1.1.26