# Diff of /trunk/escript/py_src/flows.py

revision 3074 by gross, Tue Jul 27 01:47:45 2010 UTC revision 3771 by jfenwick, Wed Jan 18 02:30:48 2012 UTC
# Line 32  Some models for flow Line 32  Some models for flow
32
33  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
34
35  from escript import *  from . import escript
36  import util  from . import util
37  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions  from .linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions
38  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES  from .pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
39
40  class DarcyFlow(object):  class DarcyFlow(object):
41     """     """
# Line 46  class DarcyFlow(object): Line 46  class DarcyFlow(object):
46
47     where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,     where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,
48
49     :note: The problem is solved in a least squares formulation.     :cvar EVAL: direct pressure gradient evaluation for flux
50       :cvar POST: global postprocessing of flux by solving the PDE *K_{ij} u_j + (w * K * l u_{k,k})_{,i}= - p_{,j} + K_{ij} g_j*
51                   where *l* is the length scale, *K* is the inverse of the permeability tensor, and *w* is a positive weighting factor.
52       :cvar SMOOTH: global smoothing by solving the PDE *K_{ij} u_j= - p_{,j} + K_{ij} g_j*
53     """     """
54         EVAL="EVAL"
55     def __init__(self, domain, useReduced=False, adaptSubTolerance=True, solveForFlux=False):     SIMPLE="EVAL"
56       POST="POST"
57       SMOOTH="SMOOTH"
58       def __init__(self, domain, useReduced=False, solver="POST", verbose=False, w=1.):
59        """        """
60        initializes the Darcy flux problem        initializes the Darcy flux problem
61        :param domain: domain of the problem        :param domain: domain of the problem
62        :type domain: `Domain`        :type domain: `Domain`
63        :param useReduced: uses reduced oreder on flux and pressure        :param useReduced: uses reduced oreder on flux and pressure
64        :type useReduced: ``bool``        :type useReduced: ``bool``
65        :param adaptSubTolerance: switches on automatic subtolerance selection        :param solver: solver method
66        :type adaptSubTolerance: ``bool``        :type solver: in [`DarcyFlow.EVAL`, `DarcyFlow.POST',  `DarcyFlow.SMOOTH' ]
67        :param solveForFlux: if True the solver solves for the flux (do not use!)        :param verbose: if ``True`` some information on the iteration progress are printed.
68        :type solveForFlux: ``bool``          :type verbose: ``bool``
69          :param w: weighting factor for `DarcyFlow.POST` solver
70          :type w: ``float``
71
72        """        """
73          if not solver in [DarcyFlow.EVAL, DarcyFlow.POST,  DarcyFlow.SMOOTH ] :
74              raise ValueError("unknown solver %d."%solver)
75
76        self.domain=domain        self.domain=domain
77        self.solveForFlux=solveForFlux        self.solver=solver
78        self.useReduced=useReduced        self.useReduced=useReduced
80        self.verbose=False        self.l=None
81                self.w=None
82        self.__pde_k=LinearPDESystem(domain)
self.__pde_k.setSymmetryOn()
if self.useReduced: self.__pde_k.setReducedOrderOn()

83        self.__pde_p=LinearSinglePDE(domain)        self.__pde_p=LinearSinglePDE(domain)
84        self.__pde_p.setSymmetryOn()        self.__pde_p.setSymmetryOn()
85        if self.useReduced: self.__pde_p.setReducedOrderOn()        if self.useReduced: self.__pde_p.setReducedOrderOn()
86
87        self.__pde_l=LinearSinglePDE(domain)   # this is here for getSolverOptionsWeighting        if self.solver  == self.EVAL:
88        # self.__pde_l.setSymmetryOn()           self.__pde_v=None
89        # if self.useReduced: self.__pde_l.setReducedOrderOn()       if self.verbose: print("DarcyFlow: simple solver is used.")
90        self.setTolerance()
91        self.setAbsoluteTolerance()        elif self.solver  == self.POST:
92        self.__f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))       if util.inf(w)<0.:
93        self.__g=Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))          raise ValueError("Weighting factor must be non-negative.")
94               if self.verbose: print("DarcyFlow: global postprocessing of flux is used.")
95             self.__pde_v=LinearPDESystem(domain)
96             self.__pde_v.setSymmetryOn()
97             if self.useReduced: self.__pde_v.setReducedOrderOn()
98         self.w=w
99             self.l=util.vol(self.domain)**(1./self.domain.getDim()) # length scale
100
101          elif self.solver  == self.SMOOTH:
102             self.__pde_v=LinearPDESystem(domain)
103             self.__pde_v.setSymmetryOn()
104             if self.useReduced: self.__pde_v.setReducedOrderOn()
105         if self.verbose: print("DarcyFlow: flux smoothing is used.")
106         self.w=0
107
108          self.__f=escript.Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))
109          self.__g=escript.Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
110          self.location_of_fixed_pressure = escript.Scalar(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
111          self.location_of_fixed_flux = escript.Vector(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
112          self.perm_scale=1.
113
114
115       def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
116          """
117          assigns values to model parameters
118
119          :param f: volumetic sources/sinks
120          :type f: scalar value on the domain (e.g. `escript.Data`)
121          :param g: flux sources/sinks
122          :type g: vector values on the domain (e.g. `escript.Data`)
123          :param location_of_fixed_pressure: mask for locations where pressure is fixed
124          :type location_of_fixed_pressure: scalar value on the domain (e.g. `escript.Data`)
125          :param location_of_fixed_flux:  mask for locations where flux is fixed.
126          :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)
127          :param permeability: permeability tensor. If scalar ``s`` is given the tensor with ``s`` on the main diagonal is used.
128          :type permeability: scalar or symmetric tensor values on the domain (e.g. `escript.Data`)
129
130          :note: the values of parameters which are not set by calling ``setValue`` are not altered.
131          :note: at any point on the boundary of the domain the pressure
132                 (``location_of_fixed_pressure`` >0) or the normal component of the
133                 flux (``location_of_fixed_flux[i]>0``) if direction of the normal
134                 is along the *x_i* axis.
135
136          """
137          if location_of_fixed_pressure!=None:
138               self.location_of_fixed_pressure=util.wherePositive(location_of_fixed_pressure)
139               self.__pde_p.setValue(q=self.location_of_fixed_pressure)
140          if location_of_fixed_flux!=None:
141              self.location_of_fixed_flux=util.wherePositive(location_of_fixed_flux)
142              if not self.__pde_v == None: self.__pde_v.setValue(q=self.location_of_fixed_flux)
143
144          if permeability!=None:
145
146         perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
147             self.perm_scale=util.Lsup(util.length(perm))
148         if self.verbose: print(("DarcyFlow: permeability scaling factor = %e."%self.perm_scale))
149             perm=perm*(1./self.perm_scale)
150
151         if perm.getRank()==0:
152
153            perm_inv=(1./perm)
154            perm_inv=perm_inv*util.kronecker(self.domain.getDim())
155            perm=perm*util.kronecker(self.domain.getDim())
156
157
158         elif perm.getRank()==2:
159            perm_inv=util.inverse(perm)
160         else:
161            raise ValueError("illegal rank of permeability.")
162
163         self.__permeability=perm
164         self.__permeability_inv=perm_inv
165
166             #====================
167         self.__pde_p.setValue(A=self.__permeability)
168             if self.solver  == self.EVAL:
169                  pass # no extra work required
170             elif self.solver  == self.POST:
171            k=util.kronecker(self.domain.getDim())
172            self.omega = self.w*util.length(perm_inv)*self.l*self.domain.getSize()
173            self.__pde_v.setValue(D=self.__permeability_inv, A=self.omega*util.outer(k,k))
174             elif self.solver  == self.SMOOTH:
175            self.__pde_v.setValue(D=self.__permeability_inv)
176
177          if g != None:
178        g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
179        if g.isEmpty():
180              g=Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
181        else:
182            if not g.getShape()==(self.domain.getDim(),): raise ValueError("illegal shape of g")
183        self.__g=g
184          if f !=None:
185         f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
186         if f.isEmpty():
187              f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
188         else:
189             if f.getRank()>0: raise ValueError("illegal rank of f.")
190         self.__f=f
191
192     def getSolverOptionsFlux(self):     def getSolverOptionsFlux(self):
193        """        """
194        Returns the solver options used to solve the flux problems        Returns the solver options used to solve the flux problems

*K^{-1} u=F*

195        :return: `SolverOptions`        :return: `SolverOptions`
196        """        """
197        return self.__pde_k.getSolverOptions()        if self.__pde_v == None:
198              return None
199          else:
200              return self.__pde_v.getSolverOptions()
201
202     def setSolverOptionsFlux(self, options=None):     def setSolverOptionsFlux(self, options=None):
203        """        """
204        Sets the solver options used to solve the flux problems        Sets the solver options used to solve the flux problems

*K^{-1}u=F*

205        If ``options`` is not present, the options are reset to default        If ``options`` is not present, the options are reset to default

206        :param options: `SolverOptions`        :param options: `SolverOptions`
:note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
207        """        """
208        return self.__pde_v.setSolverOptions(options)        if not self.__pde_v == None:
209              self.__pde_v.setSolverOptions(options)
210
211     def getSolverOptionsPressure(self):     def getSolverOptionsPressure(self):
212        """        """
213        Returns the solver options used to solve the pressure problems        Returns the solver options used to solve the pressure problems

*(Q^* K Q)p=-Q^*G*

214        :return: `SolverOptions`        :return: `SolverOptions`
215        """        """
216        return self.__pde_p.getSolverOptions()        return self.__pde_p.getSolverOptions()
# Line 119  class DarcyFlow(object): Line 218  class DarcyFlow(object):
218     def setSolverOptionsPressure(self, options=None):     def setSolverOptionsPressure(self, options=None):
219        """        """
220        Sets the solver options used to solve the pressure problems        Sets the solver options used to solve the pressure problems

*(Q^* K Q)p=-Q^*G*

221        If ``options`` is not present, the options are reset to default        If ``options`` is not present, the options are reset to default
222
223        :param options: `SolverOptions`        :param options: `SolverOptions`
224        :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.        :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
225        """        """
226        return self.__pde_p.setSolverOptions(options)        return self.__pde_p.setSolverOptions(options)

def getSolverOptionsWeighting(self):
"""
Returns the solver options used to solve the pressure problems

*(D K D^*)p=-D F*

:return: `SolverOptions`
"""
return self.__pde_l.getSolverOptions()

def setSolverOptionsWeighting(self, options=None):
"""
Sets the solver options used to solve the pressure problems

*(D K D^*)p=-D F*

If ``options`` is not present, the options are reset to default

:param options: `SolverOptions`
:note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
"""
return self.__pde_l.setSolverOptions(options)

def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
"""
assigns values to model parameters

:param f: volumetic sources/sinks
:type f: scalar value on the domain (e.g. `Data`)
:param g: flux sources/sinks
:type g: vector values on the domain (e.g. `Data`)
:param location_of_fixed_pressure: mask for locations where pressure is fixed
:type location_of_fixed_pressure: scalar value on the domain (e.g. `Data`)
:param location_of_fixed_flux:  mask for locations where flux is fixed.
:type location_of_fixed_flux: vector values on the domain (e.g. `Data`)
:param permeability: permeability tensor. If scalar ``s`` is given the tensor with
``s`` on the main diagonal is used.
:type permeability: scalar or tensor values on the domain (e.g. `Data`)
:note: the values of parameters which are not set by calling ``setValue`` are not altered.
:note: at any point on the boundary of the domain the pressure (``location_of_fixed_pressure`` >0)
or the normal component of the flux (``location_of_fixed_flux[i]>0`` if direction of the normal
is along the *x_i* axis.
"""
if f !=None:
f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("X"))
if f.isEmpty():
f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))
else:
if f.getRank()>0: raise ValueError,"illegal rank of f."
self.__f=f
if g !=None:
g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
if g.isEmpty():
g=Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
else:
if not g.getShape()==(self.domain.getDim(),):
raise ValueError,"illegal shape of g"
self.__g=g
if location_of_fixed_pressure!=None:
self.__pde_p.setValue(q=location_of_fixed_pressure)
#self.__pde_l.setValue(q=location_of_fixed_pressure)
if location_of_fixed_flux!=None:
self.__pde_k.setValue(q=location_of_fixed_flux)

if permeability!=None:
perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
if perm.getRank()==0:
perm_inv=(1./perm)*util.kronecker(self.domain.getDim())
perm=perm*util.kronecker(self.domain.getDim())
elif perm.getRank()==2:
perm_inv=util.inverse(perm)
else:
raise ValueError,"illegal rank of permeability."

self.__permeability=perm
self.__permeability_inv=perm_inv
self.__l =(util.longestEdge(self.domain)**2*util.length(self.__permeability_inv))/10
#self.__l =(self.domain.getSize()**2*util.length(self.__permeability_inv))/10
if  self.solveForFlux:
self.__pde_k.setValue(D=self.__permeability_inv)
else:
self.__pde_k.setValue(D=self.__permeability_inv, A=self.__l*util.outer(util.kronecker(self.domain),util.kronecker(self.domain)))
self.__pde_p.setValue(A=self.__permeability)
#self.__pde_l.setValue(D=1/self.__l)
#self.__pde_l.setValue(A=self.__permeability)

def __applWeight(self, v, f=None):
# solves L p = f-Dv with p = 0
if self.getSolverOptionsWeighting().isVerbose() or self.verbose: print "DarcyFlux: Applying weighting operator"
if f == None:
return -util.div(v)*self.__l
else:
return (f-util.div(v))*self.__l
# if f == None:
#      self.__pde_l.setValue(Y=-util.div(v))
# else:
#      return (f-util.div(v))/self.__l
# return self.__pde_l.getSolution()
227
228     def __getPressure(self, v, p0, g=None):     def solve(self, u0, p0):
# solves (G*KG)p = G^(g-v) with p = p0 where location_of_fixed_pressure>0
if self.getSolverOptionsPressure().isVerbose() or self.verbose: print "DarcyFlux: Pressure update"
if g == None:
self.__pde_p.setValue(X=-v, r=p0)
else:
self.__pde_p.setValue(X=g-v, r=p0)
p=self.__pde_p.getSolution()
return p

def __Aprod_v(self,dv):
# calculates: (a,b,c) = (K^{-1}(dv + KG * dp), L^{-1}Ddv, dp)  with (G*KG)dp = - G^*dv
dp=self.__getPressure(dv, p0=Data()) # dp = (G*KG)^{-1} (0-G^*dv)
b= - self.__applWeight(dv) # b = - (D K D^*)^{-1} (0-Dv)
return ArithmeticTuple(a,b,-dp)

def __Msolve_PCG_v(self,r):
# K^{-1} u = r[0] + D^*r[1] = K^{-1}(dv + KG * dp) + D^*L^{-1}Ddv
if self.getSolverOptionsFlux().isVerbose() or self.verbose: print "DarcyFlux: Applying preconditioner"
self.__pde_k.setValue(X=r[1]*util.kronecker(self.domain), Y=r[0], r=Data())
# self.__pde_p.getOperator().saveMM("prec.mm")
return self.__pde_k.getSolution()

def __inner_PCG_v(self,v,r):
return util.integrate(util.inner(v,r[0])+util.div(v)*r[1])

def __Aprod_p(self,dp):
if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
self.__pde_k.setValue(Y=-Gdp,X=Data(), r=Data())
du=self.__pde_k.getSolution()
# self.__pde_v.getOperator().saveMM("proj.mm")
return ArithmeticTuple(util.tensor_mult(self.__permeability,Gdp),-du)

def __getFlux(self,p, v0, f=None, g=None):
# solves (K^{-1}+D^*L^{-1} D) v = D^*L^{-1}f + K^{-1}g - Gp
if f!=None:
self.__pde_k.setValue(X=self.__applWeight(v0*0,self.__f)*util.kronecker(self.domain))
self.__pde_k.setValue(r=v0)
g2=util.tensor_mult(self.__permeability_inv,g)
if p == None:
self.__pde_k.setValue(Y=g2)
else:
return self.__pde_k.getSolution()

#v=self.__getFlux(p, u0, f=self.__f, g=g2)
def __Msolve_PCG_p(self,r):
if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"
self.__pde_p.setValue(X=r[0]-r[1], Y=Data(), r=Data(), y=Data())
# self.__pde_p.getOperator().saveMM("prec.mm")
return self.__pde_p.getSolution()

def __inner_PCG_p(self,p,r):

def __L2(self,v):
return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))

def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):
229        """        """
230        solves the problem.        solves the problem.
231
The iteration is terminated if the residual norm is less then self.getTolerance().

232        :param u0: initial guess for the flux. At locations in the domain marked by ``location_of_fixed_flux`` the value of ``u0`` is kept unchanged.        :param u0: initial guess for the flux. At locations in the domain marked by ``location_of_fixed_flux`` the value of ``u0`` is kept unchanged.
233        :type u0: vector value on the domain (e.g. `Data`).        :type u0: vector value on the domain (e.g. `escript.Data`).
234        :param p0: initial guess for the pressure. At locations in the domain marked by ``location_of_fixed_pressure`` the value of ``p0`` is kept unchanged.        :param p0: initial guess for the pressure. At locations in the domain marked by ``location_of_fixed_pressure`` the value of ``p0`` is kept unchanged.
235        :type p0: scalar value on the domain (e.g. `Data`).        :type p0: scalar value on the domain (e.g. `escript.Data`).
:param verbose: if set some information on iteration progress are printed
:type verbose: ``bool``
236        :return: flux and pressure        :return: flux and pressure
237        :rtype: ``tuple`` of `Data`.        :rtype: ``tuple`` of `escript.Data`.

:note: The problem is solved as a least squares form
*(K^[-1]+D^* (DKD^*)^[-1] D)u+G p=D^* (DKD^*)^[-1] f + K^[-1]g*
*G^*u+G^* K Gp=G^*g*

where *D* is the *div* operator and *(Gp)_i=p_{,i}* for the permeability *K=k_{ij}*.
"""
self.verbose=verbose
rtol=self.getTolerance()
atol=self.getAbsoluteTolerance()
self.setSubProblemTolerance()
num_corrections=0
converged=False
norm_r=None

# Eliminate the hydrostatic pressure:
if self.verbose: print "DarcyFlux: calculate hydrostatic pressure component."
self.__pde_p.setValue(X=self.__g, r=p0, y=-util.inner(self.domain.getNormal(),u0))
p0=self.__pde_p.getSolution()
norm_g2=util.integrate(util.inner(g2,util.tensor_mult(self.__permeability_inv,g2)))**0.5

p=p0*0
if self.solveForFlux:
v=u0.copy()
else:
v=self.__getFlux(p, u0, f=self.__f, g=g2)

while not converged and norm_g2 > 0:
KGp=util.tensor_mult(self.__permeability,Gp)
if self.verbose:
def_p=g2-(v+KGp)
def_v=self.__f-util.div(v)
print "DarcyFlux: L2: g-v-K*grad(p) = %e (v = %e)."%(self.__L2(def_p),self.__L2(v))
print "DarcyFlux: K^{-1}-norm of v = %e."%util.integrate(util.inner(v,util.tensor_mult(self.__permeability_inv,v)))**0.5
print "DarcyFlux: K^{-1}-norm of g2 = %e."%norm_g2
print "DarcyFlux: K-norm of grad(dp) = %e."%util.integrate(util.inner(Gp,KGp))**0.5
ATOL=atol+rtol*norm_g2
if self.verbose: print "DarcyFlux: absolute tolerance ATOL = %e."%(ATOL,)
if norm_r == None or norm_r>ATOL:
if num_corrections>max_num_corrections:
raise ValueError,"maximum number of correction steps reached."

if self.solveForFlux:
# initial residual is r=K^{-1}*(g-v-K*Gp)+D^*L^{-1}(f-Du)
v,r, norm_r=PCG(ArithmeticTuple(util.tensor_mult(self.__permeability_inv,g2-v)-Gp,self.__applWeight(v,self.__f),p),
self.__Aprod_v,
v,
self.__Msolve_PCG_v,
self.__inner_PCG_v,
atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
p=r[2]
else:
# initial residual is r=G^*(g2-KGp - v)
p,r, norm_r=PCG(ArithmeticTuple(g2-KGp,v),
self.__Aprod_p,
p,
self.__Msolve_PCG_p,
self.__inner_PCG_p,
atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
v=r[1]
if self.verbose: print "DarcyFlux: residual norm = %e."%norm_r
num_corrections+=1
else:
if self.verbose: print "DarcyFlux: stopping criterium reached."
converged=True
return v,p+p0
def setTolerance(self,rtol=1e-4):
"""
sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if

*|g-v-K gard(p)|_PCG <= atol + rtol * |K^{1/2}g2|_0*

where ``atol`` is an absolut tolerance (see `setAbsoluteTolerance`).

:param rtol: relative tolerance for the pressure
:type rtol: non-negative ``float``
"""
if rtol<0:
raise ValueError,"Relative tolerance needs to be non-negative."
self.__rtol=rtol
def getTolerance(self):
"""
returns the relative tolerance
:return: current relative tolerance
:rtype: ``float``
"""
return self.__rtol
238
def setAbsoluteTolerance(self,atol=0.):
239        """        """
240        sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if        self.__pde_p.setValue(X=self.__g * 1./self.perm_scale,
241                                Y=self.__f * 1./self.perm_scale,
242                                y= - util.inner(self.domain.getNormal(),u0 * self.location_of_fixed_flux * 1./self.perm_scale ),
243                                r=p0)
244          p=self.__pde_p.getSolution()
245          u = self.getFlux(p, u0)
246          return u,p
247
248        *|g-v-K gard(p)|_PCG <= atol + rtol * |K^{1/2}g2|_0*     def getFlux(self,p, u0=None):

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}*.

:param atol: absolute tolerance for the pressure
:type atol: non-negative ``float``
"""
if atol<0:
raise ValueError,"Absolute tolerance needs to be non-negative."
self.__atol=atol
def getAbsoluteTolerance(self):
"""
returns the absolute tolerance
:return: current absolute tolerance
:rtype: ``float``
"""
return self.__atol
def getSubProblemTolerance(self):
"""
Returns a suitable subtolerance
:type: ``float``
"""
return max(util.EPSILON**(0.5),self.getTolerance()**2)

def setSubProblemTolerance(self):
"""
Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected.
"""
sub_tol=self.getSubProblemTolerance()
self.getSolverOptionsFlux().setTolerance(sub_tol)
self.getSolverOptionsFlux().setAbsoluteTolerance(0.)
self.getSolverOptionsPressure().setTolerance(sub_tol)
self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
self.getSolverOptionsWeighting().setTolerance(sub_tol)
self.getSolverOptionsWeighting().setAbsoluteTolerance(0.)
if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol

class DarcyFlowOld(object):
"""
solves the problem

*u_i+k_{ij}*p_{,j} = g_i*
*u_{i,i} = f*

where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,

:note: The problem is solved in a least squares formulation.
"""

def __init__(self, domain, weight=None, useReduced=False, adaptSubTolerance=True):
"""
initializes the Darcy flux problem
:param domain: domain of the problem
:type domain: `Domain`
:param useReduced: uses reduced oreder on flux and pressure
:type useReduced: ``bool``
:param adaptSubTolerance: switches on automatic subtolerance selection
"""
self.domain=domain
if weight == None:
s=self.domain.getSize()
self.__l=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2
# self.__l=(3.*util.longestEdge(self.domain))**2
#self.__l=(0.1*util.longestEdge(self.domain)*s/util.sup(s))**2
else:
self.__l=weight
self.__pde_v=LinearPDESystem(domain)
if useReduced: self.__pde_v.setReducedOrderOn()
self.__pde_v.setSymmetryOn()
self.__pde_v.setValue(D=util.kronecker(domain), A=self.__l*util.outer(util.kronecker(domain),util.kronecker(domain)))
self.__pde_p=LinearSinglePDE(domain)
self.__pde_p.setSymmetryOn()
if useReduced: self.__pde_p.setReducedOrderOn()
self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
self.setTolerance()
self.setAbsoluteTolerance()
self.verbose=False
def getSolverOptionsFlux(self):
"""
Returns the solver options used to solve the flux problems

*(I+D^*D)u=F*

:return: `SolverOptions`
"""
return self.__pde_v.getSolverOptions()
def setSolverOptionsFlux(self, options=None):
"""
Sets the solver options used to solve the flux problems

*(I+D^*D)u=F*

If ``options`` is not present, the options are reset to default
:param options: `SolverOptions`
:note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
"""
return self.__pde_v.setSolverOptions(options)
def getSolverOptionsPressure(self):
"""
Returns the solver options used to solve the pressure problems

*(Q^*Q)p=Q^*G*

:return: `SolverOptions`
"""
return self.__pde_p.getSolverOptions()
def setSolverOptionsPressure(self, options=None):
"""
Sets the solver options used to solve the pressure problems

*(Q^*Q)p=Q^*G*

If ``options`` is not present, the options are reset to default
:param options: `SolverOptions`
:note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
"""
return self.__pde_p.setSolverOptions(options)

def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
"""
assigns values to model parameters

:param f: volumetic sources/sinks
:type f: scalar value on the domain (e.g. `Data`)
:param g: flux sources/sinks
:type g: vector values on the domain (e.g. `Data`)
:param location_of_fixed_pressure: mask for locations where pressure is fixed
:type location_of_fixed_pressure: scalar value on the domain (e.g. `Data`)
:param location_of_fixed_flux:  mask for locations where flux is fixed.
:type location_of_fixed_flux: vector values on the domain (e.g. `Data`)
:param permeability: permeability tensor. If scalar ``s`` is given the tensor with
``s`` on the main diagonal is used. If vector ``v`` is given the tensor with
``v`` on the main diagonal is used.
:type permeability: scalar, vector or tensor values on the domain (e.g. `Data`)

:note: the values of parameters which are not set by calling ``setValue`` are not altered.
:note: at any point on the boundary of the domain the pressure (``location_of_fixed_pressure`` >0)
or the normal component of the flux (``location_of_fixed_flux[i]>0`` if direction of the normal
is along the *x_i* axis.
"""
if f !=None:
f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
if f.isEmpty():
f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
else:
if f.getRank()>0: raise ValueError,"illegal rank of f."
self.__f=f
if g !=None:
g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
if g.isEmpty():
g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
else:
if not g.getShape()==(self.domain.getDim(),):
raise ValueError,"illegal shape of g"
self.__g=g

if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)
if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux)

if permeability!=None:
perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
if perm.getRank()==0:
perm=perm*util.kronecker(self.domain.getDim())
elif perm.getRank()==1:
perm, perm2=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A")), perm
for i in range(self.domain.getDim()): perm[i,i]=perm2[i]
elif perm.getRank()==2:
pass
else:
raise ValueError,"illegal rank of permeability."
self.__permeability=perm
self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))

def setTolerance(self,rtol=1e-4):
"""
sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if

*|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*

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}*.

:param rtol: relative tolerance for the pressure
:type rtol: non-negative ``float``
"""
if rtol<0:
raise ValueError,"Relative tolerance needs to be non-negative."
self.__rtol=rtol
def getTolerance(self):
"""
returns the relative tolerance

:return: current relative tolerance
:rtype: ``float``
"""
return self.__rtol

def setAbsoluteTolerance(self,atol=0.):
249          """          """
250          sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if          returns the flux for a given pressure ``p`` where the flux is equal to ``u0``

*|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*

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}*.

:param atol: absolute tolerance for the pressure
:type atol: non-negative ``float``
"""
if atol<0:
raise ValueError,"Absolute tolerance needs to be non-negative."
self.__atol=atol
def getAbsoluteTolerance(self):
"""
returns the absolute tolerance

:return: current absolute tolerance
:rtype: ``float``
"""
return self.__atol
def getSubProblemTolerance(self):
"""
Returns a suitable subtolerance
@type: ``float``
"""
return max(util.EPSILON**(0.75),self.getTolerance()**2)
def setSubProblemTolerance(self):
"""
Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected.
"""
sub_tol=self.getSubProblemTolerance()
self.getSolverOptionsFlux().setTolerance(sub_tol)
self.getSolverOptionsFlux().setAbsoluteTolerance(0.)
self.getSolverOptionsPressure().setTolerance(sub_tol)
self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol

def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):
"""
solves the problem.

The iteration is terminated if the residual norm is less then self.getTolerance().

: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.
:type u0: vector value on the domain (e.g. `Data`).
: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.
:type p0: scalar value on the domain (e.g. `Data`).
:param verbose: if set some information on iteration progress are printed
:type verbose: ``bool``
:return: flux and pressure
:rtype: ``tuple`` of `Data`.

:note: The problem is solved as a least squares form

*(I+D^*D)u+Qp=D^*f+g*
*Q^*u+Q^*Qp=Q^*g*

where *D* is the *div* operator and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.
We eliminate the flux form the problem by setting

*u=(I+D^*D)^{-1}(D^*f-g-Qp)* with u=u0 on location_of_fixed_flux

form the first equation. Inserted into the second equation we get

*Q^*(I-(I+D^*D)^{-1})Qp= Q^*(g-(I+D^*D)^{-1}(D^*f+g))* with p=p0  on location_of_fixed_pressure

which is solved using the PCG method (precondition is *Q^*Q*). In each iteration step
PDEs with operator *I+D^*D* and with *Q^*Q* needs to be solved using a sub iteration scheme.
"""
self.verbose=verbose
rtol=self.getTolerance()
atol=self.getAbsoluteTolerance()
self.setSubProblemTolerance()
num_corrections=0
converged=False
p=p0
norm_r=None
while not converged:
v=self.getFlux(p, fixed_flux=u0)
Qp=self.__Q(p)
norm_v=self.__L2(v)
norm_Qp=self.__L2(Qp)
if norm_v == 0.:
if norm_Qp == 0.:
return v,p
else:
fac=norm_Qp
else:
if norm_Qp == 0.:
fac=norm_v
else:
fac=2./(1./norm_v+1./norm_Qp)
ATOL=(atol+rtol*fac)
if self.verbose:
print "DarcyFlux: L2 norm of v = %e."%norm_v
print "DarcyFlux: L2 norm of k.util.grad(p) = %e."%norm_Qp
print "DarcyFlux: L2 defect u = %e."%(util.integrate(util.length(self.__g-util.interpolate(v,Function(self.domain))-Qp)**2)**(0.5),)
print "DarcyFlux: L2 defect div(v) = %e."%(util.integrate((self.__f-util.div(v))**2)**(0.5),)
print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
if norm_r == None or norm_r>ATOL:
if num_corrections>max_num_corrections:
raise ValueError,"maximum number of correction steps reached."
p,r, norm_r=PCG(self.__g-util.interpolate(v,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)
num_corrections+=1
else:
converged=True
return v,p
def __L2(self,v):
return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))

def __Q(self,p):

def __Aprod(self,dp):
if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
Qdp=self.__Q(dp)
self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())
du=self.__pde_v.getSolution()
# self.__pde_v.getOperator().saveMM("proj.mm")
return Qdp+du
def __inner_GMRES(self,r,s):
return util.integrate(util.inner(r,s))

def __inner_PCG(self,p,r):
return util.integrate(util.inner(self.__Q(p), r))

def __Msolve_PCG(self,r):
if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"
self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=Data(), r=Data())
# self.__pde_p.getOperator().saveMM("prec.mm")
return self.__pde_p.getSolution()

def getFlux(self,p=None, fixed_flux=Data()):
"""
returns the flux for a given pressure ``p`` where the flux is equal to ``fixed_flux``
251          on locations where ``location_of_fixed_flux`` is positive (see `setValue`).          on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
252          Note that ``g`` and ``f`` are used, see `setValue`.          Notice that ``g`` is used, see `setValue`.
253
254          :param p: pressure.          :param p: pressure.
255          :type p: scalar value on the domain (e.g. `Data`).          :type p: scalar value on the domain (e.g. `escript.Data`).
256          :param fixed_flux: flux on the locations of the domain marked be ``location_of_fixed_flux``.          :param u0: flux on the locations of the domain marked be ``location_of_fixed_flux``.
257          :type fixed_flux: vector values on the domain (e.g. `Data`).          :type u0: vector values on the domain (e.g. `escript.Data`) or ``None``
258          :return: flux          :return: flux
259          :rtype: `Data`          :rtype: `escript.Data`
: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}*
for the permeability *k_{ij}*
260          """          """
261      self.setSubProblemTolerance()          if self.solver  == self.EVAL:
262          g=self.__g             u = self.__g-self.perm_scale * util.tensor_mult(self.__permeability,util.grad(p))
263          f=self.__f          elif self.solver  == self.POST or self.solver  == self.SMOOTH:
264          self.__pde_v.setValue(X=self.__l*f*util.kronecker(self.domain), r=fixed_flux)              self.__pde_v.setValue(Y=util.tensor_mult(self.__permeability_inv,self.__g * 1./self.perm_scale)-util.grad(p))
265          if p == None:              if u0 == None:
266             self.__pde_v.setValue(Y=g)             self.__pde_v.setValue(r=escript.Data())
267          else:          else:
268             self.__pde_v.setValue(Y=g-self.__Q(p))                 if not isinstance(u0, escript.Data) : u0 = escript.Vector(u0, escript.Solution(self.domain))
269          return self.__pde_v.getSolution()             self.__pde_v.setValue(r=1./self.perm_scale * u0)
270                u= self.__pde_v.getSolution() * self.perm_scale
271        return u
272
274       """       """
275       solves       solves
# Line 778  class StokesProblemCartesian(Homogeneous Line 288  class StokesProblemCartesian(Homogeneous
288              sp.setTolerance()              sp.setTolerance()
289              sp.initialize(...)              sp.initialize(...)
290              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
291                sp.setStokesEquation(...) # new values for some parameters
292                v1,p1=sp.solve(v,p)
293       """       """
294       def __init__(self,domain,**kwargs):       def __init__(self,domain,**kwargs):
295           """           """
# Line 792  class StokesProblemCartesian(Homogeneous Line 304  class StokesProblemCartesian(Homogeneous
304           """           """
306           self.domain=domain           self.domain=domain
307           self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())           self.__pde_v=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
308           self.__pde_u.setSymmetryOn()           self.__pde_v.setSymmetryOn()
309
310           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
311           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
# Line 810  class StokesProblemCartesian(Homogeneous Line 322  class StokesProblemCartesian(Homogeneous
322
323       :rtype: `SolverOptions`       :rtype: `SolverOptions`
324       """       """
325       return self.__pde_u.getSolverOptions()       return self.__pde_v.getSolverOptions()
326       def setSolverOptionsVelocity(self, options=None):       def setSolverOptionsVelocity(self, options=None):
327           """           """
328       set the solver options for solving the equation for velocity.       set the solver options for solving the equation for velocity.
# Line 818  class StokesProblemCartesian(Homogeneous Line 330  class StokesProblemCartesian(Homogeneous
330       :param options: new solver  options       :param options: new solver  options
331       :type options: `SolverOptions`       :type options: `SolverOptions`
332       """       """
333           self.__pde_u.setSolverOptions(options)           self.__pde_v.setSolverOptions(options)
334       def getSolverOptionsPressure(self):       def getSolverOptionsPressure(self):
335           """           """
336       returns the solver options used  solve the equation for pressure.       returns the solver options used  solve the equation for pressure.
# Line 854  class StokesProblemCartesian(Homogeneous Line 366  class StokesProblemCartesian(Homogeneous
367           """           """
368           updates the Stokes equation to consider dependencies from ``v`` and ``p``           updates the Stokes equation to consider dependencies from ``v`` and ``p``
369           :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values.           :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values to model parameters.
370           """           """
371           pass           pass
# Line 875  class StokesProblemCartesian(Homogeneous Line 387  class StokesProblemCartesian(Homogeneous
387          if eta !=None:          if eta !=None:
388              k=util.kronecker(self.domain.getDim())              k=util.kronecker(self.domain.getDim())
389              kk=util.outer(k,k)              kk=util.outer(k,k)
390              self.eta=util.interpolate(eta, Function(self.domain))              self.eta=util.interpolate(eta, escript.Function(self.domain))
391          self.__pde_prec.setValue(D=1/self.eta)          self.__pde_prec.setValue(D=1/self.eta)
392              self.__pde_u.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))              self.__pde_v.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))
393          if restoration_factor!=None:          if restoration_factor!=None:
394              n=self.domain.getNormal()              n=self.domain.getNormal()
395              self.__pde_u.setValue(d=restoration_factor*util.outer(n,n))              self.__pde_v.setValue(d=restoration_factor*util.outer(n,n))
398          if f!=None: self.__f=f          if f!=None: self.__f=f
399          if surface_stress!=None: self.__surface_stress=surface_stress          if surface_stress!=None: self.__surface_stress=surface_stress
400          if stress!=None: self.__stress=stress          if stress!=None: self.__stress=stress
401
403          """          """
404          assigns values to the model parameters          assigns values to the model parameters
405
# Line 927  class StokesProblemCartesian(Homogeneous Line 439  class StokesProblemCartesian(Homogeneous
439           :return: inner product of element p and Bv=-div(v)           :return: inner product of element p and Bv=-div(v)
440           :rtype: ``float``           :rtype: ``float``
441           """           """
442           return util.integrate(util.interpolate(p,Function(self.domain))*util.interpolate(Bv,Function(self.domain)))           return util.integrate(util.interpolate(p,escript.Function(self.domain))*util.interpolate(Bv, escript.Function(self.domain)))
443
444       def inner_p(self,p0,p1):       def inner_p(self,p0,p1):
445           """           """
# Line 938  class StokesProblemCartesian(Homogeneous Line 450  class StokesProblemCartesian(Homogeneous
450           :return: inner product of p0 and p1           :return: inner product of p0 and p1
451           :rtype: ``float``           :rtype: ``float``
452           """           """
453           s0=util.interpolate(p0,Function(self.domain))           s0=util.interpolate(p0, escript.Function(self.domain))
454           s1=util.interpolate(p1,Function(self.domain))           s1=util.interpolate(p1, escript.Function(self.domain))
455           return util.integrate(s0*s1)           return util.integrate(s0*s1)
456
457       def norm_v(self,v):       def norm_v(self,v):
# Line 955  class StokesProblemCartesian(Homogeneous Line 467  class StokesProblemCartesian(Homogeneous
467
468       def getDV(self, p, v, tol):       def getDV(self, p, v, tol):
469           """           """
470           return the value for v for a given p (overwrite)           return the value for v for a given p
471
472           :param p: a pressure           :param p: a pressure
473           :param v: a initial guess for the value v to return.           :param v: a initial guess for the value v to return.
475           """           """
477           self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress)           self.__pde_v.setValue(Y=self.__f, y=self.__surface_stress)
478       self.getSolverOptionsVelocity().setTolerance(tol)       self.getSolverOptionsVelocity().setTolerance(tol)
479       self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)       self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
480           if self.__stress.isEmpty():           if self.__stress.isEmpty():
482           else:           else:
484           out=self.__pde_u.getSolution()           out=self.__pde_v.getSolution()
485           return  out           return  out
486
487       def norm_Bv(self,Bv):       def norm_Bv(self,Bv):
# Line 979  class StokesProblemCartesian(Homogeneous Line 491  class StokesProblemCartesian(Homogeneous
491          :rtype: equal to the type of p          :rtype: equal to the type of p
492          :note: boundary conditions on p should be zero!          :note: boundary conditions on p should be zero!
493          """          """
494          return util.sqrt(util.integrate(util.interpolate(Bv,Function(self.domain))**2))          return util.sqrt(util.integrate(util.interpolate(Bv, escript.Function(self.domain))**2))
495
496       def solve_AinvBt(self,p, tol):       def solve_AinvBt(self,p, tol):
497           """           """
# Line 989  class StokesProblemCartesian(Homogeneous Line 501  class StokesProblemCartesian(Homogeneous
501           :return: the solution of *Av=B^*p*           :return: the solution of *Av=B^*p*
502           :note: boundary conditions on v should be zero!           :note: boundary conditions on v should be zero!
503           """           """
504           self.__pde_u.setValue(Y=Data(), y=Data(), X=-p*util.kronecker(self.domain))           self.__pde_v.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))
505           out=self.__pde_u.getSolution()           out=self.__pde_v.getSolution()
506           return  out           return  out
507
508       def solve_prec(self,Bv, tol):       def solve_prec(self,Bv, tol):

Legend:
 Removed from v.3074 changed lines Added in v.3771