/[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 3462 - (show annotations)
Mon Feb 7 02:01:50 2011 UTC (8 years, 7 months ago) by caltinay
File MIME type: text/x-python
File size: 44083 byte(s)
Removed debug output.

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 :note: The problem is solved in a least squares formulation.
50 """
51
52 def __init__(self, domain, useReduced=False, adaptSubTolerance=True, solveForFlux=False, useVPIteration=True, weighting_scale=0.1):
53 """
54 initializes the Darcy flux problem
55 :param domain: domain of the problem
56 :type domain: `Domain`
57 :param useReduced: uses reduced oreder on flux and pressure
58 :type useReduced: ``bool``
59 :param adaptSubTolerance: switches on automatic subtolerance selection
60 :type adaptSubTolerance: ``bool``
61 :param solveForFlux: if True the solver solves for the flux (do not use!)
62 :type solveForFlux: ``bool``
63 :param useVPIteration: if True altenative iteration over v and p is performed. Otherwise V and P are calculated in a single PDE.
64 :type useVPIteration: ``bool``
65 """
66 self.domain=domain
67 self.useVPIteration=useVPIteration
68 self.useReduced=useReduced
69 self.weighting_scale=weighting_scale
70 if self.useVPIteration:
71 self.solveForFlux=solveForFlux
72 self.__adaptSubTolerance=adaptSubTolerance
73 self.verbose=False
74
75 self.__pde_k=LinearPDESystem(domain)
76 self.__pde_k.setSymmetryOn()
77 if self.useReduced: self.__pde_k.setReducedOrderOn()
78
79 self.__pde_p=LinearSinglePDE(domain)
80 self.__pde_p.setSymmetryOn()
81 if self.useReduced: self.__pde_p.setReducedOrderOn()
82 self.setTolerance()
83 self.setAbsoluteTolerance()
84 else:
85 self.__pde_k=LinearPDE(self.domain, numEquations=self.domain.getDim()+1)
86 self.__pde_k.setSymmetryOn()
87 if self.useReduced: self.__pde_k.setReducedOrderOn()
88 C=self.__pde_k.createCoefficient("C")
89 B=self.__pde_k.createCoefficient("B")
90 for i in range(self.domain.getDim()):
91 C[i,self.domain.getDim(),i]=1
92 B[self.domain.getDim(),i,i]=1
93 self.__pde_k.setValue(C=C, B=B)
94 self.__f=escript.Scalar(0,self.__pde_k.getFunctionSpaceForCoefficient("X"))
95 self.__g=escript.Vector(0,self.__pde_k.getFunctionSpaceForCoefficient("Y"))
96
97 def getSolverOptionsFlux(self):
98 """
99 Returns the solver options used to solve the flux problems
100
101 *K^{-1} u=F*
102
103 :return: `SolverOptions`
104 """
105 return self.__pde_k.getSolverOptions()
106
107 def setSolverOptionsFlux(self, options=None):
108 """
109 Sets the solver options used to solve the flux problems
110
111 *K^{-1}u=F*
112
113 If ``options`` is not present, the options are reset to default
114
115 :param options: `SolverOptions`
116 :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
117 """
118 return self.__pde_v.setSolverOptions(options)
119
120 def getSolverOptionsPressure(self):
121 """
122 Returns the solver options used to solve the pressure problems
123
124 *(Q^* K Q)p=-Q^*G*
125
126 :return: `SolverOptions`
127 """
128 return self.__pde_p.getSolverOptions()
129
130 def setSolverOptionsPressure(self, options=None):
131 """
132 Sets the solver options used to solve the pressure problems
133
134 *(Q^* K Q)p=-Q^*G*
135
136 If ``options`` is not present, the options are reset to default
137
138 :param options: `SolverOptions`
139 :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
140 """
141 return self.__pde_p.setSolverOptions(options)
142
143
144 def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
145 """
146 assigns values to model parameters
147
148 :param f: volumetic sources/sinks
149 :type f: scalar value on the domain (e.g. `escript.Data`)
150 :param g: flux sources/sinks
151 :type g: vector values on the domain (e.g. `escript.Data`)
152 :param location_of_fixed_pressure: mask for locations where pressure is fixed
153 :type location_of_fixed_pressure: scalar value on the domain (e.g. `escript.Data`)
154 :param location_of_fixed_flux: mask for locations where flux is fixed.
155 :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)
156 :param permeability: permeability tensor. If scalar ``s`` is given the tensor with ``s`` on the main diagonal is used.
157 :type permeability: scalar or tensor values on the domain (e.g. `escript.Data`)
158
159 :note: the values of parameters which are not set by calling ``setValue`` are not altered.
160 :note: at any point on the boundary of the domain the pressure
161 (``location_of_fixed_pressure`` >0) or the normal component of the
162 flux (``location_of_fixed_flux[i]>0``) if direction of the normal
163 is along the *x_i* axis.
164
165 """
166 if self.useVPIteration:
167 if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)
168 if location_of_fixed_flux!=None: self.__pde_k.setValue(q=location_of_fixed_flux)
169 else:
170 q=self.__pde_k.getCoefficient("q")
171 if q.isEmpty(): q=self.__pde_k.createCoefficient("q")
172 if location_of_fixed_pressure!=None: q[self.domain.getDim()]=location_of_fixed_pressure
173 if location_of_fixed_flux!=None: q[:self.domain.getDim()]=location_of_fixed_flux
174 self.__pde_k.setValue(q=q)
175
176 # flux is rescaled by the factor mean value(perm_inv)*length where length**self.domain.getDim()=vol(self.domain)
177 if permeability!=None:
178 perm=util.interpolate(permeability,self.__pde_k.getFunctionSpaceForCoefficient("A"))
179 V=util.vol(self.domain)
180 if perm.getRank()==0:
181 perm_inv=(1./perm)
182 if self.useVPIteration:
183 self.scale=1.
184 else:
185 self.scale=util.integrate(perm_inv)*V**(1./self.domain.getDim()-1.)
186
187 perm_inv=perm_inv*((1./self.scale)*util.kronecker(self.domain.getDim()))
188 perm=perm*(self.scale*util.kronecker(self.domain.getDim()))
189 elif perm.getRank()==2:
190 perm_inv=util.inverse(perm)
191 if self.useVPIteration:
192 self.scale=1.
193 else:
194 self.scale=util.sqrt(util.integrate(util.length(perm_inv)**2)*V**(2./self.domain.getDim()-1.)/self.domain.getDim())
195 perm_inv*=(1./self.scale)
196 perm=perm*self.scale
197 else:
198 raise ValueError,"illegal rank of permeability."
199
200 self.__permeability=perm
201 self.__permeability_inv=perm_inv
202
203 self.__l2 =(util.longestEdge(self.domain)**2*util.length(self.__permeability_inv))*self.weighting_scale
204 if self.useVPIteration:
205 if self.solveForFlux:
206 self.__pde_k.setValue(D=self.__permeability_inv)
207 else:
208 self.__pde_k.setValue(D=self.__permeability_inv, A=self.__l2*util.outer(util.kronecker(self.domain),util.kronecker(self.domain)))
209 self.__pde_p.setValue(A=self.__permeability)
210 else:
211 D=self.__pde_k.createCoefficient("D")
212 A=self.__pde_k.createCoefficient("A")
213 D[:self.domain.getDim(),:self.domain.getDim()]=self.__permeability_inv
214 for i in range(self.domain.getDim()):
215 for j in range(self.domain.getDim()):
216 A[i,i,j,j]=self.__l2
217 A[self.domain.getDim(),:,self.domain.getDim(),:]=self.__permeability
218 self.__pde_k.setValue(A=A, D=D)
219 if g !=None:
220 g=util.interpolate(g, self.__pde_k.getFunctionSpaceForCoefficient("Y"))
221 if g.isEmpty():
222 g=Vector(0,self.__pde_k.getFunctionSpaceForCoefficient("Y"))
223 else:
224 if not g.getShape()==(self.domain.getDim(),):
225 raise ValueError,"illegal shape of g"
226 self.__g=g
227 elif permeability!=None:
228 X
229 if f !=None:
230 f=util.interpolate(f, self.__pde_k.getFunctionSpaceForCoefficient("X"))
231 if f.isEmpty():
232 f=Scalar(0,self.__pde_k.getFunctionSpaceForCoefficient("X"))
233 else:
234 if f.getRank()>0: raise ValueError,"illegal rank of f."
235 self.__f=f
236
237
238 def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):
239 """
240 solves the problem.
241
242 The iteration is terminated if the residual norm is less then self.getTolerance().
243
244 :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.
245 :type u0: vector value on the domain (e.g. `escript.Data`).
246 :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.
247 :type p0: scalar value on the domain (e.g. `escript.Data`).
248 :param verbose: if set some information on iteration progress are printed
249 :type verbose: ``bool``
250 :return: flux and pressure
251 :rtype: ``tuple`` of `escript.Data`.
252
253 :note: The problem is solved as a least squares form
254 *(K^[-1]+D^* l2 D)u+G p=D^* l2 * f + K^[-1]g*
255 *G^*u+*G^* K Gp=G^*g*
256 where *D* is the *div* operator and *(Gp)_i=p_{,i}* for the permeability *K=k_{ij}*.
257 """
258 self.verbose=verbose
259 if self.useVPIteration:
260 return self.__solveVP(u0,p0,max_iter,max_num_corrections)
261 else:
262 X=self.__pde_k.createCoefficient("X")
263 Y=self.__pde_k.createCoefficient("Y")
264 Y[:self.domain.getDim()]=self.scale*util.tensor_mult(self.__permeability_inv,self.__g)
265 rtmp=self.__f * self.__l2 * self.scale
266 for i in range(self.domain.getDim()): X[i,i]=rtmp
267 X[self.domain.getDim(),:]=self.__g*self.scale
268 r=self.__pde_k.createCoefficient("r")
269 r[:self.domain.getDim()]=u0*self.scale
270 r[self.domain.getDim()]=p0
271 self.__pde_k.setValue(X=X, Y=Y, r=r)
272 self.__pde_k.getSolverOptions().setVerbosity(self.verbose)
273 #self.__pde_k.getSolverOptions().setPreconditioner(self.__pde_k.getSolverOptions().AMG)
274 self.__pde_k.getSolverOptions().setSolverMethod(self.__pde_k.getSolverOptions().DIRECT)
275 U=self.__pde_k.getSolution()
276 # self.__pde_k.getOperator().saveMM("k.mm")
277 u=U[:self.domain.getDim()]*(1./self.scale)
278 p=U[self.domain.getDim()]
279 if self.verbose:
280 KGp=util.tensor_mult(self.__permeability,util.grad(p)/self.scale)
281 def_p=self.__g-(u+KGp)
282 def_v=self.__f-util.div(u, self.__pde_k.getFunctionSpaceForCoefficient("X"))
283 print "DarcyFlux: L2: g-v-K*grad(p) = %e (v = %e)."%(self.__L2(def_p),self.__L2(u))
284 print "DarcyFlux: L2: f-div(v) = %e (grad(v) = %e)."%(self.__L2(def_v),self.__L2(util.grad(u)))
285 return u,p
286
287 def __solveVP(self,u0,p0, max_iter=100, max_num_corrections=10):
288 """
289 solves the problem.
290
291 The iteration is terminated if the residual norm is less than self.getTolerance().
292
293 :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.
294 :type u0: vector value on the domain (e.g. `escript.Data`).
295 :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.
296 :type p0: scalar value on the domain (e.g. `escript.Data`).
297 :return: flux and pressure
298 :rtype: ``tuple`` of `escript.Data`.
299
300 :note: The problem is solved as a least squares form
301 *(K^[-1]+D^* (DKD^*)^[-1] D)u+G p=D^* (DKD^*)^[-1] f + K^[-1]g*
302 *G^*u+*G^* K Gp=G^*g*
303 where *D* is the *div* operator and *(Gp)_i=p_{,i}* for the permeability *K=k_{ij}*.
304 """
305 rtol=self.getTolerance()
306 atol=self.getAbsoluteTolerance()
307 self.setSubProblemTolerance()
308 num_corrections=0
309 converged=False
310 norm_r=None
311
312 # Eliminate the hydrostatic pressure:
313 if self.verbose: print "DarcyFlux: calculate hydrostatic pressure component."
314 self.__pde_p.setValue(X=self.__g, r=p0, y=-util.inner(self.domain.getNormal(),u0))
315 p0=self.__pde_p.getSolution()
316 g2=self.__g - util.tensor_mult(self.__permeability, util.grad(p0))
317 norm_g2=util.integrate(util.inner(g2,util.tensor_mult(self.__permeability_inv,g2)))**0.5
318
319 p=p0*0
320 if self.solveForFlux:
321 v=u0.copy()
322 else:
323 v=self.__getFlux(p, u0, f=self.__f, g=g2)
324
325 while not converged and norm_g2 > 0:
326 Gp=util.grad(p)
327 KGp=util.tensor_mult(self.__permeability,Gp)
328 if self.verbose:
329 def_p=g2-(v+KGp)
330 def_v=self.__f-util.div(v)
331 print "DarcyFlux: L2: g-v-K*grad(p) = %e (v = %e)."%(self.__L2(def_p),self.__L2(v))
332 print "DarcyFlux: L2: f-div(v) = %e (grad(v) = %e)."%(self.__L2(def_v),self.__L2(util.grad(v)))
333 print "DarcyFlux: K^{-1}-norm of v = %e."%util.integrate(util.inner(v,util.tensor_mult(self.__permeability_inv,v)))**0.5
334 print "DarcyFlux: K^{-1}-norm of g2 = %e."%norm_g2
335 print "DarcyFlux: K-norm of grad(dp) = %e."%util.integrate(util.inner(Gp,KGp))**0.5
336 ATOL=atol+rtol*norm_g2
337 if self.verbose: print "DarcyFlux: absolute tolerance ATOL = %e."%(ATOL,)
338 if norm_r == None or norm_r>ATOL:
339 if num_corrections>max_num_corrections:
340 raise ValueError,"maximum number of correction steps reached."
341
342 if self.solveForFlux:
343 # initial residual is r=K^{-1}*(g-v-K*Gp)+D^*L^{-1}(f-Du)
344 v,r, norm_r=PCG(ArithmeticTuple(util.tensor_mult(self.__permeability_inv,g2-v)-Gp,self.__applWeight(v,self.__f),p),
345 self.__Aprod_v,
346 v,
347 self.__Msolve_PCG_v,
348 self.__inner_PCG_v,
349 atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
350 p=r[2]
351 else:
352 # initial residual is r=G^*(g2-KGp - v)
353 p,r, norm_r=PCG(ArithmeticTuple(g2-KGp,v),
354 self.__Aprod_p,
355 p,
356 self.__Msolve_PCG_p,
357 self.__inner_PCG_p,
358 atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
359 v=r[1]
360 if self.verbose: print "DarcyFlux: residual norm = %e."%norm_r
361 num_corrections+=1
362 else:
363 if self.verbose: print "DarcyFlux: stopping criterium reached."
364 converged=True
365 return v,p+p0
366
367 def __applWeight(self, v, f=None):
368 # solves L p = f-Dv with p = 0
369 if self.verbose: print "DarcyFlux: Applying weighting operator"
370 if f == None:
371 return -util.div(v)*self.__l2
372 else:
373 return (f-util.div(v))*self.__l2
374 def __getPressure(self, v, p0, g=None):
375 # solves (G*KG)p = G^(g-v) with p = p0 where location_of_fixed_pressure>0
376 if self.getSolverOptionsPressure().isVerbose() or self.verbose: print "DarcyFlux: Pressure update"
377 if g == None:
378 self.__pde_p.setValue(X=-v, r=p0)
379 else:
380 self.__pde_p.setValue(X=g-v, r=p0)
381 p=self.__pde_p.getSolution()
382 return p
383
384 def __Aprod_v(self,dv):
385 # calculates: (a,b,c) = (K^{-1}(dv + KG * dp), L^{-1}Ddv, dp) with (G*KG)dp = - G^*dv
386 dp=self.__getPressure(dv, p0=escript.Data()) # dp = (G*KG)^{-1} (0-G^*dv)
387 a=util.tensor_mult(self.__permeability_inv,dv)+util.grad(dp) # a= K^{-1}u+G*dp
388 b= - self.__applWeight(dv) # b = - (D K D^*)^{-1} (0-Dv)
389 return ArithmeticTuple(a,b,-dp)
390
391 def __Msolve_PCG_v(self,r):
392 # K^{-1} u = r[0] + D^*r[1] = K^{-1}(dv + KG * dp) + D^*L^{-1}Ddv
393 if self.getSolverOptionsFlux().isVerbose() or self.verbose: print "DarcyFlux: Applying preconditioner"
394 self.__pde_k.setValue(X=r[1]*util.kronecker(self.domain), Y=r[0], r=escript.Data())
395 return self.__pde_k.getSolution()
396
397 def __inner_PCG_v(self,v,r):
398 return util.integrate(util.inner(v,r[0])+util.div(v)*r[1])
399
400 def __Aprod_p(self,dp):
401 if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
402 Gdp=util.grad(dp)
403 self.__pde_k.setValue(Y=-Gdp,X=escript.Data(), r=escript.Data())
404 du=self.__pde_k.getSolution()
405 return ArithmeticTuple(util.tensor_mult(self.__permeability,Gdp),-du)
406
407 def __getFlux(self,p, v0, f=None, g=None):
408 # solves (K^{-1}+D^*L^{-1} D) v = D^*L^{-1}f + K^{-1}g - Gp
409 if f!=None:
410 self.__pde_k.setValue(X=self.__applWeight(v0*0,self.__f)*util.kronecker(self.domain))
411 self.__pde_k.setValue(r=v0)
412 g2=util.tensor_mult(self.__permeability_inv,g)
413 if p == None:
414 self.__pde_k.setValue(Y=g2)
415 else:
416 self.__pde_k.setValue(Y=g2-util.grad(p))
417 return self.__pde_k.getSolution()
418
419 #v=self.__getFlux(p, u0, f=self.__f, g=g2)
420 def __Msolve_PCG_p(self,r):
421 if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"
422 self.__pde_p.setValue(X=r[0]-r[1], Y=escript.Data(), r=escript.Data(), y=escript.Data())
423 return self.__pde_p.getSolution()
424
425 def __inner_PCG_p(self,p,r):
426 return util.integrate(util.inner(util.grad(p), r[0]-r[1]))
427
428 def __L2(self,v):
429 return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))
430
431 def __L2_r(self,v):
432 return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.ReducedFunction(self.domain)))**2))
433
434 def setTolerance(self,rtol=1e-4):
435 """
436 sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if
437
438 *|g-v-K gard(p)|_PCG <= atol + rtol * |K^{1/2}g2|_0*
439
440 where ``atol`` is an absolut tolerance (see `setAbsoluteTolerance`).
441
442 :param rtol: relative tolerance for the pressure
443 :type rtol: non-negative ``float``
444 """
445 if rtol<0:
446 raise ValueError,"Relative tolerance needs to be non-negative."
447 self.__rtol=rtol
448 def getTolerance(self):
449 """
450 returns the relative tolerance
451 :return: current relative tolerance
452 :rtype: ``float``
453 """
454 return self.__rtol
455
456 def setAbsoluteTolerance(self,atol=0.):
457 """
458 sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if
459
460 *|g-v-K gard(p)|_PCG <= atol + rtol * |K^{1/2}g2|_0*
461
462
463 where ``rtol`` is an absolut tolerance (see `setTolerance`), *|f|^2 = integrate(length(f)^2)* and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.
464
465 :param atol: absolute tolerance for the pressure
466 :type atol: non-negative ``float``
467 """
468 if atol<0:
469 raise ValueError,"Absolute tolerance needs to be non-negative."
470 self.__atol=atol
471 def getAbsoluteTolerance(self):
472 """
473 returns the absolute tolerance
474 :return: current absolute tolerance
475 :rtype: ``float``
476 """
477 return self.__atol
478 def getSubProblemTolerance(self):
479 """
480 Returns a suitable subtolerance
481 :type: ``float``
482 """
483 return max(util.EPSILON**(0.5),self.getTolerance()**2)
484
485 def setSubProblemTolerance(self):
486 """
487 Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected.
488 """
489 if self.__adaptSubTolerance:
490 sub_tol=self.getSubProblemTolerance()
491 self.getSolverOptionsFlux().setTolerance(sub_tol)
492 self.getSolverOptionsFlux().setAbsoluteTolerance(0.)
493 self.getSolverOptionsPressure().setTolerance(sub_tol)
494 self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
495 if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol
496
497
498 class DarcyFlowOld(object):
499 """
500 solves the problem
501
502 *u_i+k_{ij}*p_{,j} = g_i*
503 *u_{i,i} = f*
504
505 where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,
506
507 :note: The problem is solved in a least squares formulation.
508 """
509
510 def __init__(self, domain, weight=None, useReduced=False, adaptSubTolerance=True):
511 """
512 initializes the Darcy flux problem
513 :param domain: domain of the problem
514 :type domain: `Domain`
515 :param useReduced: uses reduced oreder on flux and pressure
516 :type useReduced: ``bool``
517 :param adaptSubTolerance: switches on automatic subtolerance selection
518 :type adaptSubTolerance: ``bool``
519 """
520 self.domain=domain
521 if weight == None:
522 s=self.domain.getSize()
523 self.__l2=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2
524 # self.__l2=(3.*util.longestEdge(self.domain))**2
525 #self.__l2=(0.1*util.longestEdge(self.domain)*s/util.sup(s))**2
526 else:
527 self.__l2=weight
528 self.__pde_v=LinearPDESystem(domain)
529 if useReduced: self.__pde_v.setReducedOrderOn()
530 self.__pde_v.setSymmetryOn()
531 self.__pde_v.setValue(D=util.kronecker(domain), A=self.__l2*util.outer(util.kronecker(domain),util.kronecker(domain)))
532 self.__pde_p=LinearSinglePDE(domain)
533 self.__pde_p.setSymmetryOn()
534 if useReduced: self.__pde_p.setReducedOrderOn()
535 self.__f=escript.Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
536 self.__g=escript.Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
537 self.setTolerance()
538 self.setAbsoluteTolerance()
539 self.__adaptSubTolerance=adaptSubTolerance
540 self.verbose=False
541 def getSolverOptionsFlux(self):
542 """
543 Returns the solver options used to solve the flux problems
544
545 *(I+D^*D)u=F*
546
547 :return: `SolverOptions`
548 """
549 return self.__pde_v.getSolverOptions()
550 def setSolverOptionsFlux(self, options=None):
551 """
552 Sets the solver options used to solve the flux problems
553
554 *(I+D^*D)u=F*
555
556 If ``options`` is not present, the options are reset to default
557 :param options: `SolverOptions`
558 :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
559 """
560 return self.__pde_v.setSolverOptions(options)
561 def getSolverOptionsPressure(self):
562 """
563 Returns the solver options used to solve the pressure problems
564
565 *(Q^*Q)p=Q^*G*
566
567 :return: `SolverOptions`
568 """
569 return self.__pde_p.getSolverOptions()
570 def setSolverOptionsPressure(self, options=None):
571 """
572 Sets the solver options used to solve the pressure problems
573
574 *(Q^*Q)p=Q^*G*
575
576 If ``options`` is not present, the options are reset to default
577 :param options: `SolverOptions`
578 :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
579 """
580 return self.__pde_p.setSolverOptions(options)
581
582 def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
583 """
584 assigns values to model parameters
585
586 :param f: volumetic sources/sinks
587 :type f: scalar value on the domain (e.g. `escript.Data`)
588 :param g: flux sources/sinks
589 :type g: vector values on the domain (e.g. `escript.Data`)
590 :param location_of_fixed_pressure: mask for locations where pressure is fixed
591 :type location_of_fixed_pressure: scalar value on the domain (e.g. `escript.Data`)
592 :param location_of_fixed_flux: mask for locations where flux is fixed.
593 :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)
594 :param permeability: permeability tensor. If scalar ``s`` is given the tensor with
595 ``s`` on the main diagonal is used. If vector ``v`` is given the tensor with
596 ``v`` on the main diagonal is used.
597 :type permeability: scalar, vector or tensor values on the domain (e.g. `escript.Data`)
598
599 :note: the values of parameters which are not set by calling ``setValue`` are not altered.
600 :note: at any point on the boundary of the domain the pressure (``location_of_fixed_pressure`` >0)
601 or the normal component of the flux (``location_of_fixed_flux[i]>0`` if direction of the normal
602 is along the *x_i* axis.
603 """
604 if f !=None:
605 f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
606 if f.isEmpty():
607 f=escript.Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
608 else:
609 if f.getRank()>0: raise ValueError,"illegal rank of f."
610 self.__f=f
611 if g !=None:
612 g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
613 if g.isEmpty():
614 g=escript.Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
615 else:
616 if not g.getShape()==(self.domain.getDim(),):
617 raise ValueError,"illegal shape of g"
618 self.__g=g
619
620 if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)
621 if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux)
622
623 if permeability!=None:
624 perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
625 if perm.getRank()==0:
626 perm=perm*util.kronecker(self.domain.getDim())
627 elif perm.getRank()==1:
628 perm, perm2=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A")), perm
629 for i in range(self.domain.getDim()): perm[i,i]=perm2[i]
630 elif perm.getRank()==2:
631 pass
632 else:
633 raise ValueError,"illegal rank of permeability."
634 self.__permeability=perm
635 self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))
636
637 def setTolerance(self,rtol=1e-4):
638 """
639 sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if
640
641 *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
642
643 where ``atol`` is an absolut tolerance (see `setAbsoluteTolerance`), *|f|^2 = integrate(length(f)^2)* and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.
644
645 :param rtol: relative tolerance for the pressure
646 :type rtol: non-negative ``float``
647 """
648 if rtol<0:
649 raise ValueError,"Relative tolerance needs to be non-negative."
650 self.__rtol=rtol
651 def getTolerance(self):
652 """
653 returns the relative tolerance
654
655 :return: current relative tolerance
656 :rtype: ``float``
657 """
658 return self.__rtol
659
660 def setAbsoluteTolerance(self,atol=0.):
661 """
662 sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if
663
664 *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
665
666 where ``rtol`` is an absolut tolerance (see `setTolerance`), *|f|^2 = integrate(length(f)^2)* and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.
667
668 :param atol: absolute tolerance for the pressure
669 :type atol: non-negative ``float``
670 """
671 if atol<0:
672 raise ValueError,"Absolute tolerance needs to be non-negative."
673 self.__atol=atol
674 def getAbsoluteTolerance(self):
675 """
676 returns the absolute tolerance
677
678 :return: current absolute tolerance
679 :rtype: ``float``
680 """
681 return self.__atol
682 def getSubProblemTolerance(self):
683 """
684 Returns a suitable subtolerance
685 @type: ``float``
686 """
687 return max(util.EPSILON**(0.75),self.getTolerance()**2)
688 def setSubProblemTolerance(self):
689 """
690 Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected.
691 """
692 if self.__adaptSubTolerance:
693 sub_tol=self.getSubProblemTolerance()
694 self.getSolverOptionsFlux().setTolerance(sub_tol)
695 self.getSolverOptionsFlux().setAbsoluteTolerance(0.)
696 self.getSolverOptionsPressure().setTolerance(sub_tol)
697 self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
698 if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol
699
700 def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):
701 """
702 solves the problem.
703
704 The iteration is terminated if the residual norm is less then self.getTolerance().
705
706 :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.
707 :type u0: vector value on the domain (e.g. `escript.Data`).
708 :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.
709 :type p0: scalar value on the domain (e.g. `escript.Data`).
710 :param verbose: if set some information on iteration progress are printed
711 :type verbose: ``bool``
712 :return: flux and pressure
713 :rtype: ``tuple`` of `escript.Data`.
714
715 :note: The problem is solved as a least squares form
716
717 *(I+D^*D)u+Qp=D^*f+g*
718 *Q^*u+Q^*Qp=Q^*g*
719
720 where *D* is the *div* operator and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.
721 We eliminate the flux form the problem by setting
722
723 *u=(I+D^*D)^{-1}(D^*f-g-Qp)* with u=u0 on location_of_fixed_flux
724
725 form the first equation. Inserted into the second equation we get
726
727 *Q^*(I-(I+D^*D)^{-1})Qp= Q^*(g-(I+D^*D)^{-1}(D^*f+g))* with p=p0 on location_of_fixed_pressure
728
729 which is solved using the PCG method (precondition is *Q^*Q*). In each iteration step
730 PDEs with operator *I+D^*D* and with *Q^*Q* needs to be solved using a sub iteration scheme.
731 """
732 self.verbose=verbose
733 rtol=self.getTolerance()
734 atol=self.getAbsoluteTolerance()
735 self.setSubProblemTolerance()
736 num_corrections=0
737 converged=False
738 p=p0
739 norm_r=None
740 while not converged:
741 v=self.getFlux(p, fixed_flux=u0)
742 Qp=self.__Q(p)
743 norm_v=self.__L2(v)
744 norm_Qp=self.__L2(Qp)
745 if norm_v == 0.:
746 if norm_Qp == 0.:
747 return v,p
748 else:
749 fac=norm_Qp
750 else:
751 if norm_Qp == 0.:
752 fac=norm_v
753 else:
754 fac=2./(1./norm_v+1./norm_Qp)
755 ATOL=(atol+rtol*fac)
756 if self.verbose:
757 print "DarcyFlux: L2 norm of v = %e."%norm_v
758 print "DarcyFlux: L2 norm of k.util.grad(p) = %e."%norm_Qp
759 print "DarcyFlux: L2 defect u = %e."%(util.integrate(util.length(self.__g-util.interpolate(v,escript.Function(self.domain))-Qp)**2)**(0.5),)
760 print "DarcyFlux: L2 defect div(v) = %e."%(util.integrate((self.__f-util.div(v))**2)**(0.5),)
761 print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
762 if norm_r == None or norm_r>ATOL:
763 if num_corrections>max_num_corrections:
764 raise ValueError,"maximum number of correction steps reached."
765 p,r, norm_r=PCG(self.__g-util.interpolate(v,escript.Function(self.domain))-Qp,self.__Aprod,p,self.__Msolve_PCG,self.__inner_PCG,atol=0.5*ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
766 num_corrections+=1
767 else:
768 converged=True
769 return v,p
770 def __L2(self,v):
771 return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))
772
773 def __Q(self,p):
774 return util.tensor_mult(self.__permeability,util.grad(p))
775
776 def __Aprod(self,dp):
777 if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
778 Qdp=self.__Q(dp)
779 self.__pde_v.setValue(Y=-Qdp,X=escript.Data(), r=escript.Data())
780 du=self.__pde_v.getSolution()
781 return Qdp+du
782 def __inner_GMRES(self,r,s):
783 return util.integrate(util.inner(r,s))
784
785 def __inner_PCG(self,p,r):
786 return util.integrate(util.inner(self.__Q(p), r))
787
788 def __Msolve_PCG(self,r):
789 if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"
790 self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=escript.Data(), r=escript.Data())
791 return self.__pde_p.getSolution()
792
793 def getFlux(self,p=None, fixed_flux=escript.Data()):
794 """
795 returns the flux for a given pressure ``p`` where the flux is equal to ``fixed_flux``
796 on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
797 Note that ``g`` and ``f`` are used, see `setValue`.
798
799 :param p: pressure.
800 :type p: scalar value on the domain (e.g. `escript.Data`).
801 :param fixed_flux: flux on the locations of the domain marked be ``location_of_fixed_flux``.
802 :type fixed_flux: vector values on the domain (e.g. `escript.Data`).
803 :return: flux
804 :rtype: `escript.Data`
805 :note: the method uses the least squares solution *u=(I+D^*D)^{-1}(D^*f-g-Qp)* where *D* is the *div* operator and *(Qp)_i=k_{ij}p_{,j}*
806 for the permeability *k_{ij}*
807 """
808 self.setSubProblemTolerance()
809 g=self.__g
810 f=self.__f
811 self.__pde_v.setValue(X=self.__l2*f*util.kronecker(self.domain), r=fixed_flux)
812 if p == None:
813 self.__pde_v.setValue(Y=g)
814 else:
815 self.__pde_v.setValue(Y=g-self.__Q(p))
816 return self.__pde_v.getSolution()
817
818 class StokesProblemCartesian(HomogeneousSaddlePointProblem):
819 """
820 solves
821
822 -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
823 u_{i,i}=0
824
825 u=0 where fixed_u_mask>0
826 eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j
827
828 if surface_stress is not given 0 is assumed.
829
830 typical usage:
831
832 sp=StokesProblemCartesian(domain)
833 sp.setTolerance()
834 sp.initialize(...)
835 v,p=sp.solve(v0,p0)
836 """
837 def __init__(self,domain,**kwargs):
838 """
839 initialize the Stokes Problem
840
841 The approximation spaces used for velocity (=Solution(domain)) and pressure (=ReducedSolution(domain)) must be
842 LBB complient, for instance using quadratic and linear approximation on the same element or using linear approximation
843 with macro elements for the pressure.
844
845 :param domain: domain of the problem.
846 :type domain: `Domain`
847 """
848 HomogeneousSaddlePointProblem.__init__(self,**kwargs)
849 self.domain=domain
850 self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
851 self.__pde_u.setSymmetryOn()
852
853 self.__pde_prec=LinearPDE(domain)
854 self.__pde_prec.setReducedOrderOn()
855 self.__pde_prec.setSymmetryOn()
856
857 self.__pde_proj=LinearPDE(domain)
858 self.__pde_proj.setReducedOrderOn()
859 self.__pde_proj.setValue(D=1)
860 self.__pde_proj.setSymmetryOn()
861
862 def getSolverOptionsVelocity(self):
863 """
864 returns the solver options used solve the equation for velocity.
865
866 :rtype: `SolverOptions`
867 """
868 return self.__pde_u.getSolverOptions()
869 def setSolverOptionsVelocity(self, options=None):
870 """
871 set the solver options for solving the equation for velocity.
872
873 :param options: new solver options
874 :type options: `SolverOptions`
875 """
876 self.__pde_u.setSolverOptions(options)
877 def getSolverOptionsPressure(self):
878 """
879 returns the solver options used solve the equation for pressure.
880 :rtype: `SolverOptions`
881 """
882 return self.__pde_prec.getSolverOptions()
883 def setSolverOptionsPressure(self, options=None):
884 """
885 set the solver options for solving the equation for pressure.
886 :param options: new solver options
887 :type options: `SolverOptions`
888 """
889 self.__pde_prec.setSolverOptions(options)
890
891 def setSolverOptionsDiv(self, options=None):
892 """
893 set the solver options for solving the equation to project the divergence of
894 the velocity onto the function space of presure.
895
896 :param options: new solver options
897 :type options: `SolverOptions`
898 """
899 self.__pde_proj.setSolverOptions(options)
900 def getSolverOptionsDiv(self):
901 """
902 returns the solver options for solving the equation to project the divergence of
903 the velocity onto the function space of presure.
904
905 :rtype: `SolverOptions`
906 """
907 return self.__pde_proj.getSolverOptions()
908
909 def updateStokesEquation(self, v, p):
910 """
911 updates the Stokes equation to consider dependencies from ``v`` and ``p``
912 :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values.
913 """
914 pass
915 def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None):
916 """
917 assigns new values to the model parameters.
918
919 :param f: external force
920 :type f: `Vector` object in `FunctionSpace` `Function` or similar
921 :param fixed_u_mask: mask of locations with fixed velocity.
922 :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
923 :param eta: viscosity
924 :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
925 :param surface_stress: normal surface stress
926 :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
927 :param stress: initial stress
928 :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
929 """
930 if eta !=None:
931 k=util.kronecker(self.domain.getDim())
932 kk=util.outer(k,k)
933 self.eta=util.interpolate(eta, escript.Function(self.domain))
934 self.__pde_prec.setValue(D=1/self.eta)
935 self.__pde_u.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))
936 if restoration_factor!=None:
937 n=self.domain.getNormal()
938 self.__pde_u.setValue(d=restoration_factor*util.outer(n,n))
939 if fixed_u_mask!=None:
940 self.__pde_u.setValue(q=fixed_u_mask)
941 if f!=None: self.__f=f
942 if surface_stress!=None: self.__surface_stress=surface_stress
943 if stress!=None: self.__stress=stress
944
945 def initialize(self,f=escript.Data(),fixed_u_mask=escript.Data(),eta=1,surface_stress=escript.Data(),stress=escript.Data(), restoration_factor=0):
946 """
947 assigns values to the model parameters
948
949 :param f: external force
950 :type f: `Vector` object in `FunctionSpace` `Function` or similar
951 :param fixed_u_mask: mask of locations with fixed velocity.
952 :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
953 :param eta: viscosity
954 :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
955 :param surface_stress: normal surface stress
956 :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
957 :param stress: initial stress
958 :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
959 """
960 self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
961
962 def Bv(self,v,tol):
963 """
964 returns inner product of element p and div(v)
965
966 :param v: a residual
967 :return: inner product of element p and div(v)
968 :rtype: ``float``
969 """
970 self.__pde_proj.setValue(Y=-util.div(v))
971 self.getSolverOptionsDiv().setTolerance(tol)
972 self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
973 out=self.__pde_proj.getSolution()
974 return out
975
976 def inner_pBv(self,p,Bv):
977 """
978 returns inner product of element p and Bv=-div(v)
979
980 :param p: a pressure increment
981 :param Bv: a residual
982 :return: inner product of element p and Bv=-div(v)
983 :rtype: ``float``
984 """
985 return util.integrate(util.interpolate(p,escript.Function(self.domain))*util.interpolate(Bv, escript.Function(self.domain)))
986
987 def inner_p(self,p0,p1):
988 """
989 Returns inner product of p0 and p1
990
991 :param p0: a pressure
992 :param p1: a pressure
993 :return: inner product of p0 and p1
994 :rtype: ``float``
995 """
996 s0=util.interpolate(p0, escript.Function(self.domain))
997 s1=util.interpolate(p1, escript.Function(self.domain))
998 return util.integrate(s0*s1)
999
1000 def norm_v(self,v):
1001 """
1002 returns the norm of v
1003
1004 :param v: a velovity
1005 :return: norm of v
1006 :rtype: non-negative ``float``
1007 """
1008 return util.sqrt(util.integrate(util.length(util.grad(v))**2))
1009
1010
1011 def getDV(self, p, v, tol):
1012 """
1013 return the value for v for a given p (overwrite)
1014
1015 :param p: a pressure
1016 :param v: a initial guess for the value v to return.
1017 :return: dv given as *Adv=(f-Av-B^*p)*
1018 """
1019 self.updateStokesEquation(v,p)
1020 self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress)
1021 self.getSolverOptionsVelocity().setTolerance(tol)
1022 self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
1023 if self.__stress.isEmpty():
1024 self.__pde_u.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
1025 else:
1026 self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
1027 out=self.__pde_u.getSolution()
1028 return out
1029
1030 def norm_Bv(self,Bv):
1031 """
1032 Returns Bv (overwrite).
1033
1034 :rtype: equal to the type of p
1035 :note: boundary conditions on p should be zero!
1036 """
1037 return util.sqrt(util.integrate(util.interpolate(Bv, escript.Function(self.domain))**2))
1038
1039 def solve_AinvBt(self,p, tol):
1040 """
1041 Solves *Av=B^*p* with accuracy `tol`
1042
1043 :param p: a pressure increment
1044 :return: the solution of *Av=B^*p*
1045 :note: boundary conditions on v should be zero!
1046 """
1047 self.__pde_u.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))
1048 out=self.__pde_u.getSolution()
1049 return out
1050
1051 def solve_prec(self,Bv, tol):
1052 """
1053 applies preconditioner for for *BA^{-1}B^** to *Bv*
1054 with accuracy `self.getSubProblemTolerance()`
1055
1056 :param Bv: velocity increment
1057 :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*
1058 :note: boundary conditions on p are zero.
1059 """
1060 self.__pde_prec.setValue(Y=Bv)
1061 self.getSolverOptionsPressure().setTolerance(tol)
1062 self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
1063 out=self.__pde_prec.getSolution()
1064 return out

  ViewVC Help
Powered by ViewVC 1.1.26