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

revision 3452 by caltinay, Tue Jan 25 01:53:57 2011 UTC revision 3569 by gross, Thu Sep 1 02:42:36 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
print dir(escript)

40  class DarcyFlow(object):  class DarcyFlow(object):
41     """     """
42     solves the problem     solves the problem
# Line 48  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, useVPIteration=True, weighting_scale=0.1):     SIMPLE="EVAL"
56       POST="POST"
57       SMOOTH="SMOOTH"
58       def __init__(self, domain, useReduced=False, solver="EVAL", 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 useVPIteration: if True altenative iteration over v and p is performed. Otherwise V and P are calculated in a single PDE.        :param w: weighting factor for `DarcyFlow.POST` solver
70        :type useVPIteration: ``bool``            :type w: ``float``
"""
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"))

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

*K^{-1} u=F*

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

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

*K^{-1}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^* K 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^* K Q)p=-Q^*G*

If ``options`` is not present, the options are reset to default
71
:param options: `SolverOptions`
:note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
72        """        """
73        return self.__pde_p.setSolverOptions(options)        if not solver in [DarcyFlow.EVAL, DarcyFlow.POST,  DarcyFlow.SMOOTH ] :
74              raise ValueError,"unknown solver %d."%solver
75
76          self.domain=domain
77          self.solver=solver
78          self.useReduced=useReduced
79          self.verbose=verbose
80          self.l=None
81          self.w=None
82
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.EVAL:
88             self.__pde_v=None
89         if self.verbose: print "DarcyFlow: simple solver is used."
90
91          elif self.solver  == self.POST:
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             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
113
114     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):
115        """        """
116        assigns values to model parameters        assigns values to model parameters
# Line 156  class DarcyFlow(object): Line 124  class DarcyFlow(object):
124        :param location_of_fixed_flux:  mask for locations where flux is fixed.        :param location_of_fixed_flux:  mask for locations where flux is fixed.
125        :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`)
126        :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.
127        :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`)
128
129        :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.
130        :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 165  class DarcyFlow(object): Line 133  class DarcyFlow(object):
133               is along the *x_i* axis.               is along the *x_i* axis.
134
135        """        """
136        if self.useVPIteration:        if location_of_fixed_pressure!=None:
137           if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)             self.location_of_fixed_pressure=util.wherePositive(location_of_fixed_pressure)
138           if location_of_fixed_flux!=None: self.__pde_k.setValue(q=location_of_fixed_flux)             self.__pde_p.setValue(q=self.location_of_fixed_pressure)
139        else:        if location_of_fixed_flux!=None:
140           q=self.__pde_k.getCoefficient("q")            self.location_of_fixed_flux=util.wherePositive(location_of_fixed_flux)
141           if q.isEmpty(): q=self.__pde_k.createCoefficient("q")            if not self.__pde_v == None: self.__pde_v.setValue(q=self.location_of_fixed_flux)
if location_of_fixed_pressure!=None: q[self.domain.getDim()]=location_of_fixed_pressure
if location_of_fixed_flux!=None: q[:self.domain.getDim()]=location_of_fixed_flux
self.__pde_k.setValue(q=q)
142
# flux is rescaled by the factor mean value(perm_inv)*length where length**self.domain.getDim()=vol(self.domain)
143        if permeability!=None:        if permeability!=None:
144       perm=util.interpolate(permeability,self.__pde_k.getFunctionSpaceForCoefficient("A"))
145           V=util.vol(self.domain)       perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
146
147       if perm.getRank()==0:       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.)
148
149          perm_inv=perm_inv*((1./self.scale)*util.kronecker(self.domain.getDim()))          perm_inv=(1./perm)
150          perm=perm*(self.scale*util.kronecker(self.domain.getDim()))          perm_inv=perm_inv*util.kronecker(self.domain.getDim())
151            perm=perm*util.kronecker(self.domain.getDim())
152
153
154       elif perm.getRank()==2:       elif perm.getRank()==2:
155          perm_inv=util.inverse(perm)          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
156       else:       else:
157          raise ValueError,"illegal rank of permeability."          raise ValueError,"illegal rank of permeability."
158
159       self.__permeability=perm       self.__permeability=perm
160       self.__permeability_inv=perm_inv       self.__permeability_inv=perm_inv
161
162       self.__l2 =(util.longestEdge(self.domain)**2*util.length(self.__permeability_inv))*self.weighting_scale           #====================
163           if self.useVPIteration:       self.__pde_p.setValue(A=self.__permeability)
164          if  self.solveForFlux:           if self.solver  == self.EVAL:
165             self.__pde_k.setValue(D=self.__permeability_inv)                pass # no extra work required
166          else:           elif self.solver  == self.POST:
167             self.__pde_k.setValue(D=self.__permeability_inv, A=self.__l2*util.outer(util.kronecker(self.domain),util.kronecker(self.domain)))          k=util.kronecker(self.domain.getDim())
168          self.__pde_p.setValue(A=self.__permeability)          self.omega = self.w*util.length(perm_inv)*self.l*self.domain.getSize()
169           else:          self.__pde_v.setValue(D=self.__permeability_inv, A=self.omega*util.outer(k,k))
170              D=self.__pde_k.createCoefficient("D")           elif self.solver  == self.SMOOTH:
171              A=self.__pde_k.createCoefficient("A")          self.__pde_v.setValue(D=self.__permeability_inv)
172              D[:self.domain.getDim(),:self.domain.getDim()]=self.__permeability_inv
173              for i in range(self.domain.getDim()):        if g != None:
174                 for j in range(self.domain.getDim()):      g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
175                   A[i,i,j,j]=self.__l2      if g.isEmpty():
176              A[self.domain.getDim(),:,self.domain.getDim(),:]=self.__permeability            g=Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
177              self.__pde_k.setValue(A=A, D=D)      else:
178        if g !=None:          if not g.getShape()==(self.domain.getDim(),): raise ValueError,"illegal shape of g"
179       g=util.interpolate(g, self.__pde_k.getFunctionSpaceForCoefficient("Y"))      self.__g=g
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
180        if f !=None:        if f !=None:
181       f=util.interpolate(f, self.__pde_k.getFunctionSpaceForCoefficient("X"))       f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
182       if f.isEmpty():       if f.isEmpty():
183            f=Scalar(0,self.__pde_k.getFunctionSpaceForCoefficient("X"))            f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
184       else:       else:
185           if f.getRank()>0: raise ValueError,"illegal rank of f."           if f.getRank()>0: raise ValueError,"illegal rank of f."
186           self.__f=f       self.__f=f
187
188       def getSolverOptionsFlux(self):
189          """
190          Returns the solver options used to solve the flux problems
191          :return: `SolverOptions`
192          """
193          if self.__pde_v == None:
194              return None
195          else:
196              return self.__pde_v.getSolverOptions()
197
198     def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):     def setSolverOptionsFlux(self, options=None):
199        """        """
200        solves the problem.        Sets the solver options used to solve the flux problems
201          If ``options`` is not present, the options are reset to default
202          :param options: `SolverOptions`
203          """
204          if not self.__pde_v == None:
205              self.__pde_v.setSolverOptions(options)
206
207       def getSolverOptionsPressure(self):
208          """
209          Returns the solver options used to solve the pressure problems
210          :return: `SolverOptions`
211          """
212          return self.__pde_p.getSolverOptions()
213
214        The iteration is terminated if the residual norm is less then self.getTolerance().     def setSolverOptionsPressure(self, options=None):

: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}*.
215        """        """
216        self.verbose=verbose        Sets the solver options used to solve the pressure problems
217        if self.useVPIteration:        If ``options`` is not present, the options are reset to default
218            return self.__solveVP(u0,p0,max_iter,max_num_corrections)
219        else:        :param options: `SolverOptions`
220            X=self.__pde_k.createCoefficient("X")        :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
221            Y=self.__pde_k.createCoefficient("Y")        """
222            Y[:self.domain.getDim()]=self.scale*util.tensor_mult(self.__permeability_inv,self.__g)        return self.__pde_p.setSolverOptions(options)
223            rtmp=self.__f * self.__l2 * self.scale
224            for i in range(self.domain.getDim()): X[i,i]=rtmp     def solve(self, u0, p0):
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):
225        """        """
226        solves the problem.        solves the problem.
227
The iteration is terminated if the residual norm is less than self.getTolerance().

228        :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.
229        :type u0: vector value on the domain (e.g. `escript.Data`).        :type u0: vector value on the domain (e.g. `escript.Data`).
230        :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.
# Line 299  class DarcyFlow(object): Line 232  class DarcyFlow(object):
232        :return: flux and pressure        :return: flux and pressure
233        :rtype: ``tuple`` of `escript.Data`.        :rtype: ``tuple`` of `escript.Data`.
234
: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))

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

def setAbsoluteTolerance(self,atol=0.):
235        """        """
236        sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if        self.__pde_p.setValue(X=self.__g ,
237                                Y=self.__f,
238                                y= - util.inner(self.domain.getNormal(),u0 * self.location_of_fixed_flux),
239                                r=p0)
240          p=self.__pde_p.getSolution()
241          u = self.getFlux(p, u0)
242          return u,p
243
244        *|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.)
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
245          """          """
246          self.domain=domain          returns the flux for a given pressure ``p`` where the flux is equal to ``u0``
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*

: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. `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()):
"""
returns the flux for a given pressure ``p`` where the flux is equal to ``fixed_flux``
247          on locations where ``location_of_fixed_flux`` is positive (see `setValue`).          on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
248          Note that ``g`` and ``f`` are used, see `setValue`.          Notice that ``g`` is used, see `setValue`.
249
250          :param p: pressure.          :param p: pressure.
251          :type p: scalar value on the domain (e.g. `escript.Data`).          :type p: scalar value on the domain (e.g. `escript.Data`).
252          :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``.
253          :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``
254          :return: flux          :return: flux
255          :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}*
256          """          """
258          g=self.__g          if self.solver  == self.EVAL:
260          self.__pde_v.setValue(X=self.__l2*f*util.kronecker(self.domain), r=fixed_flux)          elif self.solver  == self.POST or self.solver  == self.SMOOTH:
261          if p == None:              self.__pde_v.setValue(Y=util.tensor_mult(self.__permeability_inv,self.__g)-util.grad(p))
262             self.__pde_v.setValue(Y=g)              if u0 == None:
263          else:             self.__pde_v.setValue(r=escript.Data())
264             self.__pde_v.setValue(Y=g-self.__Q(p))          else:
265          return self.__pde_v.getSolution()             self.__pde_v.setValue(r=u0)
266                u= self.__pde_v.getSolution()
267        return u
268
270       """       """
271       solves       solves
# Line 849  class StokesProblemCartesian(Homogeneous Line 298  class StokesProblemCartesian(Homogeneous
298           """           """
300           self.domain=domain           self.domain=domain
301           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())
302           self.__pde_u.setSymmetryOn()           self.__pde_v.setSymmetryOn()
303
304           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
305           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
# Line 867  class StokesProblemCartesian(Homogeneous Line 316  class StokesProblemCartesian(Homogeneous
316
317       :rtype: `SolverOptions`       :rtype: `SolverOptions`
318       """       """
319       return self.__pde_u.getSolverOptions()       return self.__pde_v.getSolverOptions()
320       def setSolverOptionsVelocity(self, options=None):       def setSolverOptionsVelocity(self, options=None):
321           """           """
322       set the solver options for solving the equation for velocity.       set the solver options for solving the equation for velocity.
# Line 875  class StokesProblemCartesian(Homogeneous Line 324  class StokesProblemCartesian(Homogeneous
324       :param options: new solver  options       :param options: new solver  options
325       :type options: `SolverOptions`       :type options: `SolverOptions`
326       """       """
327           self.__pde_u.setSolverOptions(options)           self.__pde_v.setSolverOptions(options)
328       def getSolverOptionsPressure(self):       def getSolverOptionsPressure(self):
329           """           """
330       returns the solver options used  solve the equation for pressure.       returns the solver options used  solve the equation for pressure.
# Line 934  class StokesProblemCartesian(Homogeneous Line 383  class StokesProblemCartesian(Homogeneous
383              kk=util.outer(k,k)              kk=util.outer(k,k)
384              self.eta=util.interpolate(eta, escript.Function(self.domain))              self.eta=util.interpolate(eta, escript.Function(self.domain))
385          self.__pde_prec.setValue(D=1/self.eta)          self.__pde_prec.setValue(D=1/self.eta)
386              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)))
387          if restoration_factor!=None:          if restoration_factor!=None:
388              n=self.domain.getNormal()              n=self.domain.getNormal()
389              self.__pde_u.setValue(d=restoration_factor*util.outer(n,n))              self.__pde_v.setValue(d=restoration_factor*util.outer(n,n))
392          if f!=None: self.__f=f          if f!=None: self.__f=f
393          if surface_stress!=None: self.__surface_stress=surface_stress          if surface_stress!=None: self.__surface_stress=surface_stress
394          if stress!=None: self.__stress=stress          if stress!=None: self.__stress=stress
# Line 1019  class StokesProblemCartesian(Homogeneous Line 468  class StokesProblemCartesian(Homogeneous
469           """           """
471           self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress)           self.__pde_v.setValue(Y=self.__f, y=self.__surface_stress)
472       self.getSolverOptionsVelocity().setTolerance(tol)       self.getSolverOptionsVelocity().setTolerance(tol)
473       self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)       self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
474           if self.__stress.isEmpty():           if self.__stress.isEmpty():
476           else:           else:
478           out=self.__pde_u.getSolution()           out=self.__pde_v.getSolution()
479           return  out           return  out
480
481       def norm_Bv(self,Bv):       def norm_Bv(self,Bv):
# Line 1046  class StokesProblemCartesian(Homogeneous Line 495  class StokesProblemCartesian(Homogeneous
495           :return: the solution of *Av=B^*p*           :return: the solution of *Av=B^*p*
496           :note: boundary conditions on v should be zero!           :note: boundary conditions on v should be zero!
497           """           """
498           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))
499           out=self.__pde_u.getSolution()           out=self.__pde_v.getSolution()
500           return  out           return  out
501
502       def solve_prec(self,Bv, tol):       def solve_prec(self,Bv, tol):

Legend:
 Removed from v.3452 changed lines Added in v.3569