/[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 2549 by jfenwick, Mon Jul 20 06:43:47 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-2009 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-2009 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"""
# Line 21  __url__="https://launchpad.net/escript-f Line 22  __url__="https://launchpad.net/escript-f
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, SolverOptions  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions
38  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES  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, weight=None, useReduced=False, adaptSubTolerance=True):      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      :param useReduced: uses reduced oreder on flux and pressure
458      @type useReduced: C{bool}      :type useReduced: ``bool``
459      @param adaptSubTolerance: switches on automatic subtolerance selection      :param adaptSubTolerance: switches on automatic subtolerance selection
460      @type adaptSubTolerance: C{bool}          :type adaptSubTolerance: ``bool``  
461          """          """
462          self.domain=domain          self.domain=domain
463          if weight == None:          if weight == None:
464             s=self.domain.getSize()             s=self.domain.getSize()
465             self.__l=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2             self.__l=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2
466             # self.__l=(3.*util.longestEdge(self.domain))**2             # self.__l=(3.*util.longestEdge(self.domain))**2
467             # self.__l=(0.1*util.longestEdge(self.domain)*s/util.sup(s))**2             #self.__l=(0.1*util.longestEdge(self.domain)*s/util.sup(s))**2
468          else:          else:
469             self.__l=weight             self.__l=weight
470          self.__pde_v=LinearPDESystem(domain)          self.__pde_v=LinearPDESystem(domain)
# Line 73  class DarcyFlow(object): Line 474  class DarcyFlow(object):
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.setTolerance()          self.setTolerance()
480          self.setAbsoluteTolerance()          self.setAbsoluteTolerance()
481      self.__adaptSubTolerance=adaptSubTolerance      self.__adaptSubTolerance=adaptSubTolerance
# Line 83  class DarcyFlow(object): Line 484  class DarcyFlow(object):
484      """      """
485      Returns the solver options used to solve the flux problems      Returns the solver options used to solve the flux problems
486            
487      M{(I+D^*D)u=F}      *(I+D^*D)u=F*
488            
489      @return: L{SolverOptions}      :return: `SolverOptions`
490      """      """
491      return self.__pde_v.getSolverOptions()      return self.__pde_v.getSolverOptions()
492      def setSolverOptionsFlux(self, options=None):      def setSolverOptionsFlux(self, options=None):
493      """      """
494      Sets the solver options used to solve the flux problems      Sets the solver options used to solve the flux problems
495            
496      M{(I+D^*D)u=F}      *(I+D^*D)u=F*
497            
498      If C{options} is not present, the options are reset to default      If ``options`` is not present, the options are reset to default
499      @param options: L{SolverOptions}      :param options: `SolverOptions`
500      @note: if the adaption of subtolerance is choosen, the tolerance set by C{options} will be overwritten before the solver is called.      :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)      return self.__pde_v.setSolverOptions(options)
503      def getSolverOptionsPressure(self):      def getSolverOptionsPressure(self):
504      """      """
505      Returns the solver options used to solve the pressure problems      Returns the solver options used to solve the pressure problems
506            
507      M{(Q^*Q)p=Q^*G}      *(Q^*Q)p=Q^*G*
508            
509      @return: L{SolverOptions}      :return: `SolverOptions`
510      """      """
511      return self.__pde_p.getSolverOptions()      return self.__pde_p.getSolverOptions()
512      def setSolverOptionsPressure(self, options=None):      def setSolverOptionsPressure(self, options=None):
513      """      """
514      Sets the solver options used to solve the pressure problems      Sets the solver options used to solve the pressure problems
515            
516      M{(Q^*Q)p=Q^*G}      *(Q^*Q)p=Q^*G*
517            
518      If C{options} is not present, the options are reset to default      If ``options`` is not present, the options are reset to default
519      @param options: L{SolverOptions}      :param options: `SolverOptions`
520      @note: if the adaption of subtolerance is choosen, the tolerance set by C{options} will be overwritten before the solver is called.      :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)      return self.__pde_p.setSolverOptions(options)
523    
# Line 124  class DarcyFlow(object): Line 525  class DarcyFlow(object):
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 177  class DarcyFlow(object): Line 578  class DarcyFlow(object):
578    
579      def setTolerance(self,rtol=1e-4):      def setTolerance(self,rtol=1e-4):
580          """          """
581          sets the relative tolerance C{rtol} used to terminate the solution process. The iteration is terminated if          sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if
582    
583          M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }          *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
584    
585          where C{atol} is an absolut tolerance (see L{setAbsoluteTolerance}), M{|f|^2 = integrate(length(f)^2)} and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.          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          @param rtol: relative tolerance for the pressure          :param rtol: relative tolerance for the pressure
588          @type rtol: non-negative C{float}          :type rtol: non-negative ``float``
589          """          """
590          if rtol<0:          if rtol<0:
591              raise ValueError,"Relative tolerance needs to be non-negative."              raise ValueError,"Relative tolerance needs to be non-negative."
# Line 193  class DarcyFlow(object): Line 594  class DarcyFlow(object):
594          """          """
595          returns the relative tolerance          returns the relative tolerance
596    
597          @return: current relative tolerance          :return: current relative tolerance
598          @rtype: C{float}          :rtype: ``float``
599          """          """
600          return self.__rtol          return self.__rtol
601    
602      def setAbsoluteTolerance(self,atol=0.):      def setAbsoluteTolerance(self,atol=0.):
603          """          """
604          sets the absolute tolerance C{atol} used to terminate the solution process. The iteration is terminated if          sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if
605    
606          M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }          *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
607    
608          where C{rtol} is an absolut tolerance (see L{setTolerance}), M{|f|^2 = integrate(length(f)^2)} and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.          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          :param atol: absolute tolerance for the pressure
611          @type atol: non-negative C{float}          :type atol: non-negative ``float``
612          """          """
613          if atol<0:          if atol<0:
614              raise ValueError,"Absolute tolerance needs to be non-negative."              raise ValueError,"Absolute tolerance needs to be non-negative."
# Line 216  class DarcyFlow(object): Line 617  class DarcyFlow(object):
617         """         """
618         returns the absolute tolerance         returns the absolute tolerance
619                
620         @return: current absolute tolerance         :return: current absolute tolerance
621         @rtype: C{float}         :rtype: ``float``
622         """         """
623         return self.__atol         return self.__atol
624      def getSubProblemTolerance(self):      def getSubProblemTolerance(self):
625      """      """
626      Returns a suitable subtolerance      Returns a suitable subtolerance
627      @type: C{float}      @type: ``float``
628      """      """
629      return max(util.EPSILON**(0.75),self.getTolerance()**2)      return max(util.EPSILON**(0.75),self.getTolerance()**2)
630      def setSubProblemTolerance(self):      def setSubProblemTolerance(self):
# Line 244  class DarcyFlow(object): Line 645  class DarcyFlow(object):
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 verbose: if set some information on iteration progress are printed           :param verbose: if set some information on iteration progress are printed
653           @type verbose: C{bool}           :type verbose: ``bool``
654           @return: flux and pressure           :return: flux and pressure
655           @rtype: C{tuple} of L{Data}.           :rtype: ``tuple`` of `Data`.
656    
657           @note: The problem is solved as a least squares form           :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           rtol=self.getTolerance()           rtol=self.getTolerance()
676           atol=self.getAbsoluteTolerance()           atol=self.getAbsoluteTolerance()
677       self.setSubProblemTolerance()       self.setSubProblemTolerance()
       
678           num_corrections=0           num_corrections=0
679           converged=False           converged=False
680           p=p0           p=p0
# Line 297  class DarcyFlow(object): Line 697  class DarcyFlow(object):
697                 ATOL=(atol+rtol*fac)                 ATOL=(atol+rtol*fac)
698                 if self.verbose:                 if self.verbose:
699                      print "DarcyFlux: L2 norm of v = %e."%norm_v                      print "DarcyFlux: L2 norm of v = %e."%norm_v
700                      print "DarcyFlux: L2 norm of k.grad(p) = %e."%norm_Qp                      print "DarcyFlux: L2 norm of k.util.grad(p) = %e."%norm_Qp
701                      print "DarcyFlux: L2 defect u = %e."%(util.integrate(util.length(self.__g-util.interpolate(v,Function(self.domain))-Qp)**2)**(0.5),)                      print "DarcyFlux: L2 defect u = %e."%(util.integrate(util.length(self.__g-util.interpolate(v,escript.Function(self.domain))-Qp)**2)**(0.5),)
702                      print "DarcyFlux: L2 defect div(v) = %e."%(util.integrate((self.__f-util.div(v))**2)**(0.5),)                      print "DarcyFlux: L2 defect div(v) = %e."%(util.integrate((self.__f-util.div(v))**2)**(0.5),)
703                      print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL                      print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
704                 if norm_r == None or norm_r>ATOL:                 if norm_r == None or norm_r>ATOL:
705                     if num_corrections>max_num_corrections:                     if num_corrections>max_num_corrections:
706                           raise ValueError,"maximum number of correction steps reached."                           raise ValueError,"maximum number of correction steps reached."
707                     p,r, norm_r=PCG(self.__g-util.interpolate(v,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)                     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                     num_corrections+=1                     num_corrections+=1
709                 else:                 else:
710                     converged=True                     converged=True
711           return v,p           return v,p
712      def __L2(self,v):      def __L2(self,v):
713           return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))           return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))
714    
715      def __Q(self,p):      def __Q(self,p):
716            return util.tensor_mult(self.__permeability,util.grad(p))            return util.tensor_mult(self.__permeability,util.grad(p))
# Line 318  class DarcyFlow(object): Line 718  class DarcyFlow(object):
718      def __Aprod(self,dp):      def __Aprod(self,dp):
719            if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"            if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
720            Qdp=self.__Q(dp)            Qdp=self.__Q(dp)
721            self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())            self.__pde_v.setValue(Y=-Qdp,X=escript.Data(), r=escript.Data())
722            du=self.__pde_v.getSolution()            du=self.__pde_v.getSolution()
723            # self.__pde_v.getOperator().saveMM("proj.mm")            # self.__pde_v.getOperator().saveMM("proj.mm")
724            return Qdp+du            return Qdp+du
# Line 330  class DarcyFlow(object): Line 730  class DarcyFlow(object):
730    
731      def __Msolve_PCG(self,r):      def __Msolve_PCG(self,r):
732        if self.getSolverOptionsPressure().isVerbose(): 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), Y=Data(), r=Data())            self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=escript.Data(), r=escript.Data())
734            # self.__pde_p.getOperator().saveMM("prec.mm")            # self.__pde_p.getOperator().saveMM("prec.mm")
735            return self.__pde_p.getSolution()            return self.__pde_p.getSolution()
736    
737      def getFlux(self,p=None, fixed_flux=Data()):      def getFlux(self,p=None, fixed_flux=escript.Data()):
738          """          """
739          returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux}          returns the flux for a given pressure ``p`` where the flux is equal to ``fixed_flux``
740          on locations where C{location_of_fixed_flux} is positive (see L{setValue}).          on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
741          Note that C{g} and C{f} are used, see L{setValue}.          Note that ``g`` and ``f`` are used, see `setValue`.
742    
743          @param p: pressure.          :param p: pressure.
744          @type p: scalar value on the domain (e.g. L{Data}).          :type p: scalar value on the domain (e.g. `Data`).
745          @param fixed_flux: flux on the locations of the domain marked be C{location_of_fixed_flux}.          :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. L{Data}).          :type fixed_flux: vector values on the domain (e.g. `Data`).
747          @param tol: relative tolerance to be used.          :return: flux
748          @type tol: positive C{float}.          :rtype: `Data`
749          @return: flux          :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          @rtype: L{Data}                 for the permeability *k_{ij}*
         @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}}  
751          """          """
752      self.setSubProblemTolerance()      self.setSubProblemTolerance()
753          g=self.__g          g=self.__g
# Line 380  class StokesProblemCartesian(Homogeneous Line 778  class StokesProblemCartesian(Homogeneous
778              sp.initialize(...)              sp.initialize(...)
779              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
780       """       """
781       def __init__(self,domain,adaptSubTolerance=True, **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       @param adaptSubTolerance: If True the tolerance for subproblem is set automatically.           with macro elements for the pressure.
788       @type adaptSubTolerance: C{bool}  
789           @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.           :param domain: domain of the problem.
790             :type domain: `Domain`
791           """           """
792           HomogeneousSaddlePointProblem.__init__(self,adaptSubTolerance=adaptSubTolerance,**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            
# Line 409  class StokesProblemCartesian(Homogeneous Line 807  class StokesProblemCartesian(Homogeneous
807           """           """
808       returns the solver options used  solve the equation for velocity.       returns the solver options used  solve the equation for velocity.
809            
810       @rtype: L{SolverOptions}       :rtype: `SolverOptions`
811       """       """
812       return self.__pde_u.getSolverOptions()       return self.__pde_u.getSolverOptions()
813       def setSolverOptionsVelocity(self, options=None):       def setSolverOptionsVelocity(self, options=None):
814           """           """
815       set the solver options for solving the equation for velocity.       set the solver options for solving the equation for velocity.
816            
817       @param options: new solver  options       :param options: new solver  options
818       @type options: L{SolverOptions}       :type options: `SolverOptions`
819       """       """
820           self.__pde_u.setSolverOptions(options)           self.__pde_u.setSolverOptions(options)
821       def getSolverOptionsPressure(self):       def getSolverOptionsPressure(self):
822           """           """
823       returns the solver options used  solve the equation for pressure.       returns the solver options used  solve the equation for pressure.
824       @rtype: L{SolverOptions}       :rtype: `SolverOptions`
825       """       """
826       return self.__pde_prec.getSolverOptions()       return self.__pde_prec.getSolverOptions()
827       def setSolverOptionsPressure(self, options=None):       def setSolverOptionsPressure(self, options=None):
828           """           """
829       set the solver options for solving the equation for pressure.       set the solver options for solving the equation for pressure.
830       @param options: new solver  options       :param options: new solver  options
831       @type options: L{SolverOptions}       :type options: `SolverOptions`
832       """       """
833       self.__pde_prec.setSolverOptions(options)       self.__pde_prec.setSolverOptions(options)
834    
# Line 439  class StokesProblemCartesian(Homogeneous Line 837  class StokesProblemCartesian(Homogeneous
837       set the solver options for solving the equation to project the divergence of       set the solver options for solving the equation to project the divergence of
838       the velocity onto the function space of presure.       the velocity onto the function space of presure.
839            
840       @param options: new solver options       :param options: new solver options
841       @type options: L{SolverOptions}       :type options: `SolverOptions`
842       """       """
843       self.__pde_prec.setSolverOptions(options)       self.__pde_proj.setSolverOptions(options)
844       def getSolverOptionsDiv(self):       def getSolverOptionsDiv(self):
845           """           """
846       returns the solver options for solving the equation to project the divergence of       returns the solver options for solving the equation to project the divergence of
847       the velocity onto the function space of presure.       the velocity onto the function space of presure.
848            
849       @rtype: L{SolverOptions}       :rtype: `SolverOptions`
850       """       """
851       return self.__pde_prec.getSolverOptions()       return self.__pde_proj.getSolverOptions()
852       def setSubProblemTolerance(self):  
853         def updateStokesEquation(self, v, p):
854           """           """
855       Updates the tolerance for subproblems           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       if self.adaptSubTolerance():           pass
859               sub_tol=self.getSubProblemTolerance()       def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None):
860           self.getSolverOptionsDiv().setTolerance(sub_tol)          """
861           self.getSolverOptionsDiv().setAbsoluteTolerance(0.)          assigns new values to the model parameters.
862           self.getSolverOptionsPressure().setTolerance(sub_tol)  
863           self.getSolverOptionsPressure().setAbsoluteTolerance(0.)          :param f: external force
864           self.getSolverOptionsVelocity().setTolerance(sub_tol)          :type f: `Vector` object in `FunctionSpace` `Function` or similar
865           self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)          :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=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):       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
903          @note: All values needs to be set.          """
904          """          self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
         self.eta=eta  
         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)  
         self.__f=f  
         self.__surface_stress=surface_stress  
         self.__stress=stress  
905    
906       def Bv(self,v):       def Bv(self,v,tol):
907           """           """
908           returns inner product of element p and div(v)           returns inner product of element p and div(v)
909    
910           @param p: a pressure increment           :param v: a residual
911           @param v: a residual           :return: inner product of element p and div(v)
912           @return: inner product of element p and div(v)           :rtype: ``float``
913           @rtype: C{float}           """
914           """           self.__pde_proj.setValue(Y=-util.div(v))
915           self.__pde_proj.setValue(Y=-util.div(v))       self.getSolverOptionsDiv().setTolerance(tol)
916           return self.__pde_proj.getSolution()       self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
917             out=self.__pde_proj.getSolution()
918             return out
919    
920       def inner_pBv(self,p,Bv):       def inner_pBv(self,p,Bv):
921           """           """
922           returns inner product of element p and Bv=-div(v)           returns inner product of element p and Bv=-div(v)
923    
924           @param p: a pressure increment           :param p: a pressure increment
925           @param v: a residual           :param Bv: a residual
926           @return: inner product of element p and Bv=-div(v)           :return: inner product of element p and Bv=-div(v)
927           @rtype: C{float}           :rtype: ``float``
928           """           """
929           return util.integrate(util.interpolate(p,Function(self.domain))*util.interpolate(Bv,Function(self.domain)))           return util.integrate(util.interpolate(p,escript.Function(self.domain))*util.interpolate(Bv, escript.Function(self.domain)))
930    
931       def inner_p(self,p0,p1):       def inner_p(self,p0,p1):
932           """           """
933           Returns inner product of p0 and p1           Returns inner product of p0 and p1
934    
935           @param p0: a pressure           :param p0: a pressure
936           @param p1: a pressure           :param p1: a pressure
937           @return: inner product of p0 and p1           :return: inner product of p0 and p1
938           @rtype: C{float}           :rtype: ``float``
939           """           """
940           s0=util.interpolate(p0/self.eta,Function(self.domain))           s0=util.interpolate(p0, escript.Function(self.domain))
941           s1=util.interpolate(p1/self.eta,Function(self.domain))           s1=util.interpolate(p1, escript.Function(self.domain))
942           return util.integrate(s0*s1)           return util.integrate(s0*s1)
943    
944       def norm_v(self,v):       def norm_v(self,v):
945           """           """
946           returns the norm of v           returns the norm of v
947    
948           @param v: a velovity           :param v: a velovity
949           @return: norm of v           :return: norm of v
950           @rtype: non-negative C{float}           :rtype: non-negative ``float``
951           """           """
952           return util.sqrt(util.integrate(util.length(util.grad(v))))           return util.sqrt(util.integrate(util.length(util.grad(v))**2))
953    
954    
955       def getV(self, p, v0):       def getDV(self, p, v, tol):
956           """           """
957           return the value for v for a given p (overwrite)           return the value for v for a given p (overwrite)
958    
959           @param p: a pressure           :param p: a pressure
960           @param v0: a initial guess for the value v to return.           :param v: a initial guess for the value v to return.
961           @return: v given as M{v= A^{-1} (f-B^*p)}           :return: dv given as *Adv=(f-Av-B^*p)*
962           """           """
963           self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress, r=v0)           self.updateStokesEquation(v,p)
964             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=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+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()           out=self.__pde_u.getSolution()
972           return  out           return  out
973    
# Line 560  class StokesProblemCartesian(Homogeneous Line 975  class StokesProblemCartesian(Homogeneous
975          """          """
976          Returns Bv (overwrite).          Returns Bv (overwrite).
977    
978          @rtype: equal to the type of p          :rtype: equal to the type of p
979          @note: boundary conditions on p should be zero!          :note: boundary conditions on p should be zero!
980          """          """
981          return util.sqrt(util.integrate(util.interpolate(Bv,Function(self.domain))**2))          return util.sqrt(util.integrate(util.interpolate(Bv, escript.Function(self.domain))**2))
982    
983       def solve_AinvBt(self,p):       def solve_AinvBt(self,p, tol):
984           """           """
985           Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}           Solves *Av=B^*p* with accuracy `tol`
986    
987           @param p: a pressure increment           :param p: a pressure increment
988           @return: the solution of M{Av=B^*p}           :return: the solution of *Av=B^*p*
989           @note: boundary conditions on v should be zero!           :note: boundary conditions on v should be zero!
990           """           """
991           self.__pde_u.setValue(Y=Data(), y=Data(), r=Data(),X=-p*util.kronecker(self.domain))           self.__pde_u.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))
992           out=self.__pde_u.getSolution()           out=self.__pde_u.getSolution()
993           return  out           return  out
994    
995       def solve_prec(self,Bv):       def solve_prec(self,Bv, tol):
996           """           """
997           applies preconditioner for for M{BA^{-1}B^*} to M{Bv}           applies preconditioner for for *BA^{-1}B^** to *Bv*
998           with accuracy L{self.getSubProblemTolerance()}           with accuracy `self.getSubProblemTolerance()`
999    
1000           @param v: velocity increment           :param Bv: velocity increment
1001           @return: M{p=P(Bv)} where M{P^{-1}} is an approximation of M{BA^{-1}B^*}           :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*
1002           @note: boundary conditions on p are zero.           :note: boundary conditions on p are zero.
1003           """           """
1004           self.__pde_prec.setValue(Y=Bv)           self.__pde_prec.setValue(Y=Bv)
1005           return self.__pde_prec.getSolution()       self.getSolverOptionsPressure().setTolerance(tol)
1006         self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
1007             out=self.__pde_prec.getSolution()
1008             return out

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

  ViewVC Help
Powered by ViewVC 1.1.26