/[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 3501 - (show annotations)
Wed Apr 13 04:07:53 2011 UTC (8 years, 5 months ago) by gross
File MIME type: text/x-python
File size: 55236 byte(s)


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

  ViewVC Help
Powered by ViewVC 1.1.26