/[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 3569 - (hide annotations)
Thu Sep 1 02:42:36 2011 UTC (8 years ago) by gross
File MIME type: text/x-python
File size: 20166 byte(s)
Darcy flow solver and docu reviewed. Coupled solvers have been removed.

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 jfenwick 3432 import escript
36 gross 1414 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 3502 class DarcyFlow(object):
41 gross 3499 """
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 gross 3569 :cvar EVAL: direct pressure gradient evaluation for flux
50     :cvar POST: global postprocessing of flux by solving the PDE *K_{ij} u_j + (w * K * l u_{k,k})_{,i}= - p_{,j} + K_{ij} g_j*
51     where *l* is the length scale, *K* is the inverse of the permeability tensor, and *w* is a positive weighting factor.
52     :cvar SMOOTH: global smoothing by solving the PDE *K_{ij} u_j= - p_{,j} + K_{ij} g_j*
53 gross 3499 """
54 gross 3569 EVAL="EVAL"
55     SIMPLE="EVAL"
56 gross 3502 POST="POST"
57 gross 3569 SMOOTH="SMOOTH"
58     def __init__(self, domain, useReduced=False, solver="EVAL", verbose=False, w=1.):
59 gross 3499 """
60     initializes the Darcy flux problem
61     :param domain: domain of the problem
62     :type domain: `Domain`
63     :param useReduced: uses reduced oreder on flux and pressure
64     :type useReduced: ``bool``
65 gross 3502 :param solver: solver method
66 gross 3569 :type solver: in [`DarcyFlow.EVAL`, `DarcyFlow.POST', `DarcyFlow.SMOOTH' ]
67 gross 3502 :param verbose: if ``True`` some information on the iteration progress are printed.
68     :type verbose: ``bool``
69     :param w: weighting factor for `DarcyFlow.POST` solver
70     :type w: ``float``
71    
72 gross 3499 """
73 gross 3569 if not solver in [DarcyFlow.EVAL, DarcyFlow.POST, DarcyFlow.SMOOTH ] :
74     raise ValueError,"unknown solver %d."%solver
75    
76 gross 3499 self.domain=domain
77 gross 3502 self.solver=solver
78 gross 3499 self.useReduced=useReduced
79 gross 3502 self.verbose=verbose
80 gross 3569 self.l=None
81     self.w=None
82    
83 gross 3502 self.__pde_p=LinearSinglePDE(domain)
84     self.__pde_p.setSymmetryOn()
85     if self.useReduced: self.__pde_p.setReducedOrderOn()
86 gross 3569
87     if self.solver == self.EVAL:
88     self.__pde_v=None
89 gross 3502 if self.verbose: print "DarcyFlow: simple solver is used."
90 gross 3569
91 gross 3502 elif self.solver == self.POST:
92     if util.inf(w)<0.:
93     raise ValueError,"Weighting factor must be non-negative."
94     if self.verbose: print "DarcyFlow: global postprocessing of flux is used."
95 gross 3569 self.__pde_v=LinearPDESystem(domain)
96     self.__pde_v.setSymmetryOn()
97     if self.useReduced: self.__pde_v.setReducedOrderOn()
98     self.w=w
99     self.l=util.vol(self.domain)**(1./self.domain.getDim()) # length scale
100    
101     elif self.solver == self.SMOOTH:
102     self.__pde_v=LinearPDESystem(domain)
103     self.__pde_v.setSymmetryOn()
104     if self.useReduced: self.__pde_v.setReducedOrderOn()
105     if self.verbose: print "DarcyFlow: flux smoothing is used."
106     self.w=0
107    
108 gross 3502 self.__f=escript.Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))
109 gross 3569 self.__g=escript.Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
110 gross 3502 self.location_of_fixed_pressure = escript.Scalar(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
111 gross 3569 self.location_of_fixed_flux = escript.Vector(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
112 gross 3502
113    
114 gross 3499 def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
115     """
116     assigns values to model parameters
117    
118     :param f: volumetic sources/sinks
119     :type f: scalar value on the domain (e.g. `escript.Data`)
120     :param g: flux sources/sinks
121     :type g: vector values on the domain (e.g. `escript.Data`)
122     :param location_of_fixed_pressure: mask for locations where pressure is fixed
123     :type location_of_fixed_pressure: scalar value on the domain (e.g. `escript.Data`)
124     :param location_of_fixed_flux: mask for locations where flux is fixed.
125     :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)
126     :param permeability: permeability tensor. If scalar ``s`` is given the tensor with ``s`` on the main diagonal is used.
127 gross 3502 :type permeability: scalar or symmetric tensor values on the domain (e.g. `escript.Data`)
128 gross 3499
129     :note: the values of parameters which are not set by calling ``setValue`` are not altered.
130     :note: at any point on the boundary of the domain the pressure
131     (``location_of_fixed_pressure`` >0) or the normal component of the
132     flux (``location_of_fixed_flux[i]>0``) if direction of the normal
133     is along the *x_i* axis.
134    
135     """
136 gross 3502 if location_of_fixed_pressure!=None:
137     self.location_of_fixed_pressure=util.wherePositive(location_of_fixed_pressure)
138     self.__pde_p.setValue(q=self.location_of_fixed_pressure)
139     if location_of_fixed_flux!=None:
140     self.location_of_fixed_flux=util.wherePositive(location_of_fixed_flux)
141 gross 3569 if not self.__pde_v == None: self.__pde_v.setValue(q=self.location_of_fixed_flux)
142 gross 3499
143     if permeability!=None:
144 gross 3502
145 gross 3569 perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
146 gross 3502
147 gross 3499 if perm.getRank()==0:
148 gross 3569
149 gross 3499 perm_inv=(1./perm)
150 gross 3569 perm_inv=perm_inv*util.kronecker(self.domain.getDim())
151     perm=perm*util.kronecker(self.domain.getDim())
152 gross 3502
153    
154 gross 3499 elif perm.getRank()==2:
155     perm_inv=util.inverse(perm)
156     else:
157     raise ValueError,"illegal rank of permeability."
158    
159     self.__permeability=perm
160     self.__permeability_inv=perm_inv
161 gross 3502
162 gross 3569 #====================
163     self.__pde_p.setValue(A=self.__permeability)
164     if self.solver == self.EVAL:
165     pass # no extra work required
166     elif self.solver == self.POST:
167 gross 3502 k=util.kronecker(self.domain.getDim())
168 gross 3569 self.omega = self.w*util.length(perm_inv)*self.l*self.domain.getSize()
169     self.__pde_v.setValue(D=self.__permeability_inv, A=self.omega*util.outer(k,k))
170     elif self.solver == self.SMOOTH:
171     self.__pde_v.setValue(D=self.__permeability_inv)
172 gross 3502
173 gross 3499 if g != None:
174 gross 3569 g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
175 gross 3499 if g.isEmpty():
176 gross 3569 g=Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
177 gross 3499 else:
178     if not g.getShape()==(self.domain.getDim(),): raise ValueError,"illegal shape of g"
179 gross 3502 self.__g=g
180 gross 3499 if f !=None:
181 gross 3502 f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
182     if f.isEmpty():
183     f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
184 gross 3499 else:
185     if f.getRank()>0: raise ValueError,"illegal rank of f."
186 gross 3502 self.__f=f
187 gross 3569
188 gross 3071 def getSolverOptionsFlux(self):
189     """
190     Returns the solver options used to solve the flux problems
191     :return: `SolverOptions`
192     """
193 gross 3569 if self.__pde_v == None:
194     return None
195     else:
196     return self.__pde_v.getSolverOptions()
197 gross 3071
198     def setSolverOptionsFlux(self, options=None):
199     """
200     Sets the solver options used to solve the flux problems
201     If ``options`` is not present, the options are reset to default
202     :param options: `SolverOptions`
203     """
204 gross 3569 if not self.__pde_v == None:
205     self.__pde_v.setSolverOptions(options)
206 gross 3071
207     def getSolverOptionsPressure(self):
208     """
209     Returns the solver options used to solve the pressure problems
210     :return: `SolverOptions`
211     """
212     return self.__pde_p.getSolverOptions()
213    
214     def setSolverOptionsPressure(self, options=None):
215     """
216     Sets the solver options used to solve the pressure problems
217     If ``options`` is not present, the options are reset to default
218    
219     :param options: `SolverOptions`
220     :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
221     """
222     return self.__pde_p.setSolverOptions(options)
223    
224 gross 3569 def solve(self, u0, p0):
225 gross 3071 """
226 gross 3502 solves the problem.
227 gross 3071
228 gross 3502 :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.
229     :type u0: vector value on the domain (e.g. `escript.Data`).
230     :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.
231     :type p0: scalar value on the domain (e.g. `escript.Data`).
232     :return: flux and pressure
233     :rtype: ``tuple`` of `escript.Data`.
234 gross 2960
235 gross 3071 """
236 gross 3569 self.__pde_p.setValue(X=self.__g ,
237     Y=self.__f,
238     y= - util.inner(self.domain.getNormal(),u0 * self.location_of_fixed_flux),
239     r=p0)
240     p=self.__pde_p.getSolution()
241     u = self.getFlux(p, u0)
242 gross 3502 return u,p
243    
244     def getFlux(self,p, u0=None):
245 gross 2100 """
246 gross 3502 returns the flux for a given pressure ``p`` where the flux is equal to ``u0``
247 jfenwick 2625 on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
248 gross 3569 Notice that ``g`` is used, see `setValue`.
249 gross 2264
250 jfenwick 2625 :param p: pressure.
251 gross 3440 :type p: scalar value on the domain (e.g. `escript.Data`).
252 gross 3502 :param u0: flux on the locations of the domain marked be ``location_of_fixed_flux``.
253     :type u0: vector values on the domain (e.g. `escript.Data`) or ``None``
254 jfenwick 2625 :return: flux
255 gross 3440 :rtype: `escript.Data`
256 gross 2264 """
257 gross 3569 u_eval=self.__g-util.tensor_mult(self.__permeability,util.grad(p))
258     if self.solver == self.EVAL:
259     u = self.__g-util.tensor_mult(self.__permeability,util.grad(p))
260     elif self.solver == self.POST or self.solver == self.SMOOTH:
261     self.__pde_v.setValue(Y=util.tensor_mult(self.__permeability_inv,self.__g)-util.grad(p))
262 gross 3502 if u0 == None:
263     self.__pde_v.setValue(r=escript.Data())
264     else:
265     self.__pde_v.setValue(r=u0)
266     u= self.__pde_v.getSolution()
267     return u
268    
269 gross 1414 class StokesProblemCartesian(HomogeneousSaddlePointProblem):
270 gross 2251 """
271 gross 2264 solves
272 gross 1414
273 gross 2208 -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
274     u_{i,i}=0
275 gross 1414
276 gross 2208 u=0 where fixed_u_mask>0
277     eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j
278 gross 1414
279 gross 2264 if surface_stress is not given 0 is assumed.
280 gross 1414
281 gross 2251 typical usage:
282 gross 1414
283 gross 2208 sp=StokesProblemCartesian(domain)
284     sp.setTolerance()
285     sp.initialize(...)
286     v,p=sp.solve(v0,p0)
287 gross 2251 """
288 gross 2719 def __init__(self,domain,**kwargs):
289 gross 2100 """
290 gross 2208 initialize the Stokes Problem
291 gross 2100
292 gross 2793 The approximation spaces used for velocity (=Solution(domain)) and pressure (=ReducedSolution(domain)) must be
293     LBB complient, for instance using quadratic and linear approximation on the same element or using linear approximation
294     with macro elements for the pressure.
295    
296     :param domain: domain of the problem.
297 jfenwick 2625 :type domain: `Domain`
298 gross 2100 """
299 gross 2719 HomogeneousSaddlePointProblem.__init__(self,**kwargs)
300 gross 1414 self.domain=domain
301 gross 3502 self.__pde_v=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
302     self.__pde_v.setSymmetryOn()
303 gross 2474
304 gross 1414 self.__pde_prec=LinearPDE(domain)
305     self.__pde_prec.setReducedOrderOn()
306     self.__pde_prec.setSymmetryOn()
307    
308 gross 2415 self.__pde_proj=LinearPDE(domain)
309     self.__pde_proj.setReducedOrderOn()
310     self.__pde_proj.setValue(D=1)
311     self.__pde_proj.setSymmetryOn()
312    
313 gross 2474 def getSolverOptionsVelocity(self):
314     """
315     returns the solver options used solve the equation for velocity.
316    
317 jfenwick 2625 :rtype: `SolverOptions`
318 gross 2474 """
319 gross 3502 return self.__pde_v.getSolverOptions()
320 gross 2474 def setSolverOptionsVelocity(self, options=None):
321     """
322     set the solver options for solving the equation for velocity.
323    
324 jfenwick 2625 :param options: new solver options
325     :type options: `SolverOptions`
326 gross 2474 """
327 gross 3502 self.__pde_v.setSolverOptions(options)
328 gross 2474 def getSolverOptionsPressure(self):
329     """
330     returns the solver options used solve the equation for pressure.
331 jfenwick 2625 :rtype: `SolverOptions`
332 gross 2474 """
333     return self.__pde_prec.getSolverOptions()
334     def setSolverOptionsPressure(self, options=None):
335     """
336     set the solver options for solving the equation for pressure.
337 jfenwick 2625 :param options: new solver options
338     :type options: `SolverOptions`
339 gross 2474 """
340     self.__pde_prec.setSolverOptions(options)
341 gross 2415
342 gross 2474 def setSolverOptionsDiv(self, options=None):
343     """
344     set the solver options for solving the equation to project the divergence of
345     the velocity onto the function space of presure.
346    
347 jfenwick 2625 :param options: new solver options
348     :type options: `SolverOptions`
349 gross 2474 """
350 artak 2689 self.__pde_proj.setSolverOptions(options)
351 gross 2474 def getSolverOptionsDiv(self):
352     """
353     returns the solver options for solving the equation to project the divergence of
354     the velocity onto the function space of presure.
355    
356 jfenwick 2625 :rtype: `SolverOptions`
357 gross 2474 """
358 artak 2689 return self.__pde_proj.getSolverOptions()
359 gross 2474
360 gross 2793 def updateStokesEquation(self, v, p):
361     """
362     updates the Stokes equation to consider dependencies from ``v`` and ``p``
363     :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values.
364     """
365     pass
366     def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None):
367     """
368     assigns new values to the model parameters.
369    
370     :param f: external force
371     :type f: `Vector` object in `FunctionSpace` `Function` or similar
372     :param fixed_u_mask: mask of locations with fixed velocity.
373     :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
374     :param eta: viscosity
375     :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
376     :param surface_stress: normal surface stress
377     :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
378     :param stress: initial stress
379     :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
380     """
381     if eta !=None:
382     k=util.kronecker(self.domain.getDim())
383     kk=util.outer(k,k)
384 jfenwick 3432 self.eta=util.interpolate(eta, escript.Function(self.domain))
385 gross 2793 self.__pde_prec.setValue(D=1/self.eta)
386 gross 3502 self.__pde_v.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))
387 gross 2793 if restoration_factor!=None:
388     n=self.domain.getNormal()
389 gross 3502 self.__pde_v.setValue(d=restoration_factor*util.outer(n,n))
390 gross 2793 if fixed_u_mask!=None:
391 gross 3502 self.__pde_v.setValue(q=fixed_u_mask)
392 gross 2793 if f!=None: self.__f=f
393     if surface_stress!=None: self.__surface_stress=surface_stress
394     if stress!=None: self.__stress=stress
395    
396 jfenwick 3432 def initialize(self,f=escript.Data(),fixed_u_mask=escript.Data(),eta=1,surface_stress=escript.Data(),stress=escript.Data(), restoration_factor=0):
397 gross 2208 """
398     assigns values to the model parameters
399 gross 2100
400 jfenwick 2625 :param f: external force
401     :type f: `Vector` object in `FunctionSpace` `Function` or similar
402     :param fixed_u_mask: mask of locations with fixed velocity.
403     :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
404     :param eta: viscosity
405     :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
406     :param surface_stress: normal surface stress
407     :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
408     :param stress: initial stress
409     :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
410 gross 2208 """
411 gross 2793 self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
412 gross 1414
413 gross 2719 def Bv(self,v,tol):
414 gross 2251 """
415     returns inner product of element p and div(v)
416 gross 1414
417 jfenwick 2625 :param v: a residual
418     :return: inner product of element p and div(v)
419     :rtype: ``float``
420 gross 2100 """
421 gross 2793 self.__pde_proj.setValue(Y=-util.div(v))
422 gross 2719 self.getSolverOptionsDiv().setTolerance(tol)
423     self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
424     out=self.__pde_proj.getSolution()
425     return out
426 gross 2208
427 gross 2445 def inner_pBv(self,p,Bv):
428     """
429     returns inner product of element p and Bv=-div(v)
430    
431 jfenwick 2625 :param p: a pressure increment
432     :param Bv: a residual
433     :return: inner product of element p and Bv=-div(v)
434     :rtype: ``float``
435 gross 2445 """
436 jfenwick 3432 return util.integrate(util.interpolate(p,escript.Function(self.domain))*util.interpolate(Bv, escript.Function(self.domain)))
437 gross 2445
438 gross 2251 def inner_p(self,p0,p1):
439 gross 2100 """
440 gross 2251 Returns inner product of p0 and p1
441 gross 1414
442 jfenwick 2625 :param p0: a pressure
443     :param p1: a pressure
444     :return: inner product of p0 and p1
445     :rtype: ``float``
446 gross 2100 """
447 jfenwick 3432 s0=util.interpolate(p0, escript.Function(self.domain))
448     s1=util.interpolate(p1, escript.Function(self.domain))
449 gross 2100 return util.integrate(s0*s1)
450 artak 1550
451 gross 2251 def norm_v(self,v):
452 gross 2100 """
453 gross 2251 returns the norm of v
454 gross 2208
455 jfenwick 2625 :param v: a velovity
456     :return: norm of v
457     :rtype: non-negative ``float``
458 gross 2100 """
459 gross 2719 return util.sqrt(util.integrate(util.length(util.grad(v))**2))
460 gross 2100
461 gross 2793
462 gross 2719 def getDV(self, p, v, tol):
463 gross 1414 """
464 gross 2251 return the value for v for a given p (overwrite)
465    
466 jfenwick 2625 :param p: a pressure
467 gross 2719 :param v: a initial guess for the value v to return.
468     :return: dv given as *Adv=(f-Av-B^*p)*
469 gross 1414 """
470 gross 2793 self.updateStokesEquation(v,p)
471 gross 3502 self.__pde_v.setValue(Y=self.__f, y=self.__surface_stress)
472 gross 2719 self.getSolverOptionsVelocity().setTolerance(tol)
473     self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
474 gross 2100 if self.__stress.isEmpty():
475 gross 3502 self.__pde_v.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
476 gross 2100 else:
477 gross 3502 self.__pde_v.setValue(X=self.__stress+p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
478     out=self.__pde_v.getSolution()
479 gross 2208 return out
480 gross 1414
481 gross 2445 def norm_Bv(self,Bv):
482 gross 2251 """
483     Returns Bv (overwrite).
484    
485 jfenwick 2625 :rtype: equal to the type of p
486     :note: boundary conditions on p should be zero!
487 gross 2251 """
488 jfenwick 3432 return util.sqrt(util.integrate(util.interpolate(Bv, escript.Function(self.domain))**2))
489 gross 2251
490 gross 2719 def solve_AinvBt(self,p, tol):
491 gross 2251 """
492 gross 2719 Solves *Av=B^*p* with accuracy `tol`
493 gross 2251
494 jfenwick 2625 :param p: a pressure increment
495     :return: the solution of *Av=B^*p*
496     :note: boundary conditions on v should be zero!
497 gross 2251 """
498 gross 3502 self.__pde_v.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))
499     out=self.__pde_v.getSolution()
500 gross 2251 return out
501    
502 gross 2719 def solve_prec(self,Bv, tol):
503 gross 2251 """
504 jfenwick 2625 applies preconditioner for for *BA^{-1}B^** to *Bv*
505     with accuracy `self.getSubProblemTolerance()`
506 gross 2251
507 jfenwick 2625 :param Bv: velocity increment
508     :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*
509     :note: boundary conditions on p are zero.
510 gross 2251 """
511 gross 2445 self.__pde_prec.setValue(Y=Bv)
512 gross 2719 self.getSolverOptionsPressure().setTolerance(tol)
513     self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
514     out=self.__pde_prec.getSolution()
515     return out

  ViewVC Help
Powered by ViewVC 1.1.26