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

Diff of /trunk/escript/py_src/flows.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 2208 by gross, Mon Jan 12 06:37:07 2009 UTC revision 3432 by jfenwick, Fri Jan 7 01:32:07 2011 UTC
# Line 1  Line 1 
1    # -*- coding: utf-8 -*-
2  ########################################################  ########################################################
3  #  #
4  # Copyright (c) 2003-2008 by University of Queensland  # Copyright (c) 2003-2010 by University of Queensland
5  # Earth Systems Science Computational Center (ESSCC)  # Earth Systems Science Computational Center (ESSCC)
6  # http://www.uq.edu.au/esscc  # http://www.uq.edu.au/esscc
7  #  #
# Line 10  Line 11 
11  #  #
12  ########################################################  ########################################################
13    
14  __copyright__="""Copyright (c) 2003-2008 by University of Queensland  __copyright__="""Copyright (c) 2003-2010 by University of Queensland
15  Earth Systems Science Computational Center (ESSCC)  Earth Systems Science Computational Center (ESSCC)
16  http://www.uq.edu.au/esscc  http://www.uq.edu.au/esscc
17  Primary Business: Queensland, Australia"""  Primary Business: Queensland, Australia"""
18  __license__="""Licensed under the Open Software License version 3.0  __license__="""Licensed under the Open Software License version 3.0
19  http://www.opensource.org/licenses/osl-3.0.php"""  http://www.opensource.org/licenses/osl-3.0.php"""
20  __url__="http://www.uq.edu.au/esscc/escript-finley"  __url__="https://launchpad.net/escript-finley"
21    
22  """  """
23  Some models for flow  Some models for flow
24    
25  @var __author__: name of author  :var __author__: name of author
26  @var __copyright__: copyrights  :var __copyright__: copyrights
27  @var __license__: licence agreement  :var __license__: licence agreement
28  @var __url__: url entry point on documentation  :var __url__: url entry point on documentation
29  @var __version__: version  :var __version__: version
30  @var __date__: date of the version  :var __date__: date of the version
31  """  """
32    
33  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
34    
35  from escript import *  import escript
36  import util  import util
37  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions
38  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
39    
40  class DarcyFlow(object):  class DarcyFlow(object):
41       """
42       solves the problem
43      
44       *u_i+k_{ij}*p_{,j} = g_i*
45       *u_{i,i} = f*
46      
47       where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,
48      
49       :note: The problem is solved in a least squares formulation.
50       """
51      
52       def __init__(self, domain, useReduced=False, adaptSubTolerance=True, solveForFlux=False):
53          """
54          initializes the Darcy flux problem
55          :param domain: domain of the problem
56          :type domain: `Domain`
57          :param useReduced: uses reduced oreder on flux and pressure
58          :type useReduced: ``bool``
59          :param adaptSubTolerance: switches on automatic subtolerance selection
60          :type adaptSubTolerance: ``bool``
61          :param solveForFlux: if True the solver solves for the flux (do not use!)
62          :type solveForFlux: ``bool``  
63          """
64          self.domain=domain
65          self.solveForFlux=solveForFlux
66          self.useReduced=useReduced
67          self.__adaptSubTolerance=adaptSubTolerance
68          self.verbose=False
69          
70          self.__pde_k=LinearPDESystem(domain)
71          self.__pde_k.setSymmetryOn()
72          if self.useReduced: self.__pde_k.setReducedOrderOn()
73    
74          self.__pde_p=LinearSinglePDE(domain)
75          self.__pde_p.setSymmetryOn()
76          if self.useReduced: self.__pde_p.setReducedOrderOn()
77    
78          self.__pde_l=LinearSinglePDE(domain)   # this is here for getSolverOptionsWeighting
79          # self.__pde_l.setSymmetryOn()
80          # if self.useReduced: self.__pde_l.setReducedOrderOn()
81          self.setTolerance()
82          self.setAbsoluteTolerance()
83          self.__f=escript.Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))
84          self.__g=escript.Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
85          
86       def getSolverOptionsFlux(self):
87          """
88          Returns the solver options used to solve the flux problems
89          
90          *K^{-1} u=F*
91          
92          :return: `SolverOptions`
93          """
94          return self.__pde_k.getSolverOptions()
95          
96       def setSolverOptionsFlux(self, options=None):
97          """
98          Sets the solver options used to solve the flux problems
99          
100          *K^{-1}u=F*
101          
102          If ``options`` is not present, the options are reset to default
103          
104          :param options: `SolverOptions`
105          :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
106          """
107          return self.__pde_v.setSolverOptions(options)
108        
109       def getSolverOptionsPressure(self):
110          """
111          Returns the solver options used to solve the pressure problems
112          
113          *(Q^* K Q)p=-Q^*G*
114          
115          :return: `SolverOptions`
116          """
117          return self.__pde_p.getSolverOptions()
118          
119       def setSolverOptionsPressure(self, options=None):
120          """
121          Sets the solver options used to solve the pressure problems
122          
123          *(Q^* K Q)p=-Q^*G*
124          
125          If ``options`` is not present, the options are reset to default
126          
127          :param options: `SolverOptions`
128          :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
129          """
130          return self.__pde_p.setSolverOptions(options)
131    
132       def getSolverOptionsWeighting(self):
133          """
134          Returns the solver options used to solve the pressure problems
135    
136          *(D K D^*)p=-D F*
137    
138          :return: `SolverOptions`
139          """
140          return self.__pde_l.getSolverOptions()
141    
142       def setSolverOptionsWeighting(self, options=None):
143          """
144          Sets the solver options used to solve the pressure problems
145    
146          *(D K D^*)p=-D F*
147    
148          If ``options`` is not present, the options are reset to default
149    
150          :param options: `SolverOptions`
151          :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
152          """
153          return self.__pde_l.setSolverOptions(options)
154    
155    
156       def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
157          """
158          assigns values to model parameters
159    
160          :param f: volumetic sources/sinks
161          :type f: scalar value on the domain (e.g. `Data`)
162          :param g: flux sources/sinks
163          :type g: vector values on the domain (e.g. `Data`)
164          :param location_of_fixed_pressure: mask for locations where pressure is fixed
165          :type location_of_fixed_pressure: scalar value on the domain (e.g. `Data`)
166          :param location_of_fixed_flux:  mask for locations where flux is fixed.
167          :type location_of_fixed_flux: vector values on the domain (e.g. `Data`)
168          :param permeability: permeability tensor. If scalar ``s`` is given the tensor with ``s`` on the main diagonal is used.
169          :type permeability: scalar or tensor values on the domain (e.g. `Data`)
170    
171          :note: the values of parameters which are not set by calling ``setValue`` are not altered.
172          :note: at any point on the boundary of the domain the pressure (``location_of_fixed_pressure`` >0)
173          or the normal component of the flux (``location_of_fixed_flux[i]>0``) if direction of the normal is along the *x_i* axis.
174    
175          """
176          if f !=None:
177         f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("X"))
178         if f.isEmpty():
179            f=escript.Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))
180         else:
181            if f.getRank()>0: raise ValueError,"illegal rank of f."
182            self.__f=f
183          if g !=None:
184         g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
185         if g.isEmpty():
186            g=escript.Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
187         else:
188            if not g.getShape()==(self.domain.getDim(),):
189               raise ValueError,"illegal shape of g"
190            self.__g=g
191          if location_of_fixed_pressure!=None:
192               self.__pde_p.setValue(q=location_of_fixed_pressure)
193               #self.__pde_l.setValue(q=location_of_fixed_pressure)
194          if location_of_fixed_flux!=None:
195               self.__pde_k.setValue(q=location_of_fixed_flux)
196                
197          if permeability!=None:
198         perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
199         if perm.getRank()==0:
200            perm_inv=(1./perm)*util.kronecker(self.domain.getDim())
201            perm=perm*util.kronecker(self.domain.getDim())
202         elif perm.getRank()==2:
203            perm_inv=util.inverse(perm)
204         else:
205            raise ValueError,"illegal rank of permeability."
206    
207         self.__permeability=perm
208         self.__permeability_inv=perm_inv
209         self.__l =(util.longestEdge(self.domain)**2*util.length(self.__permeability_inv))/10
210         #self.__l =(self.domain.getSize()**2*util.length(self.__permeability_inv))/10
211         if  self.solveForFlux:
212            self.__pde_k.setValue(D=self.__permeability_inv)
213         else:
214            self.__pde_k.setValue(D=self.__permeability_inv, A=self.__l*util.outer(util.kronecker(self.domain),util.kronecker(self.domain)))
215         self.__pde_p.setValue(A=self.__permeability)
216         #self.__pde_l.setValue(D=1/self.__l)
217             #self.__pde_l.setValue(A=self.__permeability)
218    
219       def __applWeight(self, v, f=None):
220          # solves L p = f-Dv with p = 0
221          if self.getSolverOptionsWeighting().isVerbose() or self.verbose: print "DarcyFlux: Applying weighting operator"
222          if f == None:
223         return -util.div(v)*self.__l
224          else:
225         return (f-util.div(v))*self.__l
226          # if f == None:
227          #      self.__pde_l.setValue(Y=-util.div(v))  
228          # else:
229          #      return (f-util.div(v))/self.__l
230          # return self.__pde_l.getSolution()
231          
232       def __getPressure(self, v, p0, g=None):
233          # solves (G*KG)p = G^(g-v) with p = p0 where location_of_fixed_pressure>0
234          if self.getSolverOptionsPressure().isVerbose() or self.verbose: print "DarcyFlux: Pressure update"
235          if g == None:
236         self.__pde_p.setValue(X=-v, r=p0)
237          else:
238         self.__pde_p.setValue(X=g-v, r=p0)      
239          p=self.__pde_p.getSolution()
240          return p
241    
242       def __Aprod_v(self,dv):
243          # calculates: (a,b,c) = (K^{-1}(dv + KG * dp), L^{-1}Ddv, dp)  with (G*KG)dp = - G^*dv  
244          dp=self.__getPressure(dv, p0=escript.Data()) # dp = (G*KG)^{-1} (0-G^*dv)
245          a=util.tensor_mult(self.__permeability_inv,dv)+util.grad(dp) # a= K^{-1}u+G*dp
246          b= - self.__applWeight(dv) # b = - (D K D^*)^{-1} (0-Dv)
247          return ArithmeticTuple(a,b,-dp)
248    
249       def __Msolve_PCG_v(self,r):
250          # K^{-1} u = r[0] + D^*r[1] = K^{-1}(dv + KG * dp) + D^*L^{-1}Ddv
251          if self.getSolverOptionsFlux().isVerbose() or self.verbose: print "DarcyFlux: Applying preconditioner"
252          self.__pde_k.setValue(X=r[1]*util.kronecker(self.domain), Y=r[0], r=escript.Data())
253          # self.__pde_p.getOperator().saveMM("prec.mm")
254          return self.__pde_k.getSolution()
255    
256       def __inner_PCG_v(self,v,r):
257          return util.integrate(util.inner(v,r[0])+util.div(v)*r[1])
258          
259       def __Aprod_p(self,dp):
260          if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
261          Gdp=util.grad(dp)
262          self.__pde_k.setValue(Y=-Gdp,X=escript.Data(), r=escript.Data())
263          du=self.__pde_k.getSolution()
264          # self.__pde_v.getOperator().saveMM("proj.mm")
265          return ArithmeticTuple(util.tensor_mult(self.__permeability,Gdp),-du)
266    
267       def __getFlux(self,p, v0, f=None, g=None):
268          # solves (K^{-1}+D^*L^{-1} D) v = D^*L^{-1}f + K^{-1}g - Gp
269          if f!=None:
270         self.__pde_k.setValue(X=self.__applWeight(v0*0,self.__f)*util.kronecker(self.domain))
271          self.__pde_k.setValue(r=v0)
272          g2=util.tensor_mult(self.__permeability_inv,g)
273          if p == None:
274         self.__pde_k.setValue(Y=g2)
275          else:
276         self.__pde_k.setValue(Y=g2-util.grad(p))
277          return self.__pde_k.getSolution()  
278          
279          #v=self.__getFlux(p, u0, f=self.__f, g=g2)      
280       def __Msolve_PCG_p(self,r):
281          if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"
282          self.__pde_p.setValue(X=r[0]-r[1], Y=escript.Data(), r=escript.Data(), y=escript.Data())
283          # self.__pde_p.getOperator().saveMM("prec.mm")
284          return self.__pde_p.getSolution()
285            
286       def __inner_PCG_p(self,p,r):
287           return util.integrate(util.inner(util.grad(p), r[0]-r[1]))
288    
289       def __L2(self,v):
290          return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))
291    
292       def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):
293          """
294          solves the problem.
295          
296          The iteration is terminated if the residual norm is less then self.getTolerance().
297    
298          :param u0: initial guess for the flux. At locations in the domain marked by ``location_of_fixed_flux`` the value of ``u0`` is kept unchanged.
299          :type u0: vector value on the domain (e.g. `Data`).
300          :param p0: initial guess for the pressure. At locations in the domain marked by ``location_of_fixed_pressure`` the value of ``p0`` is kept unchanged.
301          :type p0: scalar value on the domain (e.g. `Data`).
302          :param verbose: if set some information on iteration progress are printed
303          :type verbose: ``bool``
304          :return: flux and pressure
305          :rtype: ``tuple`` of `Data`.
306    
307          :note: The problem is solved as a least squares form
308          *(K^[-1]+D^* (DKD^*)^[-1] D)u+G p=D^* (DKD^*)^[-1] f + K^[-1]g*
309          *G^*u+*G^* K Gp=G^*g*
310          where *D* is the *div* operator and *(Gp)_i=p_{,i}* for the permeability *K=k_{ij}*.
311          """
312          self.verbose=verbose
313          rtol=self.getTolerance()
314          atol=self.getAbsoluteTolerance()
315          self.setSubProblemTolerance()
316          num_corrections=0
317          converged=False
318          norm_r=None
319          
320          # Eliminate the hydrostatic pressure:
321          if self.verbose: print "DarcyFlux: calculate hydrostatic pressure component."
322          self.__pde_p.setValue(X=self.__g, r=p0, y=-util.inner(self.domain.getNormal(),u0))        
323          p0=self.__pde_p.getSolution()
324          g2=self.__g - util.tensor_mult(self.__permeability, util.grad(p0))
325          norm_g2=util.integrate(util.inner(g2,util.tensor_mult(self.__permeability_inv,g2)))**0.5    
326    
327          p=p0*0
328          if self.solveForFlux:
329         v=u0.copy()
330          else:
331         v=self.__getFlux(p, u0, f=self.__f, g=g2)
332    
333          while not converged and norm_g2 > 0:
334         Gp=util.grad(p)
335         KGp=util.tensor_mult(self.__permeability,Gp)
336         if self.verbose:
337            def_p=g2-(v+KGp)
338            def_v=self.__f-util.div(v)
339            print "DarcyFlux: L2: g-v-K*grad(p) = %e (v = %e)."%(self.__L2(def_p),self.__L2(v))
340            print "DarcyFlux: L2: f-div(v) = %e (grad(v) = %e)."%(self.__L2(def_v),self.__L2(util.grad(v)))
341            print "DarcyFlux: K^{-1}-norm of v = %e."%util.integrate(util.inner(v,util.tensor_mult(self.__permeability_inv,v)))**0.5
342            print "DarcyFlux: K^{-1}-norm of g2 = %e."%norm_g2
343            print "DarcyFlux: K-norm of grad(dp) = %e."%util.integrate(util.inner(Gp,KGp))**0.5
344         ATOL=atol+rtol*norm_g2
345         if self.verbose: print "DarcyFlux: absolute tolerance ATOL = %e."%(ATOL,)
346         if norm_r == None or norm_r>ATOL:
347            if num_corrections>max_num_corrections:
348               raise ValueError,"maximum number of correction steps reached."
349          
350            if self.solveForFlux:
351               # initial residual is r=K^{-1}*(g-v-K*Gp)+D^*L^{-1}(f-Du)
352               v,r, norm_r=PCG(ArithmeticTuple(util.tensor_mult(self.__permeability_inv,g2-v)-Gp,self.__applWeight(v,self.__f),p),
353                   self.__Aprod_v,
354                   v,
355                   self.__Msolve_PCG_v,
356                   self.__inner_PCG_v,
357                   atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
358               p=r[2]
359            else:
360               # initial residual is r=G^*(g2-KGp - v)
361               p,r, norm_r=PCG(ArithmeticTuple(g2-KGp,v),
362                     self.__Aprod_p,
363                     p,
364                     self.__Msolve_PCG_p,
365                     self.__inner_PCG_p,
366                     atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
367               v=r[1]
368            if self.verbose: print "DarcyFlux: residual norm = %e."%norm_r
369            num_corrections+=1
370         else:
371            if self.verbose: print "DarcyFlux: stopping criterium reached."
372            converged=True
373          return v,p+p0
374       def setTolerance(self,rtol=1e-4):
375          """
376          sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if
377    
378          *|g-v-K gard(p)|_PCG <= atol + rtol * |K^{1/2}g2|_0*
379          
380          where ``atol`` is an absolut tolerance (see `setAbsoluteTolerance`).
381          
382          :param rtol: relative tolerance for the pressure
383          :type rtol: non-negative ``float``
384          """
385          if rtol<0:
386         raise ValueError,"Relative tolerance needs to be non-negative."
387          self.__rtol=rtol
388       def getTolerance(self):
389          """
390          returns the relative tolerance
391          :return: current relative tolerance
392          :rtype: ``float``
393          """
394          return self.__rtol
395    
396       def setAbsoluteTolerance(self,atol=0.):
397          """
398          sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if
399          
400          *|g-v-K gard(p)|_PCG <= atol + rtol * |K^{1/2}g2|_0*
401    
402    
403          where ``rtol`` is an absolut tolerance (see `setTolerance`), *|f|^2 = integrate(length(f)^2)* and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.
404    
405          :param atol: absolute tolerance for the pressure
406          :type atol: non-negative ``float``
407          """
408          if atol<0:
409         raise ValueError,"Absolute tolerance needs to be non-negative."
410          self.__atol=atol
411       def getAbsoluteTolerance(self):
412          """
413          returns the absolute tolerance
414          :return: current absolute tolerance
415          :rtype: ``float``
416          """
417          return self.__atol
418       def getSubProblemTolerance(self):
419          """
420          Returns a suitable subtolerance
421          :type: ``float``
422          """
423          return max(util.EPSILON**(0.5),self.getTolerance()**2)
424    
425       def setSubProblemTolerance(self):
426          """
427          Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected.
428          """
429          if self.__adaptSubTolerance:
430         sub_tol=self.getSubProblemTolerance()
431         self.getSolverOptionsFlux().setTolerance(sub_tol)
432         self.getSolverOptionsFlux().setAbsoluteTolerance(0.)
433         self.getSolverOptionsPressure().setTolerance(sub_tol)
434         self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
435         self.getSolverOptionsWeighting().setTolerance(sub_tol)
436         self.getSolverOptionsWeighting().setAbsoluteTolerance(0.)
437         if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol
438    
439    
440    class DarcyFlowOld(object):
441      """      """
442      solves the problem      solves the problem
443    
444      M{u_i+k_{ij}*p_{,j} = g_i}      *u_i+k_{ij}*p_{,j} = g_i*
445      M{u_{i,i} = f}      *u_{i,i} = f*
446    
447      where M{p} represents the pressure and M{u} the Darcy flux. M{k} represents the permeability,      where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,
448    
449      @note: The problem is solved in a least squares formulation.      :note: The problem is solved in a least squares formulation.
450      """      """
451    
452      def __init__(self, domain,useReduced=False):      def __init__(self, domain, weight=None, useReduced=False, adaptSubTolerance=True):
453          """          """
454          initializes the Darcy flux problem          initializes the Darcy flux problem
455          @param domain: domain of the problem          :param domain: domain of the problem
456          @type domain: L{Domain}          :type domain: `Domain`
457        :param useReduced: uses reduced oreder on flux and pressure
458        :type useReduced: ``bool``
459        :param adaptSubTolerance: switches on automatic subtolerance selection
460        :type adaptSubTolerance: ``bool``  
461          """          """
462          self.domain=domain          self.domain=domain
463            if weight == None:
464               s=self.domain.getSize()
465               self.__l=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2
466               # self.__l=(3.*util.longestEdge(self.domain))**2
467               #self.__l=(0.1*util.longestEdge(self.domain)*s/util.sup(s))**2
468            else:
469               self.__l=weight
470          self.__pde_v=LinearPDESystem(domain)          self.__pde_v=LinearPDESystem(domain)
471          if useReduced: self.__pde_v.setReducedOrderOn()          if useReduced: self.__pde_v.setReducedOrderOn()
472          self.__pde_v.setSymmetryOn()          self.__pde_v.setSymmetryOn()
473          self.__pde_v.setValue(D=util.kronecker(domain), A=util.outer(util.kronecker(domain),util.kronecker(domain)))          self.__pde_v.setValue(D=util.kronecker(domain), A=self.__l*util.outer(util.kronecker(domain),util.kronecker(domain)))
474          self.__pde_p=LinearSinglePDE(domain)          self.__pde_p=LinearSinglePDE(domain)
475          self.__pde_p.setSymmetryOn()          self.__pde_p.setSymmetryOn()
476          if useReduced: self.__pde_p.setReducedOrderOn()          if useReduced: self.__pde_p.setReducedOrderOn()
477          self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))          self.__f=escript.Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
478          self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))          self.__g=escript.Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
479          self.__ATOL= None          self.setTolerance()
480            self.setAbsoluteTolerance()
481        self.__adaptSubTolerance=adaptSubTolerance
482        self.verbose=False
483        def getSolverOptionsFlux(self):
484        """
485        Returns the solver options used to solve the flux problems
486        
487        *(I+D^*D)u=F*
488        
489        :return: `SolverOptions`
490        """
491        return self.__pde_v.getSolverOptions()
492        def setSolverOptionsFlux(self, options=None):
493        """
494        Sets the solver options used to solve the flux problems
495        
496        *(I+D^*D)u=F*
497        
498        If ``options`` is not present, the options are reset to default
499        :param options: `SolverOptions`
500        :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
501        """
502        return self.__pde_v.setSolverOptions(options)
503        def getSolverOptionsPressure(self):
504        """
505        Returns the solver options used to solve the pressure problems
506        
507        *(Q^*Q)p=Q^*G*
508        
509        :return: `SolverOptions`
510        """
511        return self.__pde_p.getSolverOptions()
512        def setSolverOptionsPressure(self, options=None):
513        """
514        Sets the solver options used to solve the pressure problems
515        
516        *(Q^*Q)p=Q^*G*
517        
518        If ``options`` is not present, the options are reset to default
519        :param options: `SolverOptions`
520        :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
521        """
522        return self.__pde_p.setSolverOptions(options)
523    
524      def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):      def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
525          """          """
526          assigns values to model parameters          assigns values to model parameters
527    
528          @param f: volumetic sources/sinks          :param f: volumetic sources/sinks
529          @type f: scalar value on the domain (e.g. L{Data})          :type f: scalar value on the domain (e.g. `Data`)
530          @param g: flux sources/sinks          :param g: flux sources/sinks
531          @type g: vector values on the domain (e.g. L{Data})          :type g: vector values on the domain (e.g. `Data`)
532          @param location_of_fixed_pressure: mask for locations where pressure is fixed          :param location_of_fixed_pressure: mask for locations where pressure is fixed
533          @type location_of_fixed_pressure: scalar value on the domain (e.g. L{Data})          :type location_of_fixed_pressure: scalar value on the domain (e.g. `Data`)
534          @param location_of_fixed_flux:  mask for locations where flux is fixed.          :param location_of_fixed_flux:  mask for locations where flux is fixed.
535          @type location_of_fixed_flux: vector values on the domain (e.g. L{Data})          :type location_of_fixed_flux: vector values on the domain (e.g. `Data`)
536          @param permeability: permeability tensor. If scalar C{s} is given the tensor with          :param permeability: permeability tensor. If scalar ``s`` is given the tensor with
537                               C{s} on the main diagonal is used. If vector C{v} is given the tensor with                               ``s`` on the main diagonal is used. If vector ``v`` is given the tensor with
538                               C{v} on the main diagonal is used.                               ``v`` on the main diagonal is used.
539          @type permeability: scalar, vector or tensor values on the domain (e.g. L{Data})          :type permeability: scalar, vector or tensor values on the domain (e.g. `Data`)
540    
541          @note: the values of parameters which are not set by calling C{setValue} are not altered.          :note: the values of parameters which are not set by calling ``setValue`` are not altered.
542          @note: at any point on the boundary of the domain the pressure (C{location_of_fixed_pressure} >0)          :note: at any point on the boundary of the domain the pressure (``location_of_fixed_pressure`` >0)
543                 or the normal component of the flux (C{location_of_fixed_flux[i]>0} if direction of the normal                 or the normal component of the flux (``location_of_fixed_flux[i]>0`` if direction of the normal
544                 is along the M{x_i} axis.                 is along the *x_i* axis.
545          """          """
546          if f !=None:          if f !=None:
547             f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))             f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
548             if f.isEmpty():             if f.isEmpty():
549                 f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))                 f=escript.Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
550             else:             else:
551                 if f.getRank()>0: raise ValueError,"illegal rank of f."                 if f.getRank()>0: raise ValueError,"illegal rank of f."
552             self.f=f             self.__f=f
553          if g !=None:            if g !=None:
554             g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))             g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
555             if g.isEmpty():             if g.isEmpty():
556               g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))               g=escript.Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
557             else:             else:
558               if not g.getShape()==(self.domain.getDim(),):               if not g.getShape()==(self.domain.getDim(),):
559                 raise ValueError,"illegal shape of g"                 raise ValueError,"illegal shape of g"
# Line 121  class DarcyFlow(object): Line 576  class DarcyFlow(object):
576             self.__permeability=perm             self.__permeability=perm
577             self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))             self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))
578    
579        def setTolerance(self,rtol=1e-4):
     def getFlux(self,p=None, fixed_flux=Data(),tol=1.e-8, show_details=False):  
         """  
         returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux}  
         on locations where C{location_of_fixed_flux} is positive (see L{setValue}).  
         Note that C{g} and C{f} are used, see L{setValue}.  
           
         @param p: pressure.  
         @type p: scalar value on the domain (e.g. L{Data}).  
         @param fixed_flux: flux on the locations of the domain marked be C{location_of_fixed_flux}.  
         @type fixed_flux: vector values on the domain (e.g. L{Data}).  
         @param tol: relative tolerance to be used.  
         @type tol: positive C{float}.  
         @return: flux  
         @rtype: L{Data}  
         @note: the method uses the least squares solution M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}}  
                for the permeability M{k_{ij}}  
580          """          """
581          self.__pde_v.setTolerance(tol)          sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if
         g=self.__g  
         f=self.__f  
         self.__pde_v.setValue(X=f*util.kronecker(self.domain), r=fixed_flux)  
         if p == None:  
            self.__pde_v.setValue(Y=g)  
         else:  
            self.__pde_v.setValue(Y=g-util.tensor_mult(self.__permeability,util.grad(p)))  
         return self.__pde_v.getSolution(verbose=show_details)  
   
     def getPressure(self,v=None, fixed_pressure=Data(),tol=1.e-8, show_details=False):  
         """  
         returns the pressure for a given flux C{v} where the pressure is equal to C{fixed_pressure}  
         on locations where C{location_of_fixed_pressure} is positive (see L{setValue}).  
         Note that C{g} is used, see L{setValue}.  
           
         @param v: flux.  
         @type v: vector-valued on the domain (e.g. L{Data}).  
         @param fixed_pressure: pressure on the locations of the domain marked be C{location_of_fixed_pressure}.  
         @type fixed_pressure: vector values on the domain (e.g. L{Data}).  
         @param tol: relative tolerance to be used.  
         @type tol: positive C{float}.  
         @return: pressure  
         @rtype: L{Data}  
         @note: the method uses the least squares solution M{p=(Q^*Q)^{-1}Q^*(g-u)} where and M{(Qp)_i=k_{ij}p_{,j}}  
                for the permeability M{k_{ij}}  
         """  
         self.__pde_v.setTolerance(tol)  
         g=self.__g  
         self.__pde_p.setValue(r=fixed_pressure)  
         if v == None:  
            self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,g-v))  
         else:  
            self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,g))  
         return self.__pde_p.getSolution(verbose=show_details)  
582    
583      def setTolerance(self,atol=0,rtol=1e-8,p_ref=None,v_ref=None):          *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
         """  
         set the tolerance C{ATOL} used to terminate the solution process. It is used  
584    
585          M{ATOL = atol + rtol * max( |g-v_ref|, |Qp_ref| )}          where ``atol`` is an absolut tolerance (see `setAbsoluteTolerance`), *|f|^2 = integrate(length(f)^2)* and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.
586    
587          where M{|f|^2 = integrate(length(f)^2)} and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}. If C{v_ref} or C{p_ref} is not present zero is assumed.          :param rtol: relative tolerance for the pressure
588            :type rtol: non-negative ``float``
         The iteration is terminated if for the current approximation C{p}, flux C{v=(I+D^*D)^{-1}(D^*f-g-Qp)} and their residual  
   
         M{r=Q^*(g-Qp-v)}  
   
         the condition  
   
         M{<(Q^*Q)^{-1} r,r> <= ATOL}  
   
         holds. M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}  
   
         @param atol: absolute tolerance for the pressure  
         @type atol: non-negative C{float}  
         @param rtol: relative tolerance for the pressure  
         @type rtol: non-negative C{float}  
         @param p_ref: reference pressure. If not present zero is used. You may use physical arguments to set a resonable value for C{p_ref}, use the  
         L{getPressure} method or use  the value from a previous time step.  
         @type p_ref: scalar value on the domain (e.g. L{Data}).  
         @param v_ref: reference velocity.  If not present zero is used. You may use physical arguments to set a resonable value for C{v_ref}, use the  
         L{getFlux} method or use  the value from a previous time step.  
         @type v_ref: vector-valued on the domain (e.g. L{Data}).  
         @return: used absolute tolerance.  
         @rtype: positive C{float}  
589          """          """
590          g=self.__g          if rtol<0:
591          if not v_ref == None:              raise ValueError,"Relative tolerance needs to be non-negative."
592             f1=util.integrate(util.length(util.interpolate(g-v_ref,Function(self.domain)))**2)          self.__rtol=rtol
         else:  
            f1=util.integrate(util.length(util.interpolate(g))**2)  
         if not p_ref == None:  
            f2=util.integrate(util.length(util.tensor_mult(self.__permeability,util.grad(p_ref)))**2)  
         else:  
            f2=0  
         self.__ATOL= atol + rtol * util.sqrt(max(f1,f2))  
         if self.__ATOL<=0:  
            raise ValueError,"Positive tolerance (=%e) is expected."%self.__ATOL  
         return self.__ATOL  
           
593      def getTolerance(self):      def getTolerance(self):
594          """          """
595          returns the current tolerance.          returns the relative tolerance
596      
597          @return: used absolute tolerance.          :return: current relative tolerance
598          @rtype: positive C{float}          :rtype: ``float``
599          """          """
600          if self.__ATOL==None:          return self.__rtol
            raise ValueError,"no tolerance is defined."  
         return self.__ATOL  
601    
602      def solve(self,u0,p0, max_iter=100, verbose=False, show_details=False, sub_rtol=1.e-8):      def setAbsoluteTolerance(self,atol=0.):
603           """          """
604            sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if
605    
606            *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
607    
608            where ``rtol`` is an absolut tolerance (see `setTolerance`), *|f|^2 = integrate(length(f)^2)* and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.
609    
610            :param atol: absolute tolerance for the pressure
611            :type atol: non-negative ``float``
612            """
613            if atol<0:
614                raise ValueError,"Absolute tolerance needs to be non-negative."
615            self.__atol=atol
616        def getAbsoluteTolerance(self):
617           """
618           returns the absolute tolerance
619          
620           :return: current absolute tolerance
621           :rtype: ``float``
622           """
623           return self.__atol
624        def getSubProblemTolerance(self):
625        """
626        Returns a suitable subtolerance
627        @type: ``float``
628        """
629        return max(util.EPSILON**(0.75),self.getTolerance()**2)
630        def setSubProblemTolerance(self):
631             """
632             Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected.
633             """
634         if self.__adaptSubTolerance:
635             sub_tol=self.getSubProblemTolerance()
636                 self.getSolverOptionsFlux().setTolerance(sub_tol)
637             self.getSolverOptionsFlux().setAbsoluteTolerance(0.)
638             self.getSolverOptionsPressure().setTolerance(sub_tol)
639             self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
640             if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol
641    
642        def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):
643             """
644           solves the problem.           solves the problem.
645    
646           The iteration is terminated if the residual norm is less then self.getTolerance().           The iteration is terminated if the residual norm is less then self.getTolerance().
647    
648           @param u0: initial guess for the flux. At locations in the domain marked by C{location_of_fixed_flux} the value of C{u0} is kept unchanged.           :param u0: initial guess for the flux. At locations in the domain marked by ``location_of_fixed_flux`` the value of ``u0`` is kept unchanged.
649           @type u0: vector value on the domain (e.g. L{Data}).           :type u0: vector value on the domain (e.g. `Data`).
650           @param p0: initial guess for the pressure. At locations in the domain marked by C{location_of_fixed_pressure} the value of C{p0} is kept unchanged.           :param p0: initial guess for the pressure. At locations in the domain marked by ``location_of_fixed_pressure`` the value of ``p0`` is kept unchanged.
651           @type p0: scalar value on the domain (e.g. L{Data}).           :type p0: scalar value on the domain (e.g. `Data`).
652           @param sub_rtol: tolerance to be used in the sub iteration. It is recommended that M{sub_rtol<rtol*5.e-3}           :param verbose: if set some information on iteration progress are printed
653           @type sub_rtol: positive-negative C{float}           :type verbose: ``bool``
654           @param verbose: if set some information on iteration progress are printed           :return: flux and pressure
655           @type verbose: C{bool}           :rtype: ``tuple`` of `Data`.
656           @param show_details:  if set information on the subiteration process are printed.  
657           @type show_details: C{bool}           :note: The problem is solved as a least squares form
          @return: flux and pressure  
          @rtype: C{tuple} of L{Data}.  
   
          @note: The problem is solved as a least squares form  
658    
659           M{(I+D^*D)u+Qp=D^*f+g}           *(I+D^*D)u+Qp=D^*f+g*
660           M{Q^*u+Q^*Qp=Q^*g}           *Q^*u+Q^*Qp=Q^*g*
661    
662           where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.           where *D* is the *div* operator and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.
663           We eliminate the flux form the problem by setting           We eliminate the flux form the problem by setting
664    
665           M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} with u=u0 on location_of_fixed_flux           *u=(I+D^*D)^{-1}(D^*f-g-Qp)* with u=u0 on location_of_fixed_flux
666    
667           form the first equation. Inserted into the second equation we get           form the first equation. Inserted into the second equation we get
668    
669           M{Q^*(I-(I+D^*D)^{-1})Qp= Q^*(g-(I+D^*D)^{-1}(D^*f+g))} with p=p0  on location_of_fixed_pressure           *Q^*(I-(I+D^*D)^{-1})Qp= Q^*(g-(I+D^*D)^{-1}(D^*f+g))* with p=p0  on location_of_fixed_pressure
670            
671           which is solved using the PCG method (precondition is M{Q^*Q}). In each iteration step           which is solved using the PCG method (precondition is *Q^*Q*). In each iteration step
672           PDEs with operator M{I+D^*D} and with M{Q^*Q} needs to be solved using a sub iteration scheme.           PDEs with operator *I+D^*D* and with *Q^*Q* needs to be solved using a sub iteration scheme.
673           """           """
674           self.verbose=verbose           self.verbose=verbose
675           self.show_details= show_details and self.verbose           rtol=self.getTolerance()
676           self.__pde_v.setTolerance(sub_rtol)           atol=self.getAbsoluteTolerance()
677           self.__pde_p.setTolerance(sub_rtol)       self.setSubProblemTolerance()
678           ATOL=self.getTolerance()           num_corrections=0
679           if self.verbose: print "DarcyFlux: absolute tolerance = %e"%ATOL           converged=False
680           #########################################################################################################################           p=p0
681           #           norm_r=None
682           #   we solve:           while not converged:
683           #                   v=self.getFlux(p, fixed_flux=u0)
684           #      Q^*(I-(I+D^*D)^{-1})Q dp =  Q^* (g-u0-Qp0 - (I+D^*D)^{-1} ( D^*(f-Du0)+g-u0-Qp0) )                 Qp=self.__Q(p)
685           #                 norm_v=self.__L2(v)
686           #   residual is                 norm_Qp=self.__L2(Qp)
687           #                 if norm_v == 0.:
688           #    r=  Q^* (g-u0-Qp0 - (I+D^*D)^{-1} ( D^*(f-Du0)+g-u0-Qp0) - Q dp +(I+D^*D)^{-1})Q dp ) = Q^* (g - Qp - v)                    if norm_Qp == 0.:
689           #                       return v,p
690           #        with v = (I+D^*D)^{-1} (D^*f+g-Qp) including BC                    else:
691           #                      fac=norm_Qp
692           #    we use (g - Qp, v) to represent the residual. not that                 else:
693           #                    if norm_Qp == 0.:
694           #    dr(dp)=( -Q(dp), dv) with dv = - (I+D^*D)^{-1} Q(dp)                      fac=norm_v
695           #                    else:
696           #   while the initial residual is                      fac=2./(1./norm_v+1./norm_Qp)
697           #                 ATOL=(atol+rtol*fac)
698           #      r0=( g - Qp0, v00) with v00=(I+D^*D)^{-1} (D^*f+g-Qp0) including BC                 if self.verbose:
699           #                        print "DarcyFlux: L2 norm of v = %e."%norm_v
700           d0=self.__g-util.tensor_mult(self.__permeability,util.grad(p0))                      print "DarcyFlux: L2 norm of k.util.grad(p) = %e."%norm_Qp
701           self.__pde_v.setValue(Y=d0, X=self.__f*util.kronecker(self.domain), r=u0)                      print "DarcyFlux: L2 defect u = %e."%(util.integrate(util.length(self.__g-util.interpolate(v,escript.Function(self.domain))-Qp)**2)**(0.5),)
702           v00=self.__pde_v.getSolution(verbose=show_details)                      print "DarcyFlux: L2 defect div(v) = %e."%(util.integrate((self.__f-util.div(v))**2)**(0.5),)
703           if self.verbose: print "DarcyFlux: range of initial flux = ",util.inf(v00), util.sup(v00)                      print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
704           self.__pde_v.setValue(r=Data())                 if norm_r == None or norm_r>ATOL:
705           # start CG                     if num_corrections>max_num_corrections:
706           r=ArithmeticTuple(d0, v00)                           raise ValueError,"maximum number of correction steps reached."
707           p,r=PCG(r,self.__Aprod_PCG,p0,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)                     p,r, norm_r=PCG(self.__g-util.interpolate(v,escript.Function(self.domain))-Qp,self.__Aprod,p,self.__Msolve_PCG,self.__inner_PCG,atol=0.5*ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
708           return r[1],p                     num_corrections+=1
709                   else:
710      def __Aprod_PCG(self,dp):                     converged=True
711            if self.show_details: print "DarcyFlux: Applying operator"           return v,p
712            #  -dr(dp) = (Qdp,du) where du = (I+D^*D)^{-1} (Qdp)      def __L2(self,v):
713            mQdp=util.tensor_mult(self.__permeability,util.grad(dp))           return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))
714            self.__pde_v.setValue(Y=mQdp,X=Data(), r=Data())  
715            du=self.__pde_v.getSolution(verbose=self.show_details)      def __Q(self,p):
716            return ArithmeticTuple(mQdp,du)            return util.tensor_mult(self.__permeability,util.grad(p))
717    
718        def __Aprod(self,dp):
719              if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
720              Qdp=self.__Q(dp)
721              self.__pde_v.setValue(Y=-Qdp,X=escript.Data(), r=escript.Data())
722              du=self.__pde_v.getSolution()
723              # self.__pde_v.getOperator().saveMM("proj.mm")
724              return Qdp+du
725        def __inner_GMRES(self,r,s):
726             return util.integrate(util.inner(r,s))
727    
728      def __inner_PCG(self,p,r):      def __inner_PCG(self,p,r):
729           a=util.tensor_mult(self.__permeability,util.grad(p))           return util.integrate(util.inner(self.__Q(p), r))
          f0=util.integrate(util.inner(a,r[0]))  
          f1=util.integrate(util.inner(a,r[1]))  
          # print "__inner_PCG:",f0,f1,"->",f0-f1  
          return f0-f1  
730    
731      def __Msolve_PCG(self,r):      def __Msolve_PCG(self,r):
732            if self.show_details: print "DarcyFlux: Applying preconditioner"        if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"
733            self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r[0]-r[1]), r=Data())            self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=escript.Data(), r=escript.Data())
734            return self.__pde_p.getSolution(verbose=self.show_details)            # self.__pde_p.getOperator().saveMM("prec.mm")
735              return self.__pde_p.getSolution()
736    
737        def getFlux(self,p=None, fixed_flux=escript.Data()):
738            """
739            returns the flux for a given pressure ``p`` where the flux is equal to ``fixed_flux``
740            on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
741            Note that ``g`` and ``f`` are used, see `setValue`.
742    
743            :param p: pressure.
744            :type p: scalar value on the domain (e.g. `Data`).
745            :param fixed_flux: flux on the locations of the domain marked be ``location_of_fixed_flux``.
746            :type fixed_flux: vector values on the domain (e.g. `Data`).
747            :return: flux
748            :rtype: `Data`
749            :note: the method uses the least squares solution *u=(I+D^*D)^{-1}(D^*f-g-Qp)* where *D* is the *div* operator and *(Qp)_i=k_{ij}p_{,j}*
750                   for the permeability *k_{ij}*
751            """
752        self.setSubProblemTolerance()
753            g=self.__g
754            f=self.__f
755            self.__pde_v.setValue(X=self.__l*f*util.kronecker(self.domain), r=fixed_flux)
756            if p == None:
757               self.__pde_v.setValue(Y=g)
758            else:
759               self.__pde_v.setValue(Y=g-self.__Q(p))
760            return self.__pde_v.getSolution()
761    
762  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
763        """       """
764        solves       solves
765    
766            -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}            -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
767                  u_{i,i}=0                  u_{i,i}=0
# Line 333  class StokesProblemCartesian(Homogeneous Line 769  class StokesProblemCartesian(Homogeneous
769            u=0 where  fixed_u_mask>0            u=0 where  fixed_u_mask>0
770            eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j            eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j
771    
772        if surface_stress is not given 0 is assumed.       if surface_stress is not given 0 is assumed.
773    
774        typical usage:       typical usage:
775    
776              sp=StokesProblemCartesian(domain)              sp=StokesProblemCartesian(domain)
777              sp.setTolerance()              sp.setTolerance()
778              sp.initialize(...)              sp.initialize(...)
779              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
780        """       """
781        def __init__(self,domain,**kwargs):       def __init__(self,domain,**kwargs):
782           """           """
783           initialize the Stokes Problem           initialize the Stokes Problem
784    
785           @param domain: domain of the problem. The approximation order needs to be two.           The approximation spaces used for velocity (=Solution(domain)) and pressure (=ReducedSolution(domain)) must be
786           @type domain: L{Domain}           LBB complient, for instance using quadratic and linear approximation on the same element or using linear approximation
787           @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.           with macro elements for the pressure.
788    
789             :param domain: domain of the problem.
790             :type domain: `Domain`
791           """           """
792           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,**kwargs)
793           self.domain=domain           self.domain=domain
          self.vol=util.integrate(1.,Function(self.domain))  
794           self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())           self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
795           self.__pde_u.setSymmetryOn()           self.__pde_u.setSymmetryOn()
796           # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)      
          # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.RILU)  
               
797           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
798           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
          # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)  
799           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
800    
801           self.__pde_proj=LinearPDE(domain)           self.__pde_proj=LinearPDE(domain)
802           self.__pde_proj.setReducedOrderOn()           self.__pde_proj.setReducedOrderOn()
803         self.__pde_proj.setValue(D=1)
804           self.__pde_proj.setSymmetryOn()           self.__pde_proj.setSymmetryOn()
          self.__pde_proj.setValue(D=1.)  
805    
806        def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):       def getSolverOptionsVelocity(self):
807             """
808         returns the solver options used  solve the equation for velocity.
809        
810         :rtype: `SolverOptions`
811         """
812         return self.__pde_u.getSolverOptions()
813         def setSolverOptionsVelocity(self, options=None):
814             """
815         set the solver options for solving the equation for velocity.
816        
817         :param options: new solver  options
818         :type options: `SolverOptions`
819         """
820             self.__pde_u.setSolverOptions(options)
821         def getSolverOptionsPressure(self):
822             """
823         returns the solver options used  solve the equation for pressure.
824         :rtype: `SolverOptions`
825         """
826         return self.__pde_prec.getSolverOptions()
827         def setSolverOptionsPressure(self, options=None):
828             """
829         set the solver options for solving the equation for pressure.
830         :param options: new solver  options
831         :type options: `SolverOptions`
832         """
833         self.__pde_prec.setSolverOptions(options)
834    
835         def setSolverOptionsDiv(self, options=None):
836             """
837         set the solver options for solving the equation to project the divergence of
838         the velocity onto the function space of presure.
839        
840         :param options: new solver options
841         :type options: `SolverOptions`
842         """
843         self.__pde_proj.setSolverOptions(options)
844         def getSolverOptionsDiv(self):
845             """
846         returns the solver options for solving the equation to project the divergence of
847         the velocity onto the function space of presure.
848        
849         :rtype: `SolverOptions`
850         """
851         return self.__pde_proj.getSolverOptions()
852    
853         def updateStokesEquation(self, v, p):
854             """
855             updates the Stokes equation to consider dependencies from ``v`` and ``p``
856             :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values.
857             """
858             pass
859         def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None):
860            """
861            assigns new values to the model parameters.
862    
863            :param f: external force
864            :type f: `Vector` object in `FunctionSpace` `Function` or similar
865            :param fixed_u_mask: mask of locations with fixed velocity.
866            :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
867            :param eta: viscosity
868            :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
869            :param surface_stress: normal surface stress
870            :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
871            :param stress: initial stress
872        :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
873            """
874            if eta !=None:
875                k=util.kronecker(self.domain.getDim())
876                kk=util.outer(k,k)
877                self.eta=util.interpolate(eta, escript.Function(self.domain))
878            self.__pde_prec.setValue(D=1/self.eta)
879                self.__pde_u.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))
880            if restoration_factor!=None:
881                n=self.domain.getNormal()
882                self.__pde_u.setValue(d=restoration_factor*util.outer(n,n))
883            if fixed_u_mask!=None:
884                self.__pde_u.setValue(q=fixed_u_mask)
885            if f!=None: self.__f=f
886            if surface_stress!=None: self.__surface_stress=surface_stress
887            if stress!=None: self.__stress=stress
888    
889         def initialize(self,f=escript.Data(),fixed_u_mask=escript.Data(),eta=1,surface_stress=escript.Data(),stress=escript.Data(), restoration_factor=0):
890          """          """
891          assigns values to the model parameters          assigns values to the model parameters
892    
893          @param f: external force          :param f: external force
894          @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar          :type f: `Vector` object in `FunctionSpace` `Function` or similar
895          @param fixed_u_mask: mask of locations with fixed velocity.          :param fixed_u_mask: mask of locations with fixed velocity.
896          @type fixed_u_mask: L{Vector} object on L{FunctionSpace} L{Solution} or similar          :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
897          @param eta: viscosity          :param eta: viscosity
898          @type eta: L{Scalar} object on L{FunctionSpace} L{Function} or similar          :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
899          @param surface_stress: normal surface stress          :param surface_stress: normal surface stress
900          @type eta: L{Vector} object on L{FunctionSpace} L{FunctionOnBoundary} or similar          :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
901          @param stress: initial stress          :param stress: initial stress
902      @type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar      :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
         @note: All values needs to be set.  
   
903          """          """
904          self.eta=eta          self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
         A =self.__pde_u.createCoefficient("A")  
     self.__pde_u.setValue(A=Data())  
         for i in range(self.domain.getDim()):  
         for j in range(self.domain.getDim()):  
             A[i,j,j,i] += 1.  
             A[i,j,i,j] += 1.  
     self.__pde_prec.setValue(D=1/self.eta)  
         self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress)  
         self.__stress=stress  
   
       def B(self,v):  
         """  
         returns div(v)  
         @rtype: equal to the type of p  
   
         @note: boundary conditions on p should be zero!  
         """  
         if self.show_details: print "apply divergence:"  
         self.__pde_proj.setValue(Y=-util.div(v))  
         self.__pde_proj.setTolerance(self.getSubProblemTolerance())  
         return self.__pde_proj.getSolution(verbose=self.show_details)  
   
       def inner_pBv(self,p,Bv):  
          """  
          returns inner product of element p and Bv  (overwrite)  
           
          @type p: equal to the type of p  
          @type Bv: equal to the type of result of operator B  
          @rtype: C{float}  
905    
906           @rtype: equal to the type of p       def Bv(self,v,tol):
907           """           """
908           s0=util.interpolate(p,Function(self.domain))           returns inner product of element p and div(v)
909           s1=util.interpolate(Bv,Function(self.domain))  
910           return util.integrate(s0*s1)           :param v: a residual
911             :return: inner product of element p and div(v)
912             :rtype: ``float``
913             """
914             self.__pde_proj.setValue(Y=-util.div(v))
915         self.getSolverOptionsDiv().setTolerance(tol)
916         self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
917             out=self.__pde_proj.getSolution()
918             return out
919    
920         def inner_pBv(self,p,Bv):
921             """
922             returns inner product of element p and Bv=-div(v)
923    
924        def inner_p(self,p0,p1):           :param p: a pressure increment
925             :param Bv: a residual
926             :return: inner product of element p and Bv=-div(v)
927             :rtype: ``float``
928           """           """
929           returns inner product of element p0 and p1  (overwrite)           return util.integrate(util.interpolate(p,escript.Function(self.domain))*util.interpolate(Bv, escript.Function(self.domain)))
           
          @type p0: equal to the type of p  
          @type p1: equal to the type of p  
          @rtype: C{float}  
930    
931           @rtype: equal to the type of p       def inner_p(self,p0,p1):
932           """           """
933           s0=util.interpolate(p0/self.eta,Function(self.domain))           Returns inner product of p0 and p1
934           s1=util.interpolate(p1/self.eta,Function(self.domain))  
935             :param p0: a pressure
936             :param p1: a pressure
937             :return: inner product of p0 and p1
938             :rtype: ``float``
939             """
940             s0=util.interpolate(p0, escript.Function(self.domain))
941             s1=util.interpolate(p1, escript.Function(self.domain))
942           return util.integrate(s0*s1)           return util.integrate(s0*s1)
943    
944        def inner_v(self,v0,v1):       def norm_v(self,v):
945           """           """
946           returns inner product of two element v0 and v1  (overwrite)           returns the norm of v
           
          @type v0: equal to the type of v  
          @type v1: equal to the type of v  
          @rtype: C{float}  
947    
948           @rtype: equal to the type of v           :param v: a velovity
949             :return: norm of v
950             :rtype: non-negative ``float``
951           """           """
952       gv0=util.grad(v0)           return util.sqrt(util.integrate(util.length(util.grad(v))**2))
      gv1=util.grad(v1)  
          return util.integrate(util.inner(gv0,gv1))  
953    
954        def solve_A(self,u,p):  
955         def getDV(self, p, v, tol):
956           """           """
957           solves Av=f-Au-B^*p (v=0 on fixed_u_mask)           return the value for v for a given p (overwrite)
958    
959             :param p: a pressure
960             :param v: a initial guess for the value v to return.
961             :return: dv given as *Adv=(f-Av-B^*p)*
962           """           """
963           if self.show_details: print "solve for velocity:"           self.updateStokesEquation(v,p)
964           self.__pde_u.setTolerance(self.getSubProblemTolerance())           self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress)
965         self.getSolverOptionsVelocity().setTolerance(tol)
966         self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
967           if self.__stress.isEmpty():           if self.__stress.isEmpty():
968              self.__pde_u.setValue(X=-2*self.eta*util.symmetric(util.grad(u))+p*util.kronecker(self.domain))              self.__pde_u.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
969           else:           else:
970              self.__pde_u.setValue(X=self.__stress-2*self.eta*util.symmetric(util.grad(u))+p*util.kronecker(self.domain))              self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
971           out=self.__pde_u.getSolution(verbose=self.show_details)           out=self.__pde_u.getSolution()
972           return  out           return  out
973    
974        def solve_prec(self,p):       def norm_Bv(self,Bv):
975           if self.show_details: print "apply preconditioner:"          """
976           self.__pde_prec.setTolerance(self.getSubProblemTolerance())          Returns Bv (overwrite).
977           self.__pde_prec.setValue(Y=p)  
978           q=self.__pde_prec.getSolution(verbose=self.show_details)          :rtype: equal to the type of p
979           return q          :note: boundary conditions on p should be zero!
980            """
981            return util.sqrt(util.integrate(util.interpolate(Bv, escript.Function(self.domain))**2))
982    
983         def solve_AinvBt(self,p, tol):
984             """
985             Solves *Av=B^*p* with accuracy `tol`
986    
987             :param p: a pressure increment
988             :return: the solution of *Av=B^*p*
989             :note: boundary conditions on v should be zero!
990             """
991             self.__pde_u.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))
992             out=self.__pde_u.getSolution()
993             return  out
994    
995         def solve_prec(self,Bv, tol):
996             """
997             applies preconditioner for for *BA^{-1}B^** to *Bv*
998             with accuracy `self.getSubProblemTolerance()`
999    
1000             :param Bv: velocity increment
1001             :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*
1002             :note: boundary conditions on p are zero.
1003             """
1004             self.__pde_prec.setValue(Y=Bv)
1005         self.getSolverOptionsPressure().setTolerance(tol)
1006         self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
1007             out=self.__pde_prec.getSolution()
1008             return out

Legend:
Removed from v.2208  
changed lines
  Added in v.3432

  ViewVC Help
Powered by ViewVC 1.1.26