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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3360 - (hide annotations)
Thu Nov 18 00:20:21 2010 UTC (8 years, 10 months ago) by jfenwick
File MIME type: text/x-python
File size: 40280 byte(s)
Fix some epydoc errors
Fixed issue 564
removed minimize.py since it hasn't worked in who knows how long
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 jfenwick 3360
160 gross 3071 :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 jfenwick 3360 :param permeability: permeability tensor. If scalar ``s`` is given the tensor with ``s`` on the main diagonal is used.
169 gross 3071 :type permeability: scalar or tensor values on the domain (e.g. `Data`)
170 jfenwick 3360
171 gross 3071 :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 jfenwick 3360 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 gross 3071 """
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 jfenwick 3360 *G^*u+*G^* K Gp=G^*g*
310 gross 3071 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 gross 2960
327 gross 3071 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 gross 2960
333 gross 3071 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 gross 2960
378 gross 3071 *|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 gross 2960
396 gross 3071 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 gross 2960
402    
403 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}*.
404 gross 2960
405 gross 3071 :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 gross 3074 return max(util.EPSILON**(0.5),self.getTolerance()**2)
424 gross 2960
425 gross 3071 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 gross 2960
439    
440     class DarcyFlowOld(object):
441     """
442     solves the problem
443    
444     *u_i+k_{ij}*p_{,j} = g_i*
445     *u_{i,i} = f*
446    
447     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.
450     """
451    
452     def __init__(self, domain, weight=None, useReduced=False, adaptSubTolerance=True):
453     """
454     initializes the Darcy flux problem
455     :param domain: domain of the problem
456     :type domain: `Domain`
457     :param useReduced: uses reduced oreder on flux and pressure
458     :type useReduced: ``bool``
459     :param adaptSubTolerance: switches on automatic subtolerance selection
460     :type adaptSubTolerance: ``bool``
461     """
462     self.domain=domain
463     if weight == None:
464 gross 2354 s=self.domain.getSize()
465     self.__l=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2
466 gross 2486 # self.__l=(3.*util.longestEdge(self.domain))**2
467 artak 2689 #self.__l=(0.1*util.longestEdge(self.domain)*s/util.sup(s))**2
468 gross 2351 else:
469     self.__l=weight
470 gross 2100 self.__pde_v=LinearPDESystem(domain)
471 gross 2208 if useReduced: self.__pde_v.setReducedOrderOn()
472     self.__pde_v.setSymmetryOn()
473 gross 2267 self.__pde_v.setValue(D=util.kronecker(domain), A=self.__l*util.outer(util.kronecker(domain),util.kronecker(domain)))
474 gross 2100 self.__pde_p=LinearSinglePDE(domain)
475     self.__pde_p.setSymmetryOn()
476 gross 2208 if useReduced: self.__pde_p.setReducedOrderOn()
477 gross 2100 self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
478     self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
479 gross 2264 self.setTolerance()
480     self.setAbsoluteTolerance()
481 gross 2474 self.__adaptSubTolerance=adaptSubTolerance
482     self.verbose=False
483     def getSolverOptionsFlux(self):
484     """
485     Returns the solver options used to solve the flux problems
486    
487 jfenwick 2625 *(I+D^*D)u=F*
488 gross 2474
489 jfenwick 2625 :return: `SolverOptions`
490 gross 2474 """
491     return self.__pde_v.getSolverOptions()
492     def setSolverOptionsFlux(self, options=None):
493     """
494     Sets the solver options used to solve the flux problems
495    
496 jfenwick 2625 *(I+D^*D)u=F*
497 gross 2474
498 jfenwick 2625 If ``options`` is not present, the options are reset to default
499     :param options: `SolverOptions`
500     :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
501 gross 2474 """
502     return self.__pde_v.setSolverOptions(options)
503     def getSolverOptionsPressure(self):
504     """
505     Returns the solver options used to solve the pressure problems
506    
507 jfenwick 2625 *(Q^*Q)p=Q^*G*
508 gross 2474
509 jfenwick 2625 :return: `SolverOptions`
510 gross 2474 """
511     return self.__pde_p.getSolverOptions()
512     def setSolverOptionsPressure(self, options=None):
513     """
514     Sets the solver options used to solve the pressure problems
515    
516 jfenwick 2625 *(Q^*Q)p=Q^*G*
517 gross 2474
518 jfenwick 2625 If ``options`` is not present, the options are reset to default
519     :param options: `SolverOptions`
520     :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
521 gross 2474 """
522     return self.__pde_p.setSolverOptions(options)
523 gross 1659
524 gross 2100 def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
525     """
526 gross 2208 assigns values to model parameters
527 gross 1659
528 jfenwick 2625 :param f: volumetic sources/sinks
529     :type f: scalar value on the domain (e.g. `Data`)
530     :param g: flux sources/sinks
531     :type g: vector values on the domain (e.g. `Data`)
532     :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. `Data`)
534     :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. `Data`)
536     :param permeability: permeability tensor. If scalar ``s`` is given the tensor with
537     ``s`` on the main diagonal is used. If vector ``v`` is given the tensor with
538     ``v`` on the main diagonal is used.
539     :type permeability: scalar, vector or tensor values on the domain (e.g. `Data`)
540 gross 1659
541 jfenwick 2625 :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 (``location_of_fixed_pressure`` >0)
543     or the normal component of the flux (``location_of_fixed_flux[i]>0`` if direction of the normal
544     is along the *x_i* axis.
545 gross 2100 """
546 gross 2264 if f !=None:
547 gross 2100 f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
548     if f.isEmpty():
549     f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
550     else:
551     if f.getRank()>0: raise ValueError,"illegal rank of f."
552 gross 2267 self.__f=f
553 gross 2264 if g !=None:
554 gross 2100 g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
555     if g.isEmpty():
556     g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
557     else:
558     if not g.getShape()==(self.domain.getDim(),):
559     raise ValueError,"illegal shape of g"
560     self.__g=g
561 gross 1659
562 gross 2100 if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)
563     if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux)
564 gross 1659
565 gross 2100 if permeability!=None:
566     perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
567     if perm.getRank()==0:
568     perm=perm*util.kronecker(self.domain.getDim())
569     elif perm.getRank()==1:
570     perm, perm2=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A")), perm
571     for i in range(self.domain.getDim()): perm[i,i]=perm2[i]
572     elif perm.getRank()==2:
573     pass
574     else:
575     raise ValueError,"illegal rank of permeability."
576     self.__permeability=perm
577     self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))
578 gross 1659
579 gross 2264 def setTolerance(self,rtol=1e-4):
580     """
581 jfenwick 2625 sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if
582 gross 1659
583 jfenwick 2625 *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
584 gross 2264
585 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}*.
586 gross 2264
587 jfenwick 2625 :param rtol: relative tolerance for the pressure
588     :type rtol: non-negative ``float``
589 gross 2100 """
590 gross 2264 if rtol<0:
591     raise ValueError,"Relative tolerance needs to be non-negative."
592     self.__rtol=rtol
593     def getTolerance(self):
594 gross 2100 """
595 gross 2264 returns the relative tolerance
596 gross 1659
597 jfenwick 2625 :return: current relative tolerance
598     :rtype: ``float``
599 caltinay 2169 """
600 gross 2264 return self.__rtol
601 gross 1659
602 gross 2264 def setAbsoluteTolerance(self,atol=0.):
603 gross 2208 """
604 jfenwick 2625 sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if
605 gross 1659
606 jfenwick 2625 *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
607 gross 2208
608 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}*.
609 gross 2208
610 jfenwick 2625 :param atol: absolute tolerance for the pressure
611     :type atol: non-negative ``float``
612 gross 2264 """
613     if atol<0:
614     raise ValueError,"Absolute tolerance needs to be non-negative."
615     self.__atol=atol
616     def getAbsoluteTolerance(self):
617     """
618     returns the absolute tolerance
619    
620 jfenwick 2625 :return: current absolute tolerance
621     :rtype: ``float``
622 gross 2264 """
623     return self.__atol
624     def getSubProblemTolerance(self):
625 gross 2474 """
626     Returns a suitable subtolerance
627 jfenwick 2625 @type: ``float``
628 gross 2474 """
629     return max(util.EPSILON**(0.75),self.getTolerance()**2)
630     def setSubProblemTolerance(self):
631 gross 2264 """
632 gross 2474 Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected.
633 gross 2264 """
634 gross 2474 if self.__adaptSubTolerance:
635     sub_tol=self.getSubProblemTolerance()
636     self.getSolverOptionsFlux().setTolerance(sub_tol)
637     self.getSolverOptionsFlux().setAbsoluteTolerance(0.)
638     self.getSolverOptionsPressure().setTolerance(sub_tol)
639     self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
640     if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol
641 gross 2208
642 gross 2474 def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):
643 gross 2264 """
644     solves the problem.
645 gross 1659
646 gross 2208 The iteration is terminated if the residual norm is less then self.getTolerance().
647 gross 1659
648 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.
649     :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 ``location_of_fixed_pressure`` the value of ``p0`` is kept unchanged.
651     :type p0: scalar value on the domain (e.g. `Data`).
652     :param verbose: if set some information on iteration progress are printed
653     :type verbose: ``bool``
654     :return: flux and pressure
655     :rtype: ``tuple`` of `Data`.
656 gross 2264
657 jfenwick 2625 :note: The problem is solved as a least squares form
658 gross 2100
659 jfenwick 2625 *(I+D^*D)u+Qp=D^*f+g*
660     *Q^*u+Q^*Qp=Q^*g*
661 gross 2100
662 jfenwick 2625 where *D* is the *div* operator and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.
663 gross 2264 We eliminate the flux form the problem by setting
664 caltinay 2169
665 jfenwick 2625 *u=(I+D^*D)^{-1}(D^*f-g-Qp)* with u=u0 on location_of_fixed_flux
666 caltinay 2169
667 gross 2208 form the first equation. Inserted into the second equation we get
668 caltinay 2169
669 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
670 gross 2264
671 jfenwick 2625 which is solved using the PCG method (precondition is *Q^*Q*). In each iteration step
672     PDEs with operator *I+D^*D* and with *Q^*Q* needs to be solved using a sub iteration scheme.
673 gross 2208 """
674 gross 2386 self.verbose=verbose
675 gross 2264 rtol=self.getTolerance()
676     atol=self.getAbsoluteTolerance()
677 gross 2474 self.setSubProblemTolerance()
678 gross 2264 num_corrections=0
679     converged=False
680     p=p0
681     norm_r=None
682     while not converged:
683 gross 2474 v=self.getFlux(p, fixed_flux=u0)
684 gross 2264 Qp=self.__Q(p)
685     norm_v=self.__L2(v)
686     norm_Qp=self.__L2(Qp)
687     if norm_v == 0.:
688     if norm_Qp == 0.:
689     return v,p
690     else:
691     fac=norm_Qp
692     else:
693     if norm_Qp == 0.:
694     fac=norm_v
695     else:
696     fac=2./(1./norm_v+1./norm_Qp)
697     ATOL=(atol+rtol*fac)
698     if self.verbose:
699     print "DarcyFlux: L2 norm of v = %e."%norm_v
700 gross 2960 print "DarcyFlux: L2 norm of k.util.grad(p) = %e."%norm_Qp
701 gross 2486 print "DarcyFlux: L2 defect u = %e."%(util.integrate(util.length(self.__g-util.interpolate(v,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),)
703 gross 2264 print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
704     if norm_r == None or norm_r>ATOL:
705     if num_corrections>max_num_corrections:
706     raise ValueError,"maximum number of correction steps reached."
707 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)
708 gross 2264 num_corrections+=1
709     else:
710     converged=True
711     return v,p
712     def __L2(self,v):
713     return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))
714    
715     def __Q(self,p):
716     return util.tensor_mult(self.__permeability,util.grad(p))
717    
718     def __Aprod(self,dp):
719 gross 2474 if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
720 gross 2264 Qdp=self.__Q(dp)
721     self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())
722 gross 2474 du=self.__pde_v.getSolution()
723 gross 2370 # self.__pde_v.getOperator().saveMM("proj.mm")
724 gross 2264 return Qdp+du
725     def __inner_GMRES(self,r,s):
726     return util.integrate(util.inner(r,s))
727    
728 gross 2100 def __inner_PCG(self,p,r):
729 gross 2264 return util.integrate(util.inner(self.__Q(p), r))
730 gross 2100
731     def __Msolve_PCG(self,r):
732 gross 2474 if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"
733 gross 2264 self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=Data(), r=Data())
734 gross 2370 # self.__pde_p.getOperator().saveMM("prec.mm")
735 gross 2474 return self.__pde_p.getSolution()
736 gross 2100
737 gross 2474 def getFlux(self,p=None, fixed_flux=Data()):
738 gross 2264 """
739 jfenwick 2625 returns the flux for a given pressure ``p`` where the flux is equal to ``fixed_flux``
740     on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
741     Note that ``g`` and ``f`` are used, see `setValue`.
742 gross 2264
743 jfenwick 2625 :param p: pressure.
744     :type p: scalar value on the domain (e.g. `Data`).
745     :param fixed_flux: flux on the locations of the domain marked be ``location_of_fixed_flux``.
746     :type fixed_flux: vector values on the domain (e.g. `Data`).
747     :return: flux
748     :rtype: `Data`
749     :note: the method uses the least squares solution *u=(I+D^*D)^{-1}(D^*f-g-Qp)* where *D* is the *div* operator and *(Qp)_i=k_{ij}p_{,j}*
750     for the permeability *k_{ij}*
751 gross 2264 """
752 gross 2474 self.setSubProblemTolerance()
753 gross 2264 g=self.__g
754     f=self.__f
755 gross 2267 self.__pde_v.setValue(X=self.__l*f*util.kronecker(self.domain), r=fixed_flux)
756 gross 2264 if p == None:
757     self.__pde_v.setValue(Y=g)
758     else:
759     self.__pde_v.setValue(Y=g-self.__Q(p))
760 gross 2474 return self.__pde_v.getSolution()
761 gross 2264
762 gross 1414 class StokesProblemCartesian(HomogeneousSaddlePointProblem):
763 gross 2251 """
764 gross 2264 solves
765 gross 1414
766 gross 2208 -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
767     u_{i,i}=0
768 gross 1414
769 gross 2208 u=0 where fixed_u_mask>0
770     eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j
771 gross 1414
772 gross 2264 if surface_stress is not given 0 is assumed.
773 gross 1414
774 gross 2251 typical usage:
775 gross 1414
776 gross 2208 sp=StokesProblemCartesian(domain)
777     sp.setTolerance()
778     sp.initialize(...)
779     v,p=sp.solve(v0,p0)
780 gross 2251 """
781 gross 2719 def __init__(self,domain,**kwargs):
782 gross 2100 """
783 gross 2208 initialize the Stokes Problem
784 gross 2100
785 gross 2793 The approximation spaces used for velocity (=Solution(domain)) and pressure (=ReducedSolution(domain)) must be
786     LBB complient, for instance using quadratic and linear approximation on the same element or using linear approximation
787     with macro elements for the pressure.
788    
789     :param domain: domain of the problem.
790 jfenwick 2625 :type domain: `Domain`
791 gross 2100 """
792 gross 2719 HomogeneousSaddlePointProblem.__init__(self,**kwargs)
793 gross 1414 self.domain=domain
794     self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
795     self.__pde_u.setSymmetryOn()
796 gross 2474
797 gross 1414 self.__pde_prec=LinearPDE(domain)
798     self.__pde_prec.setReducedOrderOn()
799     self.__pde_prec.setSymmetryOn()
800    
801 gross 2415 self.__pde_proj=LinearPDE(domain)
802     self.__pde_proj.setReducedOrderOn()
803     self.__pde_proj.setValue(D=1)
804     self.__pde_proj.setSymmetryOn()
805    
806 gross 2474 def getSolverOptionsVelocity(self):
807     """
808     returns the solver options used solve the equation for velocity.
809    
810 jfenwick 2625 :rtype: `SolverOptions`
811 gross 2474 """
812     return self.__pde_u.getSolverOptions()
813     def setSolverOptionsVelocity(self, options=None):
814     """
815     set the solver options for solving the equation for velocity.
816    
817 jfenwick 2625 :param options: new solver options
818     :type options: `SolverOptions`
819 gross 2474 """
820     self.__pde_u.setSolverOptions(options)
821     def getSolverOptionsPressure(self):
822     """
823     returns the solver options used solve the equation for pressure.
824 jfenwick 2625 :rtype: `SolverOptions`
825 gross 2474 """
826     return self.__pde_prec.getSolverOptions()
827     def setSolverOptionsPressure(self, options=None):
828     """
829     set the solver options for solving the equation for pressure.
830 jfenwick 2625 :param options: new solver options
831     :type options: `SolverOptions`
832 gross 2474 """
833     self.__pde_prec.setSolverOptions(options)
834 gross 2415
835 gross 2474 def setSolverOptionsDiv(self, options=None):
836     """
837     set the solver options for solving the equation to project the divergence of
838     the velocity onto the function space of presure.
839    
840 jfenwick 2625 :param options: new solver options
841     :type options: `SolverOptions`
842 gross 2474 """
843 artak 2689 self.__pde_proj.setSolverOptions(options)
844 gross 2474 def getSolverOptionsDiv(self):
845     """
846     returns the solver options for solving the equation to project the divergence of
847     the velocity onto the function space of presure.
848    
849 jfenwick 2625 :rtype: `SolverOptions`
850 gross 2474 """
851 artak 2689 return self.__pde_proj.getSolverOptions()
852 gross 2474
853 gross 2793 def updateStokesEquation(self, v, p):
854     """
855     updates the Stokes equation to consider dependencies from ``v`` and ``p``
856     :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values.
857     """
858     pass
859     def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None):
860     """
861     assigns new values to the model parameters.
862    
863     :param f: external force
864     :type f: `Vector` object in `FunctionSpace` `Function` or similar
865     :param fixed_u_mask: mask of locations with fixed velocity.
866     :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
867     :param eta: viscosity
868     :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
869     :param surface_stress: normal surface stress
870     :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
871     :param stress: initial stress
872     :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
873     """
874     if eta !=None:
875     k=util.kronecker(self.domain.getDim())
876     kk=util.outer(k,k)
877     self.eta=util.interpolate(eta, 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 gross 2620 def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data(), restoration_factor=0):
890 gross 2208 """
891     assigns values to the model parameters
892 gross 2100
893 jfenwick 2625 :param f: external force
894     :type f: `Vector` object in `FunctionSpace` `Function` or similar
895     :param fixed_u_mask: mask of locations with fixed velocity.
896     :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
897     :param eta: viscosity
898     :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
899     :param surface_stress: normal surface stress
900     :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
901     :param stress: initial stress
902     :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
903 gross 2208 """
904 gross 2793 self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
905 gross 1414
906 gross 2719 def Bv(self,v,tol):
907 gross 2251 """
908     returns inner product of element p and div(v)
909 gross 1414
910 jfenwick 2625 :param v: a residual
911     :return: inner product of element p and div(v)
912     :rtype: ``float``
913 gross 2100 """
914 gross 2793 self.__pde_proj.setValue(Y=-util.div(v))
915 gross 2719 self.getSolverOptionsDiv().setTolerance(tol)
916     self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
917     out=self.__pde_proj.getSolution()
918     return out
919 gross 2208
920 gross 2445 def inner_pBv(self,p,Bv):
921     """
922     returns inner product of element p and Bv=-div(v)
923    
924 jfenwick 2625 :param p: a pressure increment
925     :param Bv: a residual
926     :return: inner product of element p and Bv=-div(v)
927     :rtype: ``float``
928 gross 2445 """
929     return util.integrate(util.interpolate(p,Function(self.domain))*util.interpolate(Bv,Function(self.domain)))
930    
931 gross 2251 def inner_p(self,p0,p1):
932 gross 2100 """
933 gross 2251 Returns inner product of p0 and p1
934 gross 1414
935 jfenwick 2625 :param p0: a pressure
936     :param p1: a pressure
937     :return: inner product of p0 and p1
938     :rtype: ``float``
939 gross 2100 """
940 gross 2719 s0=util.interpolate(p0,Function(self.domain))
941     s1=util.interpolate(p1,Function(self.domain))
942 gross 2100 return util.integrate(s0*s1)
943 artak 1550
944 gross 2251 def norm_v(self,v):
945 gross 2100 """
946 gross 2251 returns the norm of v
947 gross 2208
948 jfenwick 2625 :param v: a velovity
949     :return: norm of v
950     :rtype: non-negative ``float``
951 gross 2100 """
952 gross 2719 return util.sqrt(util.integrate(util.length(util.grad(v))**2))
953 gross 2100
954 gross 2793
955 gross 2719 def getDV(self, p, v, tol):
956 gross 1414 """
957 gross 2251 return the value for v for a given p (overwrite)
958    
959 jfenwick 2625 :param p: a pressure
960 gross 2719 :param v: a initial guess for the value v to return.
961     :return: dv given as *Adv=(f-Av-B^*p)*
962 gross 1414 """
963 gross 2793 self.updateStokesEquation(v,p)
964 gross 2719 self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress)
965     self.getSolverOptionsVelocity().setTolerance(tol)
966     self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
967 gross 2100 if self.__stress.isEmpty():
968 gross 2719 self.__pde_u.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
969 gross 2100 else:
970 gross 2719 self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
971 gross 2474 out=self.__pde_u.getSolution()
972 gross 2208 return out
973 gross 1414
974 gross 2445 def norm_Bv(self,Bv):
975 gross 2251 """
976     Returns Bv (overwrite).
977    
978 jfenwick 2625 :rtype: equal to the type of p
979     :note: boundary conditions on p should be zero!
980 gross 2251 """
981 gross 2445 return util.sqrt(util.integrate(util.interpolate(Bv,Function(self.domain))**2))
982 gross 2251
983 gross 2719 def solve_AinvBt(self,p, tol):
984 gross 2251 """
985 gross 2719 Solves *Av=B^*p* with accuracy `tol`
986 gross 2251
987 jfenwick 2625 :param p: a pressure increment
988     :return: the solution of *Av=B^*p*
989     :note: boundary conditions on v should be zero!
990 gross 2251 """
991 gross 2719 self.__pde_u.setValue(Y=Data(), y=Data(), X=-p*util.kronecker(self.domain))
992 gross 2474 out=self.__pde_u.getSolution()
993 gross 2251 return out
994    
995 gross 2719 def solve_prec(self,Bv, tol):
996 gross 2251 """
997 jfenwick 2625 applies preconditioner for for *BA^{-1}B^** to *Bv*
998     with accuracy `self.getSubProblemTolerance()`
999 gross 2251
1000 jfenwick 2625 :param Bv: velocity increment
1001     :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*
1002     :note: boundary conditions on p are zero.
1003 gross 2251 """
1004 gross 2445 self.__pde_prec.setValue(Y=Bv)
1005 gross 2719 self.getSolverOptionsPressure().setTolerance(tol)
1006     self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
1007     out=self.__pde_prec.getSolution()
1008     return out

  ViewVC Help
Powered by ViewVC 1.1.26