/[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 3452 by caltinay, Tue Jan 25 01:53:57 2011 UTC revision 3501 by gross, Wed Apr 13 04:07:53 2011 UTC
# Line 37  import util Line 37  import util
37  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions
38  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
39    
 print dir(escript)  
40    
41    
42    class DarcyFlowNew(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 stabelized formulation.
52       """
53       def __init__(self, domain, useReduced=False, useVPIteration=True, *args, **kargs):
54          """
55          initializes the Darcy flux problem
56          :param domain: domain of the problem
57          :type domain: `Domain`
58          :param useReduced: uses reduced oreder on flux and pressure
59          :type useReduced: ``bool``
60          :param adaptSubTolerance: switches on automatic subtolerance selection
61          :param useVPIteration: if True altenative iteration over v and p is performed. Otherwise V and P are calculated in a single PDE.
62          :type useVPIteration: ``bool``    
63          """
64          self.domain=domain
65          self.useVPIteration=useVPIteration
66          self.useReduced=useReduced
67          self.verbose=False
68    
69          if self.useVPIteration:
70         # this does not work yet
71             self.__pde_k=LinearPDESystem(domain)
72             self.__pde_k.setSymmetryOn()
73             if self.useReduced: self.__pde_k.setReducedOrderOn()
74      
75             self.__pde_p=LinearSinglePDE(domain)
76             self.__pde_p.setSymmetryOn()
77             if self.useReduced: self.__pde_p.setReducedOrderOn()
78          else:
79             self.__pde_k=LinearPDE(self.domain, numEquations=self.domain.getDim()+1)
80             self.__pde_k.setSymmetryOff()
81            
82             if self.useReduced: self.__pde_k.setReducedOrderOn()
83             C=self.__pde_k.createCoefficient("C")
84             B=self.__pde_k.createCoefficient("B")
85             for i in range(self.domain.getDim()):
86                B[i,i,self.domain.getDim()]=-1
87                C[self.domain.getDim(),i,i]=1
88                C[i,self.domain.getDim(),i]=-0.5
89                B[self.domain.getDim(),i,i]=0.5
90             self.__pde_k.setValue(C=C, B=B)
91          self.__f=escript.Scalar(0,self.__pde_k.getFunctionSpaceForCoefficient("X"))
92          self.__g=escript.Vector(0,self.__pde_k.getFunctionSpaceForCoefficient("Y"))
93          self.location_of_fixed_pressure = escript.Scalar(0, self.__pde_k.getFunctionSpaceForCoefficient("q"))
94          self.location_of_fixed_flux = escript.Vector(0, self.__pde_k.getFunctionSpaceForCoefficient("q"))
95          self.scale=1.
96          self.setTolerance()
97          
98    
99       def __L2(self,v):
100             return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))  
101       def __inner_GMRES(self,r,s):
102             return util.integrate(util.inner(r,s))
103            
104       def __Aprod_GMRES(self,p):
105          self.__pde_k.setValue(Y=-0.5*util.grad(p), X=-p*util.kronecker(self.__pde_k.getDomain()) )
106          du=self.__pde_k.getSolution()
107          self.__pde_p.setValue(Y=-util.div(du), X=0.5*(-du+util.tensor_mult(self.__permeability,util.grad(p))))
108          return self.__pde_p.getSolution()
109            
110       def getSolverOptionsFlux(self):
111          """
112          Returns the solver options used to solve the flux problems
113          
114          *K^{-1} u=F*
115          
116          :return: `SolverOptions`
117          """
118          return self.__pde_k.getSolverOptions()      
119          
120       def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
121          """
122          assigns values to model parameters
123    
124          :param f: volumetic sources/sinks
125          :type f: scalar value on the domain (e.g. `escript.Data`)
126          :param g: flux sources/sinks
127          :type g: vector values on the domain (e.g. `escript.Data`)
128          :param location_of_fixed_pressure: mask for locations where pressure is fixed
129          :type location_of_fixed_pressure: scalar value on the domain (e.g. `escript.Data`)
130          :param location_of_fixed_flux:  mask for locations where flux is fixed.
131          :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)
132          :param permeability: permeability tensor. If scalar ``s`` is given the tensor with ``s`` on the main diagonal is used.
133          :type permeability: scalar or tensor values on the domain (e.g. `escript.Data`)
134    
135          :note: the values of parameters which are not set by calling ``setValue`` are not altered.
136          :note: at any point on the boundary of the domain the pressure
137                 (``location_of_fixed_pressure`` >0) or the normal component of the
138                 flux (``location_of_fixed_flux[i]>0``) if direction of the normal
139                 is along the *x_i* axis.
140    
141          """
142          if location_of_fixed_pressure!=None: self.location_of_fixed_pressure=util.wherePositive(location_of_fixed_pressure)
143          if location_of_fixed_flux!=None: self.location_of_fixed_flux=util.wherePositive(location_of_fixed_flux)
144          
145          if self.useVPIteration:
146             if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=self.location_of_fixed_pressure)
147             if location_of_fixed_flux!=None: self.__pde_k.setValue(q=self.location_of_fixed_flux)
148          else:
149         if location_of_fixed_pressure!=None or location_of_fixed_flux!=None:
150            q=self.__pde_k.createCoefficient("q")
151            q[self.domain.getDim()]=self.location_of_fixed_pressure
152            q[:self.domain.getDim()]=self.location_of_fixed_flux
153            self.__pde_k.setValue(q=q)
154                
155          # flux is rescaled by the factor mean value(perm_inv)*length where length**self.domain.getDim()=vol(self.domain)
156          if permeability!=None:
157         perm=util.interpolate(permeability,self.__pde_k.getFunctionSpaceForCoefficient("A"))
158             V=util.vol(self.domain)
159         if perm.getRank()==0:
160            perm_inv=(1./perm)
161                self.scale=util.integrate(perm_inv)*V**(1./self.domain.getDim()-1.)
162            perm_inv=perm_inv*((1./self.scale)*util.kronecker(self.domain.getDim()))
163            perm=perm*(self.scale*util.kronecker(self.domain.getDim()))
164         elif perm.getRank()==2:
165            perm_inv=util.inverse(perm)
166                self.scale=util.sqrt(util.integrate(util.length(perm_inv)**2)*V**(2./self.domain.getDim()-1.)/self.domain.getDim())
167            perm_inv*=(1./self.scale)
168            perm=perm*self.scale
169         else:
170            raise ValueError,"illegal rank of permeability."
171            
172         self.__permeability=perm
173         self.__permeability_inv=perm_inv
174         if self.useVPIteration:
175                self.__pde_k.setValue(D=0.5*self.__permeability_inv)
176            self.__pde_p.setValue(A=0.5*self.__permeability)
177             else:
178                D=self.__pde_k.createCoefficient("D")
179                A=self.__pde_k.createCoefficient("A")
180                D[:self.domain.getDim(),:self.domain.getDim()]=0.5*self.__permeability_inv
181                A[self.domain.getDim(),:,self.domain.getDim(),:]=0.5*self.__permeability
182                self.__pde_k.setValue(A=A, D=D)
183          if g != None:
184        g=util.interpolate(g, self.__pde_k.getFunctionSpaceForCoefficient("Y"))
185        if g.isEmpty():
186              g=Vector(0,self.__pde_k.getFunctionSpaceForCoefficient("Y"))
187        else:
188            if not g.getShape()==(self.domain.getDim(),): raise ValueError,"illegal shape of g"
189            self.__g=g
190          if f !=None:
191         f=util.interpolate(f, self.__pde_k.getFunctionSpaceForCoefficient("X"))
192         if f.isEmpty():
193          
194              f=Scalar(0,self.__pde_k.getFunctionSpaceForCoefficient("X"))
195         else:
196             if f.getRank()>0: raise ValueError,"illegal rank of f."
197             self.__f=f
198            
199       def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10, iter_restart=20):
200          """
201          solves the problem.
202          
203          The iteration is terminated if the residual norm is less then self.getTolerance().
204    
205          :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.
206          :type u0: vector value on the domain (e.g. `escript.Data`).
207          :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.
208          :type p0: scalar value on the domain (e.g. `escript.Data`).
209          :param verbose: if set some information on iteration progress are printed
210          :type verbose: ``bool``
211          :return: flux and pressure
212          :rtype: ``tuple`` of `escript.Data`.
213    
214          """
215          u0_b=u0*self.location_of_fixed_flux
216          p0_b=p0*self.location_of_fixed_pressure/self.scale
217          f=self.__f-util.div(u0_b)
218          g=self.__g-u0_b - util.tensor_mult(self.__permeability,util.grad(p0_b))
219          self.verbose=verbose
220          if self.useVPIteration:
221            # get u:
222                rtol=self.getTolerance()
223                p_init=0*p0_b
224            self.__pde_k.setValue(Y=0.5*(util.tensor_mult(self.__permeability_inv,g)+util.grad(p_init)) ,
225                                      X=p_init*util.kronecker(self.__pde_k.getDomain()))
226            du=self.__pde_k.getSolution()
227            self.__pde_p.setValue(Y=f-util.div(du), X=0.5*(g-du))
228            p=GMRES(self.__pde_p.getSolution(),
229                    self.__Aprod_GMRES,
230                p_init,
231                self.__inner_GMRES,
232                atol=0,
233                rtol=rtol,
234                iter_max=max_iter,
235                iter_restart=iter_restart, verbose=self.verbose,P_R=None)
236                self.__pde_k.setValue(Y=0.5*( util.tensor_mult(self.__permeability_inv,g) + util.grad(p)) ,
237                                      X=p*util.kronecker(self.__pde_k.getDomain()))
238                #self.__pde_k.setValue(Y=0.5*util.grad(p), X=p*util.kronecker(self.__pde_k.getDomain()) )
239            u=self.__pde_k.getSolution()
240          else:
241              X=self.__pde_k.createCoefficient("X")
242              Y=self.__pde_k.createCoefficient("Y")
243              Y[:self.domain.getDim()]=0.5*util.tensor_mult(self.__permeability_inv,g)
244              Y[self.domain.getDim()]=f
245              X[self.domain.getDim(),:]=g*0.5
246              self.__pde_k.setValue(X=X, Y=Y)
247              self.__pde_k.getSolverOptions().setVerbosity(self.verbose)
248              #self.__pde_k.getSolverOptions().setPreconditioner(self.__pde_k.getSolverOptions().AMG)
249              self.__pde_k.getSolverOptions().setSolverMethod(self.__pde_k.getSolverOptions().DIRECT)
250              U=self.__pde_k.getSolution()
251              u=U[:self.domain.getDim()]
252              p=U[self.domain.getDim()]
253          # self.__pde_k.getOperator().saveMM("k.mm")
254          u=u0_b+u
255          p=(p0_b+p)*self.scale
256          if self.verbose:
257            KGp=util.tensor_mult(self.__permeability,util.grad(p)/self.scale)
258            def_p=self.__g-(u+KGp)
259            def_v=self.__f-util.div(u, self.__pde_k.getFunctionSpaceForCoefficient("X"))
260            print "DarcyFlux: L2: g-v-K*grad(p) = %e (v = %e)."%(self.__L2(def_p),self.__L2(u))
261            print "DarcyFlux: L2: f-div(v) = %e (grad(v) = %e)."%(self.__L2(def_v),self.__L2(util.grad(u)))
262          return u,p
263       def setTolerance(self,rtol=1e-4):
264          """
265          sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if
266    
267          *|g-v-K gard(p)|_PCG <= atol + rtol * |K^{1/2}g2|_0*
268          
269          where ``atol`` is an absolut tolerance (see `setAbsoluteTolerance`).
270          
271          :param rtol: relative tolerance for the pressure
272          :type rtol: non-negative ``float``
273          """
274          if rtol<0:
275         raise ValueError,"Relative tolerance needs to be non-negative."
276          self.__rtol=rtol
277       def getTolerance(self):
278          """
279          returns the relative tolerance
280          :return: current relative tolerance
281          :rtype: ``float``
282          """
283          return self.__rtol
284          
285  class DarcyFlow(object):  class DarcyFlow(object):
286     """     """
287     solves the problem     solves the problem
# Line 51  class DarcyFlow(object): Line 294  class DarcyFlow(object):
294     :note: The problem is solved in a least squares formulation.     :note: The problem is solved in a least squares formulation.
295     """     """
296        
297     def __init__(self, domain, useReduced=False, adaptSubTolerance=True, solveForFlux=False, useVPIteration=True, weighting_scale=0.1):     def __init__(self, domain, useReduced=False, adaptSubTolerance=True, solveForFlux=False, useVPIteration=True, weighting_scale=1.):
298        """        """
299        initializes the Darcy flux problem        initializes the Darcy flux problem
300        :param domain: domain of the problem        :param domain: domain of the problem
# Line 166  class DarcyFlow(object): Line 409  class DarcyFlow(object):
409    
410        """        """
411        if self.useVPIteration:        if self.useVPIteration:
412           if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)           if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=util.wherePositive(location_of_fixed_pressure))
413           if location_of_fixed_flux!=None: self.__pde_k.setValue(q=location_of_fixed_flux)           if location_of_fixed_flux!=None: self.__pde_k.setValue(q=util.wherePositive(location_of_fixed_flux))
414        else:        else:
415           q=self.__pde_k.getCoefficient("q")           q=self.__pde_k.getCoefficient("q")
416           if q.isEmpty(): q=self.__pde_k.createCoefficient("q")           if q.isEmpty(): q=self.__pde_k.createCoefficient("q")
417           if location_of_fixed_pressure!=None: q[self.domain.getDim()]=location_of_fixed_pressure           if location_of_fixed_pressure!=None: q[self.domain.getDim()]=util.wherePositive(location_of_fixed_pressure)
418           if location_of_fixed_flux!=None: q[:self.domain.getDim()]=location_of_fixed_flux           if location_of_fixed_flux!=None: q[:self.domain.getDim()]=util.wherePositive(location_of_fixed_flux)
419           self.__pde_k.setValue(q=q)           self.__pde_k.setValue(q=q)
420                            
421        # flux is rescaled by the factor mean value(perm_inv)*length where length**self.domain.getDim()=vol(self.domain)        # flux is rescaled by the factor mean value(perm_inv)*length where length**self.domain.getDim()=vol(self.domain)
422          V=util.vol(self.domain)
423        if permeability!=None:        if permeability!=None:
424       perm=util.interpolate(permeability,self.__pde_k.getFunctionSpaceForCoefficient("A"))       perm=util.interpolate(permeability,self.__pde_k.getFunctionSpaceForCoefficient("A"))
425           V=util.vol(self.domain)          
426       if perm.getRank()==0:       if perm.getRank()==0:
427          perm_inv=(1./perm)          perm_inv=(1./perm)
428              if self.useVPIteration:              if self.useVPIteration:
# Line 202  class DarcyFlow(object): Line 446  class DarcyFlow(object):
446       self.__permeability=perm       self.__permeability=perm
447       self.__permeability_inv=perm_inv       self.__permeability_inv=perm_inv
448    
449       self.__l2 =(util.longestEdge(self.domain)**2*util.length(self.__permeability_inv))*self.weighting_scale       self.__l2 = V**(1./self.domain.getDim())*util.length(self.__permeability_inv)*self.domain.getSize()*self.weighting_scale
450           if self.useVPIteration:           if self.useVPIteration:
451          if  self.solveForFlux:          if  self.solveForFlux:
452             self.__pde_k.setValue(D=self.__permeability_inv)             self.__pde_k.setValue(D=self.__permeability_inv)

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

  ViewVC Help
Powered by ViewVC 1.1.26