/[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 6651 - (hide annotations)
Wed Feb 7 02:12:08 2018 UTC (18 months, 1 week ago) by jfenwick
File MIME type: text/x-python
File size: 24238 byte(s)
Make everyone sad by touching all the files

Copyright dates update

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

  ViewVC Help
Powered by ViewVC 1.1.26