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

revision 3501 by gross, Wed Apr 13 04:07:53 2011 UTC revision 3502 by gross, Thu Apr 28 05:06:24 2011 UTC
# Line 37  import util Line 37  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 DarcyFlowNew(object):
41     """     """
42     solves the problem     solves the problem
43
# Line 48  class DarcyFlowNew(object): Line 46  class DarcyFlowNew(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 stabelized formulation.     :cvar SIMPLE: simple solver
50       :cvar POST: solver using global postprocessing of flux
51       :cvar STAB: solver uses (non-symmetric) stabilization
52       :cvar SYMSTAB: solver uses symmetric stabilization
53     """     """
54     def __init__(self, domain, useReduced=False, useVPIteration=True, *args, **kargs):     SIMPLE="SIMPLE"
55       POST="POST"
56       STAB="STAB"
57       SYMSTAB="SYMSTAB"
58       def __init__(self, domain, useReduced=False, solver="SYMSTAB", 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        :param useVPIteration: if True altenative iteration over v and p is performed. Otherwise V and P are calculated in a single PDE.        :type solver: in [`DarcyFlow.SIMPLE`, `DarcyFlow.POST', `DarcyFlow.STAB`, `DarcyFlow.SYMSTAB` ]
67        :type useVPIteration: ``bool``            :param verbose: if ``True`` some information on the iteration progress are printed.
68          :type verbose: ``bool``
69          :param w: weighting factor for `DarcyFlow.POST` solver
70          :type w: ``float``
71
72        """        """
73        self.domain=domain        self.domain=domain
74        self.useVPIteration=useVPIteration        self.solver=solver
75        self.useReduced=useReduced        self.useReduced=useReduced
76        self.verbose=False        self.verbose=verbose

if self.useVPIteration:
# this does not work yet
self.__pde_k=LinearPDESystem(domain)
self.__pde_k.setSymmetryOn()
if self.useReduced: self.__pde_k.setReducedOrderOn()

self.__pde_p=LinearSinglePDE(domain)
self.__pde_p.setSymmetryOn()
if self.useReduced: self.__pde_p.setReducedOrderOn()
else:
self.__pde_k=LinearPDE(self.domain, numEquations=self.domain.getDim()+1)
self.__pde_k.setSymmetryOff()

if self.useReduced: self.__pde_k.setReducedOrderOn()
C=self.__pde_k.createCoefficient("C")
B=self.__pde_k.createCoefficient("B")
for i in range(self.domain.getDim()):
B[i,i,self.domain.getDim()]=-1
C[self.domain.getDim(),i,i]=1
C[i,self.domain.getDim(),i]=-0.5
B[self.domain.getDim(),i,i]=0.5
self.__pde_k.setValue(C=C, B=B)
self.__f=escript.Scalar(0,self.__pde_k.getFunctionSpaceForCoefficient("X"))
self.__g=escript.Vector(0,self.__pde_k.getFunctionSpaceForCoefficient("Y"))
self.location_of_fixed_pressure = escript.Scalar(0, self.__pde_k.getFunctionSpaceForCoefficient("q"))
self.location_of_fixed_flux = escript.Vector(0, self.__pde_k.getFunctionSpaceForCoefficient("q"))
77        self.scale=1.        self.scale=1.
self.setTolerance()
78

def __L2(self,v):
return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))
def __inner_GMRES(self,r,s):
return util.integrate(util.inner(r,s))

def __Aprod_GMRES(self,p):
du=self.__pde_k.getSolution()
return self.__pde_p.getSolution()

def getSolverOptionsFlux(self):
"""
Returns the solver options used to solve the flux problems

*K^{-1} u=F*

:return: `SolverOptions`
"""
return self.__pde_k.getSolverOptions()
79
80          self.__pde_v=LinearPDESystem(domain)
81          self.__pde_v.setSymmetryOn()
82          if self.useReduced: self.__pde_v.setReducedOrderOn()
83          self.__pde_p=LinearSinglePDE(domain)
84          self.__pde_p.setSymmetryOn()
85          if self.useReduced: self.__pde_p.setReducedOrderOn()
86
87          if self.solver  == self.SIMPLE:
88         if self.verbose: print "DarcyFlow: simple solver is used."
89             self.__pde_v.setValue(D=util.kronecker(self.domain.getDim()))
90          elif self.solver  == self.POST:
91         self.w=w
92         if util.inf(w)<0.:
93            raise ValueError,"Weighting factor must be non-negative."
94         if self.verbose: print "DarcyFlow: global postprocessing of flux is used."
95          elif self.solver  == self.STAB:
96          if self.verbose: print "DarcyFlow: (non-symmetric) stabilization is used."
97          elif  self.solver  == self.SYMSTAB:
98          if self.verbose: print "DarcyFlow: symmetric stabilization is used."
99          else:
100        raise ValueError,"unknown solver %s"%self.solver
101          self.__f=escript.Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))
102          self.__g=escript.Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
103          self.location_of_fixed_pressure = escript.Scalar(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
104          self.location_of_fixed_flux = escript.Vector(0, self.__pde_v.getFunctionSpaceForCoefficient("q"))
105          self.setTolerance()
106
107
108     def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):     def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
109        """        """
110        assigns values to model parameters        assigns values to model parameters
# Line 130  class DarcyFlowNew(object): Line 118  class DarcyFlowNew(object):
118        :param location_of_fixed_flux:  mask for locations where flux is fixed.        :param location_of_fixed_flux:  mask for locations where flux is fixed.
119        :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)        :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)
120        :param permeability: permeability tensor. If scalar ``s`` is given the tensor with ``s`` on the main diagonal is used.        :param permeability: permeability tensor. If scalar ``s`` is given the tensor with ``s`` on the main diagonal is used.
121        :type permeability: scalar or tensor values on the domain (e.g. `escript.Data`)        :type permeability: scalar or symmetric tensor values on the domain (e.g. `escript.Data`)
122
123        :note: the values of parameters which are not set by calling ``setValue`` are not altered.        :note: the values of parameters which are not set by calling ``setValue`` are not altered.
124        :note: at any point on the boundary of the domain the pressure        :note: at any point on the boundary of the domain the pressure
# Line 139  class DarcyFlowNew(object): Line 127  class DarcyFlowNew(object):
127               is along the *x_i* axis.               is along the *x_i* axis.
128
129        """        """
130        if location_of_fixed_pressure!=None: self.location_of_fixed_pressure=util.wherePositive(location_of_fixed_pressure)        if location_of_fixed_pressure!=None:
131        if location_of_fixed_flux!=None: self.location_of_fixed_flux=util.wherePositive(location_of_fixed_flux)             self.location_of_fixed_pressure=util.wherePositive(location_of_fixed_pressure)
132               self.__pde_p.setValue(q=self.location_of_fixed_pressure)
133          if location_of_fixed_flux!=None:
134              self.location_of_fixed_flux=util.wherePositive(location_of_fixed_flux)
135              self.__pde_v.setValue(q=self.location_of_fixed_flux)
136
if self.useVPIteration:
if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=self.location_of_fixed_pressure)
if location_of_fixed_flux!=None: self.__pde_k.setValue(q=self.location_of_fixed_flux)
else:
if location_of_fixed_pressure!=None or location_of_fixed_flux!=None:
q=self.__pde_k.createCoefficient("q")
q[self.domain.getDim()]=self.location_of_fixed_pressure
q[:self.domain.getDim()]=self.location_of_fixed_flux
self.__pde_k.setValue(q=q)
137
138        # flux is rescaled by the factor mean value(perm_inv)*length where length**self.domain.getDim()=vol(self.domain)        # pressure  is rescaled by the factor 1/self.scale
139        if permeability!=None:        if permeability!=None:
140       perm=util.interpolate(permeability,self.__pde_k.getFunctionSpaceForCoefficient("A"))
141         perm=util.interpolate(permeability,self.__pde_v.getFunctionSpaceForCoefficient("A"))
142           V=util.vol(self.domain)           V=util.vol(self.domain)
143             l=V**(1./self.domain.getDim())
144
145       if perm.getRank()==0:       if perm.getRank()==0:
146          perm_inv=(1./perm)          perm_inv=(1./perm)
147              self.scale=util.integrate(perm_inv)*V**(1./self.domain.getDim()-1.)              self.scale=util.integrate(perm_inv)/V*l
148          perm_inv=perm_inv*((1./self.scale)*util.kronecker(self.domain.getDim()))          perm_inv=perm_inv*((1./self.scale)*util.kronecker(self.domain.getDim()))
149          perm=perm*(self.scale*util.kronecker(self.domain.getDim()))          perm=perm*(self.scale*util.kronecker(self.domain.getDim()))
150
151
152       elif perm.getRank()==2:       elif perm.getRank()==2:
153          perm_inv=util.inverse(perm)          perm_inv=util.inverse(perm)
154              self.scale=util.sqrt(util.integrate(util.length(perm_inv)**2)*V**(2./self.domain.getDim()-1.)/self.domain.getDim())              self.scale=util.sqrt(util.integrate(util.length(perm_inv)**2)/V)*l
155          perm_inv*=(1./self.scale)          perm_inv*=(1./self.scale)
156          perm=perm*self.scale          perm=perm*self.scale
157       else:       else:
# Line 171  class DarcyFlowNew(object): Line 159  class DarcyFlowNew(object):
159
160       self.__permeability=perm       self.__permeability=perm
161       self.__permeability_inv=perm_inv       self.__permeability_inv=perm_inv
162       if self.useVPIteration:       if self.verbose: print "DarcyFlow: scaling factor for pressure is %e."%self.scale
163              self.__pde_k.setValue(D=0.5*self.__permeability_inv)
164         if self.solver  == self.SIMPLE:
165            self.__pde_p.setValue(A=self.__permeability)
166         elif self.solver  == self.POST:
167            self.__pde_p.setValue(A=self.__permeability)
168            k=util.kronecker(self.domain.getDim())
169            self.lamb = self.w*util.length(perm_inv)*l
170            self.__pde_v.setValue(D=self.__permeability_inv, A=self.lamb*self.domain.getSize()*util.outer(k,k))
171         elif self.solver  == self.STAB:
172          self.__pde_p.setValue(A=0.5*self.__permeability)          self.__pde_p.setValue(A=0.5*self.__permeability)
173           else:          self.__pde_v.setValue(D=0.5*self.__permeability_inv)
174              D=self.__pde_k.createCoefficient("D")       elif  self.solver  == self.SYMSTAB:
175              A=self.__pde_k.createCoefficient("A")          self.__pde_p.setValue(A=0.5*self.__permeability)
176              D[:self.domain.getDim(),:self.domain.getDim()]=0.5*self.__permeability_inv          self.__pde_v.setValue(D=0.5*self.__permeability_inv)
177              A[self.domain.getDim(),:,self.domain.getDim(),:]=0.5*self.__permeability
self.__pde_k.setValue(A=A, D=D)
178        if g != None:        if g != None:
179      g=util.interpolate(g, self.__pde_k.getFunctionSpaceForCoefficient("Y"))      g=util.interpolate(g, self.__pde_v.getFunctionSpaceForCoefficient("Y"))
180      if g.isEmpty():      if g.isEmpty():
181            g=Vector(0,self.__pde_k.getFunctionSpaceForCoefficient("Y"))            g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
182      else:      else:
183          if not g.getShape()==(self.domain.getDim(),): raise ValueError,"illegal shape of g"          if not g.getShape()==(self.domain.getDim(),): raise ValueError,"illegal shape of g"
184          self.__g=g      self.__g=g
185        if f !=None:        if f !=None:
186       f=util.interpolate(f, self.__pde_k.getFunctionSpaceForCoefficient("X"))       f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
187       if f.isEmpty():       if f.isEmpty():
188                    f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
f=Scalar(0,self.__pde_k.getFunctionSpaceForCoefficient("X"))
189       else:       else:
190           if f.getRank()>0: raise ValueError,"illegal rank of f."           if f.getRank()>0: raise ValueError,"illegal rank of f."
191           self.__f=f       self.__f=f

def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10, iter_restart=20):
"""
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. `escript.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. `escript.Data`).
:param verbose: if set some information on iteration progress are printed
:type verbose: ``bool``
:return: flux and pressure
:rtype: ``tuple`` of `escript.Data`.

"""
u0_b=u0*self.location_of_fixed_flux
p0_b=p0*self.location_of_fixed_pressure/self.scale
f=self.__f-util.div(u0_b)
self.verbose=verbose
if self.useVPIteration:
# get u:
rtol=self.getTolerance()
p_init=0*p0_b
X=p_init*util.kronecker(self.__pde_k.getDomain()))
du=self.__pde_k.getSolution()
self.__pde_p.setValue(Y=f-util.div(du), X=0.5*(g-du))
p=GMRES(self.__pde_p.getSolution(),
self.__Aprod_GMRES,
p_init,
self.__inner_GMRES,
atol=0,
rtol=rtol,
iter_max=max_iter,
iter_restart=iter_restart, verbose=self.verbose,P_R=None)
X=p*util.kronecker(self.__pde_k.getDomain()))
u=self.__pde_k.getSolution()
else:
X=self.__pde_k.createCoefficient("X")
Y=self.__pde_k.createCoefficient("Y")
Y[:self.domain.getDim()]=0.5*util.tensor_mult(self.__permeability_inv,g)
Y[self.domain.getDim()]=f
X[self.domain.getDim(),:]=g*0.5
self.__pde_k.setValue(X=X, Y=Y)
self.__pde_k.getSolverOptions().setVerbosity(self.verbose)
#self.__pde_k.getSolverOptions().setPreconditioner(self.__pde_k.getSolverOptions().AMG)
self.__pde_k.getSolverOptions().setSolverMethod(self.__pde_k.getSolverOptions().DIRECT)
U=self.__pde_k.getSolution()
u=U[:self.domain.getDim()]
p=U[self.domain.getDim()]
# self.__pde_k.getOperator().saveMM("k.mm")
u=u0_b+u
p=(p0_b+p)*self.scale
if self.verbose:
def_p=self.__g-(u+KGp)
def_v=self.__f-util.div(u, self.__pde_k.getFunctionSpaceForCoefficient("X"))
print "DarcyFlux: L2: g-v-K*grad(p) = %e (v = %e)."%(self.__L2(def_p),self.__L2(u))
return u,p
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

class DarcyFlow(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, useReduced=False, adaptSubTolerance=True, solveForFlux=False, useVPIteration=True, weighting_scale=1.):
"""
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
:param solveForFlux: if True the solver solves for the flux (do not use!)
:type solveForFlux: ``bool``
:param useVPIteration: if True altenative iteration over v and p is performed. Otherwise V and P are calculated in a single PDE.
:type useVPIteration: ``bool``
"""
self.domain=domain
self.useVPIteration=useVPIteration
self.useReduced=useReduced
self.weighting_scale=weighting_scale
if self.useVPIteration:
self.solveForFlux=solveForFlux
self.verbose=False

self.__pde_k=LinearPDESystem(domain)
self.__pde_k.setSymmetryOn()
if self.useReduced: self.__pde_k.setReducedOrderOn()

self.__pde_p=LinearSinglePDE(domain)
self.__pde_p.setSymmetryOn()
if self.useReduced: self.__pde_p.setReducedOrderOn()
self.setTolerance()
self.setAbsoluteTolerance()
else:
self.__pde_k=LinearPDE(self.domain, numEquations=self.domain.getDim()+1)
self.__pde_k.setSymmetryOn()
if self.useReduced: self.__pde_k.setReducedOrderOn()
C=self.__pde_k.createCoefficient("C")
B=self.__pde_k.createCoefficient("B")
for i in range(self.domain.getDim()):
C[i,self.domain.getDim(),i]=1
B[self.domain.getDim(),i,i]=1
self.__pde_k.setValue(C=C, B=B)
self.__f=escript.Scalar(0,self.__pde_k.getFunctionSpaceForCoefficient("X"))
self.__g=escript.Vector(0,self.__pde_k.getFunctionSpaceForCoefficient("Y"))

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()        return self.__pde_v.getSolverOptions()
198
199     def setSolverOptionsFlux(self, options=None):     def setSolverOptionsFlux(self, options=None):
200        """        """
201        Sets the solver options used to solve the flux problems        Sets the solver options used to solve the flux problems

*K^{-1}u=F*

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

203        :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.
204        """        """
205        return self.__pde_v.setSolverOptions(options)        return self.__pde_v.setSolverOptions(options)
206
207     def getSolverOptionsPressure(self):     def getSolverOptionsPressure(self):
208        """        """
209        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*

210        :return: `SolverOptions`        :return: `SolverOptions`
211        """        """
212        return self.__pde_p.getSolverOptions()        return self.__pde_p.getSolverOptions()
# Line 375  class DarcyFlow(object): Line 214  class DarcyFlow(object):
214     def setSolverOptionsPressure(self, options=None):     def setSolverOptionsPressure(self, options=None):
215        """        """
216        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*

217        If ``options`` is not present, the options are reset to default        If ``options`` is not present, the options are reset to default
218
219        :param options: `SolverOptions`        :param options: `SolverOptions`
220        :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.
221        """        """
222        return self.__pde_p.setSolverOptions(options)        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. `escript.Data`)
:param g: flux sources/sinks
:type g: vector values on the domain (e.g. `escript.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. `escript.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. `escript.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. `escript.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 self.useVPIteration:
if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=util.wherePositive(location_of_fixed_pressure))
if location_of_fixed_flux!=None: self.__pde_k.setValue(q=util.wherePositive(location_of_fixed_flux))
else:
q=self.__pde_k.getCoefficient("q")
if q.isEmpty(): q=self.__pde_k.createCoefficient("q")
if location_of_fixed_pressure!=None: q[self.domain.getDim()]=util.wherePositive(location_of_fixed_pressure)
if location_of_fixed_flux!=None: q[:self.domain.getDim()]=util.wherePositive(location_of_fixed_flux)
self.__pde_k.setValue(q=q)

# flux is rescaled by the factor mean value(perm_inv)*length where length**self.domain.getDim()=vol(self.domain)
V=util.vol(self.domain)
if permeability!=None:
perm=util.interpolate(permeability,self.__pde_k.getFunctionSpaceForCoefficient("A"))

if perm.getRank()==0:
perm_inv=(1./perm)
if self.useVPIteration:
self.scale=1.
else:
self.scale=util.integrate(perm_inv)*V**(1./self.domain.getDim()-1.)

perm_inv=perm_inv*((1./self.scale)*util.kronecker(self.domain.getDim()))
perm=perm*(self.scale*util.kronecker(self.domain.getDim()))
elif perm.getRank()==2:
perm_inv=util.inverse(perm)
if self.useVPIteration:
self.scale=1.
else:
self.scale=util.sqrt(util.integrate(util.length(perm_inv)**2)*V**(2./self.domain.getDim()-1.)/self.domain.getDim())
perm_inv*=(1./self.scale)
perm=perm*self.scale
else:
raise ValueError,"illegal rank of permeability."

self.__permeability=perm
self.__permeability_inv=perm_inv

self.__l2 = V**(1./self.domain.getDim())*util.length(self.__permeability_inv)*self.domain.getSize()*self.weighting_scale
if self.useVPIteration:
if  self.solveForFlux:
self.__pde_k.setValue(D=self.__permeability_inv)
else:
self.__pde_k.setValue(D=self.__permeability_inv, A=self.__l2*util.outer(util.kronecker(self.domain),util.kronecker(self.domain)))
self.__pde_p.setValue(A=self.__permeability)
else:
D=self.__pde_k.createCoefficient("D")
A=self.__pde_k.createCoefficient("A")
D[:self.domain.getDim(),:self.domain.getDim()]=self.__permeability_inv
for i in range(self.domain.getDim()):
for j in range(self.domain.getDim()):
A[i,i,j,j]=self.__l2
A[self.domain.getDim(),:,self.domain.getDim(),:]=self.__permeability
self.__pde_k.setValue(A=A, D=D)
if g !=None:
g=util.interpolate(g, self.__pde_k.getFunctionSpaceForCoefficient("Y"))
if g.isEmpty():
g=Vector(0,self.__pde_k.getFunctionSpaceForCoefficient("Y"))
else:
if not g.getShape()==(self.domain.getDim(),):
raise ValueError,"illegal shape of g"
self.__g=g
elif permeability!=None:
X
if f !=None:
f=util.interpolate(f, self.__pde_k.getFunctionSpaceForCoefficient("X"))
if f.isEmpty():
f=Scalar(0,self.__pde_k.getFunctionSpaceForCoefficient("X"))
else:
if f.getRank()>0: raise ValueError,"illegal rank of f."
self.__f=f

223
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. `escript.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. `escript.Data`).
:param verbose: if set some information on iteration progress are printed
:type verbose: ``bool``
:return: flux and pressure
:rtype: ``tuple`` of `escript.Data`.

:note: The problem is solved as a least squares form
*(K^[-1]+D^* l2 D)u+G p=D^* l2 * 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
if self.useVPIteration:
return self.__solveVP(u0,p0,max_iter,max_num_corrections)
else:
X=self.__pde_k.createCoefficient("X")
Y=self.__pde_k.createCoefficient("Y")
Y[:self.domain.getDim()]=self.scale*util.tensor_mult(self.__permeability_inv,self.__g)
rtmp=self.__f * self.__l2 * self.scale
for i in range(self.domain.getDim()): X[i,i]=rtmp
X[self.domain.getDim(),:]=self.__g*self.scale
r=self.__pde_k.createCoefficient("r")
r[:self.domain.getDim()]=u0*self.scale
r[self.domain.getDim()]=p0
self.__pde_k.setValue(X=X, Y=Y, r=r)
self.__pde_k.getSolverOptions().setVerbosity(self.verbose)
#self.__pde_k.getSolverOptions().setPreconditioner(self.__pde_k.getSolverOptions().AMG)
self.__pde_k.getSolverOptions().setSolverMethod(self.__pde_k.getSolverOptions().DIRECT)
U=self.__pde_k.getSolution()
# self.__pde_k.getOperator().saveMM("k.mm")
u=U[:self.domain.getDim()]*(1./self.scale)
p=U[self.domain.getDim()]
if self.verbose:
def_p=self.__g-(u+KGp)
def_v=self.__f-util.div(u, self.__pde_k.getFunctionSpaceForCoefficient("X"))
print "DarcyFlux: L2: g-v-K*grad(p) = %e (v = %e)."%(self.__L2(def_p),self.__L2(u))
return u,p

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

The iteration is terminated if the residual norm is less than 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. `escript.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. `escript.Data`).
:return: flux and pressure
: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}*.
"""
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 __applWeight(self, v, f=None):
# solves L p = f-Dv with p = 0
if self.verbose: print "DarcyFlux: Applying weighting operator"
if f == None:
return -util.div(v)*self.__l2
else:
return (f-util.div(v))*self.__l2
def __getPressure(self, v, p0, g=None):
# 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=escript.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=escript.Data())
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=escript.Data(), r=escript.Data())
du=self.__pde_k.getSolution()
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=escript.Data(), r=escript.Data(), y=escript.Data())
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,escript.Function(self.domain)))**2))

def __L2_r(self,v):
return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.ReducedFunction(self.domain)))**2))

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

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

where ``atol`` is an absolut tolerance (see `setAbsoluteTolerance`).
227
228        :param rtol: relative tolerance for the pressure        :param rtol: relative tolerance for the pressure
229        :type rtol: non-negative ``float``        :type rtol: non-negative ``float``
# Line 691  class DarcyFlow(object): Line 231  class DarcyFlow(object):
231        if rtol<0:        if rtol<0:
232       raise ValueError,"Relative tolerance needs to be non-negative."       raise ValueError,"Relative tolerance needs to be non-negative."
233        self.__rtol=rtol        self.__rtol=rtol
234
235     def getTolerance(self):     def getTolerance(self):
236        """        """
237        returns the relative tolerance        returns the relative tolerance
# Line 698  class DarcyFlow(object): Line 239  class DarcyFlow(object):
239        :rtype: ``float``        :rtype: ``float``
240        """        """
241        return self.__rtol        return self.__rtol
242
243     def setAbsoluteTolerance(self,atol=0.):     def solve(self,u0,p0, max_iter=100, iter_restart=20):
244        """        """
245        sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if        solves the problem.
246
247        *|g-v-K gard(p)|_PCG <= atol + rtol * |K^{1/2}g2|_0*        The iteration is terminated if the residual norm is less then self.getTolerance().

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}*.
248
249        :param atol: absolute tolerance for the pressure        :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.
250        :type atol: non-negative ``float``        :type u0: vector value on the domain (e.g. `escript.Data`).
251        """        :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.
252        if atol<0:        :type p0: scalar value on the domain (e.g. `escript.Data`).
253       raise ValueError,"Absolute tolerance needs to be non-negative."        :param max_iter: maximum number of (outer) iteration steps for the stabilization solvers,
254        self.__atol=atol        :type max_iter: ``int``
255     def getAbsoluteTolerance(self):        :param iter_restart: number of steps after which the iteration is restarted. The larger ``iter_restart`` the larger the required memory.
256        """                             A small value for ``iter_restart`` may require a large number of iteration steps or may even lead to a failure
257        returns the absolute tolerance                             of the iteration. ``iter_restart`` is relevant for the stabilization solvers only.
258        :return: current absolute tolerance        :type iter_restart: ``int``
259        :rtype: ``float``        :return: flux and pressure
260        """        :rtype: ``tuple`` of `escript.Data`.
return self.__atol
def getSubProblemTolerance(self):
"""
Returns a suitable subtolerance
:type: ``float``
"""
return max(util.EPSILON**(0.5),self.getTolerance()**2)
261
def setSubProblemTolerance(self):
"""
Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected.
262        """        """
263        if self.__adaptSubTolerance:        # rescale initial guess:
264       sub_tol=self.getSubProblemTolerance()        p0=p0/self.scale
265       self.getSolverOptionsFlux().setTolerance(sub_tol)        if self.solver  == self.SIMPLE or self.solver  == self.POST :
266       self.getSolverOptionsFlux().setAbsoluteTolerance(0.)          self.__pde_p.setValue(X=self.__g ,
267       self.getSolverOptionsPressure().setTolerance(sub_tol)                                Y=self.__f,
268       self.getSolverOptionsPressure().setAbsoluteTolerance(0.)                                y=-util.inner(self.domain.getNormal(),u0 * self.location_of_fixed_flux),
269       if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol                                r=p0)
270            p=self.__pde_p.getSolution()
271            u = self.getFlux(p, u0)
272  class DarcyFlowOld(object):        elif  self.solver  == self.STAB:
273      """      u,p = self.__solve_STAB(u0,p0, max_iter, iter_restart)
274      solves the problem        elif  self.solver  == self.SYMSTAB:
275        u,p = self.__solve_SYMSTAB(u0,p0, max_iter, iter_restart)
*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.__l2=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2
# self.__l2=(3.*util.longestEdge(self.domain))**2
#self.__l2=(0.1*util.longestEdge(self.domain)*s/util.sup(s))**2
else:
self.__l2=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.__l2*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=escript.Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
self.__g=escript.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*
276
277      :return: `SolverOptions`        if self.verbose:
279      return self.__pde_v.getSolverOptions()          def_p=self.__g-(u+KGp)
280      def setSolverOptionsFlux(self, options=None):          def_v=self.__f-util.div(u, self.__pde_v.getFunctionSpaceForCoefficient("X"))
281      """          print "DarcyFlux: |g-u-K*grad(p)|_2 = %e (|u|_2 = %e)."%(self.__L2(def_p),self.__L2(u))
282      Sets the solver options used to solve the flux problems          print "DarcyFlux: |f-div(u)|_2 = %e (|grad(u)|_2 = %e)."%(self.__L2(def_v),self.__L2(util.grad(u)))
283              #rescale result
284      *(I+D^*D)u=F*        p=p*self.scale
285              return u,p
286      If ``options`` is not present, the options are reset to default
287      :param options: `SolverOptions`     def getFlux(self,p, u0=None):
: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. `escript.Data`)
:param g: flux sources/sinks
:type g: vector values on the domain (e.g. `escript.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. `escript.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. `escript.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. `escript.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=escript.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=escript.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.):
"""
sets the absolute tolerance ``atol`` 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 ``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. `escript.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. `escript.Data`).
:param verbose: if set some information on iteration progress are printed
:type verbose: ``bool``
:return: flux and pressure
:rtype: ``tuple`` of `escript.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,escript.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,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)
num_corrections+=1
else:
converged=True
return v,p
def __L2(self,v):
return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.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=escript.Data(), r=escript.Data())
du=self.__pde_v.getSolution()
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=escript.Data(), r=escript.Data())
return self.__pde_p.getSolution()

def getFlux(self,p=None, fixed_flux=escript.Data()):
288          """          """
289          returns the flux for a given pressure ``p`` where the flux is equal to ``fixed_flux``          returns the flux for a given pressure ``p`` where the flux is equal to ``u0``
290          on locations where ``location_of_fixed_flux`` is positive (see `setValue`).          on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
291          Note that ``g`` and ``f`` are used, see `setValue`.          Notice that ``g`` and ``f`` are used, see `setValue`.
292
293          :param p: pressure.          :param p: pressure.
294          :type p: scalar value on the domain (e.g. `escript.Data`).          :type p: scalar value on the domain (e.g. `escript.Data`).
295          :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``.
296          :type fixed_flux: vector values on the domain (e.g. `escript.Data`).          :type u0: vector values on the domain (e.g. `escript.Data`) or ``None``
297          :return: flux          :return: flux
298          :rtype: `escript.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}*
299          """          """
300      self.setSubProblemTolerance()          if self.solver  == self.SIMPLE or self.solver  == self.POST  :
302          f=self.__f              self.__pde_v.setValue(Y=self.__g-KGp, X=escript.Data())
303          self.__pde_v.setValue(X=self.__l2*f*util.kronecker(self.domain), r=fixed_flux)              if u0 == None:
304          if p == None:             self.__pde_v.setValue(r=escript.Data())
305             self.__pde_v.setValue(Y=g)          else:
306          else:             self.__pde_v.setValue(r=u0)
307             self.__pde_v.setValue(Y=g-self.__Q(p))              u= self.__pde_v.getSolution()
308          return self.__pde_v.getSolution()      elif self.solver  == self.POST:
310                                      X=self.lamb * self.__f * util.kronecker(self.domain.getDim()))
311                if u0 == None:
312               self.__pde_v.setValue(r=escript.Data())
313            else:
314               self.__pde_v.setValue(r=u0)
315                u= self.__pde_v.getSolution()
316        elif self.solver  == self.STAB:
318             self.__pde_v.setValue(Y=0.5*(util.tensor_mult(self.__permeability_inv,self.__g)+gp),
319                                   X= p * util.kronecker(self.domain.getDim()),
320                                   y= - p * self.domain.getNormal())
321             if u0 == None:
322               self.__pde_v.setValue(r=escript.Data())
323             else:
324               self.__pde_v.setValue(r=u0)
325             u= self.__pde_v.getSolution()
326        elif  self.solver  == self.SYMSTAB:
328             self.__pde_v.setValue(Y=0.5*(util.tensor_mult(self.__permeability_inv,self.__g)-gp),
329                                   X= escript.Data() ,
330                                   y= escript.Data() )
331             if u0 == None:
332               self.__pde_v.setValue(r=escript.Data())
333             else:
334               self.__pde_v.setValue(r=u0)
335             u= self.__pde_v.getSolution()
336        return u
337
338
339       def __solve_STAB(self, u0, p0, max_iter, iter_restart):
340              # p0 is used as an initial guess
341          u=self.getFlux(p0, u0)
342              self.__pde_p.setValue( Y=self.__f-util.div(u),
343                                     X=0.5*(self.__g - u - util.tensor_mult(self.__permeability,util.grad(p0)) ),
344                                     y= escript.Data(),
345                                     r=escript.Data())
346          dp=self.__pde_p.getSolution()
347          p=GMRES(dp,
348                  self.__STAB_Aprod,
349              p0,
350              self.__inner,
351              atol=self.__norm(p0+dp)*self.getTolerance() ,
352              rtol=0.,
353              iter_max=max_iter,
354              iter_restart=iter_restart,
355              verbose=self.verbose,P_R=None)
356
357              u=self.getFlux(p, u0)
358              return u,p
359
360       def __solve_SYMSTAB(self, u0, p0, max_iter, iter_restart):
361              # p0 is used as an initial guess
362          u=self.getFlux(p0, u0)
363              self.__pde_p.setValue( Y= self.__f,
364                                     X=  0.5*(self.__g + u - util.tensor_mult(self.__permeability,util.grad(p0)) ),
365                                     y= -   util.inner(self.domain.getNormal(), u),
366                                     r=escript.Data())
367          dp=self.__pde_p.getSolution()
368          p=GMRES(dp,
369                  self.__SYMSTAB_Aprod,
370              p0,
371              self.__inner,
372              atol=self.__norm(p0+dp)*self.getTolerance() ,
373              rtol=0.,
374              iter_max=max_iter,
375              iter_restart=iter_restart,
376              verbose=self.verbose,P_R=None)
377
378              u=self.getFlux(p, u0)
379              return u,p
380
381       def __L2(self,v):
382             return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))
383
384       def __norm(self,r):
385             return util.sqrt(self.__inner(r,r))
386
387       def __inner(self,r,s):
388             return util.integrate(util.inner(r,s), escript.Function(self.domain))
389
390       def __STAB_Aprod(self,p):
392          self.__pde_v.setValue(Y=-0.5*gp,
393                                X=-p*util.kronecker(self.__pde_v.getDomain()),
394                                y= p * self.domain.getNormal(),
395                                r=escript.Data())
396          u = -self.__pde_v.getSolution()
397          self.__pde_p.setValue(Y=util.div(u),
398                                X=0.5*(u+util.tensor_mult(self.__permeability,gp)),
399                                y=escript.Data(),
400                                r=escript.Data())
401
402          return  self.__pde_p.getSolution()
403
404       def __SYMSTAB_Aprod(self,p):
406          self.__pde_v.setValue(Y=0.5*gp ,
407                                X=escript.Data(),
408                                y=escript.Data(),
409                                r=escript.Data())
410          u = -self.__pde_v.getSolution()
411          self.__pde_p.setValue(Y=escript.Data(),
412                                X=0.5*(-u+util.tensor_mult(self.__permeability,gp)),
413                                y=   util.inner(self.domain.getNormal(), u),
414                                r=escript.Data())
415
416          return  self.__pde_p.getSolution()
417
418
420       """       """
# Line 1093  class StokesProblemCartesian(Homogeneous Line 448  class StokesProblemCartesian(Homogeneous
448           """           """
450           self.domain=domain           self.domain=domain
451           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())
452           self.__pde_u.setSymmetryOn()           self.__pde_v.setSymmetryOn()
453
454           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
455           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
# Line 1111  class StokesProblemCartesian(Homogeneous Line 466  class StokesProblemCartesian(Homogeneous
466
467       :rtype: `SolverOptions`       :rtype: `SolverOptions`
468       """       """
469       return self.__pde_u.getSolverOptions()       return self.__pde_v.getSolverOptions()
470       def setSolverOptionsVelocity(self, options=None):       def setSolverOptionsVelocity(self, options=None):
471           """           """
472       set the solver options for solving the equation for velocity.       set the solver options for solving the equation for velocity.
# Line 1119  class StokesProblemCartesian(Homogeneous Line 474  class StokesProblemCartesian(Homogeneous
474       :param options: new solver  options       :param options: new solver  options
475       :type options: `SolverOptions`       :type options: `SolverOptions`
476       """       """
477           self.__pde_u.setSolverOptions(options)           self.__pde_v.setSolverOptions(options)
478       def getSolverOptionsPressure(self):       def getSolverOptionsPressure(self):
479           """           """
480       returns the solver options used  solve the equation for pressure.       returns the solver options used  solve the equation for pressure.
# Line 1178  class StokesProblemCartesian(Homogeneous Line 533  class StokesProblemCartesian(Homogeneous
533              kk=util.outer(k,k)              kk=util.outer(k,k)
534              self.eta=util.interpolate(eta, escript.Function(self.domain))              self.eta=util.interpolate(eta, escript.Function(self.domain))
535          self.__pde_prec.setValue(D=1/self.eta)          self.__pde_prec.setValue(D=1/self.eta)
536              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)))
537          if restoration_factor!=None:          if restoration_factor!=None:
538              n=self.domain.getNormal()              n=self.domain.getNormal()
539              self.__pde_u.setValue(d=restoration_factor*util.outer(n,n))              self.__pde_v.setValue(d=restoration_factor*util.outer(n,n))
542          if f!=None: self.__f=f          if f!=None: self.__f=f
543          if surface_stress!=None: self.__surface_stress=surface_stress          if surface_stress!=None: self.__surface_stress=surface_stress
544          if stress!=None: self.__stress=stress          if stress!=None: self.__stress=stress
# Line 1263  class StokesProblemCartesian(Homogeneous Line 618  class StokesProblemCartesian(Homogeneous
619           """           """
621           self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress)           self.__pde_v.setValue(Y=self.__f, y=self.__surface_stress)
622       self.getSolverOptionsVelocity().setTolerance(tol)       self.getSolverOptionsVelocity().setTolerance(tol)
623       self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)       self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
624           if self.__stress.isEmpty():           if self.__stress.isEmpty():
626           else:           else:
628           out=self.__pde_u.getSolution()           out=self.__pde_v.getSolution()
629           return  out           return  out
630
631       def norm_Bv(self,Bv):       def norm_Bv(self,Bv):
# Line 1290  class StokesProblemCartesian(Homogeneous Line 645  class StokesProblemCartesian(Homogeneous
645           :return: the solution of *Av=B^*p*           :return: the solution of *Av=B^*p*
646           :note: boundary conditions on v should be zero!           :note: boundary conditions on v should be zero!
647           """           """
648           self.__pde_u.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))           self.__pde_v.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))
649           out=self.__pde_u.getSolution()           out=self.__pde_v.getSolution()
650           return  out           return  out
651
652       def solve_prec(self,Bv, tol):       def solve_prec(self,Bv, tol):

Legend:
 Removed from v.3501 changed lines Added in v.3502