/[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 5024 - (show annotations)
Tue Jun 10 08:46:07 2014 UTC (5 years, 2 months ago) by jfenwick
File MIME type: text/x-python
File size: 24192 byte(s)
Fixing more == None

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

  ViewVC Help
Powered by ViewVC 1.1.26