/[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 3499 by gross, Fri Apr 8 00:14:35 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    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

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

  ViewVC Help
Powered by ViewVC 1.1.26