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

  ViewVC Help
Powered by ViewVC 1.1.26