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

Annotation of /trunk/escript/py_src/rheologies.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2440 - (hide annotations)
Wed May 27 08:45:55 2009 UTC (11 years, 5 months ago) by gross
File MIME type: text/x-python
File size: 28257 byte(s)
convection works now again.
1 ksteube 1809 ########################################################
2 gross 1639 #
3 ksteube 1809 # Copyright (c) 2003-2008 by University of Queensland
4     # Earth Systems Science Computational Center (ESSCC)
5     # http://www.uq.edu.au/esscc
6 gross 1639 #
7 ksteube 1809 # 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 gross 1639 #
11 ksteube 1809 ########################################################
12 gross 1639
13 ksteube 1809 __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 jfenwick 2344 __url__="https://launchpad.net/escript-finley"
20 ksteube 1809
21 gross 1639 """
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 gross 2415 from flows import StokesProblemCartesian
37     from pdetools import MaxIterReached
38 gross 1639
39 gross 2300 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 gross 2438 M{eta=eta_N*(tau/tau_t)**(1./power-1.)}
45 gross 2300
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 gross 2301 self.__friction=None
70 gross 2300 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 gross 2440 def setEtaTolerance(self,rtol=1.e-4):
93 gross 2300 """
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 gross 2301 def setDruckerPragerLaw(self,tau_Y=None,friction=None):
112 gross 2300 """
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 gross 2438 def getEtaEff(self,gamma_dot, eta0=None, pressure=None,dt=None, iter_max=30):
228 gross 2300 """
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 gross 2301 SMALL=1./(util.DBLE_MAX/100.)
245 gross 2300 numMaterial=self.getNumMaterials()
246 gross 2438 s=[p-1. for p in self.getPower() ]
247 gross 2300 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 gross 2301 if dt<=0: raise ValueError,"time step size must be positive."
269     if tau_Y==None and fric==None:
270 gross 2300 eta_max=None
271     else:
272 gross 2301 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 gross 2300 rtol=self.getEtaTolerance()
280     iter =0
281     converged=False
282     tau=eta_eff*gamma_dot
283 gross 2415 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 gross 2300 while not converged:
285     if iter>max(iter_max,1):
286 gross 2301 raise RuntimeError,"tolerance not reached after %s steps."%max(iter_max,1)
287 gross 2300 #===========================================
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 gross 2415 if self.__verbose: print "PowerLaw: step %s: correction = %s, max eta_eff = %s, max tau= %s"%(iter, d, l,util.Lsup(tau))
302 gross 2300 converged= d<= rtol* l
303 gross 2415 if self.__verbose: print "PowerLaw: Start calculation of eta_eff finalized after %s steps."%iter
304 gross 2300 return eta_eff
305    
306     #====================================================================================================================================
307 gross 2415 class Rheology(object):
308 gross 1639 """
309 gross 2415 General framework to implement a rheology
310 gross 1639 """
311 gross 2415 def __init__(self, domain, stress=None, v=None, p=None, t=0, verbose=True):
312 gross 1639 """
313 gross 2415 Initializes the rheology
314 gross 1639
315     @param domain: problem domain
316 gross 2415 @type domain: L{Domain}
317     @param stress: initial (deviatoric) stress
318     @type stress: a tensor value/field of order 2
319 gross 1639 @param v: initial velocity field
320 gross 2415 @type stress: a vector value/field
321 gross 2100 @param p: initial pressure
322 gross 2415 @type p: a scalar value/field
323 gross 1639 @param t: initial time
324 gross 2415 @type t: C{float}
325 gross 1639 """
326     self.__domain=domain
327     self.__t=t
328 gross 2415 self.__verbose=verbose
329 gross 1639 #=======================
330 gross 2415 #
331 gross 1639 # state variables:
332     #
333 gross 2415 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 caltinay 2169
348 gross 2415 @return: the domain
349     @rtype: L{Domain}
350     """
351     return self.__domain
352 gross 1639
353 gross 2415 def getTime(self):
354     """
355     Returns current time.
356 caltinay 2169
357 gross 2415 @return: current time
358     @rtype: C{float}
359 gross 1639 """
360 gross 2415 return self.__t
361    
362     def setExternals(self, F=None, f=None, fixed_v_mask=None, v_boundary=None):
363 caltinay 2169 """
364 gross 2415 sets external forces and velocity constraints
365 gross 2100
366 gross 2415 @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 gross 1639 """
376 gross 2415 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 caltinay 2169
385 gross 2415 @return: external force
386     @rtype: L{Data}
387 gross 1639 """
388 gross 2415 return self.__F
389 gross 2100
390 gross 2415 def getSurfaceForce(self):
391 gross 1639 """
392 gross 2415 Returns the surface force
393    
394     @return: surface force
395     @rtype: L{Data}
396 gross 1639 """
397 gross 2415 return self.__f
398 gross 2100
399 gross 2415 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 gross 2100 def setTolerance(self,tol=1.e-4):
534 gross 1639 """
535 gross 2415 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 gross 1639 """
541 gross 2415 if tol<=0.:
542     raise ValueError,"tolerance must be non-negative."
543 gross 2100 self.__tol=tol
544    
545     def getTolerance(self):
546 gross 1639 """
547 gross 2415 Returns the set tolerance for terminate the iteration on a time step.
548    
549     @rtype: positive C{float}
550 gross 1639 """
551 gross 2100 return self.__tol
552    
553 gross 2415 #=======================================================================
554     def setFlowTolerance(self, tol=1.e-4):
555 gross 1639 """
556 gross 2415 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 gross 1639 """
562 gross 2415 pass
563     def getFlowTolerance(self):
564     """
565     Returns the relative tolerance for the flow solver
566 caltinay 2169
567 gross 2415 @return: tolerance of the flow solver
568     @rtype: C{float}
569     @note: Typically this method is overwritten by a subclass.
570 gross 1639 """
571 gross 2415 pass
572     def setFlowSubTolerance(self, tol=1.e-8):
573 gross 1639 """
574 gross 2415 Sets the relative tolerance for the subsolver of the flow solver
575 gross 2100
576 gross 2415 @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 caltinay 2169 """
580 gross 2415 pass
581     def getFlowSubTolerance(self):
582     """
583     Returns the relative tolerance for the subsolver of the flow solver
584 gross 1639
585 gross 2415 @return: tolerance of the flow subsolver
586     @rtype: C{float}
587     @note: Typically this method is overwritten by a subclass.
588 gross 1639 """
589 gross 2415 pass
590 gross 1639
591 gross 2100
592 gross 2415 #====================================================================================================================================
593 gross 1639
594 gross 2415 class IncompressibleIsotropicFlowCartesian(PowerLaw,Rheology):
595     """
596     This class implements the rheology of an isotropic Kelvin material.
597 caltinay 2169
598 gross 2415 Typical usage::
599 caltinay 2169
600 gross 2415 sp = IncompressibleIsotropicFlow(domain, stress=0, v=0)
601     sp.setTolerance()
602     sp.initialize(...)
603     v,p = sp.solve()
604 gross 2100
605 gross 2415 @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 gross 2438 see U{doi: 10.1111/j.1365-246X.2005.02742.x<http://www3.interscience.wiley.com/journal/118661486/abstract>}
610 gross 2415 """
611     def __init__(self, domain, stress=0, v=0, p=0, t=0, numMaterials=1, verbose=True):
612     """
613     Initializes the model.
614    
615     @param domain: problem domain
616     @type domain: L{Domain}
617     @param stress: initial (deviatoric) stress
618     @type stress: a tensor value/field of order 2
619     @param v: initial velocity field
620     @type stress: a vector value/field
621     @param p: initial pressure
622     @type p: a scalar value/field
623     @param t: initial time
624     @type t: C{float}
625     @param numMaterials: number of materials
626     @type numMaterials: C{int}
627     @param verbose: if C{True} some informations are printed.
628     @type verbose: C{bool}
629     """
630     PowerLaw. __init__(self, numMaterials,verbose)
631     Rheology. __init__(self, domain, stress, v, p, t, verbose)
632     self.__solver=StokesProblemCartesian(self.getDomain(),verbose=verbose)
633     self.__eta_eff=None
634     self.setTolerance()
635     self.setFlowTolerance()
636     self.setFlowSubTolerance()
637    
638     def update(self, dt, iter_max=100, inner_iter_max=20, verbose=False, usePCG=True):
639 gross 1639 """
640 caltinay 2169 Updates stress, velocity and pressure for time increment dt.
641 gross 1639
642     @param dt: time increment
643 caltinay 2169 @param inner_iter_max: maximum number of iteration steps in the
644     incompressible solver
645     @param verbose: prints some infos in the incompressible solver
646 gross 1639 """
647 gross 2415 if self.checkVerbose(): print "IncompressibleIsotropicFlowCartesian: start iteration for t = %s."%(self.getTime()+dt,)
648     v_last=self.getVelocity()
649     s_last=self.getDeviatoricStress()
650     F=self.getForce()
651     f=self.getSurfaceForce()
652     mask_v,v_b=self.getVelocityConstraint()
653     mu=self.getElasticShearModulus()
654     #=========================================================================
655     #
656     # we use velocity and pressure from the last time step as initial guess:
657     #
658     v=v_last
659     p=self.getPressure()
660     #
661     # calculate eta_eff if we don't have one or elasticity is present.
662     #
663     if self.__eta_eff == None or mu!=None:
664     D=self.__getDeviatoricStrain(v)
665     if mu==None:
666     gamma=util.sqrt(2.)*util.length(D)
667     else:
668     gamma=util.sqrt(2.)*util.length(D+s_last/(2*dt*mu))
669     if self.__eta_eff == None:
670     eta0=None
671     else:
672     eta0=self.__eta_eff
673     eta_eff=self.getEtaEff(gamma, pressure=p,dt=dt, eta0=eta0, iter_max=iter_max)
674     if self.checkVerbose(): print "IncompressibleIsotropicFlowCartesian: eta_eff has been initialied."
675     else:
676     eta_eff = self.__eta_eff
677     iter=0
678     converged=False
679     while not converged:
680     #
681     # intialize the solver
682     #
683     if mu==None:
684     stress0=Data()
685     else:
686     stress0=-(eta_eff/(dt*mu))*s_last
687     self.__solver.initialize(f=F,fixed_u_mask=mask_v,eta=eta_eff,surface_stress=f,stress=stress0)
688     #
689     # get a new velcocity and pressure:
690     #
691     if mask_v.isEmpty() or v_b.isEmpty():
692     v0=v
693     else:
694     v0=v_b*mask_v+v*(1.-mask_v)
695     v,p=self.__solver.solve(v0,p,show_details=False,
696     verbose=self.checkVerbose(),max_iter=inner_iter_max,usePCG=usePCG)
697     #
698     # update eta_eff:
699     #
700     D=self.__getDeviatoricStrain(v)
701     if mu==None:
702     gamma=util.sqrt(2.)*util.length(D)
703     else:
704     gamma=util.sqrt(2.)*util.length(D+s_last/(2*dt*mu))
705     eta_eff_old ,eta_eff=eta_eff, self.getEtaEff(gamma, pressure=p,dt=dt, eta0=eta_eff, iter_max=iter_max)
706     if self.checkVerbose(): print "IncompressibleIsotropicFlowCartesian: eta_eff has been updated."
707     #
708     # check the change on eta_eff:
709     #
710     diff=util.Lsup(eta_eff_old-eta_eff)
711     n=util.Lsup(eta_eff)
712 gross 2440 if self.checkVerbose(): print "IncompressibleIsotropicFlowCartesian: step %s: max. rel. change in eta_eff is %s."%(iter,diff/n)
713 gross 2415 converged = diff <= self.getTolerance()* n
714     iter+=1
715     if iter >= iter_max:
716     raise MaxIterReached,"maximum number of iteration steps on time step %e reached."%(self.getTime()+dt)
717     #
718     # finally we can update the return values:
719     #
720     self.setPressure(p)
721     self.setVelocity(v)
722     self.setDeviatoricStrain(D)
723     if mu==None:
724     stress=(2*eta_eff)*D
725     else:
726     stress=(2.*eta_eff)*(D+s_last/(2*dt*mu))
727     self.setDeviatoricStress(stress)
728     self.__eta_eff = eta_eff
729     self.setTime(self.getTime()+dt)
730     if self.checkVerbose(): print "IncompressibleIsotropicFlowCartesian: iteration on time step %s completed after %s steps."%(self.getTime(),iter)
731     return self.getVelocity(), self.getPressure()
732 gross 1639
733 gross 2440 def getCurrentEtaEff(self):
734     """
735     returns the effective viscosity
736     """
737     return self.__eta_eff
738    
739 gross 2415 def __getDeviatoricStrain(self, v):
740 gross 2100 """
741 gross 2415 Returns deviatoric strain of velocity v:
742 gross 2100 """
743 gross 2432 return util.deviatoric(util.symmetric(util.grad(v)))
744 caltinay 2169
745 gross 2440 def setFlowTolerance(self, tol=1.e-4):
746 gross 2100 """
747 gross 2415 Sets the relative tolerance for the flow solver. See L{StokesProblemCartesian.setTolerance} for details.
748 gross 1659
749 gross 2415 @param tol: desired relative tolerance for the flow solver
750     @type tol: positive C{float}
751 gross 2100 """
752 gross 2415 self.__solver.setTolerance(tol)
753     def getFlowTolerance(self):
754 gross 2100 """
755 gross 2415 Returns the relative tolerance for the flow solver
756 caltinay 2169
757 gross 2415 @return: tolerance of the flow solver
758     @rtype: C{float}
759 gross 2100 """
760 gross 2415 return self.__solver.getTolerance()
761 gross 2440 def setFlowSubTolerance(self, tol=1.e-8):
762 gross 2100 """
763 gross 2415 Sets the relative tolerance for the subsolver of the flow solver. See L{StokesProblemCartesian.setSubProblemTolerance} for details
764 gross 1639
765 gross 2415 @param tol: desired relative tolerance for the subsolver
766     @type tol: positive C{float}
767 gross 1639 """
768 gross 2415 self.__solver.setSubProblemTolerance(tol)
769     def getFlowSubTolerance(self):
770 gross 1639 """
771 gross 2415 Returns the relative tolerance for the subsolver of the flow solver
772 gross 1639
773 gross 2415 @return: tolerance of the flow subsolver
774     @rtype: C{float}
775 gross 2100 """
776 gross 2415 return self.__solver.getSubProblemTolerance()
777 gross 1639

  ViewVC Help
Powered by ViewVC 1.1.26