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

Annotation of /trunk/escriptcore/py_src/flows.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3074 - (hide annotations)
Tue Jul 27 01:47:45 2010 UTC (9 years ago) by gross
Original Path: trunk/escript/py_src/flows.py
File MIME type: text/x-python
File size: 40295 byte(s)
tolerance for PCG relaxed.
1 gross 3071 # -*- coding: utf-8 -*-
2 ksteube 1809 ########################################################
3 gross 1414 #
4 jfenwick 2881 # Copyright (c) 2003-2010 by University of Queensland
5 ksteube 1809 # Earth Systems Science Computational Center (ESSCC)
6     # http://www.uq.edu.au/esscc
7 gross 1414 #
8 ksteube 1809 # Primary Business: Queensland, Australia
9     # Licensed under the Open Software License version 3.0
10     # http://www.opensource.org/licenses/osl-3.0.php
11 gross 1414 #
12 ksteube 1809 ########################################################
13 gross 1414
14 jfenwick 2881 __copyright__="""Copyright (c) 2003-2010 by University of Queensland
15 ksteube 1809 Earth Systems Science Computational Center (ESSCC)
16     http://www.uq.edu.au/esscc
17     Primary Business: Queensland, Australia"""
18     __license__="""Licensed under the Open Software License version 3.0
19     http://www.opensource.org/licenses/osl-3.0.php"""
20 jfenwick 2344 __url__="https://launchpad.net/escript-finley"
21 ksteube 1809
22 gross 1414 """
23     Some models for flow
24    
25 jfenwick 2625 :var __author__: name of author
26     :var __copyright__: copyrights
27     :var __license__: licence agreement
28     :var __url__: url entry point on documentation
29     :var __version__: version
30     :var __date__: date of the version
31 gross 1414 """
32    
33     __author__="Lutz Gross, l.gross@uq.edu.au"
34    
35     from escript import *
36     import util
37 gross 2474 from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions
38 gross 2264 from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
39 gross 1414
40 gross 2100 class DarcyFlow(object):
41 gross 3071 """
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 gross 1659
74 gross 3071 self.__pde_p=LinearSinglePDE(domain)
75     self.__pde_p.setSymmetryOn()
76     if self.useReduced: self.__pde_p.setReducedOrderOn()
77 gross 1659
78 gross 3071 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=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))
84     self.__g=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 gross 1659
132 gross 3071 def getSolverOptionsWeighting(self):
133     """
134     Returns the solver options used to solve the pressure problems
135 gross 1659
136 gross 3071 *(D K D^*)p=-D F*
137 gross 2960
138 gross 3071 :return: `SolverOptions`
139     """
140     return self.__pde_l.getSolverOptions()
141 gross 2960
142 gross 3071 def setSolverOptionsWeighting(self, options=None):
143     """
144     Sets the solver options used to solve the pressure problems
145 gross 2960
146 gross 3071 *(D K D^*)p=-D F*
147 gross 2960
148 gross 3071 If ``options`` is not present, the options are reset to default
149 gross 2960
150 gross 3071 :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 gross 2960
155    
156 gross 3071 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
169     ``s`` on the main diagonal is used.
170     :type permeability: scalar or tensor values on the domain (e.g. `Data`)
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
174     is along the *x_i* axis.
175     """
176     if f !=None:
177     f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("X"))
178     if f.isEmpty():
179     f=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=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 gross 2960
207 gross 3071 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 gross 2960
219 gross 3071 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 gross 2960
242 gross 3071 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=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 gross 2960
249 gross 3071 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=Data())
253     # self.__pde_p.getOperator().saveMM("prec.mm")
254     return self.__pde_k.getSolution()
255 gross 2960
256 gross 3071 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=Data(), r=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 gross 2960
267 gross 3071 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=Data(), r=Data(), y=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 gross 2960
289 gross 3071 def __L2(self,v):
290     return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))
291 gross 2960
292 gross 3071 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 gross 2960
298 gross 3071 :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 gross 2960
307 gross 3071 :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 gross 2960
311 gross 3071 where *D* is the *div* operator and *(Gp)_i=p_{,i}* for the permeability *K=k_{ij}*.
312     """
313     self.verbose=verbose
314     rtol=self.getTolerance()
315     atol=self.getAbsoluteTolerance()
316     self.setSubProblemTolerance()
317     num_corrections=0
318     converged=False
319     norm_r=None
320    
321     # Eliminate the hydrostatic pressure:
322     if self.verbose: print "DarcyFlux: calculate hydrostatic pressure component."
323     self.__pde_p.setValue(X=self.__g, r=p0, y=-util.inner(self.domain.getNormal(),u0))
324     p0=self.__pde_p.getSolution()
325     g2=self.__g - util.tensor_mult(self.__permeability, util.grad(p0))
326     norm_g2=util.integrate(util.inner(g2,util.tensor_mult(self.__permeability_inv,g2)))**0.5
327 gross 2960
328 gross 3071 p=p0*0
329     if self.solveForFlux:
330     v=u0.copy()
331     else:
332     v=self.__getFlux(p, u0, f=self.__f, g=g2)
333 gross 2960
334 gross 3071 while not converged and norm_g2 > 0:
335     Gp=util.grad(p)
336     KGp=util.tensor_mult(self.__permeability,Gp)
337     if self.verbose:
338     def_p=g2-(v+KGp)
339     def_v=self.__f-util.div(v)
340     print "DarcyFlux: L2: g-v-K*grad(p) = %e (v = %e)."%(self.__L2(def_p),self.__L2(v))
341     print "DarcyFlux: L2: f-div(v) = %e (grad(v) = %e)."%(self.__L2(def_v),self.__L2(util.grad(v)))
342     print "DarcyFlux: K^{-1}-norm of v = %e."%util.integrate(util.inner(v,util.tensor_mult(self.__permeability_inv,v)))**0.5
343     print "DarcyFlux: K^{-1}-norm of g2 = %e."%norm_g2
344     print "DarcyFlux: K-norm of grad(dp) = %e."%util.integrate(util.inner(Gp,KGp))**0.5
345     ATOL=atol+rtol*norm_g2
346     if self.verbose: print "DarcyFlux: absolute tolerance ATOL = %e."%(ATOL,)
347     if norm_r == None or norm_r>ATOL:
348     if num_corrections>max_num_corrections:
349     raise ValueError,"maximum number of correction steps reached."
350    
351     if self.solveForFlux:
352     # initial residual is r=K^{-1}*(g-v-K*Gp)+D^*L^{-1}(f-Du)
353     v,r, norm_r=PCG(ArithmeticTuple(util.tensor_mult(self.__permeability_inv,g2-v)-Gp,self.__applWeight(v,self.__f),p),
354     self.__Aprod_v,
355     v,
356     self.__Msolve_PCG_v,
357     self.__inner_PCG_v,
358     atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
359     p=r[2]
360     else:
361     # initial residual is r=G^*(g2-KGp - v)
362     p,r, norm_r=PCG(ArithmeticTuple(g2-KGp,v),
363     self.__Aprod_p,
364     p,
365     self.__Msolve_PCG_p,
366     self.__inner_PCG_p,
367     atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
368     v=r[1]
369     if self.verbose: print "DarcyFlux: residual norm = %e."%norm_r
370     num_corrections+=1
371     else:
372     if self.verbose: print "DarcyFlux: stopping criterium reached."
373     converged=True
374     return v,p+p0
375     def setTolerance(self,rtol=1e-4):
376     """
377     sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if
378 gross 2960
379 gross 3071 *|g-v-K gard(p)|_PCG <= atol + rtol * |K^{1/2}g2|_0*
380    
381     where ``atol`` is an absolut tolerance (see `setAbsoluteTolerance`).
382    
383     :param rtol: relative tolerance for the pressure
384     :type rtol: non-negative ``float``
385     """
386     if rtol<0:
387     raise ValueError,"Relative tolerance needs to be non-negative."
388     self.__rtol=rtol
389     def getTolerance(self):
390     """
391     returns the relative tolerance
392     :return: current relative tolerance
393     :rtype: ``float``
394     """
395     return self.__rtol
396 gross 2960
397 gross 3071 def setAbsoluteTolerance(self,atol=0.):
398     """
399     sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if
400    
401     *|g-v-K gard(p)|_PCG <= atol + rtol * |K^{1/2}g2|_0*
402 gross 2960
403    
404 gross 3071 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}*.
405 gross 2960
406 gross 3071 :param atol: absolute tolerance for the pressure
407     :type atol: non-negative ``float``
408     """
409     if atol<0:
410     raise ValueError,"Absolute tolerance needs to be non-negative."
411     self.__atol=atol
412     def getAbsoluteTolerance(self):
413     """
414     returns the absolute tolerance
415     :return: current absolute tolerance
416     :rtype: ``float``
417     """
418     return self.__atol
419     def getSubProblemTolerance(self):
420     """
421     Returns a suitable subtolerance
422     :type: ``float``
423     """
424 gross 3074 return max(util.EPSILON**(0.5),self.getTolerance()**2)
425 gross 2960
426 gross 3071 def setSubProblemTolerance(self):
427     """
428     Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected.
429     """
430     if self.__adaptSubTolerance:
431     sub_tol=self.getSubProblemTolerance()
432     self.getSolverOptionsFlux().setTolerance(sub_tol)
433     self.getSolverOptionsFlux().setAbsoluteTolerance(0.)
434     self.getSolverOptionsPressure().setTolerance(sub_tol)
435     self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
436     self.getSolverOptionsWeighting().setTolerance(sub_tol)
437     self.getSolverOptionsWeighting().setAbsoluteTolerance(0.)
438     if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol
439 gross 2960
440    
441     class DarcyFlowOld(object):
442     """
443     solves the problem
444    
445     *u_i+k_{ij}*p_{,j} = g_i*
446     *u_{i,i} = f*
447    
448     where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,
449    
450     :note: The problem is solved in a least squares formulation.
451     """
452    
453     def __init__(self, domain, weight=None, useReduced=False, adaptSubTolerance=True):
454     """
455     initializes the Darcy flux problem
456     :param domain: domain of the problem
457     :type domain: `Domain`
458     :param useReduced: uses reduced oreder on flux and pressure
459     :type useReduced: ``bool``
460     :param adaptSubTolerance: switches on automatic subtolerance selection
461     :type adaptSubTolerance: ``bool``
462     """
463     self.domain=domain
464     if weight == None:
465 gross 2354 s=self.domain.getSize()
466     self.__l=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2
467 gross 2486 # self.__l=(3.*util.longestEdge(self.domain))**2
468 artak 2689 #self.__l=(0.1*util.longestEdge(self.domain)*s/util.sup(s))**2
469 gross 2351 else:
470     self.__l=weight
471 gross 2100 self.__pde_v=LinearPDESystem(domain)
472 gross 2208 if useReduced: self.__pde_v.setReducedOrderOn()
473     self.__pde_v.setSymmetryOn()
474 gross 2267 self.__pde_v.setValue(D=util.kronecker(domain), A=self.__l*util.outer(util.kronecker(domain),util.kronecker(domain)))
475 gross 2100 self.__pde_p=LinearSinglePDE(domain)
476     self.__pde_p.setSymmetryOn()
477 gross 2208 if useReduced: self.__pde_p.setReducedOrderOn()
478 gross 2100 self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
479     self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
480 gross 2264 self.setTolerance()
481     self.setAbsoluteTolerance()
482 gross 2474 self.__adaptSubTolerance=adaptSubTolerance
483     self.verbose=False
484     def getSolverOptionsFlux(self):
485     """
486     Returns the solver options used to solve the flux problems
487    
488 jfenwick 2625 *(I+D^*D)u=F*
489 gross 2474
490 jfenwick 2625 :return: `SolverOptions`
491 gross 2474 """
492     return self.__pde_v.getSolverOptions()
493     def setSolverOptionsFlux(self, options=None):
494     """
495     Sets the solver options used to solve the flux problems
496    
497 jfenwick 2625 *(I+D^*D)u=F*
498 gross 2474
499 jfenwick 2625 If ``options`` is not present, the options are reset to default
500     :param options: `SolverOptions`
501     :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
502 gross 2474 """
503     return self.__pde_v.setSolverOptions(options)
504     def getSolverOptionsPressure(self):
505     """
506     Returns the solver options used to solve the pressure problems
507    
508 jfenwick 2625 *(Q^*Q)p=Q^*G*
509 gross 2474
510 jfenwick 2625 :return: `SolverOptions`
511 gross 2474 """
512     return self.__pde_p.getSolverOptions()
513     def setSolverOptionsPressure(self, options=None):
514     """
515     Sets the solver options used to solve the pressure problems
516    
517 jfenwick 2625 *(Q^*Q)p=Q^*G*
518 gross 2474
519 jfenwick 2625 If ``options`` is not present, the options are reset to default
520     :param options: `SolverOptions`
521     :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
522 gross 2474 """
523     return self.__pde_p.setSolverOptions(options)
524 gross 1659
525 gross 2100 def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
526     """
527 gross 2208 assigns values to model parameters
528 gross 1659
529 jfenwick 2625 :param f: volumetic sources/sinks
530     :type f: scalar value on the domain (e.g. `Data`)
531     :param g: flux sources/sinks
532     :type g: vector values on the domain (e.g. `Data`)
533     :param location_of_fixed_pressure: mask for locations where pressure is fixed
534     :type location_of_fixed_pressure: scalar value on the domain (e.g. `Data`)
535     :param location_of_fixed_flux: mask for locations where flux is fixed.
536     :type location_of_fixed_flux: vector values on the domain (e.g. `Data`)
537     :param permeability: permeability tensor. If scalar ``s`` is given the tensor with
538     ``s`` on the main diagonal is used. If vector ``v`` is given the tensor with
539     ``v`` on the main diagonal is used.
540     :type permeability: scalar, vector or tensor values on the domain (e.g. `Data`)
541 gross 1659
542 jfenwick 2625 :note: the values of parameters which are not set by calling ``setValue`` are not altered.
543     :note: at any point on the boundary of the domain the pressure (``location_of_fixed_pressure`` >0)
544     or the normal component of the flux (``location_of_fixed_flux[i]>0`` if direction of the normal
545     is along the *x_i* axis.
546 gross 2100 """
547 gross 2264 if f !=None:
548 gross 2100 f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
549     if f.isEmpty():
550     f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
551     else:
552     if f.getRank()>0: raise ValueError,"illegal rank of f."
553 gross 2267 self.__f=f
554 gross 2264 if g !=None:
555 gross 2100 g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
556     if g.isEmpty():
557     g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
558     else:
559     if not g.getShape()==(self.domain.getDim(),):
560     raise ValueError,"illegal shape of g"
561     self.__g=g
562 gross 1659
563 gross 2100 if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)
564     if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux)
565 gross 1659
566 gross 2100 if permeability!=None:
567     perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
568     if perm.getRank()==0:
569     perm=perm*util.kronecker(self.domain.getDim())
570     elif perm.getRank()==1:
571     perm, perm2=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A")), perm
572     for i in range(self.domain.getDim()): perm[i,i]=perm2[i]
573     elif perm.getRank()==2:
574     pass
575     else:
576     raise ValueError,"illegal rank of permeability."
577     self.__permeability=perm
578     self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))
579 gross 1659
580 gross 2264 def setTolerance(self,rtol=1e-4):
581     """
582 jfenwick 2625 sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if
583 gross 1659
584 jfenwick 2625 *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
585 gross 2264
586 jfenwick 2625 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}*.
587 gross 2264
588 jfenwick 2625 :param rtol: relative tolerance for the pressure
589     :type rtol: non-negative ``float``
590 gross 2100 """
591 gross 2264 if rtol<0:
592     raise ValueError,"Relative tolerance needs to be non-negative."
593     self.__rtol=rtol
594     def getTolerance(self):
595 gross 2100 """
596 gross 2264 returns the relative tolerance
597 gross 1659
598 jfenwick 2625 :return: current relative tolerance
599     :rtype: ``float``
600 caltinay 2169 """
601 gross 2264 return self.__rtol
602 gross 1659
603 gross 2264 def setAbsoluteTolerance(self,atol=0.):
604 gross 2208 """
605 jfenwick 2625 sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if
606 gross 1659
607 jfenwick 2625 *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
608 gross 2208
609 jfenwick 2625 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}*.
610 gross 2208
611 jfenwick 2625 :param atol: absolute tolerance for the pressure
612     :type atol: non-negative ``float``
613 gross 2264 """
614     if atol<0:
615     raise ValueError,"Absolute tolerance needs to be non-negative."
616     self.__atol=atol
617     def getAbsoluteTolerance(self):
618     """
619     returns the absolute tolerance
620    
621 jfenwick 2625 :return: current absolute tolerance
622     :rtype: ``float``
623 gross 2264 """
624     return self.__atol
625     def getSubProblemTolerance(self):
626 gross 2474 """
627     Returns a suitable subtolerance
628 jfenwick 2625 @type: ``float``
629 gross 2474 """
630     return max(util.EPSILON**(0.75),self.getTolerance()**2)
631     def setSubProblemTolerance(self):
632 gross 2264 """
633 gross 2474 Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected.
634 gross 2264 """
635 gross 2474 if self.__adaptSubTolerance:
636     sub_tol=self.getSubProblemTolerance()
637     self.getSolverOptionsFlux().setTolerance(sub_tol)
638     self.getSolverOptionsFlux().setAbsoluteTolerance(0.)
639     self.getSolverOptionsPressure().setTolerance(sub_tol)
640     self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
641     if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol
642 gross 2208
643 gross 2474 def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):
644 gross 2264 """
645     solves the problem.
646 gross 1659
647 gross 2208 The iteration is terminated if the residual norm is less then self.getTolerance().
648 gross 1659
649 jfenwick 2625 :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.
650     :type u0: vector value on the domain (e.g. `Data`).
651     :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.
652     :type p0: scalar value on the domain (e.g. `Data`).
653     :param verbose: if set some information on iteration progress are printed
654     :type verbose: ``bool``
655     :return: flux and pressure
656     :rtype: ``tuple`` of `Data`.
657 gross 2264
658 jfenwick 2625 :note: The problem is solved as a least squares form
659 gross 2100
660 jfenwick 2625 *(I+D^*D)u+Qp=D^*f+g*
661     *Q^*u+Q^*Qp=Q^*g*
662 gross 2100
663 jfenwick 2625 where *D* is the *div* operator and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.
664 gross 2264 We eliminate the flux form the problem by setting
665 caltinay 2169
666 jfenwick 2625 *u=(I+D^*D)^{-1}(D^*f-g-Qp)* with u=u0 on location_of_fixed_flux
667 caltinay 2169
668 gross 2208 form the first equation. Inserted into the second equation we get
669 caltinay 2169
670 jfenwick 2625 *Q^*(I-(I+D^*D)^{-1})Qp= Q^*(g-(I+D^*D)^{-1}(D^*f+g))* with p=p0 on location_of_fixed_pressure
671 gross 2264
672 jfenwick 2625 which is solved using the PCG method (precondition is *Q^*Q*). In each iteration step
673     PDEs with operator *I+D^*D* and with *Q^*Q* needs to be solved using a sub iteration scheme.
674 gross 2208 """
675 gross 2386 self.verbose=verbose
676 gross 2264 rtol=self.getTolerance()
677     atol=self.getAbsoluteTolerance()
678 gross 2474 self.setSubProblemTolerance()
679 gross 2264 num_corrections=0
680     converged=False
681     p=p0
682     norm_r=None
683     while not converged:
684 gross 2474 v=self.getFlux(p, fixed_flux=u0)
685 gross 2264 Qp=self.__Q(p)
686     norm_v=self.__L2(v)
687     norm_Qp=self.__L2(Qp)
688     if norm_v == 0.:
689     if norm_Qp == 0.:
690     return v,p
691     else:
692     fac=norm_Qp
693     else:
694     if norm_Qp == 0.:
695     fac=norm_v
696     else:
697     fac=2./(1./norm_v+1./norm_Qp)
698     ATOL=(atol+rtol*fac)
699     if self.verbose:
700     print "DarcyFlux: L2 norm of v = %e."%norm_v
701 gross 2960 print "DarcyFlux: L2 norm of k.util.grad(p) = %e."%norm_Qp
702 gross 2486 print "DarcyFlux: L2 defect u = %e."%(util.integrate(util.length(self.__g-util.interpolate(v,Function(self.domain))-Qp)**2)**(0.5),)
703     print "DarcyFlux: L2 defect div(v) = %e."%(util.integrate((self.__f-util.div(v))**2)**(0.5),)
704 gross 2264 print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
705     if norm_r == None or norm_r>ATOL:
706     if num_corrections>max_num_corrections:
707     raise ValueError,"maximum number of correction steps reached."
708 gross 2354 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)
709 gross 2264 num_corrections+=1
710     else:
711     converged=True
712     return v,p
713     def __L2(self,v):
714     return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))
715    
716     def __Q(self,p):
717     return util.tensor_mult(self.__permeability,util.grad(p))
718    
719     def __Aprod(self,dp):
720 gross 2474 if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
721 gross 2264 Qdp=self.__Q(dp)
722     self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())
723 gross 2474 du=self.__pde_v.getSolution()
724 gross 2370 # self.__pde_v.getOperator().saveMM("proj.mm")
725 gross 2264 return Qdp+du
726     def __inner_GMRES(self,r,s):
727     return util.integrate(util.inner(r,s))
728    
729 gross 2100 def __inner_PCG(self,p,r):
730 gross 2264 return util.integrate(util.inner(self.__Q(p), r))
731 gross 2100
732     def __Msolve_PCG(self,r):
733 gross 2474 if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"
734 gross 2264 self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=Data(), r=Data())
735 gross 2370 # self.__pde_p.getOperator().saveMM("prec.mm")
736 gross 2474 return self.__pde_p.getSolution()
737 gross 2100
738 gross 2474 def getFlux(self,p=None, fixed_flux=Data()):
739 gross 2264 """
740 jfenwick 2625 returns the flux for a given pressure ``p`` where the flux is equal to ``fixed_flux``
741     on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
742     Note that ``g`` and ``f`` are used, see `setValue`.
743 gross 2264
744 jfenwick 2625 :param p: pressure.
745     :type p: scalar value on the domain (e.g. `Data`).
746     :param fixed_flux: flux on the locations of the domain marked be ``location_of_fixed_flux``.
747     :type fixed_flux: vector values on the domain (e.g. `Data`).
748     :return: flux
749     :rtype: `Data`
750     :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}*
751     for the permeability *k_{ij}*
752 gross 2264 """
753 gross 2474 self.setSubProblemTolerance()
754 gross 2264 g=self.__g
755     f=self.__f
756 gross 2267 self.__pde_v.setValue(X=self.__l*f*util.kronecker(self.domain), r=fixed_flux)
757 gross 2264 if p == None:
758     self.__pde_v.setValue(Y=g)
759     else:
760     self.__pde_v.setValue(Y=g-self.__Q(p))
761 gross 2474 return self.__pde_v.getSolution()
762 gross 2264
763 gross 1414 class StokesProblemCartesian(HomogeneousSaddlePointProblem):
764 gross 2251 """
765 gross 2264 solves
766 gross 1414
767 gross 2208 -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
768     u_{i,i}=0
769 gross 1414
770 gross 2208 u=0 where fixed_u_mask>0
771     eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j
772 gross 1414
773 gross 2264 if surface_stress is not given 0 is assumed.
774 gross 1414
775 gross 2251 typical usage:
776 gross 1414
777 gross 2208 sp=StokesProblemCartesian(domain)
778     sp.setTolerance()
779     sp.initialize(...)
780     v,p=sp.solve(v0,p0)
781 gross 2251 """
782 gross 2719 def __init__(self,domain,**kwargs):
783 gross 2100 """
784 gross 2208 initialize the Stokes Problem
785 gross 2100
786 gross 2793 The approximation spaces used for velocity (=Solution(domain)) and pressure (=ReducedSolution(domain)) must be
787     LBB complient, for instance using quadratic and linear approximation on the same element or using linear approximation
788     with macro elements for the pressure.
789    
790     :param domain: domain of the problem.
791 jfenwick 2625 :type domain: `Domain`
792 gross 2100 """
793 gross 2719 HomogeneousSaddlePointProblem.__init__(self,**kwargs)
794 gross 1414 self.domain=domain
795     self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
796     self.__pde_u.setSymmetryOn()
797 gross 2474
798 gross 1414 self.__pde_prec=LinearPDE(domain)
799     self.__pde_prec.setReducedOrderOn()
800     self.__pde_prec.setSymmetryOn()
801    
802 gross 2415 self.__pde_proj=LinearPDE(domain)
803     self.__pde_proj.setReducedOrderOn()
804     self.__pde_proj.setValue(D=1)
805     self.__pde_proj.setSymmetryOn()
806    
807 gross 2474 def getSolverOptionsVelocity(self):
808     """
809     returns the solver options used solve the equation for velocity.
810    
811 jfenwick 2625 :rtype: `SolverOptions`
812 gross 2474 """
813     return self.__pde_u.getSolverOptions()
814     def setSolverOptionsVelocity(self, options=None):
815     """
816     set the solver options for solving the equation for velocity.
817    
818 jfenwick 2625 :param options: new solver options
819     :type options: `SolverOptions`
820 gross 2474 """
821     self.__pde_u.setSolverOptions(options)
822     def getSolverOptionsPressure(self):
823     """
824     returns the solver options used solve the equation for pressure.
825 jfenwick 2625 :rtype: `SolverOptions`
826 gross 2474 """
827     return self.__pde_prec.getSolverOptions()
828     def setSolverOptionsPressure(self, options=None):
829     """
830     set the solver options for solving the equation for pressure.
831 jfenwick 2625 :param options: new solver options
832     :type options: `SolverOptions`
833 gross 2474 """
834     self.__pde_prec.setSolverOptions(options)
835 gross 2415
836 gross 2474 def setSolverOptionsDiv(self, options=None):
837     """
838     set the solver options for solving the equation to project the divergence of
839     the velocity onto the function space of presure.
840    
841 jfenwick 2625 :param options: new solver options
842     :type options: `SolverOptions`
843 gross 2474 """
844 artak 2689 self.__pde_proj.setSolverOptions(options)
845 gross 2474 def getSolverOptionsDiv(self):
846     """
847     returns the solver options for solving the equation to project the divergence of
848     the velocity onto the function space of presure.
849    
850 jfenwick 2625 :rtype: `SolverOptions`
851 gross 2474 """
852 artak 2689 return self.__pde_proj.getSolverOptions()
853 gross 2474
854 gross 2793 def updateStokesEquation(self, v, p):
855     """
856     updates the Stokes equation to consider dependencies from ``v`` and ``p``
857     :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values.
858     """
859     pass
860     def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None):
861     """
862     assigns new values to the model parameters.
863    
864     :param f: external force
865     :type f: `Vector` object in `FunctionSpace` `Function` or similar
866     :param fixed_u_mask: mask of locations with fixed velocity.
867     :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
868     :param eta: viscosity
869     :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
870     :param surface_stress: normal surface stress
871     :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
872     :param stress: initial stress
873     :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
874     """
875     if eta !=None:
876     k=util.kronecker(self.domain.getDim())
877     kk=util.outer(k,k)
878     self.eta=util.interpolate(eta, Function(self.domain))
879     self.__pde_prec.setValue(D=1/self.eta)
880     self.__pde_u.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))
881     if restoration_factor!=None:
882     n=self.domain.getNormal()
883     self.__pde_u.setValue(d=restoration_factor*util.outer(n,n))
884     if fixed_u_mask!=None:
885     self.__pde_u.setValue(q=fixed_u_mask)
886     if f!=None: self.__f=f
887     if surface_stress!=None: self.__surface_stress=surface_stress
888     if stress!=None: self.__stress=stress
889    
890 gross 2620 def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data(), restoration_factor=0):
891 gross 2208 """
892     assigns values to the model parameters
893 gross 2100
894 jfenwick 2625 :param f: external force
895     :type f: `Vector` object in `FunctionSpace` `Function` or similar
896     :param fixed_u_mask: mask of locations with fixed velocity.
897     :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
898     :param eta: viscosity
899     :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
900     :param surface_stress: normal surface stress
901     :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
902     :param stress: initial stress
903     :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
904 gross 2208 """
905 gross 2793 self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
906 gross 1414
907 gross 2719 def Bv(self,v,tol):
908 gross 2251 """
909     returns inner product of element p and div(v)
910 gross 1414
911 jfenwick 2625 :param v: a residual
912     :return: inner product of element p and div(v)
913     :rtype: ``float``
914 gross 2100 """
915 gross 2793 self.__pde_proj.setValue(Y=-util.div(v))
916 gross 2719 self.getSolverOptionsDiv().setTolerance(tol)
917     self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
918     out=self.__pde_proj.getSolution()
919     return out
920 gross 2208
921 gross 2445 def inner_pBv(self,p,Bv):
922     """
923     returns inner product of element p and Bv=-div(v)
924    
925 jfenwick 2625 :param p: a pressure increment
926     :param Bv: a residual
927     :return: inner product of element p and Bv=-div(v)
928     :rtype: ``float``
929 gross 2445 """
930     return util.integrate(util.interpolate(p,Function(self.domain))*util.interpolate(Bv,Function(self.domain)))
931    
932 gross 2251 def inner_p(self,p0,p1):
933 gross 2100 """
934 gross 2251 Returns inner product of p0 and p1
935 gross 1414
936 jfenwick 2625 :param p0: a pressure
937     :param p1: a pressure
938     :return: inner product of p0 and p1
939     :rtype: ``float``
940 gross 2100 """
941 gross 2719 s0=util.interpolate(p0,Function(self.domain))
942     s1=util.interpolate(p1,Function(self.domain))
943 gross 2100 return util.integrate(s0*s1)
944 artak 1550
945 gross 2251 def norm_v(self,v):
946 gross 2100 """
947 gross 2251 returns the norm of v
948 gross 2208
949 jfenwick 2625 :param v: a velovity
950     :return: norm of v
951     :rtype: non-negative ``float``
952 gross 2100 """
953 gross 2719 return util.sqrt(util.integrate(util.length(util.grad(v))**2))
954 gross 2100
955 gross 2793
956 gross 2719 def getDV(self, p, v, tol):
957 gross 1414 """
958 gross 2251 return the value for v for a given p (overwrite)
959    
960 jfenwick 2625 :param p: a pressure
961 gross 2719 :param v: a initial guess for the value v to return.
962     :return: dv given as *Adv=(f-Av-B^*p)*
963 gross 1414 """
964 gross 2793 self.updateStokesEquation(v,p)
965 gross 2719 self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress)
966     self.getSolverOptionsVelocity().setTolerance(tol)
967     self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
968 gross 2100 if self.__stress.isEmpty():
969 gross 2719 self.__pde_u.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
970 gross 2100 else:
971 gross 2719 self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
972 gross 2474 out=self.__pde_u.getSolution()
973 gross 2208 return out
974 gross 1414
975 gross 2445 def norm_Bv(self,Bv):
976 gross 2251 """
977     Returns Bv (overwrite).
978    
979 jfenwick 2625 :rtype: equal to the type of p
980     :note: boundary conditions on p should be zero!
981 gross 2251 """
982 gross 2445 return util.sqrt(util.integrate(util.interpolate(Bv,Function(self.domain))**2))
983 gross 2251
984 gross 2719 def solve_AinvBt(self,p, tol):
985 gross 2251 """
986 gross 2719 Solves *Av=B^*p* with accuracy `tol`
987 gross 2251
988 jfenwick 2625 :param p: a pressure increment
989     :return: the solution of *Av=B^*p*
990     :note: boundary conditions on v should be zero!
991 gross 2251 """
992 gross 2719 self.__pde_u.setValue(Y=Data(), y=Data(), X=-p*util.kronecker(self.domain))
993 gross 2474 out=self.__pde_u.getSolution()
994 gross 2251 return out
995    
996 gross 2719 def solve_prec(self,Bv, tol):
997 gross 2251 """
998 jfenwick 2625 applies preconditioner for for *BA^{-1}B^** to *Bv*
999     with accuracy `self.getSubProblemTolerance()`
1000 gross 2251
1001 jfenwick 2625 :param Bv: velocity increment
1002     :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*
1003     :note: boundary conditions on p are zero.
1004 gross 2251 """
1005 gross 2445 self.__pde_prec.setValue(Y=Bv)
1006 gross 2719 self.getSolverOptionsPressure().setTolerance(tol)
1007     self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
1008     out=self.__pde_prec.getSolution()
1009     return out

  ViewVC Help
Powered by ViewVC 1.1.26