/[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 4326 - (hide annotations)
Wed Mar 20 02:50:46 2013 UTC (6 years, 5 months ago) by caltinay
Original Path: trunk/escript/py_src/flows.py
File MIME type: text/x-python
File size: 24005 byte(s)
Removed stray print statement.

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

  ViewVC Help
Powered by ViewVC 1.1.26