/[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 2793 by gross, Tue Dec 1 06:10:10 2009 UTC revision 3452 by caltinay, Tue Jan 25 01:53:57 2011 UTC
# Line 1  Line 1 
1    # -*- coding: utf-8 -*-
2  ########################################################  ########################################################
3  #  #
4  # Copyright (c) 2003-2009 by University of Queensland  # Copyright (c) 2003-2010 by University of Queensland
5  # Earth Systems Science Computational Center (ESSCC)  # Earth Systems Science Computational Center (ESSCC)
6  # http://www.uq.edu.au/esscc  # http://www.uq.edu.au/esscc
7  #  #
# Line 10  Line 11 
11  #  #
12  ########################################################  ########################################################
13    
14  __copyright__="""Copyright (c) 2003-2009 by University of Queensland  __copyright__="""Copyright (c) 2003-2010 by University of Queensland
15  Earth Systems Science Computational Center (ESSCC)  Earth Systems Science Computational Center (ESSCC)
16  http://www.uq.edu.au/esscc  http://www.uq.edu.au/esscc
17  Primary Business: Queensland, Australia"""  Primary Business: Queensland, Australia"""
# 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    print dir(escript)
41    
42  class DarcyFlow(object):  class DarcyFlow(object):
43       """
44       solves the problem
45      
46       *u_i+k_{ij}*p_{,j} = g_i*
47       *u_{i,i} = f*
48      
49       where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,
50      
51       :note: The problem is solved in a least squares formulation.
52       """
53      
54       def __init__(self, domain, useReduced=False, adaptSubTolerance=True, solveForFlux=False, useVPIteration=True, weighting_scale=0.1):
55          """
56          initializes the Darcy flux problem
57          :param domain: domain of the problem
58          :type domain: `Domain`
59          :param useReduced: uses reduced oreder on flux and pressure
60          :type useReduced: ``bool``
61          :param adaptSubTolerance: switches on automatic subtolerance selection
62          :type adaptSubTolerance: ``bool``
63          :param solveForFlux: if True the solver solves for the flux (do not use!)
64          :type solveForFlux: ``bool``  
65          :param useVPIteration: if True altenative iteration over v and p is performed. Otherwise V and P are calculated in a single PDE.
66          :type useVPIteration: ``bool``    
67          """
68          self.domain=domain
69          self.useVPIteration=useVPIteration
70          self.useReduced=useReduced
71          self.weighting_scale=weighting_scale
72          if self.useVPIteration:
73             self.solveForFlux=solveForFlux
74             self.__adaptSubTolerance=adaptSubTolerance
75             self.verbose=False
76            
77             self.__pde_k=LinearPDESystem(domain)
78             self.__pde_k.setSymmetryOn()
79             if self.useReduced: self.__pde_k.setReducedOrderOn()
80      
81             self.__pde_p=LinearSinglePDE(domain)
82             self.__pde_p.setSymmetryOn()
83             if self.useReduced: self.__pde_p.setReducedOrderOn()
84             self.setTolerance()
85             self.setAbsoluteTolerance()
86          else:
87             self.__pde_k=LinearPDE(self.domain, numEquations=self.domain.getDim()+1)
88             self.__pde_k.setSymmetryOn()
89             if self.useReduced: self.__pde_k.setReducedOrderOn()
90             C=self.__pde_k.createCoefficient("C")
91             B=self.__pde_k.createCoefficient("B")
92             for i in range(self.domain.getDim()):
93                C[i,self.domain.getDim(),i]=1
94                B[self.domain.getDim(),i,i]=1
95             self.__pde_k.setValue(C=C, B=B)
96          self.__f=escript.Scalar(0,self.__pde_k.getFunctionSpaceForCoefficient("X"))
97          self.__g=escript.Vector(0,self.__pde_k.getFunctionSpaceForCoefficient("Y"))
98          
99       def getSolverOptionsFlux(self):
100          """
101          Returns the solver options used to solve the flux problems
102          
103          *K^{-1} u=F*
104          
105          :return: `SolverOptions`
106          """
107          return self.__pde_k.getSolverOptions()
108          
109       def setSolverOptionsFlux(self, options=None):
110          """
111          Sets the solver options used to solve the flux problems
112          
113          *K^{-1}u=F*
114          
115          If ``options`` is not present, the options are reset to default
116          
117          :param options: `SolverOptions`
118          :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
119          """
120          return self.__pde_v.setSolverOptions(options)
121        
122       def getSolverOptionsPressure(self):
123          """
124          Returns the solver options used to solve the pressure problems
125          
126          *(Q^* K Q)p=-Q^*G*
127          
128          :return: `SolverOptions`
129          """
130          return self.__pde_p.getSolverOptions()
131          
132       def setSolverOptionsPressure(self, options=None):
133          """
134          Sets the solver options used to solve the pressure problems
135          
136          *(Q^* K Q)p=-Q^*G*
137          
138          If ``options`` is not present, the options are reset to default
139          
140          :param options: `SolverOptions`
141          :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
142          """
143          return self.__pde_p.setSolverOptions(options)
144    
145    
146       def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
147          """
148          assigns values to model parameters
149    
150          :param f: volumetic sources/sinks
151          :type f: scalar value on the domain (e.g. `escript.Data`)
152          :param g: flux sources/sinks
153          :type g: vector values on the domain (e.g. `escript.Data`)
154          :param location_of_fixed_pressure: mask for locations where pressure is fixed
155          :type location_of_fixed_pressure: scalar value on the domain (e.g. `escript.Data`)
156          :param location_of_fixed_flux:  mask for locations where flux is fixed.
157          :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)
158          :param permeability: permeability tensor. If scalar ``s`` is given the tensor with ``s`` on the main diagonal is used.
159          :type permeability: scalar or tensor values on the domain (e.g. `escript.Data`)
160    
161          :note: the values of parameters which are not set by calling ``setValue`` are not altered.
162          :note: at any point on the boundary of the domain the pressure
163                 (``location_of_fixed_pressure`` >0) or the normal component of the
164                 flux (``location_of_fixed_flux[i]>0``) if direction of the normal
165                 is along the *x_i* axis.
166    
167          """
168          if self.useVPIteration:
169             if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)
170             if location_of_fixed_flux!=None: self.__pde_k.setValue(q=location_of_fixed_flux)
171          else:
172             q=self.__pde_k.getCoefficient("q")
173             if q.isEmpty(): q=self.__pde_k.createCoefficient("q")
174             if location_of_fixed_pressure!=None: q[self.domain.getDim()]=location_of_fixed_pressure
175             if location_of_fixed_flux!=None: q[:self.domain.getDim()]=location_of_fixed_flux
176             self.__pde_k.setValue(q=q)
177                
178          # flux is rescaled by the factor mean value(perm_inv)*length where length**self.domain.getDim()=vol(self.domain)
179          if permeability!=None:
180         perm=util.interpolate(permeability,self.__pde_k.getFunctionSpaceForCoefficient("A"))
181             V=util.vol(self.domain)
182         if perm.getRank()==0:
183            perm_inv=(1./perm)
184                if self.useVPIteration:
185                  self.scale=1.
186                else:
187                  self.scale=util.integrate(perm_inv)*V**(1./self.domain.getDim()-1.)
188    
189            perm_inv=perm_inv*((1./self.scale)*util.kronecker(self.domain.getDim()))
190            perm=perm*(self.scale*util.kronecker(self.domain.getDim()))
191         elif perm.getRank()==2:
192            perm_inv=util.inverse(perm)
193                if self.useVPIteration:
194                  self.scale=1.
195                else:
196                  self.scale=util.sqrt(util.integrate(util.length(perm_inv)**2)*V**(2./self.domain.getDim()-1.)/self.domain.getDim())
197              perm_inv*=(1./self.scale)
198              perm=perm*self.scale
199         else:
200            raise ValueError,"illegal rank of permeability."
201    
202         self.__permeability=perm
203         self.__permeability_inv=perm_inv
204    
205         self.__l2 =(util.longestEdge(self.domain)**2*util.length(self.__permeability_inv))*self.weighting_scale
206             if self.useVPIteration:
207            if  self.solveForFlux:
208               self.__pde_k.setValue(D=self.__permeability_inv)
209            else:
210               self.__pde_k.setValue(D=self.__permeability_inv, A=self.__l2*util.outer(util.kronecker(self.domain),util.kronecker(self.domain)))
211            self.__pde_p.setValue(A=self.__permeability)
212             else:
213                D=self.__pde_k.createCoefficient("D")
214                A=self.__pde_k.createCoefficient("A")
215                D[:self.domain.getDim(),:self.domain.getDim()]=self.__permeability_inv
216                for i in range(self.domain.getDim()):
217                   for j in range(self.domain.getDim()):
218                     A[i,i,j,j]=self.__l2
219                A[self.domain.getDim(),:,self.domain.getDim(),:]=self.__permeability
220                self.__pde_k.setValue(A=A, D=D)
221          if g !=None:
222         g=util.interpolate(g, self.__pde_k.getFunctionSpaceForCoefficient("Y"))
223         if g.isEmpty():
224              g=Vector(0,self.__pde_k.getFunctionSpaceForCoefficient("Y"))
225         else:
226            if not g.getShape()==(self.domain.getDim(),):
227                  raise ValueError,"illegal shape of g"
228            self.__g=g
229          elif permeability!=None:
230                 X
231          if f !=None:
232         f=util.interpolate(f, self.__pde_k.getFunctionSpaceForCoefficient("X"))
233         if f.isEmpty():
234              f=Scalar(0,self.__pde_k.getFunctionSpaceForCoefficient("X"))
235         else:
236             if f.getRank()>0: raise ValueError,"illegal rank of f."
237             self.__f=f
238    
239          
240       def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):
241          """
242          solves the problem.
243          
244          The iteration is terminated if the residual norm is less then self.getTolerance().
245    
246          :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.
247          :type u0: vector value on the domain (e.g. `escript.Data`).
248          :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.
249          :type p0: scalar value on the domain (e.g. `escript.Data`).
250          :param verbose: if set some information on iteration progress are printed
251          :type verbose: ``bool``
252          :return: flux and pressure
253          :rtype: ``tuple`` of `escript.Data`.
254    
255          :note: The problem is solved as a least squares form
256                 *(K^[-1]+D^* l2 D)u+G p=D^* l2 * f + K^[-1]g*
257                 *G^*u+*G^* K Gp=G^*g*
258                 where *D* is the *div* operator and *(Gp)_i=p_{,i}* for the permeability *K=k_{ij}*.
259          """
260          self.verbose=verbose
261          if self.useVPIteration:
262              return self.__solveVP(u0,p0,max_iter,max_num_corrections)
263          else:
264              X=self.__pde_k.createCoefficient("X")
265              Y=self.__pde_k.createCoefficient("Y")
266              Y[:self.domain.getDim()]=self.scale*util.tensor_mult(self.__permeability_inv,self.__g)
267              rtmp=self.__f * self.__l2 * self.scale
268              for i in range(self.domain.getDim()): X[i,i]=rtmp
269              X[self.domain.getDim(),:]=self.__g*self.scale
270              r=self.__pde_k.createCoefficient("r")
271              r[:self.domain.getDim()]=u0*self.scale
272              r[self.domain.getDim()]=p0
273              self.__pde_k.setValue(X=X, Y=Y, r=r)
274              self.__pde_k.getSolverOptions().setVerbosity(self.verbose)
275              #self.__pde_k.getSolverOptions().setPreconditioner(self.__pde_k.getSolverOptions().AMG)
276              self.__pde_k.getSolverOptions().setSolverMethod(self.__pde_k.getSolverOptions().DIRECT)
277              U=self.__pde_k.getSolution()
278              # self.__pde_k.getOperator().saveMM("k.mm")
279              u=U[:self.domain.getDim()]*(1./self.scale)
280              p=U[self.domain.getDim()]
281              if self.verbose:
282            KGp=util.tensor_mult(self.__permeability,util.grad(p)/self.scale)
283            def_p=self.__g-(u+KGp)
284            def_v=self.__f-util.div(u, self.__pde_k.getFunctionSpaceForCoefficient("X"))
285            print "DarcyFlux: L2: g-v-K*grad(p) = %e (v = %e)."%(self.__L2(def_p),self.__L2(u))
286            print "DarcyFlux: L2: f-div(v) = %e (grad(v) = %e)."%(self.__L2(def_v),self.__L2(util.grad(u)))
287              return u,p
288    
289       def __solveVP(self,u0,p0, max_iter=100, max_num_corrections=10):
290          """
291          solves the problem.
292          
293          The iteration is terminated if the residual norm is less than self.getTolerance().
294    
295          :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.
296          :type u0: vector value on the domain (e.g. `escript.Data`).
297          :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.
298          :type p0: scalar value on the domain (e.g. `escript.Data`).
299          :return: flux and pressure
300          :rtype: ``tuple`` of `escript.Data`.
301    
302          :note: The problem is solved as a least squares form
303                 *(K^[-1]+D^* (DKD^*)^[-1] D)u+G p=D^* (DKD^*)^[-1] f + K^[-1]g*
304                 *G^*u+*G^* K Gp=G^*g*
305                 where *D* is the *div* operator and *(Gp)_i=p_{,i}* for the permeability *K=k_{ij}*.
306          """
307          rtol=self.getTolerance()
308          atol=self.getAbsoluteTolerance()
309          self.setSubProblemTolerance()
310          num_corrections=0
311          converged=False
312          norm_r=None
313          
314          # Eliminate the hydrostatic pressure:
315          if self.verbose: print "DarcyFlux: calculate hydrostatic pressure component."
316          self.__pde_p.setValue(X=self.__g, r=p0, y=-util.inner(self.domain.getNormal(),u0))        
317          p0=self.__pde_p.getSolution()
318          g2=self.__g - util.tensor_mult(self.__permeability, util.grad(p0))
319          norm_g2=util.integrate(util.inner(g2,util.tensor_mult(self.__permeability_inv,g2)))**0.5    
320    
321          p=p0*0
322          if self.solveForFlux:
323         v=u0.copy()
324          else:
325         v=self.__getFlux(p, u0, f=self.__f, g=g2)
326    
327          while not converged and norm_g2 > 0:
328         Gp=util.grad(p)
329         KGp=util.tensor_mult(self.__permeability,Gp)
330         if self.verbose:
331            def_p=g2-(v+KGp)
332            def_v=self.__f-util.div(v)
333            print "DarcyFlux: L2: g-v-K*grad(p) = %e (v = %e)."%(self.__L2(def_p),self.__L2(v))
334            print "DarcyFlux: L2: f-div(v) = %e (grad(v) = %e)."%(self.__L2(def_v),self.__L2(util.grad(v)))
335            print "DarcyFlux: K^{-1}-norm of v = %e."%util.integrate(util.inner(v,util.tensor_mult(self.__permeability_inv,v)))**0.5
336            print "DarcyFlux: K^{-1}-norm of g2 = %e."%norm_g2
337            print "DarcyFlux: K-norm of grad(dp) = %e."%util.integrate(util.inner(Gp,KGp))**0.5
338         ATOL=atol+rtol*norm_g2
339         if self.verbose: print "DarcyFlux: absolute tolerance ATOL = %e."%(ATOL,)
340         if norm_r == None or norm_r>ATOL:
341            if num_corrections>max_num_corrections:
342               raise ValueError,"maximum number of correction steps reached."
343          
344            if self.solveForFlux:
345               # initial residual is r=K^{-1}*(g-v-K*Gp)+D^*L^{-1}(f-Du)
346               v,r, norm_r=PCG(ArithmeticTuple(util.tensor_mult(self.__permeability_inv,g2-v)-Gp,self.__applWeight(v,self.__f),p),
347                   self.__Aprod_v,
348                   v,
349                   self.__Msolve_PCG_v,
350                   self.__inner_PCG_v,
351                   atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
352               p=r[2]
353            else:
354               # initial residual is r=G^*(g2-KGp - v)
355               p,r, norm_r=PCG(ArithmeticTuple(g2-KGp,v),
356                     self.__Aprod_p,
357                     p,
358                     self.__Msolve_PCG_p,
359                     self.__inner_PCG_p,
360                     atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
361               v=r[1]
362            if self.verbose: print "DarcyFlux: residual norm = %e."%norm_r
363            num_corrections+=1
364         else:
365            if self.verbose: print "DarcyFlux: stopping criterium reached."
366            converged=True
367          return v,p+p0
368    
369       def __applWeight(self, v, f=None):
370          # solves L p = f-Dv with p = 0
371          if self.verbose: print "DarcyFlux: Applying weighting operator"
372          if f == None:
373         return -util.div(v)*self.__l2
374          else:
375         return (f-util.div(v))*self.__l2
376       def __getPressure(self, v, p0, g=None):
377          # solves (G*KG)p = G^(g-v) with p = p0 where location_of_fixed_pressure>0
378          if self.getSolverOptionsPressure().isVerbose() or self.verbose: print "DarcyFlux: Pressure update"
379          if g == None:
380         self.__pde_p.setValue(X=-v, r=p0)
381          else:
382         self.__pde_p.setValue(X=g-v, r=p0)      
383          p=self.__pde_p.getSolution()
384          return p
385    
386       def __Aprod_v(self,dv):
387          # calculates: (a,b,c) = (K^{-1}(dv + KG * dp), L^{-1}Ddv, dp)  with (G*KG)dp = - G^*dv  
388          dp=self.__getPressure(dv, p0=escript.Data()) # dp = (G*KG)^{-1} (0-G^*dv)
389          a=util.tensor_mult(self.__permeability_inv,dv)+util.grad(dp) # a= K^{-1}u+G*dp
390          b= - self.__applWeight(dv) # b = - (D K D^*)^{-1} (0-Dv)
391          return ArithmeticTuple(a,b,-dp)
392    
393       def __Msolve_PCG_v(self,r):
394          # K^{-1} u = r[0] + D^*r[1] = K^{-1}(dv + KG * dp) + D^*L^{-1}Ddv
395          if self.getSolverOptionsFlux().isVerbose() or self.verbose: print "DarcyFlux: Applying preconditioner"
396          self.__pde_k.setValue(X=r[1]*util.kronecker(self.domain), Y=r[0], r=escript.Data())
397          return self.__pde_k.getSolution()
398    
399       def __inner_PCG_v(self,v,r):
400          return util.integrate(util.inner(v,r[0])+util.div(v)*r[1])
401          
402       def __Aprod_p(self,dp):
403          if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
404          Gdp=util.grad(dp)
405          self.__pde_k.setValue(Y=-Gdp,X=escript.Data(), r=escript.Data())
406          du=self.__pde_k.getSolution()
407          return ArithmeticTuple(util.tensor_mult(self.__permeability,Gdp),-du)
408    
409       def __getFlux(self,p, v0, f=None, g=None):
410          # solves (K^{-1}+D^*L^{-1} D) v = D^*L^{-1}f + K^{-1}g - Gp
411          if f!=None:
412         self.__pde_k.setValue(X=self.__applWeight(v0*0,self.__f)*util.kronecker(self.domain))
413          self.__pde_k.setValue(r=v0)
414          g2=util.tensor_mult(self.__permeability_inv,g)
415          if p == None:
416         self.__pde_k.setValue(Y=g2)
417          else:
418         self.__pde_k.setValue(Y=g2-util.grad(p))
419          return self.__pde_k.getSolution()  
420          
421          #v=self.__getFlux(p, u0, f=self.__f, g=g2)      
422       def __Msolve_PCG_p(self,r):
423          if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"
424          self.__pde_p.setValue(X=r[0]-r[1], Y=escript.Data(), r=escript.Data(), y=escript.Data())
425          return self.__pde_p.getSolution()
426            
427       def __inner_PCG_p(self,p,r):
428           return util.integrate(util.inner(util.grad(p), r[0]-r[1]))
429    
430       def __L2(self,v):
431          return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))
432    
433       def __L2_r(self,v):
434          return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.ReducedFunction(self.domain)))**2))
435    
436       def setTolerance(self,rtol=1e-4):
437          """
438          sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if
439    
440          *|g-v-K gard(p)|_PCG <= atol + rtol * |K^{1/2}g2|_0*
441          
442          where ``atol`` is an absolut tolerance (see `setAbsoluteTolerance`).
443          
444          :param rtol: relative tolerance for the pressure
445          :type rtol: non-negative ``float``
446          """
447          if rtol<0:
448         raise ValueError,"Relative tolerance needs to be non-negative."
449          self.__rtol=rtol
450       def getTolerance(self):
451          """
452          returns the relative tolerance
453          :return: current relative tolerance
454          :rtype: ``float``
455          """
456          return self.__rtol
457    
458       def setAbsoluteTolerance(self,atol=0.):
459          """
460          sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if
461          
462          *|g-v-K gard(p)|_PCG <= atol + rtol * |K^{1/2}g2|_0*
463    
464    
465          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}*.
466    
467          :param atol: absolute tolerance for the pressure
468          :type atol: non-negative ``float``
469          """
470          if atol<0:
471         raise ValueError,"Absolute tolerance needs to be non-negative."
472          self.__atol=atol
473       def getAbsoluteTolerance(self):
474          """
475          returns the absolute tolerance
476          :return: current absolute tolerance
477          :rtype: ``float``
478          """
479          return self.__atol
480       def getSubProblemTolerance(self):
481          """
482          Returns a suitable subtolerance
483          :type: ``float``
484          """
485          return max(util.EPSILON**(0.5),self.getTolerance()**2)
486    
487       def setSubProblemTolerance(self):
488          """
489          Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected.
490          """
491          if self.__adaptSubTolerance:
492         sub_tol=self.getSubProblemTolerance()
493         self.getSolverOptionsFlux().setTolerance(sub_tol)
494         self.getSolverOptionsFlux().setAbsoluteTolerance(0.)
495         self.getSolverOptionsPressure().setTolerance(sub_tol)
496         self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
497         if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol
498    
499    
500    class DarcyFlowOld(object):
501      """      """
502      solves the problem      solves the problem
503    
# Line 61  class DarcyFlow(object): Line 522  class DarcyFlow(object):
522          self.domain=domain          self.domain=domain
523          if weight == None:          if weight == None:
524             s=self.domain.getSize()             s=self.domain.getSize()
525             self.__l=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2             self.__l2=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2
526             # self.__l=(3.*util.longestEdge(self.domain))**2             # self.__l2=(3.*util.longestEdge(self.domain))**2
527             #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
528          else:          else:
529             self.__l=weight             self.__l2=weight
530          self.__pde_v=LinearPDESystem(domain)          self.__pde_v=LinearPDESystem(domain)
531          if useReduced: self.__pde_v.setReducedOrderOn()          if useReduced: self.__pde_v.setReducedOrderOn()
532          self.__pde_v.setSymmetryOn()          self.__pde_v.setSymmetryOn()
533          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)))
534          self.__pde_p=LinearSinglePDE(domain)          self.__pde_p=LinearSinglePDE(domain)
535          self.__pde_p.setSymmetryOn()          self.__pde_p.setSymmetryOn()
536          if useReduced: self.__pde_p.setReducedOrderOn()          if useReduced: self.__pde_p.setReducedOrderOn()
537          self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))          self.__f=escript.Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
538          self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))          self.__g=escript.Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
539          self.setTolerance()          self.setTolerance()
540          self.setAbsoluteTolerance()          self.setAbsoluteTolerance()
541      self.__adaptSubTolerance=adaptSubTolerance      self.__adaptSubTolerance=adaptSubTolerance
# Line 125  class DarcyFlow(object): Line 586  class DarcyFlow(object):
586          assigns values to model parameters          assigns values to model parameters
587    
588          :param f: volumetic sources/sinks          :param f: volumetic sources/sinks
589          :type f: scalar value on the domain (e.g. `Data`)          :type f: scalar value on the domain (e.g. `escript.Data`)
590          :param g: flux sources/sinks          :param g: flux sources/sinks
591          :type g: vector values on the domain (e.g. `Data`)          :type g: vector values on the domain (e.g. `escript.Data`)
592          :param location_of_fixed_pressure: mask for locations where pressure is fixed          :param location_of_fixed_pressure: mask for locations where pressure is fixed
593          :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`)
594          :param location_of_fixed_flux:  mask for locations where flux is fixed.          :param location_of_fixed_flux:  mask for locations where flux is fixed.
595          :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`)
596          :param permeability: permeability tensor. If scalar ``s`` is given the tensor with          :param permeability: permeability tensor. If scalar ``s`` is given the tensor with
597                               ``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
598                               ``v`` on the main diagonal is used.                               ``v`` on the main diagonal is used.
599          :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`)
600    
601          :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.
602          :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 145  class DarcyFlow(object): Line 606  class DarcyFlow(object):
606          if f !=None:          if f !=None:
607             f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))             f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
608             if f.isEmpty():             if f.isEmpty():
609                 f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))                 f=escript.Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
610             else:             else:
611                 if f.getRank()>0: raise ValueError,"illegal rank of f."                 if f.getRank()>0: raise ValueError,"illegal rank of f."
612             self.__f=f             self.__f=f
613          if g !=None:          if g !=None:
614             g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))             g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
615             if g.isEmpty():             if g.isEmpty():
616               g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))               g=escript.Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
617             else:             else:
618               if not g.getShape()==(self.domain.getDim(),):               if not g.getShape()==(self.domain.getDim(),):
619                 raise ValueError,"illegal shape of g"                 raise ValueError,"illegal shape of g"
# Line 245  class DarcyFlow(object): Line 706  class DarcyFlow(object):
706           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().
707    
708           :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.
709           :type u0: vector value on the domain (e.g. `Data`).           :type u0: vector value on the domain (e.g. `escript.Data`).
710           :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.
711           :type p0: scalar value on the domain (e.g. `Data`).           :type p0: scalar value on the domain (e.g. `escript.Data`).
712           :param verbose: if set some information on iteration progress are printed           :param verbose: if set some information on iteration progress are printed
713           :type verbose: ``bool``           :type verbose: ``bool``
714           :return: flux and pressure           :return: flux and pressure
715           :rtype: ``tuple`` of `Data`.           :rtype: ``tuple`` of `escript.Data`.
716    
717           :note: The problem is solved as a least squares form           :note: The problem is solved as a least squares form
718    
# Line 296  class DarcyFlow(object): Line 757  class DarcyFlow(object):
757                 ATOL=(atol+rtol*fac)                 ATOL=(atol+rtol*fac)
758                 if self.verbose:                 if self.verbose:
759                      print "DarcyFlux: L2 norm of v = %e."%norm_v                      print "DarcyFlux: L2 norm of v = %e."%norm_v
760                      print "DarcyFlux: L2 norm of k.grad(p) = %e."%norm_Qp                      print "DarcyFlux: L2 norm of k.util.grad(p) = %e."%norm_Qp
761                      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),)
762                      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),)
763                      print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL                      print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
764                 if norm_r == None or norm_r>ATOL:                 if norm_r == None or norm_r>ATOL:
765                     if num_corrections>max_num_corrections:                     if num_corrections>max_num_corrections:
766                           raise ValueError,"maximum number of correction steps reached."                           raise ValueError,"maximum number of correction steps reached."
767                     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)
768                     num_corrections+=1                     num_corrections+=1
769                 else:                 else:
770                     converged=True                     converged=True
771           return v,p           return v,p
772      def __L2(self,v):      def __L2(self,v):
773           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))
774    
775      def __Q(self,p):      def __Q(self,p):
776            return util.tensor_mult(self.__permeability,util.grad(p))            return util.tensor_mult(self.__permeability,util.grad(p))
# Line 317  class DarcyFlow(object): Line 778  class DarcyFlow(object):
778      def __Aprod(self,dp):      def __Aprod(self,dp):
779            if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"            if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
780            Qdp=self.__Q(dp)            Qdp=self.__Q(dp)
781            self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())            self.__pde_v.setValue(Y=-Qdp,X=escript.Data(), r=escript.Data())
782            du=self.__pde_v.getSolution()            du=self.__pde_v.getSolution()
           # self.__pde_v.getOperator().saveMM("proj.mm")  
783            return Qdp+du            return Qdp+du
784      def __inner_GMRES(self,r,s):      def __inner_GMRES(self,r,s):
785           return util.integrate(util.inner(r,s))           return util.integrate(util.inner(r,s))
# Line 329  class DarcyFlow(object): Line 789  class DarcyFlow(object):
789    
790      def __Msolve_PCG(self,r):      def __Msolve_PCG(self,r):
791        if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"        if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"
792            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")  
793            return self.__pde_p.getSolution()            return self.__pde_p.getSolution()
794    
795      def getFlux(self,p=None, fixed_flux=Data()):      def getFlux(self,p=None, fixed_flux=escript.Data()):
796          """          """
797          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``
798          on locations where ``location_of_fixed_flux`` is positive (see `setValue`).          on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
799          Note that ``g`` and ``f`` are used, see `setValue`.          Note that ``g`` and ``f`` are used, see `setValue`.
800    
801          :param p: pressure.          :param p: pressure.
802          :type p: scalar value on the domain (e.g. `Data`).          :type p: scalar value on the domain (e.g. `escript.Data`).
803          :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``.
804          :type fixed_flux: vector values on the domain (e.g. `Data`).          :type fixed_flux: vector values on the domain (e.g. `escript.Data`).
805          :return: flux          :return: flux
806          :rtype: `Data`          :rtype: `escript.Data`
807          :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}*
808                 for the permeability *k_{ij}*                 for the permeability *k_{ij}*
809          """          """
810      self.setSubProblemTolerance()      self.setSubProblemTolerance()
811          g=self.__g          g=self.__g
812          f=self.__f          f=self.__f
813          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)
814          if p == None:          if p == None:
815             self.__pde_v.setValue(Y=g)             self.__pde_v.setValue(Y=g)
816          else:          else:
# Line 473  class StokesProblemCartesian(Homogeneous Line 932  class StokesProblemCartesian(Homogeneous
932          if eta !=None:          if eta !=None:
933              k=util.kronecker(self.domain.getDim())              k=util.kronecker(self.domain.getDim())
934              kk=util.outer(k,k)              kk=util.outer(k,k)
935              self.eta=util.interpolate(eta, Function(self.domain))              self.eta=util.interpolate(eta, escript.Function(self.domain))
936          self.__pde_prec.setValue(D=1/self.eta)          self.__pde_prec.setValue(D=1/self.eta)
937              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)))
938          if restoration_factor!=None:          if restoration_factor!=None:
# Line 485  class StokesProblemCartesian(Homogeneous Line 944  class StokesProblemCartesian(Homogeneous
944          if surface_stress!=None: self.__surface_stress=surface_stress          if surface_stress!=None: self.__surface_stress=surface_stress
945          if stress!=None: self.__stress=stress          if stress!=None: self.__stress=stress
946    
947       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):
948          """          """
949          assigns values to the model parameters          assigns values to the model parameters
950    
# Line 525  class StokesProblemCartesian(Homogeneous Line 984  class StokesProblemCartesian(Homogeneous
984           :return: inner product of element p and Bv=-div(v)           :return: inner product of element p and Bv=-div(v)
985           :rtype: ``float``           :rtype: ``float``
986           """           """
987           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)))
988    
989       def inner_p(self,p0,p1):       def inner_p(self,p0,p1):
990           """           """
# Line 536  class StokesProblemCartesian(Homogeneous Line 995  class StokesProblemCartesian(Homogeneous
995           :return: inner product of p0 and p1           :return: inner product of p0 and p1
996           :rtype: ``float``           :rtype: ``float``
997           """           """
998           s0=util.interpolate(p0,Function(self.domain))           s0=util.interpolate(p0, escript.Function(self.domain))
999           s1=util.interpolate(p1,Function(self.domain))           s1=util.interpolate(p1, escript.Function(self.domain))
1000           return util.integrate(s0*s1)           return util.integrate(s0*s1)
1001    
1002       def norm_v(self,v):       def norm_v(self,v):
# Line 577  class StokesProblemCartesian(Homogeneous Line 1036  class StokesProblemCartesian(Homogeneous
1036          :rtype: equal to the type of p          :rtype: equal to the type of p
1037          :note: boundary conditions on p should be zero!          :note: boundary conditions on p should be zero!
1038          """          """
1039          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))
1040    
1041       def solve_AinvBt(self,p, tol):       def solve_AinvBt(self,p, tol):
1042           """           """
# Line 587  class StokesProblemCartesian(Homogeneous Line 1046  class StokesProblemCartesian(Homogeneous
1046           :return: the solution of *Av=B^*p*           :return: the solution of *Av=B^*p*
1047           :note: boundary conditions on v should be zero!           :note: boundary conditions on v should be zero!
1048           """           """
1049           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))
1050           out=self.__pde_u.getSolution()           out=self.__pde_u.getSolution()
1051           return  out           return  out
1052    

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

  ViewVC Help
Powered by ViewVC 1.1.26