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

revision 2881 by jfenwick, Thu Jan 28 02:03:15 2010 UTC revision 3515 by gross, Thu May 19 08:20:57 2011 UTC
# Line 1  Line 1
1    # -*- coding: utf-8 -*-
2  ########################################################  ########################################################
3  #  #
4  # Copyright (c) 2003-2010 by University of Queensland  # Copyright (c) 2003-2010 by University of Queensland
# Line 31  Some models for flow Line 32  Some models for flow
32
33  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
34
35  from escript import *  import escript
36  import util  import util
37  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions
38  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
39
40  class DarcyFlow(object):  class DarcyFlow(object):
41      """     """
42      solves the problem     solves the problem
43
44      *u_i+k_{ij}*p_{,j} = g_i*     *u_i+k_{ij}*p_{,j} = g_i*
45      *u_{i,i} = f*     *u_{i,i} = f*
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 SIMPLE: simple solver
50      """     :cvar POST: solver using global postprocessing of flux
51       :cvar STAB: solver uses (non-symmetric) stabilization
52      def __init__(self, domain, weight=None, useReduced=False, adaptSubTolerance=True):     :cvar SYMSTAB: solver uses symmetric stabilization
53          """     """
54          initializes the Darcy flux problem     SIMPLE="SIMPLE"
55          :param domain: domain of the problem     POST="POST"
56          :type domain: `Domain`     STAB="STAB"
57      :param useReduced: uses reduced oreder on flux and pressure     SYMSTAB="SYMSTAB"
58      :type useReduced: ``bool``     def __init__(self, domain, useReduced=False, solver="SYMSTAB", verbose=False, w=1.):
59      :param adaptSubTolerance: switches on automatic subtolerance selection        """
60      :type adaptSubTolerance: ``bool``          initializes the Darcy flux problem
61          """        :param domain: domain of the problem
62          self.domain=domain        :type domain: `Domain`
63          if weight == None:        :param useReduced: uses reduced oreder on flux and pressure
64             s=self.domain.getSize()        :type useReduced: ``bool``
65             self.__l=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2        :param solver: solver method
66             # self.__l=(3.*util.longestEdge(self.domain))**2        :type solver: in [`DarcyFlow.SIMPLE`, `DarcyFlow.POST', `DarcyFlow.STAB`, `DarcyFlow.SYMSTAB` ]
67             #self.__l=(0.1*util.longestEdge(self.domain)*s/util.sup(s))**2        :param verbose: if ``True`` some information on the iteration progress are printed.
68          else:        :type verbose: ``bool``
69             self.__l=weight        :param w: weighting factor for `DarcyFlow.POST` solver
70          self.__pde_v=LinearPDESystem(domain)        :type w: ``float``
71          if useReduced: self.__pde_v.setReducedOrderOn()
72          self.__pde_v.setSymmetryOn()        """
73          self.__pde_v.setValue(D=util.kronecker(domain), A=self.__l*util.outer(util.kronecker(domain),util.kronecker(domain)))        self.domain=domain
74          self.__pde_p=LinearSinglePDE(domain)        self.solver=solver
75          self.__pde_p.setSymmetryOn()        self.useReduced=useReduced
76          if useReduced: self.__pde_p.setReducedOrderOn()        self.verbose=verbose
77          self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))        self.scale=1.
78          self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
79          self.setTolerance()
80          self.setAbsoluteTolerance()        self.__pde_v=LinearPDESystem(domain)
82      self.verbose=False        if self.useReduced: self.__pde_v.setReducedOrderOn()
83      def getSolverOptionsFlux(self):        self.__pde_p=LinearSinglePDE(domain)
84      """        self.__pde_p.setSymmetryOn()
85      Returns the solver options used to solve the flux problems        if self.useReduced: self.__pde_p.setReducedOrderOn()
86
87      *(I+D^*D)u=F*        if self.solver  == self.SIMPLE:
88             if self.verbose: print "DarcyFlow: simple solver is used."
89      :return: `SolverOptions`           self.__pde_v.setValue(D=util.kronecker(self.domain.getDim()))
90      """        elif self.solver  == self.POST:
91      return self.__pde_v.getSolverOptions()       self.w=w
92      def setSolverOptionsFlux(self, options=None):       if util.inf(w)<0.:
93      """          raise ValueError,"Weighting factor must be non-negative."
94      Sets the solver options used to solve the flux problems       if self.verbose: print "DarcyFlow: global postprocessing of flux is used."
95              elif self.solver  == self.STAB:
96      *(I+D^*D)u=F*        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):
109          """
110          assigns values to model parameters
111
112          :param f: volumetic sources/sinks
113          :type f: scalar value on the domain (e.g. `escript.Data`)
114          :param g: flux sources/sinks
115          :type g: vector values on the domain (e.g. `escript.Data`)
116          :param location_of_fixed_pressure: mask for locations where pressure is fixed
117          :type location_of_fixed_pressure: scalar value on the domain (e.g. `escript.Data`)
118          :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`)
120          :param permeability: permeability tensor. If scalar ``s`` is given the tensor with ``s`` on the main diagonal is used.
121          :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.
124          :note: at any point on the boundary of the domain the pressure
125                 (``location_of_fixed_pressure`` >0) or the normal component of the
126                 flux (``location_of_fixed_flux[i]>0``) if direction of the normal
127                 is along the *x_i* axis.
128
129          """
130          if location_of_fixed_pressure!=None:
131               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
137
138          # pressure  is rescaled by the factor 1/self.scale
139          if permeability!=None:
140
141      If ``options`` is not present, the options are reset to default       perm=util.interpolate(permeability,self.__pde_v.getFunctionSpaceForCoefficient("A"))
142      :param options: `SolverOptions`           V=util.vol(self.domain)
143      :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.           l=V**(1./self.domain.getDim())
144      """
145      return self.__pde_v.setSolverOptions(options)       if perm.getRank()==0:
146      def getSolverOptionsPressure(self):          perm_inv=(1./perm)
147      """              self.scale=util.integrate(perm_inv)/V*l
148      Returns the solver options used to solve the pressure problems          perm_inv=perm_inv*((1./self.scale)*util.kronecker(self.domain.getDim()))
149                perm=perm*(self.scale*util.kronecker(self.domain.getDim()))
150      *(Q^*Q)p=Q^*G*
151
152      :return: `SolverOptions`       elif perm.getRank()==2:
153      """          perm_inv=util.inverse(perm)
154      return self.__pde_p.getSolverOptions()              self.scale=util.sqrt(util.integrate(util.length(perm_inv)**2)/V)*l
155      def setSolverOptionsPressure(self, options=None):          perm_inv*=(1./self.scale)
156      """          perm=perm*self.scale
157      Sets the solver options used to solve the pressure problems       else:
158                raise ValueError,"illegal rank of permeability."
159      *(Q^*Q)p=Q^*G*
160         self.__permeability=perm
161         self.__permeability_inv=perm_inv
162         if self.verbose: print "DarcyFlow: scaling factor for pressure is %e."%self.scale
163
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)
173            self.__pde_v.setValue(D=0.5*self.__permeability_inv)
174         elif  self.solver  == self.SYMSTAB:
175            self.__pde_p.setValue(A=0.5*self.__permeability)
176            self.__pde_v.setValue(D=0.5*self.__permeability_inv)
177
178          if g != None:
179        g=util.interpolate(g, self.__pde_v.getFunctionSpaceForCoefficient("Y"))
180        if g.isEmpty():
181              g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
182        else:
183            if not g.getShape()==(self.domain.getDim(),): raise ValueError,"illegal shape of g"
184        self.__g=g
185          if f !=None:
186         f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
187         if f.isEmpty():
188              f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
189         else:
190             if f.getRank()>0: raise ValueError,"illegal rank of f."
191         self.__f=f
192       def getSolverOptionsFlux(self):
193          """
194          Returns the solver options used to solve the flux problems
195          :return: `SolverOptions`
196          """
197          return self.__pde_v.getSolverOptions()
198
199       def setSolverOptionsFlux(self, options=None):
200          """
201          Sets the solver options used to solve the flux problems
202          If ``options`` is not present, the options are reset to default
203          :param options: `SolverOptions`
204          """
205          return 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       def setSolverOptionsPressure(self, options=None):
215          """
216          Sets the solver options used to solve the pressure problems
217          If ``options`` is not present, the options are reset to default
218
219          :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.
221          """
222          return self.__pde_p.setSolverOptions(options)
223
224       def setTolerance(self,rtol=1e-4):
225          """
226          sets the relative tolerance ``rtol`` for the pressure for the stabelized solvers.
227
228          :param rtol: relative tolerance for the pressure
229          :type rtol: non-negative ``float``
230          """
231          if rtol<0:
232         raise ValueError,"Relative tolerance needs to be non-negative."
233          self.__rtol=rtol
234
235       def getTolerance(self):
236          """
237          returns the relative tolerance
238          :return: current relative tolerance
239          :rtype: ``float``
240          """
241          return self.__rtol
242
243       def solve(self,u0,p0, max_iter=100, iter_restart=20):
244          """
245          solves the problem.
246
247          The iteration is terminated if the residual norm is less then self.getTolerance().
248
249          :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 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          :type p0: scalar value on the domain (e.g. `escript.Data`).
253          :param max_iter: maximum number of (outer) iteration steps for the stabilization solvers,
254          :type max_iter: ``int``
255          :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                               of the iteration. ``iter_restart`` is relevant for the stabilization solvers only.
258          :type iter_restart: ``int``
259          :return: flux and pressure
260          :rtype: ``tuple`` of `escript.Data`.
261
262          """
263          # rescale initial guess:
264          p0=p0/self.scale
265          if self.solver  == self.SIMPLE or self.solver  == self.POST :
266            self.__pde_p.setValue(X=self.__g ,
267                                  Y=self.__f,
268                                  y= - util.inner(self.domain.getNormal(),u0 * self.location_of_fixed_flux),
269                                  r=p0)
270            p=self.__pde_p.getSolution()
271            u = self.getFlux(p, u0)
272          elif  self.solver  == self.STAB:
273        u,p = self.__solve_STAB(u0,p0, max_iter, iter_restart)
274          elif  self.solver  == self.SYMSTAB:
275        u,p = self.__solve_SYMSTAB(u0,p0, max_iter, iter_restart)
276
277      If ``options`` is not present, the options are reset to default        if self.verbose:
279      :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.          def_p=self.__g-(u+KGp)
280      """          def_v=self.__f-util.div(u, self.__pde_v.getFunctionSpaceForCoefficient("X"))
281      return self.__pde_p.setSolverOptions(options)          print "DarcyFlux: |g-u-K*grad(p)|_2 = %e (|u|_2 = %e)."%(self.__L2(def_p),self.__L2(u))
283      def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):        #rescale result
284          """        p=p*self.scale
285          assigns values to model parameters        return u,p
286
287          :param f: volumetic sources/sinks     def getFlux(self,p, u0=None):
:type f: scalar value on the domain (e.g. `Data`)
:param g: flux sources/sinks
:type g: vector values on the domain (e.g. `Data`)
:param location_of_fixed_pressure: mask for locations where pressure is fixed
:type location_of_fixed_pressure: scalar value on the domain (e.g. `Data`)
:param location_of_fixed_flux:  mask for locations where flux is fixed.
:type location_of_fixed_flux: vector values on the domain (e.g. `Data`)
:param permeability: permeability tensor. If scalar ``s`` is given the tensor with
``s`` on the main diagonal is used. If vector ``v`` is given the tensor with
``v`` on the main diagonal is used.
:type permeability: scalar, vector or tensor values on the domain (e.g. `Data`)

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

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

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

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

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

where ``atol`` is an absolut tolerance (see `setAbsoluteTolerance`), *|f|^2 = integrate(length(f)^2)* and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.

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

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

def setAbsoluteTolerance(self,atol=0.):
"""
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. `Data`).
:param p0: initial guess for the pressure. At locations in the domain marked by ``location_of_fixed_pressure`` the value of ``p0`` is kept unchanged.
:type p0: scalar value on the domain (e.g. `Data`).
:param verbose: if set some information on iteration progress are printed
:type verbose: ``bool``
:return: flux and pressure
:rtype: ``tuple`` of `Data`.

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

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

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

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

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

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

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

def __Q(self,p):

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

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

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

def getFlux(self,p=None, fixed_flux=Data()):
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. `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. `Data`).          :type u0: vector values on the domain (e.g. `escript.Data`) or ``None``
297          :return: flux          :return: flux
298          :rtype: `Data`          :rtype: `escript.Data`
:note: the method uses the least squares solution *u=(I+D^*D)^{-1}(D^*f-g-Qp)* where *D* is the *div* operator and *(Qp)_i=k_{ij}p_{,j}*
for the permeability *k_{ij}*
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.__l*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
347          dp=self.__pde_p.getSolution()
348          p=GMRES(dp,
349                  self.__STAB_Aprod,
350              p0,
351              self.__inner,
352              atol=self.__norm(p0+dp)*self.getTolerance() ,
353              rtol=0.,
354              iter_max=max_iter,
355              iter_restart=iter_restart,
356              verbose=self.verbose,P_R=None)
357
358              u=self.getFlux(p, u0)
359              return u,p
360
361       def __solve_SYMSTAB(self, u0, p0, max_iter, iter_restart):
362              # p0 is used as an initial guess
363          u=self.getFlux(p0, u0)
364              self.__pde_p.setValue( Y= self.__f,
365                                     X=  0.5*(self.__g + u - util.tensor_mult(self.__permeability,util.grad(p0)) ),
366                                     y=  -  util.inner(self.domain.getNormal(), u),
367                                     r=escript.Data())
368          dp=self.__pde_p.getSolution()
369
370          p=GMRES(dp,
371                  self.__SYMSTAB_Aprod,
372              p0,
373              self.__inner,
374              atol=self.__norm(p0+dp)*self.getTolerance() ,
375              rtol=0.,
376              iter_max=max_iter,
377              iter_restart=iter_restart,
378              verbose=self.verbose,P_R=None)
379
380              u=self.getFlux(p, u0)
381              return u,p
382
383       def __L2(self,v):
384             return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))
385
386       def __norm(self,r):
387             return util.sqrt(self.__inner(r,r))
388
389       def __inner(self,r,s):
390             return util.integrate(util.inner(r,s), escript.Function(self.domain))
391
392       def __STAB_Aprod(self,p):
394          self.__pde_v.setValue(Y=-0.5*gp,
395                                X=-p*util.kronecker(self.__pde_v.getDomain()),
396                                y= p * self.domain.getNormal(),
397                                r=escript.Data())
398          u = -self.__pde_v.getSolution()
399          self.__pde_p.setValue(Y=util.div(u),
400                                X=0.5*(u+util.tensor_mult(self.__permeability,gp)),
401                                y=escript.Data(),
402                                r=escript.Data())
403
404          return  self.__pde_p.getSolution()
405
406       def __SYMSTAB_Aprod(self,p):
408          self.__pde_v.setValue(Y=0.5*gp ,
409                                X=escript.Data(),
410                                y=escript.Data(),
411                                r=escript.Data())
412          u = -self.__pde_v.getSolution()
413          self.__pde_p.setValue(Y=escript.Data(),
414                                X=0.5*(-u+util.tensor_mult(self.__permeability,gp)),
415                                y=escript.Data(),
416                                r=escript.Data())
417
418          return  self.__pde_p.getSolution()
419
420
422       """       """
# Line 390  class StokesProblemCartesian(Homogeneous Line 450  class StokesProblemCartesian(Homogeneous
450           """           """
452           self.domain=domain           self.domain=domain
453           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())
454           self.__pde_u.setSymmetryOn()           self.__pde_v.setSymmetryOn()
455
456           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
457           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
# Line 408  class StokesProblemCartesian(Homogeneous Line 468  class StokesProblemCartesian(Homogeneous
468
469       :rtype: `SolverOptions`       :rtype: `SolverOptions`
470       """       """
471       return self.__pde_u.getSolverOptions()       return self.__pde_v.getSolverOptions()
472       def setSolverOptionsVelocity(self, options=None):       def setSolverOptionsVelocity(self, options=None):
473           """           """
474       set the solver options for solving the equation for velocity.       set the solver options for solving the equation for velocity.
# Line 416  class StokesProblemCartesian(Homogeneous Line 476  class StokesProblemCartesian(Homogeneous
476       :param options: new solver  options       :param options: new solver  options
477       :type options: `SolverOptions`       :type options: `SolverOptions`
478       """       """
479           self.__pde_u.setSolverOptions(options)           self.__pde_v.setSolverOptions(options)
480       def getSolverOptionsPressure(self):       def getSolverOptionsPressure(self):
481           """           """
482       returns the solver options used  solve the equation for pressure.       returns the solver options used  solve the equation for pressure.
# Line 473  class StokesProblemCartesian(Homogeneous Line 533  class StokesProblemCartesian(Homogeneous
533          if eta !=None:          if eta !=None:
534              k=util.kronecker(self.domain.getDim())              k=util.kronecker(self.domain.getDim())
535              kk=util.outer(k,k)              kk=util.outer(k,k)
536              self.eta=util.interpolate(eta, Function(self.domain))              self.eta=util.interpolate(eta, escript.Function(self.domain))
537          self.__pde_prec.setValue(D=1/self.eta)          self.__pde_prec.setValue(D=1/self.eta)
538              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)))
539          if restoration_factor!=None:          if restoration_factor!=None:
540              n=self.domain.getNormal()              n=self.domain.getNormal()
541              self.__pde_u.setValue(d=restoration_factor*util.outer(n,n))              self.__pde_v.setValue(d=restoration_factor*util.outer(n,n))
544          if f!=None: self.__f=f          if f!=None: self.__f=f
545          if surface_stress!=None: self.__surface_stress=surface_stress          if surface_stress!=None: self.__surface_stress=surface_stress
546          if stress!=None: self.__stress=stress          if stress!=None: self.__stress=stress
547
549          """          """
550          assigns values to the model parameters          assigns values to the model parameters
551
# Line 525  class StokesProblemCartesian(Homogeneous Line 585  class StokesProblemCartesian(Homogeneous
585           :return: inner product of element p and Bv=-div(v)           :return: inner product of element p and Bv=-div(v)
586           :rtype: ``float``           :rtype: ``float``
587           """           """
588           return util.integrate(util.interpolate(p,Function(self.domain))*util.interpolate(Bv,Function(self.domain)))           return util.integrate(util.interpolate(p,escript.Function(self.domain))*util.interpolate(Bv, escript.Function(self.domain)))
589
590       def inner_p(self,p0,p1):       def inner_p(self,p0,p1):
591           """           """
# Line 536  class StokesProblemCartesian(Homogeneous Line 596  class StokesProblemCartesian(Homogeneous
596           :return: inner product of p0 and p1           :return: inner product of p0 and p1
597           :rtype: ``float``           :rtype: ``float``
598           """           """
599           s0=util.interpolate(p0,Function(self.domain))           s0=util.interpolate(p0, escript.Function(self.domain))
600           s1=util.interpolate(p1,Function(self.domain))           s1=util.interpolate(p1, escript.Function(self.domain))
601           return util.integrate(s0*s1)           return util.integrate(s0*s1)
602
603       def norm_v(self,v):       def norm_v(self,v):
# Line 560  class StokesProblemCartesian(Homogeneous Line 620  class StokesProblemCartesian(Homogeneous
621           """           """
623           self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress)           self.__pde_v.setValue(Y=self.__f, y=self.__surface_stress)
624       self.getSolverOptionsVelocity().setTolerance(tol)       self.getSolverOptionsVelocity().setTolerance(tol)
625       self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)       self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
626           if self.__stress.isEmpty():           if self.__stress.isEmpty():
628           else:           else:
630           out=self.__pde_u.getSolution()           out=self.__pde_v.getSolution()
631           return  out           return  out
632
633       def norm_Bv(self,Bv):       def norm_Bv(self,Bv):
# Line 577  class StokesProblemCartesian(Homogeneous Line 637  class StokesProblemCartesian(Homogeneous
637          :rtype: equal to the type of p          :rtype: equal to the type of p
638          :note: boundary conditions on p should be zero!          :note: boundary conditions on p should be zero!
639          """          """
640          return util.sqrt(util.integrate(util.interpolate(Bv,Function(self.domain))**2))          return util.sqrt(util.integrate(util.interpolate(Bv, escript.Function(self.domain))**2))
641
642       def solve_AinvBt(self,p, tol):       def solve_AinvBt(self,p, tol):
643           """           """
# Line 587  class StokesProblemCartesian(Homogeneous Line 647  class StokesProblemCartesian(Homogeneous
647           :return: the solution of *Av=B^*p*           :return: the solution of *Av=B^*p*
648           :note: boundary conditions on v should be zero!           :note: boundary conditions on v should be zero!
649           """           """
650           self.__pde_u.setValue(Y=Data(), y=Data(), X=-p*util.kronecker(self.domain))           self.__pde_v.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))
651           out=self.__pde_u.getSolution()           out=self.__pde_v.getSolution()
652           return  out           return  out
653
654       def solve_prec(self,Bv, tol):       def solve_prec(self,Bv, tol):

Legend:
 Removed from v.2881 changed lines Added in v.3515