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

  ViewVC Help
Powered by ViewVC 1.1.26