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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3499 - (hide annotations)
Fri Apr 8 00:14:35 2011 UTC (8 years, 5 months ago) by gross
File MIME type: text/x-python
File size: 54972 byte(s)


1 gross 3071 # -*- coding: utf-8 -*-
2 ksteube 1809 ########################################################
3 gross 1414 #
4 jfenwick 2881 # Copyright (c) 2003-2010 by University of Queensland
5 ksteube 1809 # Earth Systems Science Computational Center (ESSCC)
6     # http://www.uq.edu.au/esscc
7 gross 1414 #
8 ksteube 1809 # 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 gross 1414 #
12 ksteube 1809 ########################################################
13 gross 1414
14 jfenwick 2881 __copyright__="""Copyright (c) 2003-2010 by University of Queensland
15 ksteube 1809 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 jfenwick 2344 __url__="https://launchpad.net/escript-finley"
21 ksteube 1809
22 gross 1414 """
23     Some models for flow
24    
25 jfenwick 2625 :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 gross 1414 """
32    
33     __author__="Lutz Gross, l.gross@uq.edu.au"
34    
35 jfenwick 3432 import escript
36 gross 1414 import util
37 gross 2474 from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions
38 gross 2264 from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
39 gross 1414
40 gross 3499
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 gross 2100 class DarcyFlow(object):
282 gross 3071 """
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 gross 3440 def __init__(self, domain, useReduced=False, adaptSubTolerance=True, solveForFlux=False, useVPIteration=True, weighting_scale=0.1):
294 gross 3071 """
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 gross 3440 :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 gross 3071 """
307     self.domain=domain
308 gross 3440 self.useVPIteration=useVPIteration
309 gross 3071 self.useReduced=useReduced
310 gross 3440 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 gross 3071
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 gross 1659
384    
385 gross 3071 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 jfenwick 3360
389 gross 3071 :param f: volumetic sources/sinks
390 gross 3440 :type f: scalar value on the domain (e.g. `escript.Data`)
391 gross 3071 :param g: flux sources/sinks
392 gross 3440 :type g: vector values on the domain (e.g. `escript.Data`)
393 gross 3071 :param location_of_fixed_pressure: mask for locations where pressure is fixed
394 gross 3440 :type location_of_fixed_pressure: scalar value on the domain (e.g. `escript.Data`)
395 gross 3071 :param location_of_fixed_flux: mask for locations where flux is fixed.
396 gross 3440 :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)
397 jfenwick 3360 :param permeability: permeability tensor. If scalar ``s`` is given the tensor with ``s`` on the main diagonal is used.
398 gross 3440 :type permeability: scalar or tensor values on the domain (e.g. `escript.Data`)
399 jfenwick 3360
400 gross 3071 :note: the values of parameters which are not set by calling ``setValue`` are not altered.
401 caltinay 3452 :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 jfenwick 3360
406 gross 3071 """
407 gross 3440 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 gross 3071
417 gross 3440 # flux is rescaled by the factor mean value(perm_inv)*length where length**self.domain.getDim()=vol(self.domain)
418 gross 3071 if permeability!=None:
419 gross 3440 perm=util.interpolate(permeability,self.__pde_k.getFunctionSpaceForCoefficient("A"))
420     V=util.vol(self.domain)
421 gross 3071 if perm.getRank()==0:
422 gross 3440 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 gross 3071 elif perm.getRank()==2:
431     perm_inv=util.inverse(perm)
432 gross 3440 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 gross 3071 else:
439     raise ValueError,"illegal rank of permeability."
440 gross 2960
441 gross 3071 self.__permeability=perm
442     self.__permeability_inv=perm_inv
443 gross 3440
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 gross 3071 else:
465 gross 3440 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 gross 2960
478 gross 3071
479 gross 3440 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 gross 2960
485 gross 3440 :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 gross 2960
494 gross 3440 :note: The problem is solved as a least squares form
495 caltinay 3452 *(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 gross 3440 """
499     self.verbose=verbose
500     if self.useVPIteration:
501     return self.__solveVP(u0,p0,max_iter,max_num_corrections)
502 gross 3071 else:
503 gross 3440 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 gross 2960
528 gross 3440 def __solveVP(self,u0,p0, max_iter=100, max_num_corrections=10):
529 gross 3071 """
530     solves the problem.
531    
532 caltinay 3452 The iteration is terminated if the residual norm is less than self.getTolerance().
533 gross 2960
534 gross 3071 :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 gross 3440 :type u0: vector value on the domain (e.g. `escript.Data`).
536 gross 3071 :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 gross 3440 :type p0: scalar value on the domain (e.g. `escript.Data`).
538 gross 3071 :return: flux and pressure
539 gross 3440 :rtype: ``tuple`` of `escript.Data`.
540 gross 2960
541 gross 3071 :note: The problem is solved as a least squares form
542 caltinay 3452 *(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 gross 3071 """
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 gross 2960
560 gross 3071 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 gross 2960
566 gross 3071 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 gross 3440
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 gross 3071 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 gross 2960
679 gross 3071 *|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 gross 2960
697 gross 3071 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 gross 2960
703    
704 gross 3071 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 gross 2960
706 gross 3071 :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 gross 3074 return max(util.EPSILON**(0.5),self.getTolerance()**2)
725 gross 2960
726 gross 3071 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 gross 2960
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 gross 2354 s=self.domain.getSize()
764 gross 3440 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 gross 2351 else:
768 gross 3440 self.__l2=weight
769 gross 2100 self.__pde_v=LinearPDESystem(domain)
770 gross 2208 if useReduced: self.__pde_v.setReducedOrderOn()
771     self.__pde_v.setSymmetryOn()
772 gross 3440 self.__pde_v.setValue(D=util.kronecker(domain), A=self.__l2*util.outer(util.kronecker(domain),util.kronecker(domain)))
773 gross 2100 self.__pde_p=LinearSinglePDE(domain)
774     self.__pde_p.setSymmetryOn()
775 gross 2208 if useReduced: self.__pde_p.setReducedOrderOn()
776 jfenwick 3432 self.__f=escript.Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
777     self.__g=escript.Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
778 gross 2264 self.setTolerance()
779     self.setAbsoluteTolerance()
780 gross 2474 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 jfenwick 2625 *(I+D^*D)u=F*
787 gross 2474
788 jfenwick 2625 :return: `SolverOptions`
789 gross 2474 """
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 jfenwick 2625 *(I+D^*D)u=F*
796 gross 2474
797 jfenwick 2625 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 gross 2474 """
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 jfenwick 2625 *(Q^*Q)p=Q^*G*
807 gross 2474
808 jfenwick 2625 :return: `SolverOptions`
809 gross 2474 """
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 jfenwick 2625 *(Q^*Q)p=Q^*G*
816 gross 2474
817 jfenwick 2625 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 gross 2474 """
821     return self.__pde_p.setSolverOptions(options)
822 gross 1659
823 gross 2100 def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
824     """
825 gross 2208 assigns values to model parameters
826 gross 1659
827 jfenwick 2625 :param f: volumetic sources/sinks
828 gross 3440 :type f: scalar value on the domain (e.g. `escript.Data`)
829 jfenwick 2625 :param g: flux sources/sinks
830 gross 3440 :type g: vector values on the domain (e.g. `escript.Data`)
831 jfenwick 2625 :param location_of_fixed_pressure: mask for locations where pressure is fixed
832 gross 3440 :type location_of_fixed_pressure: scalar value on the domain (e.g. `escript.Data`)
833 jfenwick 2625 :param location_of_fixed_flux: mask for locations where flux is fixed.
834 gross 3440 :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)
835 jfenwick 2625 :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 gross 3440 :type permeability: scalar, vector or tensor values on the domain (e.g. `escript.Data`)
839 gross 1659
840 jfenwick 2625 :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 gross 2100 """
845 gross 2264 if f !=None:
846 gross 2100 f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
847     if f.isEmpty():
848 jfenwick 3432 f=escript.Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
849 gross 2100 else:
850     if f.getRank()>0: raise ValueError,"illegal rank of f."
851 gross 2267 self.__f=f
852 gross 2264 if g !=None:
853 gross 2100 g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
854     if g.isEmpty():
855 jfenwick 3432 g=escript.Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
856 gross 2100 else:
857     if not g.getShape()==(self.domain.getDim(),):
858     raise ValueError,"illegal shape of g"
859     self.__g=g
860 gross 1659
861 gross 2100 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 gross 1659
864 gross 2100 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 gross 1659
878 gross 2264 def setTolerance(self,rtol=1e-4):
879     """
880 jfenwick 2625 sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if
881 gross 1659
882 jfenwick 2625 *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
883 gross 2264
884 jfenwick 2625 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 gross 2264
886 jfenwick 2625 :param rtol: relative tolerance for the pressure
887     :type rtol: non-negative ``float``
888 gross 2100 """
889 gross 2264 if rtol<0:
890     raise ValueError,"Relative tolerance needs to be non-negative."
891     self.__rtol=rtol
892     def getTolerance(self):
893 gross 2100 """
894 gross 2264 returns the relative tolerance
895 gross 1659
896 jfenwick 2625 :return: current relative tolerance
897     :rtype: ``float``
898 caltinay 2169 """
899 gross 2264 return self.__rtol
900 gross 1659
901 gross 2264 def setAbsoluteTolerance(self,atol=0.):
902 gross 2208 """
903 jfenwick 2625 sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if
904 gross 1659
905 jfenwick 2625 *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
906 gross 2208
907 jfenwick 2625 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 gross 2208
909 jfenwick 2625 :param atol: absolute tolerance for the pressure
910     :type atol: non-negative ``float``
911 gross 2264 """
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 jfenwick 2625 :return: current absolute tolerance
920     :rtype: ``float``
921 gross 2264 """
922     return self.__atol
923     def getSubProblemTolerance(self):
924 gross 2474 """
925     Returns a suitable subtolerance
926 jfenwick 2625 @type: ``float``
927 gross 2474 """
928     return max(util.EPSILON**(0.75),self.getTolerance()**2)
929     def setSubProblemTolerance(self):
930 gross 2264 """
931 gross 2474 Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected.
932 gross 2264 """
933 gross 2474 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 gross 2208
941 gross 2474 def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):
942 gross 2264 """
943     solves the problem.
944 gross 1659
945 gross 2208 The iteration is terminated if the residual norm is less then self.getTolerance().
946 gross 1659
947 jfenwick 2625 :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 gross 3440 :type u0: vector value on the domain (e.g. `escript.Data`).
949 jfenwick 2625 :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 gross 3440 :type p0: scalar value on the domain (e.g. `escript.Data`).
951 jfenwick 2625 :param verbose: if set some information on iteration progress are printed
952     :type verbose: ``bool``
953     :return: flux and pressure
954 gross 3440 :rtype: ``tuple`` of `escript.Data`.
955 gross 2264
956 jfenwick 2625 :note: The problem is solved as a least squares form
957 gross 2100
958 jfenwick 2625 *(I+D^*D)u+Qp=D^*f+g*
959     *Q^*u+Q^*Qp=Q^*g*
960 gross 2100
961 jfenwick 2625 where *D* is the *div* operator and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.
962 gross 2264 We eliminate the flux form the problem by setting
963 caltinay 2169
964 jfenwick 2625 *u=(I+D^*D)^{-1}(D^*f-g-Qp)* with u=u0 on location_of_fixed_flux
965 caltinay 2169
966 gross 2208 form the first equation. Inserted into the second equation we get
967 caltinay 2169
968 jfenwick 2625 *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 gross 2264
970 jfenwick 2625 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 gross 2208 """
973 gross 2386 self.verbose=verbose
974 gross 2264 rtol=self.getTolerance()
975     atol=self.getAbsoluteTolerance()
976 gross 2474 self.setSubProblemTolerance()
977 gross 2264 num_corrections=0
978     converged=False
979     p=p0
980     norm_r=None
981     while not converged:
982 gross 2474 v=self.getFlux(p, fixed_flux=u0)
983 gross 2264 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 gross 2960 print "DarcyFlux: L2 norm of k.util.grad(p) = %e."%norm_Qp
1000 jfenwick 3432 print "DarcyFlux: L2 defect u = %e."%(util.integrate(util.length(self.__g-util.interpolate(v,escript.Function(self.domain))-Qp)**2)**(0.5),)
1001 gross 2486 print "DarcyFlux: L2 defect div(v) = %e."%(util.integrate((self.__f-util.div(v))**2)**(0.5),)
1002 gross 2264 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 jfenwick 3432 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 gross 2264 num_corrections+=1
1008     else:
1009     converged=True
1010     return v,p
1011     def __L2(self,v):
1012 jfenwick 3432 return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))
1013 gross 2264
1014     def __Q(self,p):
1015     return util.tensor_mult(self.__permeability,util.grad(p))
1016    
1017     def __Aprod(self,dp):
1018 gross 2474 if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
1019 gross 2264 Qdp=self.__Q(dp)
1020 jfenwick 3432 self.__pde_v.setValue(Y=-Qdp,X=escript.Data(), r=escript.Data())
1021 gross 2474 du=self.__pde_v.getSolution()
1022 gross 2264 return Qdp+du
1023     def __inner_GMRES(self,r,s):
1024     return util.integrate(util.inner(r,s))
1025    
1026 gross 2100 def __inner_PCG(self,p,r):
1027 gross 2264 return util.integrate(util.inner(self.__Q(p), r))
1028 gross 2100
1029     def __Msolve_PCG(self,r):
1030 gross 2474 if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"
1031 jfenwick 3432 self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=escript.Data(), r=escript.Data())
1032 gross 2474 return self.__pde_p.getSolution()
1033 gross 2100
1034 jfenwick 3432 def getFlux(self,p=None, fixed_flux=escript.Data()):
1035 gross 2264 """
1036 jfenwick 2625 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 gross 2264
1040 jfenwick 2625 :param p: pressure.
1041 gross 3440 :type p: scalar value on the domain (e.g. `escript.Data`).
1042 jfenwick 2625 :param fixed_flux: flux on the locations of the domain marked be ``location_of_fixed_flux``.
1043 gross 3440 :type fixed_flux: vector values on the domain (e.g. `escript.Data`).
1044 jfenwick 2625 :return: flux
1045 gross 3440 :rtype: `escript.Data`
1046 jfenwick 2625 :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 gross 2264 """
1049 gross 2474 self.setSubProblemTolerance()
1050 gross 2264 g=self.__g
1051     f=self.__f
1052 gross 3440 self.__pde_v.setValue(X=self.__l2*f*util.kronecker(self.domain), r=fixed_flux)
1053 gross 2264 if p == None:
1054     self.__pde_v.setValue(Y=g)
1055     else:
1056     self.__pde_v.setValue(Y=g-self.__Q(p))
1057 gross 2474 return self.__pde_v.getSolution()
1058 gross 2264
1059 gross 1414 class StokesProblemCartesian(HomogeneousSaddlePointProblem):
1060 gross 2251 """
1061 gross 2264 solves
1062 gross 1414
1063 gross 2208 -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
1064     u_{i,i}=0
1065 gross 1414
1066 gross 2208 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 gross 1414
1069 gross 2264 if surface_stress is not given 0 is assumed.
1070 gross 1414
1071 gross 2251 typical usage:
1072 gross 1414
1073 gross 2208 sp=StokesProblemCartesian(domain)
1074     sp.setTolerance()
1075     sp.initialize(...)
1076     v,p=sp.solve(v0,p0)
1077 gross 2251 """
1078 gross 2719 def __init__(self,domain,**kwargs):
1079 gross 2100 """
1080 gross 2208 initialize the Stokes Problem
1081 gross 2100
1082 gross 2793 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 jfenwick 2625 :type domain: `Domain`
1088 gross 2100 """
1089 gross 2719 HomogeneousSaddlePointProblem.__init__(self,**kwargs)
1090 gross 1414 self.domain=domain
1091     self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
1092     self.__pde_u.setSymmetryOn()
1093 gross 2474
1094 gross 1414 self.__pde_prec=LinearPDE(domain)
1095     self.__pde_prec.setReducedOrderOn()
1096     self.__pde_prec.setSymmetryOn()
1097    
1098 gross 2415 self.__pde_proj=LinearPDE(domain)
1099     self.__pde_proj.setReducedOrderOn()
1100     self.__pde_proj.setValue(D=1)
1101     self.__pde_proj.setSymmetryOn()
1102    
1103 gross 2474 def getSolverOptionsVelocity(self):
1104     """
1105     returns the solver options used solve the equation for velocity.
1106    
1107 jfenwick 2625 :rtype: `SolverOptions`
1108 gross 2474 """
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 jfenwick 2625 :param options: new solver options
1115     :type options: `SolverOptions`
1116 gross 2474 """
1117     self.__pde_u.setSolverOptions(options)
1118     def getSolverOptionsPressure(self):
1119     """
1120     returns the solver options used solve the equation for pressure.
1121 jfenwick 2625 :rtype: `SolverOptions`
1122 gross 2474 """
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 jfenwick 2625 :param options: new solver options
1128     :type options: `SolverOptions`
1129 gross 2474 """
1130     self.__pde_prec.setSolverOptions(options)
1131 gross 2415
1132 gross 2474 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 jfenwick 2625 :param options: new solver options
1138     :type options: `SolverOptions`
1139 gross 2474 """
1140 artak 2689 self.__pde_proj.setSolverOptions(options)
1141 gross 2474 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 jfenwick 2625 :rtype: `SolverOptions`
1147 gross 2474 """
1148 artak 2689 return self.__pde_proj.getSolverOptions()
1149 gross 2474
1150 gross 2793 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 jfenwick 3432 self.eta=util.interpolate(eta, escript.Function(self.domain))
1175 gross 2793 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 jfenwick 3432 def initialize(self,f=escript.Data(),fixed_u_mask=escript.Data(),eta=1,surface_stress=escript.Data(),stress=escript.Data(), restoration_factor=0):
1187 gross 2208 """
1188     assigns values to the model parameters
1189 gross 2100
1190 jfenwick 2625 :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 gross 2208 """
1201 gross 2793 self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
1202 gross 1414
1203 gross 2719 def Bv(self,v,tol):
1204 gross 2251 """
1205     returns inner product of element p and div(v)
1206 gross 1414
1207 jfenwick 2625 :param v: a residual
1208     :return: inner product of element p and div(v)
1209     :rtype: ``float``
1210 gross 2100 """
1211 gross 2793 self.__pde_proj.setValue(Y=-util.div(v))
1212 gross 2719 self.getSolverOptionsDiv().setTolerance(tol)
1213     self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
1214     out=self.__pde_proj.getSolution()
1215     return out
1216 gross 2208
1217 gross 2445 def inner_pBv(self,p,Bv):
1218     """
1219     returns inner product of element p and Bv=-div(v)
1220    
1221 jfenwick 2625 :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 gross 2445 """
1226 jfenwick 3432 return util.integrate(util.interpolate(p,escript.Function(self.domain))*util.interpolate(Bv, escript.Function(self.domain)))
1227 gross 2445
1228 gross 2251 def inner_p(self,p0,p1):
1229 gross 2100 """
1230 gross 2251 Returns inner product of p0 and p1
1231 gross 1414
1232 jfenwick 2625 :param p0: a pressure
1233     :param p1: a pressure
1234     :return: inner product of p0 and p1
1235     :rtype: ``float``
1236 gross 2100 """
1237 jfenwick 3432 s0=util.interpolate(p0, escript.Function(self.domain))
1238     s1=util.interpolate(p1, escript.Function(self.domain))
1239 gross 2100 return util.integrate(s0*s1)
1240 artak 1550
1241 gross 2251 def norm_v(self,v):
1242 gross 2100 """
1243 gross 2251 returns the norm of v
1244 gross 2208
1245 jfenwick 2625 :param v: a velovity
1246     :return: norm of v
1247     :rtype: non-negative ``float``
1248 gross 2100 """
1249 gross 2719 return util.sqrt(util.integrate(util.length(util.grad(v))**2))
1250 gross 2100
1251 gross 2793
1252 gross 2719 def getDV(self, p, v, tol):
1253 gross 1414 """
1254 gross 2251 return the value for v for a given p (overwrite)
1255    
1256 jfenwick 2625 :param p: a pressure
1257 gross 2719 :param v: a initial guess for the value v to return.
1258     :return: dv given as *Adv=(f-Av-B^*p)*
1259 gross 1414 """
1260 gross 2793 self.updateStokesEquation(v,p)
1261 gross 2719 self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress)
1262     self.getSolverOptionsVelocity().setTolerance(tol)
1263     self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
1264 gross 2100 if self.__stress.isEmpty():
1265 gross 2719 self.__pde_u.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
1266 gross 2100 else:
1267 gross 2719 self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
1268 gross 2474 out=self.__pde_u.getSolution()
1269 gross 2208 return out
1270 gross 1414
1271 gross 2445 def norm_Bv(self,Bv):
1272 gross 2251 """
1273     Returns Bv (overwrite).
1274    
1275 jfenwick 2625 :rtype: equal to the type of p
1276     :note: boundary conditions on p should be zero!
1277 gross 2251 """
1278 jfenwick 3432 return util.sqrt(util.integrate(util.interpolate(Bv, escript.Function(self.domain))**2))
1279 gross 2251
1280 gross 2719 def solve_AinvBt(self,p, tol):
1281 gross 2251 """
1282 gross 2719 Solves *Av=B^*p* with accuracy `tol`
1283 gross 2251
1284 jfenwick 2625 :param p: a pressure increment
1285     :return: the solution of *Av=B^*p*
1286     :note: boundary conditions on v should be zero!
1287 gross 2251 """
1288 jfenwick 3432 self.__pde_u.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))
1289 gross 2474 out=self.__pde_u.getSolution()
1290 gross 2251 return out
1291    
1292 gross 2719 def solve_prec(self,Bv, tol):
1293 gross 2251 """
1294 jfenwick 2625 applies preconditioner for for *BA^{-1}B^** to *Bv*
1295     with accuracy `self.getSubProblemTolerance()`
1296 gross 2251
1297 jfenwick 2625 :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 gross 2251 """
1301 gross 2445 self.__pde_prec.setValue(Y=Bv)
1302 gross 2719 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