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

  ViewVC Help
Powered by ViewVC 1.1.26