/[escript]/trunk/escript/py_src/flows.py
ViewVC logotype

Contents of /trunk/escript/py_src/flows.py

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26