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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26