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

  ViewVC Help
Powered by ViewVC 1.1.26