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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3620 - (hide annotations)
Wed Oct 5 06:06:14 2011 UTC (7 years, 10 months ago) by gross
Original Path: trunk/escript/py_src/flows.py
File MIME type: text/x-python
File size: 20546 byte(s)
bug fixed
1 gross 3071 # -*- coding: utf-8 -*-
2 ksteube 1809 ########################################################
3 gross 1414 #
4 jfenwick 2881 # Copyright (c) 2003-2010 by University of Queensland
5 ksteube 1809 # Earth Systems Science Computational Center (ESSCC)
6     # http://www.uq.edu.au/esscc
7 gross 1414 #
8 ksteube 1809 # Primary Business: Queensland, Australia
9     # Licensed under the Open Software License version 3.0
10     # http://www.opensource.org/licenses/osl-3.0.php
11 gross 1414 #
12 ksteube 1809 ########################################################
13 gross 1414
14 jfenwick 2881 __copyright__="""Copyright (c) 2003-2010 by University of Queensland
15 ksteube 1809 Earth Systems Science Computational Center (ESSCC)
16     http://www.uq.edu.au/esscc
17     Primary Business: Queensland, Australia"""
18     __license__="""Licensed under the Open Software License version 3.0
19     http://www.opensource.org/licenses/osl-3.0.php"""
20 jfenwick 2344 __url__="https://launchpad.net/escript-finley"
21 ksteube 1809
22 gross 1414 """
23     Some models for flow
24    
25 jfenwick 2625 :var __author__: name of author
26     :var __copyright__: copyrights
27     :var __license__: licence agreement
28     :var __url__: url entry point on documentation
29     :var __version__: version
30     :var __date__: date of the version
31 gross 1414 """
32    
33     __author__="Lutz Gross, l.gross@uq.edu.au"
34    
35 jfenwick 3432 import escript
36 gross 1414 import util
37 gross 2474 from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions
38 gross 2264 from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
39 gross 1414
40 gross 3502 class DarcyFlow(object):
41 gross 3499 """
42     solves the problem
43    
44     *u_i+k_{ij}*p_{,j} = g_i*
45     *u_{i,i} = f*
46    
47     where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,
48    
49 gross 3569 :cvar EVAL: direct pressure gradient evaluation for flux
50     :cvar POST: global postprocessing of flux by solving the PDE *K_{ij} u_j + (w * K * l u_{k,k})_{,i}= - p_{,j} + K_{ij} g_j*
51     where *l* is the length scale, *K* is the inverse of the permeability tensor, and *w* is a positive weighting factor.
52     :cvar SMOOTH: global smoothing by solving the PDE *K_{ij} u_j= - p_{,j} + K_{ij} g_j*
53 gross 3499 """
54 gross 3569 EVAL="EVAL"
55     SIMPLE="EVAL"
56 gross 3502 POST="POST"
57 gross 3569 SMOOTH="SMOOTH"
58     def __init__(self, domain, useReduced=False, solver="EVAL", verbose=False, w=1.):
59 gross 3499 """
60     initializes the Darcy flux problem
61     :param domain: domain of the problem
62     :type domain: `Domain`
63     :param useReduced: uses reduced oreder on flux and pressure
64     :type useReduced: ``bool``
65 gross 3502 :param solver: solver method
66 gross 3569 :type solver: in [`DarcyFlow.EVAL`, `DarcyFlow.POST', `DarcyFlow.SMOOTH' ]
67 gross 3502 :param verbose: if ``True`` some information on the iteration progress are printed.
68     :type verbose: ``bool``
69     :param w: weighting factor for `DarcyFlow.POST` solver
70     :type w: ``float``
71    
72 gross 3499 """
73 gross 3569 if not solver in [DarcyFlow.EVAL, DarcyFlow.POST, DarcyFlow.SMOOTH ] :
74     raise ValueError,"unknown solver %d."%solver
75    
76 gross 3499 self.domain=domain
77 gross 3502 self.solver=solver
78 gross 3499 self.useReduced=useReduced
79 gross 3502 self.verbose=verbose
80 gross 3569 self.l=None
81     self.w=None
82    
83 gross 3502 self.__pde_p=LinearSinglePDE(domain)
84     self.__pde_p.setSymmetryOn()
85     if self.useReduced: self.__pde_p.setReducedOrderOn()
86 gross 3569
87     if self.solver == self.EVAL:
88     self.__pde_v=None
89 gross 3502 if self.verbose: print "DarcyFlow: simple solver is used."
90 gross 3569
91 gross 3502 elif self.solver == self.POST:
92     if util.inf(w)<0.:
93     raise ValueError,"Weighting factor must be non-negative."
94     if self.verbose: print "DarcyFlow: global postprocessing of flux is used."
95 gross 3569 self.__pde_v=LinearPDESystem(domain)
96     self.__pde_v.setSymmetryOn()
97     if self.useReduced: self.__pde_v.setReducedOrderOn()
98     self.w=w
99     self.l=util.vol(self.domain)**(1./self.domain.getDim()) # length scale
100    
101     elif self.solver == self.SMOOTH:
102     self.__pde_v=LinearPDESystem(domain)
103     self.__pde_v.setSymmetryOn()
104     if self.useReduced: self.__pde_v.setReducedOrderOn()
105     if self.verbose: print "DarcyFlow: flux smoothing is used."
106     self.w=0
107    
108 gross 3502 self.__f=escript.Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))
109 gross 3569 self.__g=escript.Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
110 gross 3502 self.location_of_fixed_pressure = escript.Scalar(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
111 gross 3569 self.location_of_fixed_flux = escript.Vector(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
112 gross 3619 self.perm_scale=1.
113 gross 3502
114    
115 gross 3499 def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
116     """
117     assigns values to model parameters
118    
119     :param f: volumetic sources/sinks
120     :type f: scalar value on the domain (e.g. `escript.Data`)
121     :param g: flux sources/sinks
122     :type g: vector values on the domain (e.g. `escript.Data`)
123     :param location_of_fixed_pressure: mask for locations where pressure is fixed
124     :type location_of_fixed_pressure: scalar value on the domain (e.g. `escript.Data`)
125     :param location_of_fixed_flux: mask for locations where flux is fixed.
126     :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)
127     :param permeability: permeability tensor. If scalar ``s`` is given the tensor with ``s`` on the main diagonal is used.
128 gross 3502 :type permeability: scalar or symmetric tensor values on the domain (e.g. `escript.Data`)
129 gross 3499
130     :note: the values of parameters which are not set by calling ``setValue`` are not altered.
131     :note: at any point on the boundary of the domain the pressure
132     (``location_of_fixed_pressure`` >0) or the normal component of the
133     flux (``location_of_fixed_flux[i]>0``) if direction of the normal
134     is along the *x_i* axis.
135    
136     """
137 gross 3502 if location_of_fixed_pressure!=None:
138     self.location_of_fixed_pressure=util.wherePositive(location_of_fixed_pressure)
139     self.__pde_p.setValue(q=self.location_of_fixed_pressure)
140     if location_of_fixed_flux!=None:
141     self.location_of_fixed_flux=util.wherePositive(location_of_fixed_flux)
142 gross 3569 if not self.__pde_v == None: self.__pde_v.setValue(q=self.location_of_fixed_flux)
143 gross 3499
144     if permeability!=None:
145 gross 3502
146 gross 3569 perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
147 gross 3619 self.perm_scale=util.Lsup(util.length(perm))
148 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 3620 self.__pde_v.setValue(r=u0/self.perm_scale)
269 gross 3619 u= self.__pde_v.getSolution() * self.perm_scale
270 gross 3502 return u
271    
272 gross 1414 class StokesProblemCartesian(HomogeneousSaddlePointProblem):
273 gross 2251 """
274 gross 2264 solves
275 gross 1414
276 gross 2208 -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
277     u_{i,i}=0
278 gross 1414
279 gross 2208 u=0 where fixed_u_mask>0
280     eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j
281 gross 1414
282 gross 2264 if surface_stress is not given 0 is assumed.
283 gross 1414
284 gross 2251 typical usage:
285 gross 1414
286 gross 2208 sp=StokesProblemCartesian(domain)
287     sp.setTolerance()
288     sp.initialize(...)
289     v,p=sp.solve(v0,p0)
290 gross 3582 sp.setStokesEquation(...) # new values for some parameters
291     v1,p1=sp.solve(v,p)
292 gross 2251 """
293 gross 2719 def __init__(self,domain,**kwargs):
294 gross 2100 """
295 gross 2208 initialize the Stokes Problem
296 gross 2100
297 gross 2793 The approximation spaces used for velocity (=Solution(domain)) and pressure (=ReducedSolution(domain)) must be
298     LBB complient, for instance using quadratic and linear approximation on the same element or using linear approximation
299     with macro elements for the pressure.
300    
301     :param domain: domain of the problem.
302 jfenwick 2625 :type domain: `Domain`
303 gross 2100 """
304 gross 2719 HomogeneousSaddlePointProblem.__init__(self,**kwargs)
305 gross 1414 self.domain=domain
306 gross 3502 self.__pde_v=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
307     self.__pde_v.setSymmetryOn()
308 gross 2474
309 gross 1414 self.__pde_prec=LinearPDE(domain)
310     self.__pde_prec.setReducedOrderOn()
311     self.__pde_prec.setSymmetryOn()
312    
313 gross 2415 self.__pde_proj=LinearPDE(domain)
314     self.__pde_proj.setReducedOrderOn()
315     self.__pde_proj.setValue(D=1)
316     self.__pde_proj.setSymmetryOn()
317    
318 gross 2474 def getSolverOptionsVelocity(self):
319     """
320     returns the solver options used solve the equation for velocity.
321    
322 jfenwick 2625 :rtype: `SolverOptions`
323 gross 2474 """
324 gross 3502 return self.__pde_v.getSolverOptions()
325 gross 2474 def setSolverOptionsVelocity(self, options=None):
326     """
327     set the solver options for solving the equation for velocity.
328    
329 jfenwick 2625 :param options: new solver options
330     :type options: `SolverOptions`
331 gross 2474 """
332 gross 3502 self.__pde_v.setSolverOptions(options)
333 gross 2474 def getSolverOptionsPressure(self):
334     """
335     returns the solver options used solve the equation for pressure.
336 jfenwick 2625 :rtype: `SolverOptions`
337 gross 2474 """
338     return self.__pde_prec.getSolverOptions()
339     def setSolverOptionsPressure(self, options=None):
340     """
341     set the solver options for solving the equation for pressure.
342 jfenwick 2625 :param options: new solver options
343     :type options: `SolverOptions`
344 gross 2474 """
345     self.__pde_prec.setSolverOptions(options)
346 gross 2415
347 gross 2474 def setSolverOptionsDiv(self, options=None):
348     """
349     set the solver options for solving the equation to project the divergence of
350     the velocity onto the function space of presure.
351    
352 jfenwick 2625 :param options: new solver options
353     :type options: `SolverOptions`
354 gross 2474 """
355 artak 2689 self.__pde_proj.setSolverOptions(options)
356 gross 2474 def getSolverOptionsDiv(self):
357     """
358     returns the solver options for solving the equation to project the divergence of
359     the velocity onto the function space of presure.
360    
361 jfenwick 2625 :rtype: `SolverOptions`
362 gross 2474 """
363 artak 2689 return self.__pde_proj.getSolverOptions()
364 gross 2474
365 gross 2793 def updateStokesEquation(self, v, p):
366     """
367     updates the Stokes equation to consider dependencies from ``v`` and ``p``
368 gross 3582 :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values to model parameters.
369 gross 2793 """
370     pass
371     def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None):
372     """
373     assigns new values to the model parameters.
374    
375     :param f: external force
376     :type f: `Vector` object in `FunctionSpace` `Function` or similar
377     :param fixed_u_mask: mask of locations with fixed velocity.
378     :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
379     :param eta: viscosity
380     :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
381     :param surface_stress: normal surface stress
382     :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
383     :param stress: initial stress
384     :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
385     """
386     if eta !=None:
387     k=util.kronecker(self.domain.getDim())
388     kk=util.outer(k,k)
389 jfenwick 3432 self.eta=util.interpolate(eta, escript.Function(self.domain))
390 gross 2793 self.__pde_prec.setValue(D=1/self.eta)
391 gross 3502 self.__pde_v.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))
392 gross 2793 if restoration_factor!=None:
393     n=self.domain.getNormal()
394 gross 3502 self.__pde_v.setValue(d=restoration_factor*util.outer(n,n))
395 gross 2793 if fixed_u_mask!=None:
396 gross 3502 self.__pde_v.setValue(q=fixed_u_mask)
397 gross 2793 if f!=None: self.__f=f
398     if surface_stress!=None: self.__surface_stress=surface_stress
399     if stress!=None: self.__stress=stress
400    
401 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):
402 gross 2208 """
403     assigns values to the model parameters
404 gross 2100
405 jfenwick 2625 :param f: external force
406     :type f: `Vector` object in `FunctionSpace` `Function` or similar
407     :param fixed_u_mask: mask of locations with fixed velocity.
408     :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
409     :param eta: viscosity
410     :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
411     :param surface_stress: normal surface stress
412     :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
413     :param stress: initial stress
414     :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
415 gross 2208 """
416 gross 2793 self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
417 gross 1414
418 gross 2719 def Bv(self,v,tol):
419 gross 2251 """
420     returns inner product of element p and div(v)
421 gross 1414
422 jfenwick 2625 :param v: a residual
423     :return: inner product of element p and div(v)
424     :rtype: ``float``
425 gross 2100 """
426 gross 2793 self.__pde_proj.setValue(Y=-util.div(v))
427 gross 2719 self.getSolverOptionsDiv().setTolerance(tol)
428     self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
429     out=self.__pde_proj.getSolution()
430     return out
431 gross 2208
432 gross 2445 def inner_pBv(self,p,Bv):
433     """
434     returns inner product of element p and Bv=-div(v)
435    
436 jfenwick 2625 :param p: a pressure increment
437     :param Bv: a residual
438     :return: inner product of element p and Bv=-div(v)
439     :rtype: ``float``
440 gross 2445 """
441 jfenwick 3432 return util.integrate(util.interpolate(p,escript.Function(self.domain))*util.interpolate(Bv, escript.Function(self.domain)))
442 gross 2445
443 gross 2251 def inner_p(self,p0,p1):
444 gross 2100 """
445 gross 2251 Returns inner product of p0 and p1
446 gross 1414
447 jfenwick 2625 :param p0: a pressure
448     :param p1: a pressure
449     :return: inner product of p0 and p1
450     :rtype: ``float``
451 gross 2100 """
452 jfenwick 3432 s0=util.interpolate(p0, escript.Function(self.domain))
453     s1=util.interpolate(p1, escript.Function(self.domain))
454 gross 2100 return util.integrate(s0*s1)
455 artak 1550
456 gross 2251 def norm_v(self,v):
457 gross 2100 """
458 gross 2251 returns the norm of v
459 gross 2208
460 jfenwick 2625 :param v: a velovity
461     :return: norm of v
462     :rtype: non-negative ``float``
463 gross 2100 """
464 gross 2719 return util.sqrt(util.integrate(util.length(util.grad(v))**2))
465 gross 2100
466 gross 2793
467 gross 2719 def getDV(self, p, v, tol):
468 gross 1414 """
469 gross 3582 return the value for v for a given p
470 gross 2251
471 jfenwick 2625 :param p: a pressure
472 gross 2719 :param v: a initial guess for the value v to return.
473     :return: dv given as *Adv=(f-Av-B^*p)*
474 gross 1414 """
475 gross 2793 self.updateStokesEquation(v,p)
476 gross 3502 self.__pde_v.setValue(Y=self.__f, y=self.__surface_stress)
477 gross 2719 self.getSolverOptionsVelocity().setTolerance(tol)
478     self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
479 gross 2100 if self.__stress.isEmpty():
480 gross 3502 self.__pde_v.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
481 gross 2100 else:
482 gross 3502 self.__pde_v.setValue(X=self.__stress+p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
483     out=self.__pde_v.getSolution()
484 gross 2208 return out
485 gross 1414
486 gross 2445 def norm_Bv(self,Bv):
487 gross 2251 """
488     Returns Bv (overwrite).
489    
490 jfenwick 2625 :rtype: equal to the type of p
491     :note: boundary conditions on p should be zero!
492 gross 2251 """
493 jfenwick 3432 return util.sqrt(util.integrate(util.interpolate(Bv, escript.Function(self.domain))**2))
494 gross 2251
495 gross 2719 def solve_AinvBt(self,p, tol):
496 gross 2251 """
497 gross 2719 Solves *Av=B^*p* with accuracy `tol`
498 gross 2251
499 jfenwick 2625 :param p: a pressure increment
500     :return: the solution of *Av=B^*p*
501     :note: boundary conditions on v should be zero!
502 gross 2251 """
503 gross 3502 self.__pde_v.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))
504     out=self.__pde_v.getSolution()
505 gross 2251 return out
506    
507 gross 2719 def solve_prec(self,Bv, tol):
508 gross 2251 """
509 jfenwick 2625 applies preconditioner for for *BA^{-1}B^** to *Bv*
510     with accuracy `self.getSubProblemTolerance()`
511 gross 2251
512 jfenwick 2625 :param Bv: velocity increment
513     :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*
514     :note: boundary conditions on p are zero.
515 gross 2251 """
516 gross 2445 self.__pde_prec.setValue(Y=Bv)
517 gross 2719 self.getSolverOptionsPressure().setTolerance(tol)
518     self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
519     out=self.__pde_prec.getSolution()
520     return out

  ViewVC Help
Powered by ViewVC 1.1.26