/[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 3990 - (show annotations)
Tue Sep 25 05:03:20 2012 UTC (6 years, 10 months ago) by caltinay
Original Path: trunk/escript/py_src/flows.py
File MIME type: text/x-python
File size: 24005 byte(s)
First set of assorted epydoc fixes/additions.

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

  ViewVC Help
Powered by ViewVC 1.1.26