/[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 3981 - (hide annotations)
Fri Sep 21 02:47:54 2012 UTC (6 years, 11 months ago) by jfenwick
Original Path: trunk/escript/py_src/flows.py
File MIME type: text/x-python
File size: 24000 byte(s)
First pass of updating copyright notices
1 gross 3071 # -*- coding: utf-8 -*-
2 jfenwick 3981 ##############################################################################
3 gross 1414 #
4 jfenwick 3911 # Copyright (c) 2003-2012 by University of Queensland
5 jfenwick 3981 # http://www.uq.edu.au
6 gross 1414 #
7 ksteube 1809 # Primary Business: Queensland, Australia
8     # Licensed under the Open Software License version 3.0
9     # http://www.opensource.org/licenses/osl-3.0.php
10 gross 1414 #
11 jfenwick 3981 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12     # Development since 2012 by School of Earth Sciences
13     #
14     ##############################################################################
15 gross 1414
16 jfenwick 3911 __copyright__="""Copyright (c) 2003-2012 by University of Queensland
17 jfenwick 3981 http://www.uq.edu.au
18 ksteube 1809 Primary Business: Queensland, Australia"""
19     __license__="""Licensed under the Open Software License version 3.0
20     http://www.opensource.org/licenses/osl-3.0.php"""
21 jfenwick 2344 __url__="https://launchpad.net/escript-finley"
22 ksteube 1809
23 gross 1414 """
24     Some models for flow
25    
26 jfenwick 2625 :var __author__: name of author
27     :var __copyright__: copyrights
28     :var __license__: licence agreement
29     :var __url__: url entry point on documentation
30     :var __version__: version
31     :var __date__: date of the version
32 gross 1414 """
33    
34     __author__="Lutz Gross, l.gross@uq.edu.au"
35    
36 jfenwick 3771 from . import escript
37     from . import util
38     from .linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions
39     from .pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
40 gross 1414
41 gross 3502 class DarcyFlow(object):
42 gross 3499 """
43     solves the problem
44    
45     *u_i+k_{ij}*p_{,j} = g_i*
46     *u_{i,i} = f*
47    
48     where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,
49    
50 gross 3569 :cvar EVAL: direct pressure gradient evaluation for flux
51     :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*
52     where *l* is the length scale, *K* is the inverse of the permeability tensor, and *w* is a positive weighting factor.
53     :cvar SMOOTH: global smoothing by solving the PDE *K_{ij} u_j= - p_{,j} + K_{ij} g_j*
54 gross 3499 """
55 gross 3569 EVAL="EVAL"
56     SIMPLE="EVAL"
57 gross 3502 POST="POST"
58 gross 3569 SMOOTH="SMOOTH"
59 gross 3630 def __init__(self, domain, useReduced=False, solver="POST", verbose=False, w=1.):
60 gross 3499 """
61     initializes the Darcy flux problem
62     :param domain: domain of the problem
63     :type domain: `Domain`
64     :param useReduced: uses reduced oreder on flux and pressure
65     :type useReduced: ``bool``
66 gross 3502 :param solver: solver method
67 gross 3569 :type solver: in [`DarcyFlow.EVAL`, `DarcyFlow.POST', `DarcyFlow.SMOOTH' ]
68 gross 3502 :param verbose: if ``True`` some information on the iteration progress are printed.
69     :type verbose: ``bool``
70     :param w: weighting factor for `DarcyFlow.POST` solver
71     :type w: ``float``
72    
73 gross 3499 """
74 gross 3569 if not solver in [DarcyFlow.EVAL, DarcyFlow.POST, DarcyFlow.SMOOTH ] :
75 jfenwick 3771 raise ValueError("unknown solver %d."%solver)
76 gross 3569
77 gross 3499 self.domain=domain
78 gross 3502 self.solver=solver
79 gross 3499 self.useReduced=useReduced
80 gross 3502 self.verbose=verbose
81 gross 3569 self.l=None
82     self.w=None
83    
84 gross 3502 self.__pde_p=LinearSinglePDE(domain)
85     self.__pde_p.setSymmetryOn()
86     if self.useReduced: self.__pde_p.setReducedOrderOn()
87 gross 3569
88     if self.solver == self.EVAL:
89     self.__pde_v=None
90 jfenwick 3852 if self.verbose: print("DarcyFlow: simple solver is used.")
91 gross 3569
92 gross 3502 elif self.solver == self.POST:
93 jfenwick 3852 if util.inf(w)<0.:
94     raise ValueError("Weighting factor must be non-negative.")
95     if self.verbose: print("DarcyFlow: global postprocessing of flux is used.")
96 gross 3569 self.__pde_v=LinearPDESystem(domain)
97     self.__pde_v.setSymmetryOn()
98     if self.useReduced: self.__pde_v.setReducedOrderOn()
99 jfenwick 3852 self.w=w
100 gross 3905 x=self.domain.getX()
101     self.l=min( [util.sup(x[i])-util.inf(x[i]) for i in xrange(self.domain.getDim()) ] )
102     #self.l=util.vol(self.domain)**(1./self.domain.getDim()) # length scale
103 gross 3569
104     elif self.solver == self.SMOOTH:
105     self.__pde_v=LinearPDESystem(domain)
106     self.__pde_v.setSymmetryOn()
107     if self.useReduced: self.__pde_v.setReducedOrderOn()
108 jfenwick 3852 if self.verbose: print("DarcyFlow: flux smoothing is used.")
109     self.w=0
110 gross 3569
111 gross 3502 self.__f=escript.Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))
112 gross 3569 self.__g=escript.Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
113 gross 3909 self.__permeability_invXg=escript.Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
114     self.__permeability_invXg_ref=util.numpy.zeros((self.domain.getDim()),util.numpy.float64)
115     self.ref_point_id=None
116     self.ref_point=util.numpy.zeros((self.domain.getDim()),util.numpy.float64)
117 gross 3502 self.location_of_fixed_pressure = escript.Scalar(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
118 gross 3569 self.location_of_fixed_flux = escript.Vector(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
119 gross 3619 self.perm_scale=1.
120 gross 3502
121    
122 gross 3499 def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
123     """
124     assigns values to model parameters
125    
126     :param f: volumetic sources/sinks
127     :type f: scalar value on the domain (e.g. `escript.Data`)
128     :param g: flux sources/sinks
129     :type g: vector values on the domain (e.g. `escript.Data`)
130     :param location_of_fixed_pressure: mask for locations where pressure is fixed
131     :type location_of_fixed_pressure: scalar value on the domain (e.g. `escript.Data`)
132     :param location_of_fixed_flux: mask for locations where flux is fixed.
133     :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)
134     :param permeability: permeability tensor. If scalar ``s`` is given the tensor with ``s`` on the main diagonal is used.
135 gross 3502 :type permeability: scalar or symmetric tensor values on the domain (e.g. `escript.Data`)
136 gross 3499
137     :note: the values of parameters which are not set by calling ``setValue`` are not altered.
138     :note: at any point on the boundary of the domain the pressure
139     (``location_of_fixed_pressure`` >0) or the normal component of the
140     flux (``location_of_fixed_flux[i]>0``) if direction of the normal
141     is along the *x_i* axis.
142    
143     """
144 gross 3502 if location_of_fixed_pressure!=None:
145 gross 3909 self.location_of_fixed_pressure=util.wherePositive(util.interpolate(location_of_fixed_pressure, self.__pde_p.getFunctionSpaceForCoefficient("q")))
146     self.ref_point_id=self.location_of_fixed_pressure.maxGlobalDataPoint()
147     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.")
148     self.ref_point=self.__pde_p.getFunctionSpaceForCoefficient("q").getX().getTupleForGlobalDataPoint(*self.ref_point_id)
149     if self.verbose: print(("DarcyFlow: reference point at %s."%(self.ref_point,)))
150 gross 3502 self.__pde_p.setValue(q=self.location_of_fixed_pressure)
151     if location_of_fixed_flux!=None:
152     self.location_of_fixed_flux=util.wherePositive(location_of_fixed_flux)
153 gross 3910 if not self.__pde_v == None:
154     self.__pde_v.setValue(q=self.location_of_fixed_flux)
155 gross 3499
156     if permeability!=None:
157 gross 3502
158 jfenwick 3852 perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
159 gross 3619 self.perm_scale=util.Lsup(util.length(perm))
160 jfenwick 3852 if self.verbose: print(("DarcyFlow: permeability scaling factor = %e."%self.perm_scale))
161 gross 3619 perm=perm*(1./self.perm_scale)
162 gross 3502
163 jfenwick 3852 if perm.getRank()==0:
164 gross 3569
165 jfenwick 3852 perm_inv=(1./perm)
166     perm_inv=perm_inv*util.kronecker(self.domain.getDim())
167     perm=perm*util.kronecker(self.domain.getDim())
168 gross 3502
169    
170 jfenwick 3852 elif perm.getRank()==2:
171     perm_inv=util.inverse(perm)
172     else:
173     raise ValueError("illegal rank of permeability.")
174 gross 3499
175 jfenwick 3852 self.__permeability=perm
176     self.__permeability_inv=perm_inv
177 gross 3502
178 gross 3569 #====================
179 jfenwick 3852 self.__pde_p.setValue(A=self.__permeability)
180 gross 3569 if self.solver == self.EVAL:
181     pass # no extra work required
182     elif self.solver == self.POST:
183 jfenwick 3852 k=util.kronecker(self.domain.getDim())
184     self.omega = self.w*util.length(perm_inv)*self.l*self.domain.getSize()
185 gross 3885 #self.__pde_v.setValue(D=self.__permeability_inv, A=self.omega*util.outer(k,k))
186     self.__pde_v.setValue(D=self.__permeability_inv, A_reduced=self.omega*util.outer(k,k))
187 gross 3569 elif self.solver == self.SMOOTH:
188 jfenwick 3852 self.__pde_v.setValue(D=self.__permeability_inv)
189 gross 3502
190 gross 3499 if g != None:
191 jfenwick 3852 g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
192     if g.isEmpty():
193     g=Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
194     else:
195     if not g.getShape()==(self.domain.getDim(),): raise ValueError("illegal shape of g")
196     self.__g=g
197 gross 3910 self.__permeability_invXg=util.tensor_mult(self.__permeability_inv,self.__g * (1./self.perm_scale ))
198 gross 3909 self.__permeability_invXg_ref=util.integrate(self.__permeability_invXg)/util.vol(self.domain)
199 gross 3499 if f !=None:
200 jfenwick 3852 f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
201     if f.isEmpty():
202     f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
203     else:
204     if f.getRank()>0: raise ValueError("illegal rank of f.")
205     self.__f=f
206 gross 3569
207 gross 3071 def getSolverOptionsFlux(self):
208     """
209     Returns the solver options used to solve the flux problems
210     :return: `SolverOptions`
211     """
212 gross 3569 if self.__pde_v == None:
213     return None
214     else:
215     return self.__pde_v.getSolverOptions()
216 gross 3071
217     def setSolverOptionsFlux(self, options=None):
218     """
219     Sets the solver options used to solve the flux problems
220     If ``options`` is not present, the options are reset to default
221     :param options: `SolverOptions`
222     """
223 gross 3569 if not self.__pde_v == None:
224     self.__pde_v.setSolverOptions(options)
225 gross 3071
226     def getSolverOptionsPressure(self):
227     """
228     Returns the solver options used to solve the pressure problems
229     :return: `SolverOptions`
230     """
231     return self.__pde_p.getSolverOptions()
232    
233     def setSolverOptionsPressure(self, options=None):
234     """
235     Sets the solver options used to solve the pressure problems
236     If ``options`` is not present, the options are reset to default
237    
238     :param options: `SolverOptions`
239     :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
240     """
241     return self.__pde_p.setSolverOptions(options)
242    
243 gross 3569 def solve(self, u0, p0):
244 gross 3071 """
245 gross 3502 solves the problem.
246 gross 3071
247 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.
248     :type u0: vector value on the domain (e.g. `escript.Data`).
249     :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.
250     :type p0: scalar value on the domain (e.g. `escript.Data`).
251     :return: flux and pressure
252     :rtype: ``tuple`` of `escript.Data`.
253 gross 2960
254 gross 3071 """
255 gross 3909 p0=util.interpolate(p0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
256     if self.ref_point_id == None:
257     p_ref=0
258     else:
259     p_ref=p0.getTupleForGlobalDataPoint(*self.ref_point_id)[0]
260     p0_hydrostatic=p_ref+util.inner(self.__permeability_invXg_ref, self.__pde_p.getFunctionSpaceForCoefficient("q").getX() - self.ref_point)
261     g_2=self.__g - util.tensor_mult(self.__permeability, self.__permeability_invXg_ref * self.perm_scale)
262     self.__pde_p.setValue(X=g_2 * 1./self.perm_scale,
263 gross 3620 Y=self.__f * 1./self.perm_scale,
264     y= - util.inner(self.domain.getNormal(),u0 * self.location_of_fixed_flux * 1./self.perm_scale ),
265 gross 3909 r=p0 - p0_hydrostatic)
266     pp=self.__pde_p.getSolution()
267     u = self._getFlux(pp, u0)
268     return u,pp + p0_hydrostatic
269 gross 3502
270     def getFlux(self,p, u0=None):
271 gross 2100 """
272 gross 3502 returns the flux for a given pressure ``p`` where the flux is equal to ``u0``
273 jfenwick 2625 on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
274 gross 3569 Notice that ``g`` is used, see `setValue`.
275 gross 2264
276 jfenwick 2625 :param p: pressure.
277 gross 3440 :type p: scalar value on the domain (e.g. `escript.Data`).
278 gross 3502 :param u0: flux on the locations of the domain marked be ``location_of_fixed_flux``.
279     :type u0: vector values on the domain (e.g. `escript.Data`) or ``None``
280 jfenwick 2625 :return: flux
281 gross 3440 :rtype: `escript.Data`
282 gross 2264 """
283 gross 3909 p=util.interpolate(p, self.__pde_p.getFunctionSpaceForCoefficient("q"))
284     if self.ref_point_id == None:
285     p_ref=0
286     else:
287     p_ref=p.getTupleForGlobalDataPoint(*self.ref_point_id)[0]
288     p_hydrostatic=p_ref+util.inner(self.__permeability_invXg_ref, self.__pde_p.getFunctionSpaceForCoefficient("q").getX() - self.ref_point)
289     return self._getFlux(p-p_hydrostatic, u0)
290    
291     def _getFlux(self,pp, u0=None):
292     """
293     returns the flux for a given pressure ``p`` where the flux is equal to ``u0``
294     on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
295     Notice that ``g`` is used, see `setValue`.
296    
297     :param p: pressure.
298     :type p: scalar value on the domain (e.g. `escript.Data`).
299     :param u0: flux on the locations of the domain marked be ``location_of_fixed_flux``.
300     :type u0: vector values on the domain (e.g. `escript.Data`) or ``None``
301     :return: flux
302     :rtype: `escript.Data`
303     """
304 gross 3569 if self.solver == self.EVAL:
305 gross 3909 u = self.__g - util.tensor_mult(self.__permeability, self.perm_scale * (util.grad(pp) + self.__permeability_invXg_ref))
306 gross 3569 elif self.solver == self.POST or self.solver == self.SMOOTH:
307 gross 3909 self.__pde_v.setValue(Y= self.__permeability_invXg - (util.grad(pp) + self.__permeability_invXg_ref))
308 gross 3910 print
309 gross 3502 if u0 == None:
310 jfenwick 3852 self.__pde_v.setValue(r=escript.Data())
311     else:
312 gross 3630 if not isinstance(u0, escript.Data) : u0 = escript.Vector(u0, escript.Solution(self.domain))
313 jfenwick 3852 self.__pde_v.setValue(r=1./self.perm_scale * u0)
314 gross 3910 u= self.__pde_v.getSolution() * self.perm_scale
315 jfenwick 3852 return u
316 gross 3502
317 gross 1414 class StokesProblemCartesian(HomogeneousSaddlePointProblem):
318 gross 2251 """
319 gross 2264 solves
320 gross 1414
321 gross 2208 -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
322     u_{i,i}=0
323 gross 1414
324 gross 2208 u=0 where fixed_u_mask>0
325     eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j
326 gross 1414
327 gross 2264 if surface_stress is not given 0 is assumed.
328 gross 1414
329 gross 2251 typical usage:
330 gross 1414
331 gross 2208 sp=StokesProblemCartesian(domain)
332     sp.setTolerance()
333     sp.initialize(...)
334     v,p=sp.solve(v0,p0)
335 gross 3582 sp.setStokesEquation(...) # new values for some parameters
336     v1,p1=sp.solve(v,p)
337 gross 2251 """
338 gross 2719 def __init__(self,domain,**kwargs):
339 gross 2100 """
340 gross 2208 initialize the Stokes Problem
341 gross 2100
342 gross 2793 The approximation spaces used for velocity (=Solution(domain)) and pressure (=ReducedSolution(domain)) must be
343     LBB complient, for instance using quadratic and linear approximation on the same element or using linear approximation
344     with macro elements for the pressure.
345    
346     :param domain: domain of the problem.
347 jfenwick 2625 :type domain: `Domain`
348 gross 2100 """
349 gross 2719 HomogeneousSaddlePointProblem.__init__(self,**kwargs)
350 gross 1414 self.domain=domain
351 gross 3502 self.__pde_v=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
352     self.__pde_v.setSymmetryOn()
353 gross 2474
354 gross 1414 self.__pde_prec=LinearPDE(domain)
355     self.__pde_prec.setReducedOrderOn()
356     self.__pde_prec.setSymmetryOn()
357    
358 gross 2415 self.__pde_proj=LinearPDE(domain)
359     self.__pde_proj.setReducedOrderOn()
360 jfenwick 3852 self.__pde_proj.setValue(D=1)
361 gross 2415 self.__pde_proj.setSymmetryOn()
362    
363 gross 2474 def getSolverOptionsVelocity(self):
364     """
365     returns the solver options used solve the equation for velocity.
366    
367 jfenwick 2625 :rtype: `SolverOptions`
368 gross 2474 """
369 jfenwick 3852 return self.__pde_v.getSolverOptions()
370 gross 2474 def setSolverOptionsVelocity(self, options=None):
371     """
372     set the solver options for solving the equation for velocity.
373    
374 jfenwick 2625 :param options: new solver options
375     :type options: `SolverOptions`
376 gross 2474 """
377 gross 3502 self.__pde_v.setSolverOptions(options)
378 gross 2474 def getSolverOptionsPressure(self):
379     """
380     returns the solver options used solve the equation for pressure.
381 jfenwick 2625 :rtype: `SolverOptions`
382 gross 2474 """
383 jfenwick 3852 return self.__pde_prec.getSolverOptions()
384 gross 2474 def setSolverOptionsPressure(self, options=None):
385     """
386     set the solver options for solving the equation for pressure.
387 jfenwick 2625 :param options: new solver options
388     :type options: `SolverOptions`
389 gross 2474 """
390 jfenwick 3852 self.__pde_prec.setSolverOptions(options)
391 gross 2415
392 gross 2474 def setSolverOptionsDiv(self, options=None):
393     """
394     set the solver options for solving the equation to project the divergence of
395     the velocity onto the function space of presure.
396    
397 jfenwick 2625 :param options: new solver options
398     :type options: `SolverOptions`
399 gross 2474 """
400 jfenwick 3852 self.__pde_proj.setSolverOptions(options)
401 gross 2474 def getSolverOptionsDiv(self):
402     """
403     returns the solver options for solving the equation to project the divergence of
404     the velocity onto the function space of presure.
405    
406 jfenwick 2625 :rtype: `SolverOptions`
407 gross 2474 """
408 jfenwick 3852 return self.__pde_proj.getSolverOptions()
409 gross 2474
410 gross 2793 def updateStokesEquation(self, v, p):
411     """
412     updates the Stokes equation to consider dependencies from ``v`` and ``p``
413 gross 3582 :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values to model parameters.
414 gross 2793 """
415     pass
416     def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None):
417     """
418     assigns new values to the model parameters.
419    
420     :param f: external force
421     :type f: `Vector` object in `FunctionSpace` `Function` or similar
422     :param fixed_u_mask: mask of locations with fixed velocity.
423     :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
424     :param eta: viscosity
425     :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
426     :param surface_stress: normal surface stress
427     :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
428     :param stress: initial stress
429     :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
430     """
431     if eta !=None:
432     k=util.kronecker(self.domain.getDim())
433     kk=util.outer(k,k)
434 jfenwick 3432 self.eta=util.interpolate(eta, escript.Function(self.domain))
435 jfenwick 3852 self.__pde_prec.setValue(D=1/self.eta)
436 gross 3502 self.__pde_v.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))
437 gross 2793 if restoration_factor!=None:
438     n=self.domain.getNormal()
439 gross 3502 self.__pde_v.setValue(d=restoration_factor*util.outer(n,n))
440 gross 2793 if fixed_u_mask!=None:
441 gross 3502 self.__pde_v.setValue(q=fixed_u_mask)
442 gross 2793 if f!=None: self.__f=f
443     if surface_stress!=None: self.__surface_stress=surface_stress
444     if stress!=None: self.__stress=stress
445    
446 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):
447 gross 2208 """
448     assigns values to the model parameters
449 gross 2100
450 jfenwick 2625 :param f: external force
451     :type f: `Vector` object in `FunctionSpace` `Function` or similar
452     :param fixed_u_mask: mask of locations with fixed velocity.
453     :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
454     :param eta: viscosity
455     :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
456     :param surface_stress: normal surface stress
457     :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
458     :param stress: initial stress
459 jfenwick 3852 :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
460 gross 2208 """
461 gross 2793 self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
462 gross 1414
463 gross 2719 def Bv(self,v,tol):
464 gross 2251 """
465     returns inner product of element p and div(v)
466 gross 1414
467 jfenwick 2625 :param v: a residual
468     :return: inner product of element p and div(v)
469     :rtype: ``float``
470 gross 2100 """
471 gross 2793 self.__pde_proj.setValue(Y=-util.div(v))
472 jfenwick 3852 self.getSolverOptionsDiv().setTolerance(tol)
473     self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
474 gross 2719 out=self.__pde_proj.getSolution()
475     return out
476 gross 2208
477 gross 2445 def inner_pBv(self,p,Bv):
478     """
479     returns inner product of element p and Bv=-div(v)
480    
481 jfenwick 2625 :param p: a pressure increment
482     :param Bv: a residual
483     :return: inner product of element p and Bv=-div(v)
484     :rtype: ``float``
485 gross 2445 """
486 jfenwick 3432 return util.integrate(util.interpolate(p,escript.Function(self.domain))*util.interpolate(Bv, escript.Function(self.domain)))
487 gross 2445
488 gross 2251 def inner_p(self,p0,p1):
489 gross 2100 """
490 gross 2251 Returns inner product of p0 and p1
491 gross 1414
492 jfenwick 2625 :param p0: a pressure
493     :param p1: a pressure
494     :return: inner product of p0 and p1
495     :rtype: ``float``
496 gross 2100 """
497 jfenwick 3432 s0=util.interpolate(p0, escript.Function(self.domain))
498     s1=util.interpolate(p1, escript.Function(self.domain))
499 gross 2100 return util.integrate(s0*s1)
500 artak 1550
501 gross 2251 def norm_v(self,v):
502 gross 2100 """
503 gross 2251 returns the norm of v
504 gross 2208
505 jfenwick 2625 :param v: a velovity
506     :return: norm of v
507     :rtype: non-negative ``float``
508 gross 2100 """
509 gross 2719 return util.sqrt(util.integrate(util.length(util.grad(v))**2))
510 gross 2100
511 gross 2793
512 gross 2719 def getDV(self, p, v, tol):
513 gross 1414 """
514 gross 3582 return the value for v for a given p
515 gross 2251
516 jfenwick 2625 :param p: a pressure
517 gross 2719 :param v: a initial guess for the value v to return.
518     :return: dv given as *Adv=(f-Av-B^*p)*
519 gross 1414 """
520 gross 2793 self.updateStokesEquation(v,p)
521 gross 3502 self.__pde_v.setValue(Y=self.__f, y=self.__surface_stress)
522 jfenwick 3852 self.getSolverOptionsVelocity().setTolerance(tol)
523     self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
524 gross 2100 if self.__stress.isEmpty():
525 gross 3502 self.__pde_v.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
526 gross 2100 else:
527 gross 3502 self.__pde_v.setValue(X=self.__stress+p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
528     out=self.__pde_v.getSolution()
529 gross 2208 return out
530 gross 1414
531 gross 2445 def norm_Bv(self,Bv):
532 gross 2251 """
533     Returns Bv (overwrite).
534    
535 jfenwick 2625 :rtype: equal to the type of p
536     :note: boundary conditions on p should be zero!
537 gross 2251 """
538 jfenwick 3432 return util.sqrt(util.integrate(util.interpolate(Bv, escript.Function(self.domain))**2))
539 gross 2251
540 gross 2719 def solve_AinvBt(self,p, tol):
541 gross 2251 """
542 gross 2719 Solves *Av=B^*p* with accuracy `tol`
543 gross 2251
544 jfenwick 2625 :param p: a pressure increment
545     :return: the solution of *Av=B^*p*
546     :note: boundary conditions on v should be zero!
547 gross 2251 """
548 gross 3502 self.__pde_v.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))
549     out=self.__pde_v.getSolution()
550 gross 2251 return out
551    
552 gross 2719 def solve_prec(self,Bv, tol):
553 gross 2251 """
554 jfenwick 2625 applies preconditioner for for *BA^{-1}B^** to *Bv*
555     with accuracy `self.getSubProblemTolerance()`
556 gross 2251
557 jfenwick 2625 :param Bv: velocity increment
558     :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*
559     :note: boundary conditions on p are zero.
560 gross 2251 """
561 gross 2445 self.__pde_prec.setValue(Y=Bv)
562 jfenwick 3852 self.getSolverOptionsPressure().setTolerance(tol)
563     self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
564 gross 2719 out=self.__pde_prec.getSolution()
565     return out

  ViewVC Help
Powered by ViewVC 1.1.26