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

  ViewVC Help
Powered by ViewVC 1.1.26