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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3909 - (show annotations)
Wed Jun 13 07:03:06 2012 UTC (7 years, 2 months ago) by gross
Original Path: trunk/escript/py_src/flows.py
File MIME type: text/x-python
File size: 23906 byte(s)
darcy solver removes hydrostatic pressure from problem now
1 # -*- coding: utf-8 -*-
2 ########################################################
3 #
4 # Copyright (c) 2003-2010 by University of Queensland
5 # Earth Systems Science Computational Center (ESSCC)
6 # http://www.uq.edu.au/esscc
7 #
8 # 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 #
12 ########################################################
13
14 __copyright__="""Copyright (c) 2003-2010 by University of Queensland
15 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 __url__="https://launchpad.net/escript-finley"
21
22 """
23 Some models for flow
24
25 :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 """
32
33 __author__="Lutz Gross, l.gross@uq.edu.au"
34
35 from . import escript
36 from . import util
37 from .linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions
38 from .pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
39
40 class DarcyFlow(object):
41 """
42 solves the problem
43
44 *u_i+k_{ij}*p_{,j} = g_i*
45 *u_{i,i} = f*
46
47 where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,
48
49 :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 """
54 EVAL="EVAL"
55 SIMPLE="EVAL"
56 POST="POST"
57 SMOOTH="SMOOTH"
58 def __init__(self, domain, useReduced=False, solver="POST", verbose=False, w=1.):
59 """
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 :param solver: solver method
66 :type solver: in [`DarcyFlow.EVAL`, `DarcyFlow.POST', `DarcyFlow.SMOOTH' ]
67 :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 """
73 if not solver in [DarcyFlow.EVAL, DarcyFlow.POST, DarcyFlow.SMOOTH ] :
74 raise ValueError("unknown solver %d."%solver)
75
76 self.domain=domain
77 self.solver=solver
78 self.useReduced=useReduced
79 self.verbose=verbose
80 self.l=None
81 self.w=None
82
83 self.__pde_p=LinearSinglePDE(domain)
84 self.__pde_p.setSymmetryOn()
85 if self.useReduced: self.__pde_p.setReducedOrderOn()
86
87 if self.solver == self.EVAL:
88 self.__pde_v=None
89 if self.verbose: print("DarcyFlow: simple solver is used.")
90
91 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 self.__pde_v=LinearPDESystem(domain)
96 self.__pde_v.setSymmetryOn()
97 if self.useReduced: self.__pde_v.setReducedOrderOn()
98 self.w=w
99 x=self.domain.getX()
100 self.l=min( [util.sup(x[i])-util.inf(x[i]) for i in xrange(self.domain.getDim()) ] )
101 #self.l=util.vol(self.domain)**(1./self.domain.getDim()) # length scale
102
103 elif self.solver == self.SMOOTH:
104 self.__pde_v=LinearPDESystem(domain)
105 self.__pde_v.setSymmetryOn()
106 if self.useReduced: self.__pde_v.setReducedOrderOn()
107 if self.verbose: print("DarcyFlow: flux smoothing is used.")
108 self.w=0
109
110 self.__f=escript.Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))
111 self.__g=escript.Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
112 self.__permeability_invXg=escript.Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
113 self.__permeability_invXg_ref=util.numpy.zeros((self.domain.getDim()),util.numpy.float64)
114 self.ref_point_id=None
115 self.ref_point=util.numpy.zeros((self.domain.getDim()),util.numpy.float64)
116 self.location_of_fixed_pressure = escript.Scalar(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
117 self.location_of_fixed_flux = escript.Vector(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
118 self.perm_scale=1.
119
120
121 def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
122 """
123 assigns values to model parameters
124
125 :param f: volumetic sources/sinks
126 :type f: scalar value on the domain (e.g. `escript.Data`)
127 :param g: flux sources/sinks
128 :type g: vector values on the domain (e.g. `escript.Data`)
129 :param location_of_fixed_pressure: mask for locations where pressure is fixed
130 :type location_of_fixed_pressure: scalar value on the domain (e.g. `escript.Data`)
131 :param location_of_fixed_flux: mask for locations where flux is fixed.
132 :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)
133 :param permeability: permeability tensor. If scalar ``s`` is given the tensor with ``s`` on the main diagonal is used.
134 :type permeability: scalar or symmetric tensor values on the domain (e.g. `escript.Data`)
135
136 :note: the values of parameters which are not set by calling ``setValue`` are not altered.
137 :note: at any point on the boundary of the domain the pressure
138 (``location_of_fixed_pressure`` >0) or the normal component of the
139 flux (``location_of_fixed_flux[i]>0``) if direction of the normal
140 is along the *x_i* axis.
141
142 """
143 if location_of_fixed_pressure!=None:
144 self.location_of_fixed_pressure=util.wherePositive(util.interpolate(location_of_fixed_pressure, self.__pde_p.getFunctionSpaceForCoefficient("q")))
145 self.ref_point_id=self.location_of_fixed_pressure.maxGlobalDataPoint()
146 if not self.location_of_fixed_pressure.getTupleForGlobalDataPoint(*self.ref_point_id)[0] > 0: raise ValueError("pressure needs to be fixed at least one point.")
147 self.ref_point=self.__pde_p.getFunctionSpaceForCoefficient("q").getX().getTupleForGlobalDataPoint(*self.ref_point_id)
148 if self.verbose: print(("DarcyFlow: reference point at %s."%(self.ref_point,)))
149 self.__pde_p.setValue(q=self.location_of_fixed_pressure)
150 if location_of_fixed_flux!=None:
151 self.location_of_fixed_flux=util.wherePositive(location_of_fixed_flux)
152 if not self.__pde_v == None: self.__pde_v.setValue(q=self.location_of_fixed_flux)
153
154 if permeability!=None:
155
156 perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
157 self.perm_scale=util.Lsup(util.length(perm))
158 if self.verbose: print(("DarcyFlow: permeability scaling factor = %e."%self.perm_scale))
159 perm=perm*(1./self.perm_scale)
160
161 if perm.getRank()==0:
162
163 perm_inv=(1./perm)
164 perm_inv=perm_inv*util.kronecker(self.domain.getDim())
165 perm=perm*util.kronecker(self.domain.getDim())
166
167
168 elif perm.getRank()==2:
169 perm_inv=util.inverse(perm)
170 else:
171 raise ValueError("illegal rank of permeability.")
172
173 self.__permeability=perm
174 self.__permeability_inv=perm_inv
175
176 #====================
177 self.__pde_p.setValue(A=self.__permeability)
178 if self.solver == self.EVAL:
179 pass # no extra work required
180 elif self.solver == self.POST:
181 k=util.kronecker(self.domain.getDim())
182 self.omega = self.w*util.length(perm_inv)*self.l*self.domain.getSize()
183 #self.__pde_v.setValue(D=self.__permeability_inv, A=self.omega*util.outer(k,k))
184 self.__pde_v.setValue(D=self.__permeability_inv, A_reduced=self.omega*util.outer(k,k))
185 elif self.solver == self.SMOOTH:
186 self.__pde_v.setValue(D=self.__permeability_inv)
187
188 if g != None:
189 g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
190 if g.isEmpty():
191 g=Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
192 else:
193 if not g.getShape()==(self.domain.getDim(),): raise ValueError("illegal shape of g")
194 self.__g=g
195 self.__permeability_invXg=util.tensor_mult(self.__permeability_inv,self.__g * (1./self.perm_scale ))
196 self.__permeability_invXg_ref=util.integrate(self.__permeability_invXg)/util.vol(self.domain)
197 if f !=None:
198 f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
199 if f.isEmpty():
200 f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
201 else:
202 if f.getRank()>0: raise ValueError("illegal rank of f.")
203 self.__f=f
204
205 def getSolverOptionsFlux(self):
206 """
207 Returns the solver options used to solve the flux problems
208 :return: `SolverOptions`
209 """
210 if self.__pde_v == None:
211 return None
212 else:
213 return self.__pde_v.getSolverOptions()
214
215 def setSolverOptionsFlux(self, options=None):
216 """
217 Sets the solver options used to solve the flux problems
218 If ``options`` is not present, the options are reset to default
219 :param options: `SolverOptions`
220 """
221 if not self.__pde_v == None:
222 self.__pde_v.setSolverOptions(options)
223
224 def getSolverOptionsPressure(self):
225 """
226 Returns the solver options used to solve the pressure problems
227 :return: `SolverOptions`
228 """
229 return self.__pde_p.getSolverOptions()
230
231 def setSolverOptionsPressure(self, options=None):
232 """
233 Sets the solver options used to solve the pressure problems
234 If ``options`` is not present, the options are reset to default
235
236 :param options: `SolverOptions`
237 :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
238 """
239 return self.__pde_p.setSolverOptions(options)
240
241 def solve(self, u0, p0):
242 """
243 solves the problem.
244
245 :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.
246 :type u0: vector value on the domain (e.g. `escript.Data`).
247 :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.
248 :type p0: scalar value on the domain (e.g. `escript.Data`).
249 :return: flux and pressure
250 :rtype: ``tuple`` of `escript.Data`.
251
252 """
253 p0=util.interpolate(p0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
254 if self.ref_point_id == None:
255 p_ref=0
256 else:
257 p_ref=p0.getTupleForGlobalDataPoint(*self.ref_point_id)[0]
258 p0_hydrostatic=p_ref+util.inner(self.__permeability_invXg_ref, self.__pde_p.getFunctionSpaceForCoefficient("q").getX() - self.ref_point)
259 g_2=self.__g - util.tensor_mult(self.__permeability, self.__permeability_invXg_ref * self.perm_scale)
260 self.__pde_p.setValue(X=g_2 * 1./self.perm_scale,
261 Y=self.__f * 1./self.perm_scale,
262 y= - util.inner(self.domain.getNormal(),u0 * self.location_of_fixed_flux * 1./self.perm_scale ),
263 r=p0 - p0_hydrostatic)
264 pp=self.__pde_p.getSolution()
265 u = self._getFlux(pp, u0)
266 return u,pp + p0_hydrostatic
267
268 def getFlux(self,p, u0=None):
269 """
270 returns the flux for a given pressure ``p`` where the flux is equal to ``u0``
271 on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
272 Notice that ``g`` is used, see `setValue`.
273
274 :param p: pressure.
275 :type p: scalar value on the domain (e.g. `escript.Data`).
276 :param u0: flux on the locations of the domain marked be ``location_of_fixed_flux``.
277 :type u0: vector values on the domain (e.g. `escript.Data`) or ``None``
278 :return: flux
279 :rtype: `escript.Data`
280 """
281 p=util.interpolate(p, self.__pde_p.getFunctionSpaceForCoefficient("q"))
282 if self.ref_point_id == None:
283 p_ref=0
284 else:
285 p_ref=p.getTupleForGlobalDataPoint(*self.ref_point_id)[0]
286 p_hydrostatic=p_ref+util.inner(self.__permeability_invXg_ref, self.__pde_p.getFunctionSpaceForCoefficient("q").getX() - self.ref_point)
287 return self._getFlux(p-p_hydrostatic, u0)
288
289 def _getFlux(self,pp, u0=None):
290 """
291 returns the flux for a given pressure ``p`` where the flux is equal to ``u0``
292 on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
293 Notice that ``g`` is used, see `setValue`.
294
295 :param p: pressure.
296 :type p: scalar value on the domain (e.g. `escript.Data`).
297 :param u0: flux on the locations of the domain marked be ``location_of_fixed_flux``.
298 :type u0: vector values on the domain (e.g. `escript.Data`) or ``None``
299 :return: flux
300 :rtype: `escript.Data`
301 """
302 if self.solver == self.EVAL:
303 u = self.__g - util.tensor_mult(self.__permeability, self.perm_scale * (util.grad(pp) + self.__permeability_invXg_ref))
304 elif self.solver == self.POST or self.solver == self.SMOOTH:
305 self.__pde_v.setValue(Y= self.__permeability_invXg - (util.grad(pp) + self.__permeability_invXg_ref))
306 if u0 == None:
307 self.__pde_v.setValue(r=escript.Data())
308 else:
309 if not isinstance(u0, escript.Data) : u0 = escript.Vector(u0, escript.Solution(self.domain))
310 self.__pde_v.setValue(r=1./self.perm_scale * u0)
311 u= self.__pde_v.getSolution() * self.perm_scale
312 return u
313
314 class StokesProblemCartesian(HomogeneousSaddlePointProblem):
315 """
316 solves
317
318 -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
319 u_{i,i}=0
320
321 u=0 where fixed_u_mask>0
322 eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j
323
324 if surface_stress is not given 0 is assumed.
325
326 typical usage:
327
328 sp=StokesProblemCartesian(domain)
329 sp.setTolerance()
330 sp.initialize(...)
331 v,p=sp.solve(v0,p0)
332 sp.setStokesEquation(...) # new values for some parameters
333 v1,p1=sp.solve(v,p)
334 """
335 def __init__(self,domain,**kwargs):
336 """
337 initialize the Stokes Problem
338
339 The approximation spaces used for velocity (=Solution(domain)) and pressure (=ReducedSolution(domain)) must be
340 LBB complient, for instance using quadratic and linear approximation on the same element or using linear approximation
341 with macro elements for the pressure.
342
343 :param domain: domain of the problem.
344 :type domain: `Domain`
345 """
346 HomogeneousSaddlePointProblem.__init__(self,**kwargs)
347 self.domain=domain
348 self.__pde_v=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
349 self.__pde_v.setSymmetryOn()
350
351 self.__pde_prec=LinearPDE(domain)
352 self.__pde_prec.setReducedOrderOn()
353 self.__pde_prec.setSymmetryOn()
354
355 self.__pde_proj=LinearPDE(domain)
356 self.__pde_proj.setReducedOrderOn()
357 self.__pde_proj.setValue(D=1)
358 self.__pde_proj.setSymmetryOn()
359
360 def getSolverOptionsVelocity(self):
361 """
362 returns the solver options used solve the equation for velocity.
363
364 :rtype: `SolverOptions`
365 """
366 return self.__pde_v.getSolverOptions()
367 def setSolverOptionsVelocity(self, options=None):
368 """
369 set the solver options for solving the equation for velocity.
370
371 :param options: new solver options
372 :type options: `SolverOptions`
373 """
374 self.__pde_v.setSolverOptions(options)
375 def getSolverOptionsPressure(self):
376 """
377 returns the solver options used solve the equation for pressure.
378 :rtype: `SolverOptions`
379 """
380 return self.__pde_prec.getSolverOptions()
381 def setSolverOptionsPressure(self, options=None):
382 """
383 set the solver options for solving the equation for pressure.
384 :param options: new solver options
385 :type options: `SolverOptions`
386 """
387 self.__pde_prec.setSolverOptions(options)
388
389 def setSolverOptionsDiv(self, options=None):
390 """
391 set the solver options for solving the equation to project the divergence of
392 the velocity onto the function space of presure.
393
394 :param options: new solver options
395 :type options: `SolverOptions`
396 """
397 self.__pde_proj.setSolverOptions(options)
398 def getSolverOptionsDiv(self):
399 """
400 returns the solver options for solving the equation to project the divergence of
401 the velocity onto the function space of presure.
402
403 :rtype: `SolverOptions`
404 """
405 return self.__pde_proj.getSolverOptions()
406
407 def updateStokesEquation(self, v, p):
408 """
409 updates the Stokes equation to consider dependencies from ``v`` and ``p``
410 :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values to model parameters.
411 """
412 pass
413 def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None):
414 """
415 assigns new values to the model parameters.
416
417 :param f: external force
418 :type f: `Vector` object in `FunctionSpace` `Function` or similar
419 :param fixed_u_mask: mask of locations with fixed velocity.
420 :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
421 :param eta: viscosity
422 :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
423 :param surface_stress: normal surface stress
424 :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
425 :param stress: initial stress
426 :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
427 """
428 if eta !=None:
429 k=util.kronecker(self.domain.getDim())
430 kk=util.outer(k,k)
431 self.eta=util.interpolate(eta, escript.Function(self.domain))
432 self.__pde_prec.setValue(D=1/self.eta)
433 self.__pde_v.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))
434 if restoration_factor!=None:
435 n=self.domain.getNormal()
436 self.__pde_v.setValue(d=restoration_factor*util.outer(n,n))
437 if fixed_u_mask!=None:
438 self.__pde_v.setValue(q=fixed_u_mask)
439 if f!=None: self.__f=f
440 if surface_stress!=None: self.__surface_stress=surface_stress
441 if stress!=None: self.__stress=stress
442
443 def initialize(self,f=escript.Data(),fixed_u_mask=escript.Data(),eta=1,surface_stress=escript.Data(),stress=escript.Data(), restoration_factor=0):
444 """
445 assigns values to the model parameters
446
447 :param f: external force
448 :type f: `Vector` object in `FunctionSpace` `Function` or similar
449 :param fixed_u_mask: mask of locations with fixed velocity.
450 :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
451 :param eta: viscosity
452 :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
453 :param surface_stress: normal surface stress
454 :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
455 :param stress: initial stress
456 :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
457 """
458 self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
459
460 def Bv(self,v,tol):
461 """
462 returns inner product of element p and div(v)
463
464 :param v: a residual
465 :return: inner product of element p and div(v)
466 :rtype: ``float``
467 """
468 self.__pde_proj.setValue(Y=-util.div(v))
469 self.getSolverOptionsDiv().setTolerance(tol)
470 self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
471 out=self.__pde_proj.getSolution()
472 return out
473
474 def inner_pBv(self,p,Bv):
475 """
476 returns inner product of element p and Bv=-div(v)
477
478 :param p: a pressure increment
479 :param Bv: a residual
480 :return: inner product of element p and Bv=-div(v)
481 :rtype: ``float``
482 """
483 return util.integrate(util.interpolate(p,escript.Function(self.domain))*util.interpolate(Bv, escript.Function(self.domain)))
484
485 def inner_p(self,p0,p1):
486 """
487 Returns inner product of p0 and p1
488
489 :param p0: a pressure
490 :param p1: a pressure
491 :return: inner product of p0 and p1
492 :rtype: ``float``
493 """
494 s0=util.interpolate(p0, escript.Function(self.domain))
495 s1=util.interpolate(p1, escript.Function(self.domain))
496 return util.integrate(s0*s1)
497
498 def norm_v(self,v):
499 """
500 returns the norm of v
501
502 :param v: a velovity
503 :return: norm of v
504 :rtype: non-negative ``float``
505 """
506 return util.sqrt(util.integrate(util.length(util.grad(v))**2))
507
508
509 def getDV(self, p, v, tol):
510 """
511 return the value for v for a given p
512
513 :param p: a pressure
514 :param v: a initial guess for the value v to return.
515 :return: dv given as *Adv=(f-Av-B^*p)*
516 """
517 self.updateStokesEquation(v,p)
518 self.__pde_v.setValue(Y=self.__f, y=self.__surface_stress)
519 self.getSolverOptionsVelocity().setTolerance(tol)
520 self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
521 if self.__stress.isEmpty():
522 self.__pde_v.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
523 else:
524 self.__pde_v.setValue(X=self.__stress+p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
525 out=self.__pde_v.getSolution()
526 return out
527
528 def norm_Bv(self,Bv):
529 """
530 Returns Bv (overwrite).
531
532 :rtype: equal to the type of p
533 :note: boundary conditions on p should be zero!
534 """
535 return util.sqrt(util.integrate(util.interpolate(Bv, escript.Function(self.domain))**2))
536
537 def solve_AinvBt(self,p, tol):
538 """
539 Solves *Av=B^*p* with accuracy `tol`
540
541 :param p: a pressure increment
542 :return: the solution of *Av=B^*p*
543 :note: boundary conditions on v should be zero!
544 """
545 self.__pde_v.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))
546 out=self.__pde_v.getSolution()
547 return out
548
549 def solve_prec(self,Bv, tol):
550 """
551 applies preconditioner for for *BA^{-1}B^** to *Bv*
552 with accuracy `self.getSubProblemTolerance()`
553
554 :param Bv: velocity increment
555 :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*
556 :note: boundary conditions on p are zero.
557 """
558 self.__pde_prec.setValue(Y=Bv)
559 self.getSolverOptionsPressure().setTolerance(tol)
560 self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
561 out=self.__pde_prec.getSolution()
562 return out

  ViewVC Help
Powered by ViewVC 1.1.26