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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3074 by gross, Tue Jul 27 01:47:45 2010 UTC revision 3499 by gross, Fri Apr 8 00:14:35 2011 UTC
# Line 32  Some models for flow Line 32  Some models for flow
32    
33  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
34    
35  from escript import *  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    
41    class DarcyFlowVeryNew(object):
42       """
43       solves the problem
44      
45       *u_i+k_{ij}*p_{,j} = g_i*
46       *u_{i,i} = f*
47      
48       where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,
49      
50       :note: The problem is solved in a stabelized formulation.
51       """
52       def __init__(self, domain, useReduced=False, *args, **kargs):
53          """
54          initializes the Darcy flux problem
55          :param domain: domain of the problem
56          :type domain: `Domain`
57          :param useReduced: uses reduced oreder on flux and pressure
58          :type useReduced: ``bool``
59          :param adaptSubTolerance: switches on automatic subtolerance selection
60          :type adaptSubTolerance: ``bool``
61          :param solveForFlux: if True the solver solves for the flux (do not use!)
62          :type solveForFlux: ``bool``  
63          :param useVPIteration: if True altenative iteration over v and p is performed. Otherwise V and P are calculated in a single PDE.
64          :type useVPIteration: ``bool``    
65          """
66          self.domain=domain
67          useVPIteration=False
68          self.useVPIteration=useVPIteration or True
69          #self.useVPIteration=False
70          self.useReduced=useReduced
71          self.verbose=False
72    
73          if self.useVPIteration:
74         # this does not work yet
75             self.__pde_k=LinearPDESystem(domain)
76             self.__pde_k.setSymmetryOn()
77             if self.useReduced: self.__pde_k.setReducedOrderOn()
78      
79             self.__pde_p=LinearSinglePDE(domain)
80             self.__pde_p.setSymmetryOn()
81             if self.useReduced: self.__pde_p.setReducedOrderOn()
82          else:
83             self.__pde_k=LinearPDE(self.domain, numEquations=self.domain.getDim()+1)
84             self.__pde_k.setSymmetryOff()
85            
86             if self.useReduced: self.__pde_k.setReducedOrderOn()
87             C=self.__pde_k.createCoefficient("C")
88             B=self.__pde_k.createCoefficient("B")
89             for i in range(self.domain.getDim()):
90                B[i,i,self.domain.getDim()]=-1
91                C[self.domain.getDim(),i,i]=1
92                C[i,self.domain.getDim(),i]=-0.5
93                B[self.domain.getDim(),i,i]=0.5
94             self.__pde_k.setValue(C=C, B=B)
95          self.__f=escript.Scalar(0,self.__pde_k.getFunctionSpaceForCoefficient("X"))
96          self.__g=escript.Vector(0,self.__pde_k.getFunctionSpaceForCoefficient("Y"))
97          self.location_of_fixed_pressure = escript.Scalar(0, self.__pde_k.getFunctionSpaceForCoefficient("q"))
98          self.location_of_fixed_flux = escript.Vector(0, self.__pde_k.getFunctionSpaceForCoefficient("q"))
99          self.scale=1.
100       def __L2(self,v):
101             return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))  
102       def __inner_GMRES(self,r,s):
103             return util.integrate(util.inner(r,s))
104            
105       def __Aprod_GMRES(self,p):
106          self.__pde_k.setValue(Y=0.5*util.grad(p), X=p*util.kronecker(self.__pde_k.getDomain()) )
107          du=self.__pde_k.getSolution()
108          self.__pde_p.setValue(Y=util.div(du), X=0.5*(du+util.tensor_mult(self.__permeability,util.grad(p))))
109          return self.__pde_p.getSolution()
110            
111       def getSolverOptionsFlux(self):
112          """
113          Returns the solver options used to solve the flux problems
114          
115          *K^{-1} u=F*
116          
117          :return: `SolverOptions`
118          """
119          return self.__pde_k.getSolverOptions()      
120          
121       def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
122          """
123          assigns values to model parameters
124    
125          :param f: volumetic sources/sinks
126          :type f: scalar value on the domain (e.g. `escript.Data`)
127          :param g: flux sources/sinks
128          :type g: vector values on the domain (e.g. `escript.Data`)
129          :param location_of_fixed_pressure: mask for locations where pressure is fixed
130          :type location_of_fixed_pressure: scalar value on the domain (e.g. `escript.Data`)
131          :param location_of_fixed_flux:  mask for locations where flux is fixed.
132          :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)
133          :param permeability: permeability tensor. If scalar ``s`` is given the tensor with ``s`` on the main diagonal is used.
134          :type permeability: scalar or tensor values on the domain (e.g. `escript.Data`)
135    
136          :note: the values of parameters which are not set by calling ``setValue`` are not altered.
137          :note: at any point on the boundary of the domain the pressure
138                 (``location_of_fixed_pressure`` >0) or the normal component of the
139                 flux (``location_of_fixed_flux[i]>0``) if direction of the normal
140                 is along the *x_i* axis.
141    
142          """
143          if location_of_fixed_pressure!=None: self.location_of_fixed_pressure=location_of_fixed_pressure
144          if location_of_fixed_flux!=None: self.location_of_fixed_flux=location_of_fixed_flux
145          
146          if self.useVPIteration:
147             if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=self.location_of_fixed_pressure)
148             if location_of_fixed_flux!=None: self.__pde_k.setValue(q=self.location_of_fixed_flux)
149          else:
150         if location_of_fixed_pressure!=None or location_of_fixed_flux!=None:
151            q=self.__pde_k.createCoefficient("q")
152            q[self.domain.getDim()]=self.location_of_fixed_pressure
153            q[:self.domain.getDim()]=self.location_of_fixed_flux
154            self.__pde_k.setValue(q=q)
155                
156          # flux is rescaled by the factor mean value(perm_inv)*length where length**self.domain.getDim()=vol(self.domain)
157          if permeability!=None:
158         perm=util.interpolate(permeability,self.__pde_k.getFunctionSpaceForCoefficient("A"))
159             V=util.vol(self.domain)
160         if perm.getRank()==0:
161            perm_inv=(1./perm)
162                self.scale=util.integrate(perm_inv)*V**(1./self.domain.getDim()-1.)
163            perm_inv=perm_inv*((1./self.scale)*util.kronecker(self.domain.getDim()))
164            perm=perm*(self.scale*util.kronecker(self.domain.getDim()))
165         elif perm.getRank()==2:
166            perm_inv=util.inverse(perm)
167                self.scale=util.sqrt(util.integrate(util.length(perm_inv)**2)*V**(2./self.domain.getDim()-1.)/self.domain.getDim())
168            perm_inv*=(1./self.scale)
169            perm=perm*self.scale
170         else:
171            raise ValueError,"illegal rank of permeability."
172            
173         self.__permeability=perm
174         self.__permeability_inv=perm_inv
175         if self.useVPIteration:
176                self.__pde_k.setValue(D=0.5*self.__permeability_inv)
177            self.__pde_p.setValue(A=0.5*self.__permeability)
178             else:
179                D=self.__pde_k.createCoefficient("D")
180                A=self.__pde_k.createCoefficient("A")
181                D[:self.domain.getDim(),:self.domain.getDim()]=0.5*self.__permeability_inv
182                A[self.domain.getDim(),:,self.domain.getDim(),:]=0.5*self.__permeability
183                self.__pde_k.setValue(A=A, D=D)
184          if g != None:
185        g=util.interpolate(g, self.__pde_k.getFunctionSpaceForCoefficient("Y"))
186        if g.isEmpty():
187              g=Vector(0,self.__pde_k.getFunctionSpaceForCoefficient("Y"))
188        else:
189            if not g.getShape()==(self.domain.getDim(),): raise ValueError,"illegal shape of g"
190            self.__g=g
191          if f !=None:
192         f=util.interpolate(f, self.__pde_k.getFunctionSpaceForCoefficient("X"))
193         if f.isEmpty():
194          
195              f=Scalar(0,self.__pde_k.getFunctionSpaceForCoefficient("X"))
196         else:
197             if f.getRank()>0: raise ValueError,"illegal rank of f."
198             self.__f=f
199            
200       def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):
201          """
202          solves the problem.
203          
204          The iteration is terminated if the residual norm is less then self.getTolerance().
205    
206          :param u0: initial guess for the flux. At locations in the domain marked by ``location_of_fixed_flux`` the value of ``u0`` is kept unchanged.
207          :type u0: vector value on the domain (e.g. `escript.Data`).
208          :param p0: initial guess for the pressure. At locations in the domain marked by ``location_of_fixed_pressure`` the value of ``p0`` is kept unchanged.
209          :type p0: scalar value on the domain (e.g. `escript.Data`).
210          :param verbose: if set some information on iteration progress are printed
211          :type verbose: ``bool``
212          :return: flux and pressure
213          :rtype: ``tuple`` of `escript.Data`.
214    
215          """
216          u0_b=u0*self.location_of_fixed_flux
217          p0_b=p0*self.location_of_fixed_pressure/self.scale
218          f=self.__f-util.div(u0_b)
219          g=self.__g-u0_b - util.tensor_mult(self.__permeability,util.grad(p0_b))
220          self.verbose=verbose
221          if self.useVPIteration:
222            # get u:
223            self.__pde_k.setValue(Y=0.5*util.tensor_mult(self.__permeability_inv,g),X=escript.Data())
224            du=self.__pde_k.getSolution()
225            self.__pde_p.setValue(Y=f-util.div(du), X=0.5*(g-du))
226            p=GMRES(self.__pde_p.getSolution(),
227                    self.__Aprod_GMRES,
228                p0_b*0,
229                self.__inner_GMRES,
230                atol=0,
231                rtol=1.e-4,
232                iter_max=100,
233                iter_restart=20, verbose=self.verbose,P_R=None)
234                self.__pde_k.setValue(Y=0.5*( util.tensor_mult(self.__permeability_inv,g) + util.grad(p)) ,
235                                      X=p*util.kronecker(self.__pde_k.getDomain()))
236            u=self.__pde_k.getSolution()
237          else:
238              X=self.__pde_k.createCoefficient("X")
239              Y=self.__pde_k.createCoefficient("Y")
240              Y[:self.domain.getDim()]=0.5*util.tensor_mult(self.__permeability_inv,g)
241              Y[self.domain.getDim()]=f
242              X[self.domain.getDim(),:]=g*0.5
243              self.__pde_k.setValue(X=X, Y=Y)
244              self.__pde_k.getSolverOptions().setVerbosity(self.verbose)
245              #self.__pde_k.getSolverOptions().setPreconditioner(self.__pde_k.getSolverOptions().AMG)
246              self.__pde_k.getSolverOptions().setSolverMethod(self.__pde_k.getSolverOptions().DIRECT)
247              U=self.__pde_k.getSolution()
248              u=U[:self.domain.getDim()]
249              p=U[self.domain.getDim()]
250          # self.__pde_k.getOperator().saveMM("k.mm")
251          u=u0_b+u
252          p=(p0_b+p)*self.scale
253          if self.verbose:
254            KGp=util.tensor_mult(self.__permeability,util.grad(p)/self.scale)
255            def_p=self.__g-(u+KGp)
256            def_v=self.__f-util.div(u, self.__pde_k.getFunctionSpaceForCoefficient("X"))
257            print "DarcyFlux: L2: g-v-K*grad(p) = %e (v = %e)."%(self.__L2(def_p),self.__L2(u))
258            print "DarcyFlux: L2: f-div(v) = %e (grad(v) = %e)."%(self.__L2(def_v),self.__L2(util.grad(u)))
259          return u,p
260       def setTolerance(self,rtol=1e-4):
261          """
262          sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if
263    
264          *|g-v-K gard(p)|_PCG <= atol + rtol * |K^{1/2}g2|_0*
265          
266          where ``atol`` is an absolut tolerance (see `setAbsoluteTolerance`).
267          
268          :param rtol: relative tolerance for the pressure
269          :type rtol: non-negative ``float``
270          """
271          if rtol<0:
272         raise ValueError,"Relative tolerance needs to be non-negative."
273          self.__rtol=rtol
274       def getTolerance(self):
275          """
276          returns the relative tolerance
277          :return: current relative tolerance
278          :rtype: ``float``
279          """
280          return self.__rtol
281  class DarcyFlow(object):  class DarcyFlow(object):
282     """     """
283     solves the problem     solves the problem
# Line 49  class DarcyFlow(object): Line 290  class DarcyFlow(object):
290     :note: The problem is solved in a least squares formulation.     :note: The problem is solved in a least squares formulation.
291     """     """
292        
293     def __init__(self, domain, useReduced=False, adaptSubTolerance=True, solveForFlux=False):     def __init__(self, domain, useReduced=False, adaptSubTolerance=True, solveForFlux=False, useVPIteration=True, weighting_scale=0.1):
294        """        """
295        initializes the Darcy flux problem        initializes the Darcy flux problem
296        :param domain: domain of the problem        :param domain: domain of the problem
# Line 60  class DarcyFlow(object): Line 301  class DarcyFlow(object):
301        :type adaptSubTolerance: ``bool``        :type adaptSubTolerance: ``bool``
302        :param solveForFlux: if True the solver solves for the flux (do not use!)        :param solveForFlux: if True the solver solves for the flux (do not use!)
303        :type solveForFlux: ``bool``          :type solveForFlux: ``bool``  
304          :param useVPIteration: if True altenative iteration over v and p is performed. Otherwise V and P are calculated in a single PDE.
305          :type useVPIteration: ``bool``    
306        """        """
307        self.domain=domain        self.domain=domain
308        self.solveForFlux=solveForFlux        self.useVPIteration=useVPIteration
309        self.useReduced=useReduced        self.useReduced=useReduced
310        self.__adaptSubTolerance=adaptSubTolerance        self.weighting_scale=weighting_scale
311        self.verbose=False        if self.useVPIteration:
312                   self.solveForFlux=solveForFlux
313        self.__pde_k=LinearPDESystem(domain)           self.__adaptSubTolerance=adaptSubTolerance
314        self.__pde_k.setSymmetryOn()           self.verbose=False
315        if self.useReduced: self.__pde_k.setReducedOrderOn()          
316             self.__pde_k=LinearPDESystem(domain)
317        self.__pde_p=LinearSinglePDE(domain)           self.__pde_k.setSymmetryOn()
318        self.__pde_p.setSymmetryOn()           if self.useReduced: self.__pde_k.setReducedOrderOn()
319        if self.useReduced: self.__pde_p.setReducedOrderOn()    
320             self.__pde_p=LinearSinglePDE(domain)
321        self.__pde_l=LinearSinglePDE(domain)   # this is here for getSolverOptionsWeighting           self.__pde_p.setSymmetryOn()
322        # self.__pde_l.setSymmetryOn()           if self.useReduced: self.__pde_p.setReducedOrderOn()
323        # if self.useReduced: self.__pde_l.setReducedOrderOn()           self.setTolerance()
324        self.setTolerance()           self.setAbsoluteTolerance()
325        self.setAbsoluteTolerance()        else:
326        self.__f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))           self.__pde_k=LinearPDE(self.domain, numEquations=self.domain.getDim()+1)
327        self.__g=Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))           self.__pde_k.setSymmetryOn()
328             if self.useReduced: self.__pde_k.setReducedOrderOn()
329             C=self.__pde_k.createCoefficient("C")
330             B=self.__pde_k.createCoefficient("B")
331             for i in range(self.domain.getDim()):
332                C[i,self.domain.getDim(),i]=1
333                B[self.domain.getDim(),i,i]=1
334             self.__pde_k.setValue(C=C, B=B)
335          self.__f=escript.Scalar(0,self.__pde_k.getFunctionSpaceForCoefficient("X"))
336          self.__g=escript.Vector(0,self.__pde_k.getFunctionSpaceForCoefficient("Y"))
337                
338     def getSolverOptionsFlux(self):     def getSolverOptionsFlux(self):
339        """        """
# Line 129  class DarcyFlow(object): Line 381  class DarcyFlow(object):
381        """        """
382        return self.__pde_p.setSolverOptions(options)        return self.__pde_p.setSolverOptions(options)
383    
    def getSolverOptionsWeighting(self):  
       """  
       Returns the solver options used to solve the pressure problems  
   
       *(D K D^*)p=-D F*  
   
       :return: `SolverOptions`  
       """  
       return self.__pde_l.getSolverOptions()  
   
    def setSolverOptionsWeighting(self, options=None):  
       """  
       Sets the solver options used to solve the pressure problems  
   
       *(D K D^*)p=-D F*  
   
       If ``options`` is not present, the options are reset to default  
   
       :param options: `SolverOptions`  
       :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.  
       """  
       return self.__pde_l.setSolverOptions(options)  
   
384    
385     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):
386        """        """
387        assigns values to model parameters        assigns values to model parameters
388          
389        :param f: volumetic sources/sinks        :param f: volumetic sources/sinks
390        :type f: scalar value on the domain (e.g. `Data`)        :type f: scalar value on the domain (e.g. `escript.Data`)
391        :param g: flux sources/sinks        :param g: flux sources/sinks
392        :type g: vector values on the domain (e.g. `Data`)        :type g: vector values on the domain (e.g. `escript.Data`)
393        :param location_of_fixed_pressure: mask for locations where pressure is fixed        :param location_of_fixed_pressure: mask for locations where pressure is fixed
394        :type location_of_fixed_pressure: scalar value on the domain (e.g. `Data`)        :type location_of_fixed_pressure: scalar value on the domain (e.g. `escript.Data`)
395        :param location_of_fixed_flux:  mask for locations where flux is fixed.        :param location_of_fixed_flux:  mask for locations where flux is fixed.
396        :type location_of_fixed_flux: vector values on the domain (e.g. `Data`)        :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)
397        :param permeability: permeability tensor. If scalar ``s`` is given the tensor with        :param permeability: permeability tensor. If scalar ``s`` is given the tensor with ``s`` on the main diagonal is used.
398        ``s`` on the main diagonal is used.        :type permeability: scalar or tensor values on the domain (e.g. `escript.Data`)
399        :type permeability: scalar or tensor values on the domain (e.g. `Data`)  
400        :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.
401        :note: at any point on the boundary of the domain the pressure (``location_of_fixed_pressure`` >0)        :note: at any point on the boundary of the domain the pressure
402        or the normal component of the flux (``location_of_fixed_flux[i]>0`` if direction of the normal               (``location_of_fixed_pressure`` >0) or the normal component of the
403        is along the *x_i* axis.               flux (``location_of_fixed_flux[i]>0``) if direction of the normal
404                 is along the *x_i* axis.
405    
406        """        """
407        if f !=None:        if self.useVPIteration:
408       f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("X"))           if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)
409       if f.isEmpty():           if location_of_fixed_flux!=None: self.__pde_k.setValue(q=location_of_fixed_flux)
410          f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))        else:
411       else:           q=self.__pde_k.getCoefficient("q")
412          if f.getRank()>0: raise ValueError,"illegal rank of f."           if q.isEmpty(): q=self.__pde_k.createCoefficient("q")
413          self.__f=f           if location_of_fixed_pressure!=None: q[self.domain.getDim()]=location_of_fixed_pressure
414        if g !=None:           if location_of_fixed_flux!=None: q[:self.domain.getDim()]=location_of_fixed_flux
415       g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))           self.__pde_k.setValue(q=q)
      if g.isEmpty():  
         g=Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))  
      else:  
         if not g.getShape()==(self.domain.getDim(),):  
            raise ValueError,"illegal shape of g"  
         self.__g=g  
       if location_of_fixed_pressure!=None:  
            self.__pde_p.setValue(q=location_of_fixed_pressure)  
            #self.__pde_l.setValue(q=location_of_fixed_pressure)  
       if location_of_fixed_flux!=None:  
            self.__pde_k.setValue(q=location_of_fixed_flux)  
416                            
417          # flux is rescaled by the factor mean value(perm_inv)*length where length**self.domain.getDim()=vol(self.domain)
418        if permeability!=None:        if permeability!=None:
419       perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))       perm=util.interpolate(permeability,self.__pde_k.getFunctionSpaceForCoefficient("A"))
420             V=util.vol(self.domain)
421       if perm.getRank()==0:       if perm.getRank()==0:
422          perm_inv=(1./perm)*util.kronecker(self.domain.getDim())          perm_inv=(1./perm)
423          perm=perm*util.kronecker(self.domain.getDim())              if self.useVPIteration:
424                  self.scale=1.
425                else:
426                  self.scale=util.integrate(perm_inv)*V**(1./self.domain.getDim()-1.)
427    
428            perm_inv=perm_inv*((1./self.scale)*util.kronecker(self.domain.getDim()))
429            perm=perm*(self.scale*util.kronecker(self.domain.getDim()))
430       elif perm.getRank()==2:       elif perm.getRank()==2:
431          perm_inv=util.inverse(perm)          perm_inv=util.inverse(perm)
432                if self.useVPIteration:
433                  self.scale=1.
434                else:
435                  self.scale=util.sqrt(util.integrate(util.length(perm_inv)**2)*V**(2./self.domain.getDim()-1.)/self.domain.getDim())
436              perm_inv*=(1./self.scale)
437              perm=perm*self.scale
438       else:       else:
439          raise ValueError,"illegal rank of permeability."          raise ValueError,"illegal rank of permeability."
440    
441       self.__permeability=perm       self.__permeability=perm
442       self.__permeability_inv=perm_inv       self.__permeability_inv=perm_inv
      self.__l =(util.longestEdge(self.domain)**2*util.length(self.__permeability_inv))/10  
      #self.__l =(self.domain.getSize()**2*util.length(self.__permeability_inv))/10  
      if  self.solveForFlux:  
         self.__pde_k.setValue(D=self.__permeability_inv)  
      else:  
         self.__pde_k.setValue(D=self.__permeability_inv, A=self.__l*util.outer(util.kronecker(self.domain),util.kronecker(self.domain)))  
      self.__pde_p.setValue(A=self.__permeability)  
      #self.__pde_l.setValue(D=1/self.__l)  
          #self.__pde_l.setValue(A=self.__permeability)  
   
    def __applWeight(self, v, f=None):  
       # solves L p = f-Dv with p = 0  
       if self.getSolverOptionsWeighting().isVerbose() or self.verbose: print "DarcyFlux: Applying weighting operator"  
       if f == None:  
      return -util.div(v)*self.__l  
       else:  
      return (f-util.div(v))*self.__l  
       # if f == None:  
       #      self.__pde_l.setValue(Y=-util.div(v))    
       # else:  
       #      return (f-util.div(v))/self.__l  
       # return self.__pde_l.getSolution()  
         
    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  
443    
444     def __Aprod_v(self,dv):       self.__l2 =(util.longestEdge(self.domain)**2*util.length(self.__permeability_inv))*self.weighting_scale
445        # calculates: (a,b,c) = (K^{-1}(dv + KG * dp), L^{-1}Ddv, dp)  with (G*KG)dp = - G^*dv             if self.useVPIteration:
446        dp=self.__getPressure(dv, p0=Data()) # dp = (G*KG)^{-1} (0-G^*dv)          if  self.solveForFlux:
447        a=util.tensor_mult(self.__permeability_inv,dv)+util.grad(dp) # a= K^{-1}u+G*dp             self.__pde_k.setValue(D=self.__permeability_inv)
448        b= - self.__applWeight(dv) # b = - (D K D^*)^{-1} (0-Dv)          else:
449        return ArithmeticTuple(a,b,-dp)             self.__pde_k.setValue(D=self.__permeability_inv, A=self.__l2*util.outer(util.kronecker(self.domain),util.kronecker(self.domain)))
450            self.__pde_p.setValue(A=self.__permeability)
451     def __Msolve_PCG_v(self,r):           else:
452        # K^{-1} u = r[0] + D^*r[1] = K^{-1}(dv + KG * dp) + D^*L^{-1}Ddv              D=self.__pde_k.createCoefficient("D")
453        if self.getSolverOptionsFlux().isVerbose() or self.verbose: print "DarcyFlux: Applying preconditioner"              A=self.__pde_k.createCoefficient("A")
454        self.__pde_k.setValue(X=r[1]*util.kronecker(self.domain), Y=r[0], r=Data())              D[:self.domain.getDim(),:self.domain.getDim()]=self.__permeability_inv
455        # self.__pde_p.getOperator().saveMM("prec.mm")              for i in range(self.domain.getDim()):
456        return self.__pde_k.getSolution()                 for j in range(self.domain.getDim()):
457                     A[i,i,j,j]=self.__l2
458     def __inner_PCG_v(self,v,r):              A[self.domain.getDim(),:,self.domain.getDim(),:]=self.__permeability
459        return util.integrate(util.inner(v,r[0])+util.div(v)*r[1])              self.__pde_k.setValue(A=A, D=D)
460                if g !=None:
461     def __Aprod_p(self,dp):       g=util.interpolate(g, self.__pde_k.getFunctionSpaceForCoefficient("Y"))
462        if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"       if g.isEmpty():
463        Gdp=util.grad(dp)            g=Vector(0,self.__pde_k.getFunctionSpaceForCoefficient("Y"))
464        self.__pde_k.setValue(Y=-Gdp,X=Data(), r=Data())       else:
465        du=self.__pde_k.getSolution()          if not g.getShape()==(self.domain.getDim(),):
466        # self.__pde_v.getOperator().saveMM("proj.mm")                raise ValueError,"illegal shape of g"
467        return ArithmeticTuple(util.tensor_mult(self.__permeability,Gdp),-du)          self.__g=g
468          elif permeability!=None:
469                 X
470          if f !=None:
471         f=util.interpolate(f, self.__pde_k.getFunctionSpaceForCoefficient("X"))
472         if f.isEmpty():
473              f=Scalar(0,self.__pde_k.getFunctionSpaceForCoefficient("X"))
474         else:
475             if f.getRank()>0: raise ValueError,"illegal rank of f."
476             self.__f=f
477    
    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:  
      self.__pde_k.setValue(Y=g2-util.grad(p))  
       return self.__pde_k.getSolution()    
478                
       #v=self.__getFlux(p, u0, f=self.__f, g=g2)        
    def __Msolve_PCG_p(self,r):  
       if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"  
       self.__pde_p.setValue(X=r[0]-r[1], Y=Data(), r=Data(), y=Data())  
       # self.__pde_p.getOperator().saveMM("prec.mm")  
       return self.__pde_p.getSolution()  
           
    def __inner_PCG_p(self,p,r):  
        return util.integrate(util.inner(util.grad(p), r[0]-r[1]))  
   
    def __L2(self,v):  
       return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))  
   
479     def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):     def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):
480        """        """
481        solves the problem.        solves the problem.
# Line 296  class DarcyFlow(object): Line 483  class DarcyFlow(object):
483        The iteration is terminated if the residual norm is less then self.getTolerance().        The iteration is terminated if the residual norm is less then self.getTolerance().
484    
485        :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.
486        :type u0: vector value on the domain (e.g. `Data`).        :type u0: vector value on the domain (e.g. `escript.Data`).
487        :param p0: initial guess for the pressure. At locations in the domain marked by ``location_of_fixed_pressure`` the value of ``p0`` is kept unchanged.        :param p0: initial guess for the pressure. At locations in the domain marked by ``location_of_fixed_pressure`` the value of ``p0`` is kept unchanged.
488        :type p0: scalar value on the domain (e.g. `Data`).        :type p0: scalar value on the domain (e.g. `escript.Data`).
489        :param verbose: if set some information on iteration progress are printed        :param verbose: if set some information on iteration progress are printed
490        :type verbose: ``bool``        :type verbose: ``bool``
491        :return: flux and pressure        :return: flux and pressure
492        :rtype: ``tuple`` of `Data`.        :rtype: ``tuple`` of `escript.Data`.
493    
494        :note: The problem is solved as a least squares form        :note: The problem is solved as a least squares form
495        *(K^[-1]+D^* (DKD^*)^[-1] D)u+G p=D^* (DKD^*)^[-1] f + K^[-1]g*               *(K^[-1]+D^* l2 D)u+G p=D^* l2 * f + K^[-1]g*
496        *G^*u+G^* K Gp=G^*g*               *G^*u+*G^* K Gp=G^*g*
497                 where *D* is the *div* operator and *(Gp)_i=p_{,i}* for the permeability *K=k_{ij}*.
       where *D* is the *div* operator and *(Gp)_i=p_{,i}* for the permeability *K=k_{ij}*.  
498        """        """
499        self.verbose=verbose        self.verbose=verbose
500          if self.useVPIteration:
501              return self.__solveVP(u0,p0,max_iter,max_num_corrections)
502          else:
503              X=self.__pde_k.createCoefficient("X")
504              Y=self.__pde_k.createCoefficient("Y")
505              Y[:self.domain.getDim()]=self.scale*util.tensor_mult(self.__permeability_inv,self.__g)
506              rtmp=self.__f * self.__l2 * self.scale
507              for i in range(self.domain.getDim()): X[i,i]=rtmp
508              X[self.domain.getDim(),:]=self.__g*self.scale
509              r=self.__pde_k.createCoefficient("r")
510              r[:self.domain.getDim()]=u0*self.scale
511              r[self.domain.getDim()]=p0
512              self.__pde_k.setValue(X=X, Y=Y, r=r)
513              self.__pde_k.getSolverOptions().setVerbosity(self.verbose)
514              #self.__pde_k.getSolverOptions().setPreconditioner(self.__pde_k.getSolverOptions().AMG)
515              self.__pde_k.getSolverOptions().setSolverMethod(self.__pde_k.getSolverOptions().DIRECT)
516              U=self.__pde_k.getSolution()
517              # self.__pde_k.getOperator().saveMM("k.mm")
518              u=U[:self.domain.getDim()]*(1./self.scale)
519              p=U[self.domain.getDim()]
520              if self.verbose:
521            KGp=util.tensor_mult(self.__permeability,util.grad(p)/self.scale)
522            def_p=self.__g-(u+KGp)
523            def_v=self.__f-util.div(u, self.__pde_k.getFunctionSpaceForCoefficient("X"))
524            print "DarcyFlux: L2: g-v-K*grad(p) = %e (v = %e)."%(self.__L2(def_p),self.__L2(u))
525            print "DarcyFlux: L2: f-div(v) = %e (grad(v) = %e)."%(self.__L2(def_v),self.__L2(util.grad(u)))
526              return u,p
527    
528       def __solveVP(self,u0,p0, max_iter=100, max_num_corrections=10):
529          """
530          solves the problem.
531          
532          The iteration is terminated if the residual norm is less than self.getTolerance().
533    
534          :param u0: initial guess for the flux. At locations in the domain marked by ``location_of_fixed_flux`` the value of ``u0`` is kept unchanged.
535          :type u0: vector value on the domain (e.g. `escript.Data`).
536          :param p0: initial guess for the pressure. At locations in the domain marked by ``location_of_fixed_pressure`` the value of ``p0`` is kept unchanged.
537          :type p0: scalar value on the domain (e.g. `escript.Data`).
538          :return: flux and pressure
539          :rtype: ``tuple`` of `escript.Data`.
540    
541          :note: The problem is solved as a least squares form
542                 *(K^[-1]+D^* (DKD^*)^[-1] D)u+G p=D^* (DKD^*)^[-1] f + K^[-1]g*
543                 *G^*u+*G^* K Gp=G^*g*
544                 where *D* is the *div* operator and *(Gp)_i=p_{,i}* for the permeability *K=k_{ij}*.
545          """
546        rtol=self.getTolerance()        rtol=self.getTolerance()
547        atol=self.getAbsoluteTolerance()        atol=self.getAbsoluteTolerance()
548        self.setSubProblemTolerance()        self.setSubProblemTolerance()
# Line 372  class DarcyFlow(object): Line 604  class DarcyFlow(object):
604          if self.verbose: print "DarcyFlux: stopping criterium reached."          if self.verbose: print "DarcyFlux: stopping criterium reached."
605          converged=True          converged=True
606        return v,p+p0        return v,p+p0
607    
608       def __applWeight(self, v, f=None):
609          # solves L p = f-Dv with p = 0
610          if self.verbose: print "DarcyFlux: Applying weighting operator"
611          if f == None:
612         return -util.div(v)*self.__l2
613          else:
614         return (f-util.div(v))*self.__l2
615       def __getPressure(self, v, p0, g=None):
616          # solves (G*KG)p = G^(g-v) with p = p0 where location_of_fixed_pressure>0
617          if self.getSolverOptionsPressure().isVerbose() or self.verbose: print "DarcyFlux: Pressure update"
618          if g == None:
619         self.__pde_p.setValue(X=-v, r=p0)
620          else:
621         self.__pde_p.setValue(X=g-v, r=p0)      
622          p=self.__pde_p.getSolution()
623          return p
624    
625       def __Aprod_v(self,dv):
626          # calculates: (a,b,c) = (K^{-1}(dv + KG * dp), L^{-1}Ddv, dp)  with (G*KG)dp = - G^*dv  
627          dp=self.__getPressure(dv, p0=escript.Data()) # dp = (G*KG)^{-1} (0-G^*dv)
628          a=util.tensor_mult(self.__permeability_inv,dv)+util.grad(dp) # a= K^{-1}u+G*dp
629          b= - self.__applWeight(dv) # b = - (D K D^*)^{-1} (0-Dv)
630          return ArithmeticTuple(a,b,-dp)
631    
632       def __Msolve_PCG_v(self,r):
633          # K^{-1} u = r[0] + D^*r[1] = K^{-1}(dv + KG * dp) + D^*L^{-1}Ddv
634          if self.getSolverOptionsFlux().isVerbose() or self.verbose: print "DarcyFlux: Applying preconditioner"
635          self.__pde_k.setValue(X=r[1]*util.kronecker(self.domain), Y=r[0], r=escript.Data())
636          return self.__pde_k.getSolution()
637    
638       def __inner_PCG_v(self,v,r):
639          return util.integrate(util.inner(v,r[0])+util.div(v)*r[1])
640          
641       def __Aprod_p(self,dp):
642          if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
643          Gdp=util.grad(dp)
644          self.__pde_k.setValue(Y=-Gdp,X=escript.Data(), r=escript.Data())
645          du=self.__pde_k.getSolution()
646          return ArithmeticTuple(util.tensor_mult(self.__permeability,Gdp),-du)
647    
648       def __getFlux(self,p, v0, f=None, g=None):
649          # solves (K^{-1}+D^*L^{-1} D) v = D^*L^{-1}f + K^{-1}g - Gp
650          if f!=None:
651         self.__pde_k.setValue(X=self.__applWeight(v0*0,self.__f)*util.kronecker(self.domain))
652          self.__pde_k.setValue(r=v0)
653          g2=util.tensor_mult(self.__permeability_inv,g)
654          if p == None:
655         self.__pde_k.setValue(Y=g2)
656          else:
657         self.__pde_k.setValue(Y=g2-util.grad(p))
658          return self.__pde_k.getSolution()  
659          
660          #v=self.__getFlux(p, u0, f=self.__f, g=g2)      
661       def __Msolve_PCG_p(self,r):
662          if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"
663          self.__pde_p.setValue(X=r[0]-r[1], Y=escript.Data(), r=escript.Data(), y=escript.Data())
664          return self.__pde_p.getSolution()
665            
666       def __inner_PCG_p(self,p,r):
667           return util.integrate(util.inner(util.grad(p), r[0]-r[1]))
668    
669       def __L2(self,v):
670          return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))
671    
672       def __L2_r(self,v):
673          return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.ReducedFunction(self.domain)))**2))
674    
675     def setTolerance(self,rtol=1e-4):     def setTolerance(self,rtol=1e-4):
676        """        """
677        sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if        sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if
# Line 433  class DarcyFlow(object): Line 733  class DarcyFlow(object):
733       self.getSolverOptionsFlux().setAbsoluteTolerance(0.)       self.getSolverOptionsFlux().setAbsoluteTolerance(0.)
734       self.getSolverOptionsPressure().setTolerance(sub_tol)       self.getSolverOptionsPressure().setTolerance(sub_tol)
735       self.getSolverOptionsPressure().setAbsoluteTolerance(0.)       self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
      self.getSolverOptionsWeighting().setTolerance(sub_tol)  
      self.getSolverOptionsWeighting().setAbsoluteTolerance(0.)  
736       if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol       if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol
737    
738    
# Line 463  class DarcyFlowOld(object): Line 761  class DarcyFlowOld(object):
761          self.domain=domain          self.domain=domain
762          if weight == None:          if weight == None:
763             s=self.domain.getSize()             s=self.domain.getSize()
764             self.__l=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2             self.__l2=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2
765             # self.__l=(3.*util.longestEdge(self.domain))**2             # self.__l2=(3.*util.longestEdge(self.domain))**2
766             #self.__l=(0.1*util.longestEdge(self.domain)*s/util.sup(s))**2             #self.__l2=(0.1*util.longestEdge(self.domain)*s/util.sup(s))**2
767          else:          else:
768             self.__l=weight             self.__l2=weight
769          self.__pde_v=LinearPDESystem(domain)          self.__pde_v=LinearPDESystem(domain)
770          if useReduced: self.__pde_v.setReducedOrderOn()          if useReduced: self.__pde_v.setReducedOrderOn()
771          self.__pde_v.setSymmetryOn()          self.__pde_v.setSymmetryOn()
772          self.__pde_v.setValue(D=util.kronecker(domain), A=self.__l*util.outer(util.kronecker(domain),util.kronecker(domain)))          self.__pde_v.setValue(D=util.kronecker(domain), A=self.__l2*util.outer(util.kronecker(domain),util.kronecker(domain)))
773          self.__pde_p=LinearSinglePDE(domain)          self.__pde_p=LinearSinglePDE(domain)
774          self.__pde_p.setSymmetryOn()          self.__pde_p.setSymmetryOn()
775          if useReduced: self.__pde_p.setReducedOrderOn()          if useReduced: self.__pde_p.setReducedOrderOn()
776          self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))          self.__f=escript.Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
777          self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))          self.__g=escript.Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
778          self.setTolerance()          self.setTolerance()
779          self.setAbsoluteTolerance()          self.setAbsoluteTolerance()
780      self.__adaptSubTolerance=adaptSubTolerance      self.__adaptSubTolerance=adaptSubTolerance
# Line 527  class DarcyFlowOld(object): Line 825  class DarcyFlowOld(object):
825          assigns values to model parameters          assigns values to model parameters
826    
827          :param f: volumetic sources/sinks          :param f: volumetic sources/sinks
828          :type f: scalar value on the domain (e.g. `Data`)          :type f: scalar value on the domain (e.g. `escript.Data`)
829          :param g: flux sources/sinks          :param g: flux sources/sinks
830          :type g: vector values on the domain (e.g. `Data`)          :type g: vector values on the domain (e.g. `escript.Data`)
831          :param location_of_fixed_pressure: mask for locations where pressure is fixed          :param location_of_fixed_pressure: mask for locations where pressure is fixed
832          :type location_of_fixed_pressure: scalar value on the domain (e.g. `Data`)          :type location_of_fixed_pressure: scalar value on the domain (e.g. `escript.Data`)
833          :param location_of_fixed_flux:  mask for locations where flux is fixed.          :param location_of_fixed_flux:  mask for locations where flux is fixed.
834          :type location_of_fixed_flux: vector values on the domain (e.g. `Data`)          :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)
835          :param permeability: permeability tensor. If scalar ``s`` is given the tensor with          :param permeability: permeability tensor. If scalar ``s`` is given the tensor with
836                               ``s`` on the main diagonal is used. If vector ``v`` is given the tensor with                               ``s`` on the main diagonal is used. If vector ``v`` is given the tensor with
837                               ``v`` on the main diagonal is used.                               ``v`` on the main diagonal is used.
838          :type permeability: scalar, vector or tensor values on the domain (e.g. `Data`)          :type permeability: scalar, vector or tensor values on the domain (e.g. `escript.Data`)
839    
840          :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.
841          :note: at any point on the boundary of the domain the pressure (``location_of_fixed_pressure`` >0)          :note: at any point on the boundary of the domain the pressure (``location_of_fixed_pressure`` >0)
# Line 547  class DarcyFlowOld(object): Line 845  class DarcyFlowOld(object):
845          if f !=None:          if f !=None:
846             f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))             f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
847             if f.isEmpty():             if f.isEmpty():
848                 f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))                 f=escript.Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
849             else:             else:
850                 if f.getRank()>0: raise ValueError,"illegal rank of f."                 if f.getRank()>0: raise ValueError,"illegal rank of f."
851             self.__f=f             self.__f=f
852          if g !=None:          if g !=None:
853             g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))             g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
854             if g.isEmpty():             if g.isEmpty():
855               g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))               g=escript.Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
856             else:             else:
857               if not g.getShape()==(self.domain.getDim(),):               if not g.getShape()==(self.domain.getDim(),):
858                 raise ValueError,"illegal shape of g"                 raise ValueError,"illegal shape of g"
# Line 647  class DarcyFlowOld(object): Line 945  class DarcyFlowOld(object):
945           The iteration is terminated if the residual norm is less then self.getTolerance().           The iteration is terminated if the residual norm is less then self.getTolerance().
946    
947           :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.
948           :type u0: vector value on the domain (e.g. `Data`).           :type u0: vector value on the domain (e.g. `escript.Data`).
949           :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.
950           :type p0: scalar value on the domain (e.g. `Data`).           :type p0: scalar value on the domain (e.g. `escript.Data`).
951           :param verbose: if set some information on iteration progress are printed           :param verbose: if set some information on iteration progress are printed
952           :type verbose: ``bool``           :type verbose: ``bool``
953           :return: flux and pressure           :return: flux and pressure
954           :rtype: ``tuple`` of `Data`.           :rtype: ``tuple`` of `escript.Data`.
955    
956           :note: The problem is solved as a least squares form           :note: The problem is solved as a least squares form
957    
# Line 699  class DarcyFlowOld(object): Line 997  class DarcyFlowOld(object):
997                 if self.verbose:                 if self.verbose:
998                      print "DarcyFlux: L2 norm of v = %e."%norm_v                      print "DarcyFlux: L2 norm of v = %e."%norm_v
999                      print "DarcyFlux: L2 norm of k.util.grad(p) = %e."%norm_Qp                      print "DarcyFlux: L2 norm of k.util.grad(p) = %e."%norm_Qp
1000                      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 u = %e."%(util.integrate(util.length(self.__g-util.interpolate(v,escript.Function(self.domain))-Qp)**2)**(0.5),)
1001                      print "DarcyFlux: L2 defect div(v) = %e."%(util.integrate((self.__f-util.div(v))**2)**(0.5),)                      print "DarcyFlux: L2 defect div(v) = %e."%(util.integrate((self.__f-util.div(v))**2)**(0.5),)
1002                      print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL                      print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
1003                 if norm_r == None or norm_r>ATOL:                 if norm_r == None or norm_r>ATOL:
1004                     if num_corrections>max_num_corrections:                     if num_corrections>max_num_corrections:
1005                           raise ValueError,"maximum number of correction steps reached."                           raise ValueError,"maximum number of correction steps reached."
1006                     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)                     p,r, norm_r=PCG(self.__g-util.interpolate(v,escript.Function(self.domain))-Qp,self.__Aprod,p,self.__Msolve_PCG,self.__inner_PCG,atol=0.5*ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
1007                     num_corrections+=1                     num_corrections+=1
1008                 else:                 else:
1009                     converged=True                     converged=True
1010           return v,p           return v,p
1011      def __L2(self,v):      def __L2(self,v):
1012           return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))           return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))
1013    
1014      def __Q(self,p):      def __Q(self,p):
1015            return util.tensor_mult(self.__permeability,util.grad(p))            return util.tensor_mult(self.__permeability,util.grad(p))
# Line 719  class DarcyFlowOld(object): Line 1017  class DarcyFlowOld(object):
1017      def __Aprod(self,dp):      def __Aprod(self,dp):
1018            if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"            if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
1019            Qdp=self.__Q(dp)            Qdp=self.__Q(dp)
1020            self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())            self.__pde_v.setValue(Y=-Qdp,X=escript.Data(), r=escript.Data())
1021            du=self.__pde_v.getSolution()            du=self.__pde_v.getSolution()
           # self.__pde_v.getOperator().saveMM("proj.mm")  
1022            return Qdp+du            return Qdp+du
1023      def __inner_GMRES(self,r,s):      def __inner_GMRES(self,r,s):
1024           return util.integrate(util.inner(r,s))           return util.integrate(util.inner(r,s))
# Line 731  class DarcyFlowOld(object): Line 1028  class DarcyFlowOld(object):
1028    
1029      def __Msolve_PCG(self,r):      def __Msolve_PCG(self,r):
1030        if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"        if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"
1031            self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=Data(), r=Data())            self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=escript.Data(), r=escript.Data())
           # self.__pde_p.getOperator().saveMM("prec.mm")  
1032            return self.__pde_p.getSolution()            return self.__pde_p.getSolution()
1033    
1034      def getFlux(self,p=None, fixed_flux=Data()):      def getFlux(self,p=None, fixed_flux=escript.Data()):
1035          """          """
1036          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 ``fixed_flux``
1037          on locations where ``location_of_fixed_flux`` is positive (see `setValue`).          on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
1038          Note that ``g`` and ``f`` are used, see `setValue`.          Note that ``g`` and ``f`` are used, see `setValue`.
1039    
1040          :param p: pressure.          :param p: pressure.
1041          :type p: scalar value on the domain (e.g. `Data`).          :type p: scalar value on the domain (e.g. `escript.Data`).
1042          :param fixed_flux: flux on the locations of the domain marked be ``location_of_fixed_flux``.          :param fixed_flux: flux on the locations of the domain marked be ``location_of_fixed_flux``.
1043          :type fixed_flux: vector values on the domain (e.g. `Data`).          :type fixed_flux: vector values on the domain (e.g. `escript.Data`).
1044          :return: flux          :return: flux
1045          :rtype: `Data`          :rtype: `escript.Data`
1046          :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}*          :note: the method uses the least squares solution *u=(I+D^*D)^{-1}(D^*f-g-Qp)* where *D* is the *div* operator and *(Qp)_i=k_{ij}p_{,j}*
1047                 for the permeability *k_{ij}*                 for the permeability *k_{ij}*
1048          """          """
1049      self.setSubProblemTolerance()      self.setSubProblemTolerance()
1050          g=self.__g          g=self.__g
1051          f=self.__f          f=self.__f
1052          self.__pde_v.setValue(X=self.__l*f*util.kronecker(self.domain), r=fixed_flux)          self.__pde_v.setValue(X=self.__l2*f*util.kronecker(self.domain), r=fixed_flux)
1053          if p == None:          if p == None:
1054             self.__pde_v.setValue(Y=g)             self.__pde_v.setValue(Y=g)
1055          else:          else:
# Line 875  class StokesProblemCartesian(Homogeneous Line 1171  class StokesProblemCartesian(Homogeneous
1171          if eta !=None:          if eta !=None:
1172              k=util.kronecker(self.domain.getDim())              k=util.kronecker(self.domain.getDim())
1173              kk=util.outer(k,k)              kk=util.outer(k,k)
1174              self.eta=util.interpolate(eta, Function(self.domain))              self.eta=util.interpolate(eta, escript.Function(self.domain))
1175          self.__pde_prec.setValue(D=1/self.eta)          self.__pde_prec.setValue(D=1/self.eta)
1176              self.__pde_u.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))              self.__pde_u.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))
1177          if restoration_factor!=None:          if restoration_factor!=None:
# Line 887  class StokesProblemCartesian(Homogeneous Line 1183  class StokesProblemCartesian(Homogeneous
1183          if surface_stress!=None: self.__surface_stress=surface_stress          if surface_stress!=None: self.__surface_stress=surface_stress
1184          if stress!=None: self.__stress=stress          if stress!=None: self.__stress=stress
1185    
1186       def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data(), restoration_factor=0):       def initialize(self,f=escript.Data(),fixed_u_mask=escript.Data(),eta=1,surface_stress=escript.Data(),stress=escript.Data(), restoration_factor=0):
1187          """          """
1188          assigns values to the model parameters          assigns values to the model parameters
1189    
# Line 927  class StokesProblemCartesian(Homogeneous Line 1223  class StokesProblemCartesian(Homogeneous
1223           :return: inner product of element p and Bv=-div(v)           :return: inner product of element p and Bv=-div(v)
1224           :rtype: ``float``           :rtype: ``float``
1225           """           """
1226           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)))
1227    
1228       def inner_p(self,p0,p1):       def inner_p(self,p0,p1):
1229           """           """
# Line 938  class StokesProblemCartesian(Homogeneous Line 1234  class StokesProblemCartesian(Homogeneous
1234           :return: inner product of p0 and p1           :return: inner product of p0 and p1
1235           :rtype: ``float``           :rtype: ``float``
1236           """           """
1237           s0=util.interpolate(p0,Function(self.domain))           s0=util.interpolate(p0, escript.Function(self.domain))
1238           s1=util.interpolate(p1,Function(self.domain))           s1=util.interpolate(p1, escript.Function(self.domain))
1239           return util.integrate(s0*s1)           return util.integrate(s0*s1)
1240    
1241       def norm_v(self,v):       def norm_v(self,v):
# Line 979  class StokesProblemCartesian(Homogeneous Line 1275  class StokesProblemCartesian(Homogeneous
1275          :rtype: equal to the type of p          :rtype: equal to the type of p
1276          :note: boundary conditions on p should be zero!          :note: boundary conditions on p should be zero!
1277          """          """
1278          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))
1279    
1280       def solve_AinvBt(self,p, tol):       def solve_AinvBt(self,p, tol):
1281           """           """
# Line 989  class StokesProblemCartesian(Homogeneous Line 1285  class StokesProblemCartesian(Homogeneous
1285           :return: the solution of *Av=B^*p*           :return: the solution of *Av=B^*p*
1286           :note: boundary conditions on v should be zero!           :note: boundary conditions on v should be zero!
1287           """           """
1288           self.__pde_u.setValue(Y=Data(), y=Data(), X=-p*util.kronecker(self.domain))           self.__pde_u.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))
1289           out=self.__pde_u.getSolution()           out=self.__pde_u.getSolution()
1290           return  out           return  out
1291    

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

  ViewVC Help
Powered by ViewVC 1.1.26