/[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 3502 - (show annotations)
Thu Apr 28 05:06:24 2011 UTC (8 years, 3 months ago) by gross
Original Path: trunk/escript/py_src/flows.py
File MIME type: text/x-python
File size: 26407 byte(s)
Darcy flow solver is now supporting 4 solvers
1 # -*- coding: utf-8 -*-
2 ########################################################
3 #
4 # Copyright (c) 2003-2010 by University of Queensland
5 # Earth Systems Science Computational Center (ESSCC)
6 # http://www.uq.edu.au/esscc
7 #
8 # 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 #
12 ########################################################
13
14 __copyright__="""Copyright (c) 2003-2010 by University of Queensland
15 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 __url__="https://launchpad.net/escript-finley"
21
22 """
23 Some models for flow
24
25 :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 """
32
33 __author__="Lutz Gross, l.gross@uq.edu.au"
34
35 import escript
36 import util
37 from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions
38 from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
39
40 class DarcyFlow(object):
41 """
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 :cvar SIMPLE: simple solver
50 :cvar POST: solver using global postprocessing of flux
51 :cvar STAB: solver uses (non-symmetric) stabilization
52 :cvar SYMSTAB: solver uses symmetric stabilization
53 """
54 SIMPLE="SIMPLE"
55 POST="POST"
56 STAB="STAB"
57 SYMSTAB="SYMSTAB"
58 def __init__(self, domain, useReduced=False, solver="SYMSTAB", verbose=False, w=1.):
59 """
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 :param solver: solver method
66 :type solver: in [`DarcyFlow.SIMPLE`, `DarcyFlow.POST', `DarcyFlow.STAB`, `DarcyFlow.SYMSTAB` ]
67 :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 """
73 self.domain=domain
74 self.solver=solver
75 self.useReduced=useReduced
76 self.verbose=verbose
77 self.scale=1.
78
79
80 self.__pde_v=LinearPDESystem(domain)
81 self.__pde_v.setSymmetryOn()
82 if self.useReduced: self.__pde_v.setReducedOrderOn()
83 self.__pde_p=LinearSinglePDE(domain)
84 self.__pde_p.setSymmetryOn()
85 if self.useReduced: self.__pde_p.setReducedOrderOn()
86
87 if self.solver == self.SIMPLE:
88 if self.verbose: print "DarcyFlow: simple solver is used."
89 self.__pde_v.setValue(D=util.kronecker(self.domain.getDim()))
90 elif self.solver == self.POST:
91 self.w=w
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 elif self.solver == self.STAB:
96 if self.verbose: print "DarcyFlow: (non-symmetric) stabilization is used."
97 elif self.solver == self.SYMSTAB:
98 if self.verbose: print "DarcyFlow: symmetric stabilization is used."
99 else:
100 raise ValueError,"unknown solver %s"%self.solver
101 self.__f=escript.Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))
102 self.__g=escript.Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
103 self.location_of_fixed_pressure = escript.Scalar(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
104 self.location_of_fixed_flux = escript.Vector(0, self.__pde_v.getFunctionSpaceForCoefficient("q"))
105 self.setTolerance()
106
107
108 def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
109 """
110 assigns values to model parameters
111
112 :param f: volumetic sources/sinks
113 :type f: scalar value on the domain (e.g. `escript.Data`)
114 :param g: flux sources/sinks
115 :type g: vector values on the domain (e.g. `escript.Data`)
116 :param location_of_fixed_pressure: mask for locations where pressure is fixed
117 :type location_of_fixed_pressure: scalar value on the domain (e.g. `escript.Data`)
118 :param location_of_fixed_flux: mask for locations where flux is fixed.
119 :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)
120 :param permeability: permeability tensor. If scalar ``s`` is given the tensor with ``s`` on the main diagonal is used.
121 :type permeability: scalar or symmetric tensor values on the domain (e.g. `escript.Data`)
122
123 :note: the values of parameters which are not set by calling ``setValue`` are not altered.
124 :note: at any point on the boundary of the domain the pressure
125 (``location_of_fixed_pressure`` >0) or the normal component of the
126 flux (``location_of_fixed_flux[i]>0``) if direction of the normal
127 is along the *x_i* axis.
128
129 """
130 if location_of_fixed_pressure!=None:
131 self.location_of_fixed_pressure=util.wherePositive(location_of_fixed_pressure)
132 self.__pde_p.setValue(q=self.location_of_fixed_pressure)
133 if location_of_fixed_flux!=None:
134 self.location_of_fixed_flux=util.wherePositive(location_of_fixed_flux)
135 self.__pde_v.setValue(q=self.location_of_fixed_flux)
136
137
138 # pressure is rescaled by the factor 1/self.scale
139 if permeability!=None:
140
141 perm=util.interpolate(permeability,self.__pde_v.getFunctionSpaceForCoefficient("A"))
142 V=util.vol(self.domain)
143 l=V**(1./self.domain.getDim())
144
145 if perm.getRank()==0:
146 perm_inv=(1./perm)
147 self.scale=util.integrate(perm_inv)/V*l
148 perm_inv=perm_inv*((1./self.scale)*util.kronecker(self.domain.getDim()))
149 perm=perm*(self.scale*util.kronecker(self.domain.getDim()))
150
151
152 elif perm.getRank()==2:
153 perm_inv=util.inverse(perm)
154 self.scale=util.sqrt(util.integrate(util.length(perm_inv)**2)/V)*l
155 perm_inv*=(1./self.scale)
156 perm=perm*self.scale
157 else:
158 raise ValueError,"illegal rank of permeability."
159
160 self.__permeability=perm
161 self.__permeability_inv=perm_inv
162 if self.verbose: print "DarcyFlow: scaling factor for pressure is %e."%self.scale
163
164 if self.solver == self.SIMPLE:
165 self.__pde_p.setValue(A=self.__permeability)
166 elif self.solver == self.POST:
167 self.__pde_p.setValue(A=self.__permeability)
168 k=util.kronecker(self.domain.getDim())
169 self.lamb = self.w*util.length(perm_inv)*l
170 self.__pde_v.setValue(D=self.__permeability_inv, A=self.lamb*self.domain.getSize()*util.outer(k,k))
171 elif self.solver == self.STAB:
172 self.__pde_p.setValue(A=0.5*self.__permeability)
173 self.__pde_v.setValue(D=0.5*self.__permeability_inv)
174 elif self.solver == self.SYMSTAB:
175 self.__pde_p.setValue(A=0.5*self.__permeability)
176 self.__pde_v.setValue(D=0.5*self.__permeability_inv)
177
178 if g != None:
179 g=util.interpolate(g, self.__pde_v.getFunctionSpaceForCoefficient("Y"))
180 if g.isEmpty():
181 g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
182 else:
183 if not g.getShape()==(self.domain.getDim(),): raise ValueError,"illegal shape of g"
184 self.__g=g
185 if f !=None:
186 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 def getSolverOptionsFlux(self):
193 """
194 Returns the solver options used to solve the flux problems
195 :return: `SolverOptions`
196 """
197 return self.__pde_v.getSolverOptions()
198
199 def setSolverOptionsFlux(self, options=None):
200 """
201 Sets the solver options used to solve the flux problems
202 If ``options`` is not present, the options are reset to default
203 :param options: `SolverOptions`
204 """
205 return self.__pde_v.setSolverOptions(options)
206
207 def getSolverOptionsPressure(self):
208 """
209 Returns the solver options used to solve the pressure problems
210 :return: `SolverOptions`
211 """
212 return self.__pde_p.getSolverOptions()
213
214 def setSolverOptionsPressure(self, options=None):
215 """
216 Sets the solver options used to solve the pressure problems
217 If ``options`` is not present, the options are reset to default
218
219 :param options: `SolverOptions`
220 :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
221 """
222 return self.__pde_p.setSolverOptions(options)
223
224 def setTolerance(self,rtol=1e-4):
225 """
226 sets the relative tolerance ``rtol`` for the pressure for the stabelized solvers.
227
228 :param rtol: relative tolerance for the pressure
229 :type rtol: non-negative ``float``
230 """
231 if rtol<0:
232 raise ValueError,"Relative tolerance needs to be non-negative."
233 self.__rtol=rtol
234
235 def getTolerance(self):
236 """
237 returns the relative tolerance
238 :return: current relative tolerance
239 :rtype: ``float``
240 """
241 return self.__rtol
242
243 def solve(self,u0,p0, max_iter=100, iter_restart=20):
244 """
245 solves the problem.
246
247 The iteration is terminated if the residual norm is less then self.getTolerance().
248
249 :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.
250 :type u0: vector value on the domain (e.g. `escript.Data`).
251 :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.
252 :type p0: scalar value on the domain (e.g. `escript.Data`).
253 :param max_iter: maximum number of (outer) iteration steps for the stabilization solvers,
254 :type max_iter: ``int``
255 :param iter_restart: number of steps after which the iteration is restarted. The larger ``iter_restart`` the larger the required memory.
256 A small value for ``iter_restart`` may require a large number of iteration steps or may even lead to a failure
257 of the iteration. ``iter_restart`` is relevant for the stabilization solvers only.
258 :type iter_restart: ``int``
259 :return: flux and pressure
260 :rtype: ``tuple`` of `escript.Data`.
261
262 """
263 # rescale initial guess:
264 p0=p0/self.scale
265 if self.solver == self.SIMPLE or self.solver == self.POST :
266 self.__pde_p.setValue(X=self.__g ,
267 Y=self.__f,
268 y=-util.inner(self.domain.getNormal(),u0 * self.location_of_fixed_flux),
269 r=p0)
270 p=self.__pde_p.getSolution()
271 u = self.getFlux(p, u0)
272 elif self.solver == self.STAB:
273 u,p = self.__solve_STAB(u0,p0, max_iter, iter_restart)
274 elif self.solver == self.SYMSTAB:
275 u,p = self.__solve_SYMSTAB(u0,p0, max_iter, iter_restart)
276
277 if self.verbose:
278 KGp=util.tensor_mult(self.__permeability,util.grad(p))
279 def_p=self.__g-(u+KGp)
280 def_v=self.__f-util.div(u, self.__pde_v.getFunctionSpaceForCoefficient("X"))
281 print "DarcyFlux: |g-u-K*grad(p)|_2 = %e (|u|_2 = %e)."%(self.__L2(def_p),self.__L2(u))
282 print "DarcyFlux: |f-div(u)|_2 = %e (|grad(u)|_2 = %e)."%(self.__L2(def_v),self.__L2(util.grad(u)))
283 #rescale result
284 p=p*self.scale
285 return u,p
286
287 def getFlux(self,p, u0=None):
288 """
289 returns the flux for a given pressure ``p`` where the flux is equal to ``u0``
290 on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
291 Notice that ``g`` and ``f`` are used, see `setValue`.
292
293 :param p: pressure.
294 :type p: scalar value on the domain (e.g. `escript.Data`).
295 :param u0: flux on the locations of the domain marked be ``location_of_fixed_flux``.
296 :type u0: vector values on the domain (e.g. `escript.Data`) or ``None``
297 :return: flux
298 :rtype: `escript.Data`
299 """
300 if self.solver == self.SIMPLE or self.solver == self.POST :
301 KGp=util.tensor_mult(self.__permeability,util.grad(p))
302 self.__pde_v.setValue(Y=self.__g-KGp, X=escript.Data())
303 if u0 == None:
304 self.__pde_v.setValue(r=escript.Data())
305 else:
306 self.__pde_v.setValue(r=u0)
307 u= self.__pde_v.getSolution()
308 elif self.solver == self.POST:
309 self.__pde_v.setValue(Y=util.tensor_mult(self.__permeability_inv,self.__g)-util.grad(p),
310 X=self.lamb * self.__f * util.kronecker(self.domain.getDim()))
311 if u0 == None:
312 self.__pde_v.setValue(r=escript.Data())
313 else:
314 self.__pde_v.setValue(r=u0)
315 u= self.__pde_v.getSolution()
316 elif self.solver == self.STAB:
317 gp=util.grad(p)
318 self.__pde_v.setValue(Y=0.5*(util.tensor_mult(self.__permeability_inv,self.__g)+gp),
319 X= p * util.kronecker(self.domain.getDim()),
320 y= - p * self.domain.getNormal())
321 if u0 == None:
322 self.__pde_v.setValue(r=escript.Data())
323 else:
324 self.__pde_v.setValue(r=u0)
325 u= self.__pde_v.getSolution()
326 elif self.solver == self.SYMSTAB:
327 gp=util.grad(p)
328 self.__pde_v.setValue(Y=0.5*(util.tensor_mult(self.__permeability_inv,self.__g)-gp),
329 X= escript.Data() ,
330 y= escript.Data() )
331 if u0 == None:
332 self.__pde_v.setValue(r=escript.Data())
333 else:
334 self.__pde_v.setValue(r=u0)
335 u= self.__pde_v.getSolution()
336 return u
337
338
339 def __solve_STAB(self, u0, p0, max_iter, iter_restart):
340 # p0 is used as an initial guess
341 u=self.getFlux(p0, u0)
342 self.__pde_p.setValue( Y=self.__f-util.div(u),
343 X=0.5*(self.__g - u - util.tensor_mult(self.__permeability,util.grad(p0)) ),
344 y= escript.Data(),
345 r=escript.Data())
346 dp=self.__pde_p.getSolution()
347 p=GMRES(dp,
348 self.__STAB_Aprod,
349 p0,
350 self.__inner,
351 atol=self.__norm(p0+dp)*self.getTolerance() ,
352 rtol=0.,
353 iter_max=max_iter,
354 iter_restart=iter_restart,
355 verbose=self.verbose,P_R=None)
356
357 u=self.getFlux(p, u0)
358 return u,p
359
360 def __solve_SYMSTAB(self, u0, p0, max_iter, iter_restart):
361 # p0 is used as an initial guess
362 u=self.getFlux(p0, u0)
363 self.__pde_p.setValue( Y= self.__f,
364 X= 0.5*(self.__g + u - util.tensor_mult(self.__permeability,util.grad(p0)) ),
365 y= - util.inner(self.domain.getNormal(), u),
366 r=escript.Data())
367 dp=self.__pde_p.getSolution()
368 p=GMRES(dp,
369 self.__SYMSTAB_Aprod,
370 p0,
371 self.__inner,
372 atol=self.__norm(p0+dp)*self.getTolerance() ,
373 rtol=0.,
374 iter_max=max_iter,
375 iter_restart=iter_restart,
376 verbose=self.verbose,P_R=None)
377
378 u=self.getFlux(p, u0)
379 return u,p
380
381 def __L2(self,v):
382 return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))
383
384 def __norm(self,r):
385 return util.sqrt(self.__inner(r,r))
386
387 def __inner(self,r,s):
388 return util.integrate(util.inner(r,s), escript.Function(self.domain))
389
390 def __STAB_Aprod(self,p):
391 gp=util.grad(p)
392 self.__pde_v.setValue(Y=-0.5*gp,
393 X=-p*util.kronecker(self.__pde_v.getDomain()),
394 y= p * self.domain.getNormal(),
395 r=escript.Data())
396 u = -self.__pde_v.getSolution()
397 self.__pde_p.setValue(Y=util.div(u),
398 X=0.5*(u+util.tensor_mult(self.__permeability,gp)),
399 y=escript.Data(),
400 r=escript.Data())
401
402 return self.__pde_p.getSolution()
403
404 def __SYMSTAB_Aprod(self,p):
405 gp=util.grad(p)
406 self.__pde_v.setValue(Y=0.5*gp ,
407 X=escript.Data(),
408 y=escript.Data(),
409 r=escript.Data())
410 u = -self.__pde_v.getSolution()
411 self.__pde_p.setValue(Y=escript.Data(),
412 X=0.5*(-u+util.tensor_mult(self.__permeability,gp)),
413 y= util.inner(self.domain.getNormal(), u),
414 r=escript.Data())
415
416 return self.__pde_p.getSolution()
417
418
419 class StokesProblemCartesian(HomogeneousSaddlePointProblem):
420 """
421 solves
422
423 -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
424 u_{i,i}=0
425
426 u=0 where fixed_u_mask>0
427 eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j
428
429 if surface_stress is not given 0 is assumed.
430
431 typical usage:
432
433 sp=StokesProblemCartesian(domain)
434 sp.setTolerance()
435 sp.initialize(...)
436 v,p=sp.solve(v0,p0)
437 """
438 def __init__(self,domain,**kwargs):
439 """
440 initialize the Stokes Problem
441
442 The approximation spaces used for velocity (=Solution(domain)) and pressure (=ReducedSolution(domain)) must be
443 LBB complient, for instance using quadratic and linear approximation on the same element or using linear approximation
444 with macro elements for the pressure.
445
446 :param domain: domain of the problem.
447 :type domain: `Domain`
448 """
449 HomogeneousSaddlePointProblem.__init__(self,**kwargs)
450 self.domain=domain
451 self.__pde_v=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
452 self.__pde_v.setSymmetryOn()
453
454 self.__pde_prec=LinearPDE(domain)
455 self.__pde_prec.setReducedOrderOn()
456 self.__pde_prec.setSymmetryOn()
457
458 self.__pde_proj=LinearPDE(domain)
459 self.__pde_proj.setReducedOrderOn()
460 self.__pde_proj.setValue(D=1)
461 self.__pde_proj.setSymmetryOn()
462
463 def getSolverOptionsVelocity(self):
464 """
465 returns the solver options used solve the equation for velocity.
466
467 :rtype: `SolverOptions`
468 """
469 return self.__pde_v.getSolverOptions()
470 def setSolverOptionsVelocity(self, options=None):
471 """
472 set the solver options for solving the equation for velocity.
473
474 :param options: new solver options
475 :type options: `SolverOptions`
476 """
477 self.__pde_v.setSolverOptions(options)
478 def getSolverOptionsPressure(self):
479 """
480 returns the solver options used solve the equation for pressure.
481 :rtype: `SolverOptions`
482 """
483 return self.__pde_prec.getSolverOptions()
484 def setSolverOptionsPressure(self, options=None):
485 """
486 set the solver options for solving the equation for pressure.
487 :param options: new solver options
488 :type options: `SolverOptions`
489 """
490 self.__pde_prec.setSolverOptions(options)
491
492 def setSolverOptionsDiv(self, options=None):
493 """
494 set the solver options for solving the equation to project the divergence of
495 the velocity onto the function space of presure.
496
497 :param options: new solver options
498 :type options: `SolverOptions`
499 """
500 self.__pde_proj.setSolverOptions(options)
501 def getSolverOptionsDiv(self):
502 """
503 returns the solver options for solving the equation to project the divergence of
504 the velocity onto the function space of presure.
505
506 :rtype: `SolverOptions`
507 """
508 return self.__pde_proj.getSolverOptions()
509
510 def updateStokesEquation(self, v, p):
511 """
512 updates the Stokes equation to consider dependencies from ``v`` and ``p``
513 :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values.
514 """
515 pass
516 def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None):
517 """
518 assigns new values to the model parameters.
519
520 :param f: external force
521 :type f: `Vector` object in `FunctionSpace` `Function` or similar
522 :param fixed_u_mask: mask of locations with fixed velocity.
523 :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
524 :param eta: viscosity
525 :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
526 :param surface_stress: normal surface stress
527 :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
528 :param stress: initial stress
529 :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
530 """
531 if eta !=None:
532 k=util.kronecker(self.domain.getDim())
533 kk=util.outer(k,k)
534 self.eta=util.interpolate(eta, escript.Function(self.domain))
535 self.__pde_prec.setValue(D=1/self.eta)
536 self.__pde_v.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))
537 if restoration_factor!=None:
538 n=self.domain.getNormal()
539 self.__pde_v.setValue(d=restoration_factor*util.outer(n,n))
540 if fixed_u_mask!=None:
541 self.__pde_v.setValue(q=fixed_u_mask)
542 if f!=None: self.__f=f
543 if surface_stress!=None: self.__surface_stress=surface_stress
544 if stress!=None: self.__stress=stress
545
546 def initialize(self,f=escript.Data(),fixed_u_mask=escript.Data(),eta=1,surface_stress=escript.Data(),stress=escript.Data(), restoration_factor=0):
547 """
548 assigns values to the model parameters
549
550 :param f: external force
551 :type f: `Vector` object in `FunctionSpace` `Function` or similar
552 :param fixed_u_mask: mask of locations with fixed velocity.
553 :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
554 :param eta: viscosity
555 :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
556 :param surface_stress: normal surface stress
557 :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
558 :param stress: initial stress
559 :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
560 """
561 self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
562
563 def Bv(self,v,tol):
564 """
565 returns inner product of element p and div(v)
566
567 :param v: a residual
568 :return: inner product of element p and div(v)
569 :rtype: ``float``
570 """
571 self.__pde_proj.setValue(Y=-util.div(v))
572 self.getSolverOptionsDiv().setTolerance(tol)
573 self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
574 out=self.__pde_proj.getSolution()
575 return out
576
577 def inner_pBv(self,p,Bv):
578 """
579 returns inner product of element p and Bv=-div(v)
580
581 :param p: a pressure increment
582 :param Bv: a residual
583 :return: inner product of element p and Bv=-div(v)
584 :rtype: ``float``
585 """
586 return util.integrate(util.interpolate(p,escript.Function(self.domain))*util.interpolate(Bv, escript.Function(self.domain)))
587
588 def inner_p(self,p0,p1):
589 """
590 Returns inner product of p0 and p1
591
592 :param p0: a pressure
593 :param p1: a pressure
594 :return: inner product of p0 and p1
595 :rtype: ``float``
596 """
597 s0=util.interpolate(p0, escript.Function(self.domain))
598 s1=util.interpolate(p1, escript.Function(self.domain))
599 return util.integrate(s0*s1)
600
601 def norm_v(self,v):
602 """
603 returns the norm of v
604
605 :param v: a velovity
606 :return: norm of v
607 :rtype: non-negative ``float``
608 """
609 return util.sqrt(util.integrate(util.length(util.grad(v))**2))
610
611
612 def getDV(self, p, v, tol):
613 """
614 return the value for v for a given p (overwrite)
615
616 :param p: a pressure
617 :param v: a initial guess for the value v to return.
618 :return: dv given as *Adv=(f-Av-B^*p)*
619 """
620 self.updateStokesEquation(v,p)
621 self.__pde_v.setValue(Y=self.__f, y=self.__surface_stress)
622 self.getSolverOptionsVelocity().setTolerance(tol)
623 self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
624 if self.__stress.isEmpty():
625 self.__pde_v.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
626 else:
627 self.__pde_v.setValue(X=self.__stress+p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
628 out=self.__pde_v.getSolution()
629 return out
630
631 def norm_Bv(self,Bv):
632 """
633 Returns Bv (overwrite).
634
635 :rtype: equal to the type of p
636 :note: boundary conditions on p should be zero!
637 """
638 return util.sqrt(util.integrate(util.interpolate(Bv, escript.Function(self.domain))**2))
639
640 def solve_AinvBt(self,p, tol):
641 """
642 Solves *Av=B^*p* with accuracy `tol`
643
644 :param p: a pressure increment
645 :return: the solution of *Av=B^*p*
646 :note: boundary conditions on v should be zero!
647 """
648 self.__pde_v.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))
649 out=self.__pde_v.getSolution()
650 return out
651
652 def solve_prec(self,Bv, tol):
653 """
654 applies preconditioner for for *BA^{-1}B^** to *Bv*
655 with accuracy `self.getSubProblemTolerance()`
656
657 :param Bv: velocity increment
658 :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*
659 :note: boundary conditions on p are zero.
660 """
661 self.__pde_prec.setValue(Y=Bv)
662 self.getSolverOptionsPressure().setTolerance(tol)
663 self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
664 out=self.__pde_prec.getSolution()
665 return out

  ViewVC Help
Powered by ViewVC 1.1.26