/[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 3500 by gross, Mon Apr 11 03:09:06 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, useVPIteration=True, *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          :param useVPIteration: if True altenative iteration over v and p is performed. Otherwise V and P are calculated in a single PDE.
61          :type useVPIteration: ``bool``    
62          """
63          self.domain=domain
64          self.useVPIteration=useVPIteration
65          self.useReduced=useReduced
66          self.verbose=False
67    
68          if self.useVPIteration:
69         # this does not work yet
70             self.__pde_k=LinearPDESystem(domain)
71             self.__pde_k.setSymmetryOn()
72             if self.useReduced: self.__pde_k.setReducedOrderOn()
73      
74             self.__pde_p=LinearSinglePDE(domain)
75             self.__pde_p.setSymmetryOn()
76             if self.useReduced: self.__pde_p.setReducedOrderOn()
77          else:
78             self.__pde_k=LinearPDE(self.domain, numEquations=self.domain.getDim()+1)
79             self.__pde_k.setSymmetryOff()
80            
81             if self.useReduced: self.__pde_k.setReducedOrderOn()
82             C=self.__pde_k.createCoefficient("C")
83             B=self.__pde_k.createCoefficient("B")
84             for i in range(self.domain.getDim()):
85                B[i,i,self.domain.getDim()]=-1
86                C[self.domain.getDim(),i,i]=1
87                C[i,self.domain.getDim(),i]=-0.5
88                B[self.domain.getDim(),i,i]=0.5
89             self.__pde_k.setValue(C=C, B=B)
90          self.__f=escript.Scalar(0,self.__pde_k.getFunctionSpaceForCoefficient("X"))
91          self.__g=escript.Vector(0,self.__pde_k.getFunctionSpaceForCoefficient("Y"))
92          self.location_of_fixed_pressure = escript.Scalar(0, self.__pde_k.getFunctionSpaceForCoefficient("q"))
93          self.location_of_fixed_flux = escript.Vector(0, self.__pde_k.getFunctionSpaceForCoefficient("q"))
94          self.scale=1.
95          self.setTolerance()
96          
97    
98       def __L2(self,v):
99             return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))  
100       def __inner_GMRES(self,r,s):
101             return util.integrate(util.inner(r,s))
102            
103       def __Aprod_GMRES(self,p):
104          self.__pde_k.setValue(Y=-0.5*util.grad(p), X=-p*util.kronecker(self.__pde_k.getDomain()) )
105          du=self.__pde_k.getSolution()
106          self.__pde_p.setValue(Y=-util.div(du), X=0.5*(-du+util.tensor_mult(self.__permeability,util.grad(p))))
107          return self.__pde_p.getSolution()
108            
109       def getSolverOptionsFlux(self):
110          """
111          Returns the solver options used to solve the flux problems
112          
113          *K^{-1} u=F*
114          
115          :return: `SolverOptions`
116          """
117          return self.__pde_k.getSolverOptions()      
118          
119       def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
120          """
121          assigns values to model parameters
122    
123          :param f: volumetic sources/sinks
124          :type f: scalar value on the domain (e.g. `escript.Data`)
125          :param g: flux sources/sinks
126          :type g: vector values on the domain (e.g. `escript.Data`)
127          :param location_of_fixed_pressure: mask for locations where pressure is fixed
128          :type location_of_fixed_pressure: scalar value on the domain (e.g. `escript.Data`)
129          :param location_of_fixed_flux:  mask for locations where flux is fixed.
130          :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)
131          :param permeability: permeability tensor. If scalar ``s`` is given the tensor with ``s`` on the main diagonal is used.
132          :type permeability: scalar or tensor values on the domain (e.g. `escript.Data`)
133    
134          :note: the values of parameters which are not set by calling ``setValue`` are not altered.
135          :note: at any point on the boundary of the domain the pressure
136                 (``location_of_fixed_pressure`` >0) or the normal component of the
137                 flux (``location_of_fixed_flux[i]>0``) if direction of the normal
138                 is along the *x_i* axis.
139    
140          """
141          if location_of_fixed_pressure!=None: self.location_of_fixed_pressure=location_of_fixed_pressure
142          if location_of_fixed_flux!=None: self.location_of_fixed_flux=location_of_fixed_flux
143          
144          if self.useVPIteration:
145             if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=self.location_of_fixed_pressure)
146             if location_of_fixed_flux!=None: self.__pde_k.setValue(q=self.location_of_fixed_flux)
147          else:
148         if location_of_fixed_pressure!=None or location_of_fixed_flux!=None:
149            q=self.__pde_k.createCoefficient("q")
150            q[self.domain.getDim()]=self.location_of_fixed_pressure
151            q[:self.domain.getDim()]=self.location_of_fixed_flux
152            self.__pde_k.setValue(q=q)
153                
154          # flux is rescaled by the factor mean value(perm_inv)*length where length**self.domain.getDim()=vol(self.domain)
155          if permeability!=None:
156         perm=util.interpolate(permeability,self.__pde_k.getFunctionSpaceForCoefficient("A"))
157             V=util.vol(self.domain)
158         if perm.getRank()==0:
159            perm_inv=(1./perm)
160                self.scale=util.integrate(perm_inv)*V**(1./self.domain.getDim()-1.)
161            perm_inv=perm_inv*((1./self.scale)*util.kronecker(self.domain.getDim()))
162            perm=perm*(self.scale*util.kronecker(self.domain.getDim()))
163         elif perm.getRank()==2:
164            perm_inv=util.inverse(perm)
165                self.scale=util.sqrt(util.integrate(util.length(perm_inv)**2)*V**(2./self.domain.getDim()-1.)/self.domain.getDim())
166            perm_inv*=(1./self.scale)
167            perm=perm*self.scale
168         else:
169            raise ValueError,"illegal rank of permeability."
170            
171         self.__permeability=perm
172         self.__permeability_inv=perm_inv
173         if self.useVPIteration:
174                self.__pde_k.setValue(D=0.5*self.__permeability_inv)
175            self.__pde_p.setValue(A=0.5*self.__permeability)
176             else:
177                D=self.__pde_k.createCoefficient("D")
178                A=self.__pde_k.createCoefficient("A")
179                D[:self.domain.getDim(),:self.domain.getDim()]=0.5*self.__permeability_inv
180                A[self.domain.getDim(),:,self.domain.getDim(),:]=0.5*self.__permeability
181                self.__pde_k.setValue(A=A, D=D)
182          if g != None:
183        g=util.interpolate(g, self.__pde_k.getFunctionSpaceForCoefficient("Y"))
184        if g.isEmpty():
185              g=Vector(0,self.__pde_k.getFunctionSpaceForCoefficient("Y"))
186        else:
187            if not g.getShape()==(self.domain.getDim(),): raise ValueError,"illegal shape of g"
188            self.__g=g
189          if f !=None:
190         f=util.interpolate(f, self.__pde_k.getFunctionSpaceForCoefficient("X"))
191         if f.isEmpty():
192          
193              f=Scalar(0,self.__pde_k.getFunctionSpaceForCoefficient("X"))
194         else:
195             if f.getRank()>0: raise ValueError,"illegal rank of f."
196             self.__f=f
197            
198       def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10, iter_restart=20):
199          """
200          solves the problem.
201          
202          The iteration is terminated if the residual norm is less then self.getTolerance().
203    
204          :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.
205          :type u0: vector value on the domain (e.g. `escript.Data`).
206          :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.
207          :type p0: scalar value on the domain (e.g. `escript.Data`).
208          :param verbose: if set some information on iteration progress are printed
209          :type verbose: ``bool``
210          :return: flux and pressure
211          :rtype: ``tuple`` of `escript.Data`.
212    
213          """
214          u0_b=u0*self.location_of_fixed_flux
215          p0_b=p0*self.location_of_fixed_pressure/self.scale
216          f=self.__f-util.div(u0_b)
217          g=self.__g-u0_b - util.tensor_mult(self.__permeability,util.grad(p0_b))
218          self.verbose=verbose
219          if self.useVPIteration:
220            # get u:
221                rtol=self.getTolerance()
222            self.__pde_k.setValue(Y=0.5*util.tensor_mult(self.__permeability_inv,g),X=escript.Data())
223            du=self.__pde_k.getSolution()
224            self.__pde_p.setValue(Y=f-util.div(du), X=0.5*(g-du))
225            p=GMRES(self.__pde_p.getSolution(),
226                    self.__Aprod_GMRES,
227                p0_b*0,
228                self.__inner_GMRES,
229                atol=0,
230                rtol=rtol,
231                iter_max=max_iter,
232                iter_restart=iter_restart, verbose=self.verbose,P_R=None)
233                self.__pde_k.setValue(Y=0.5*( util.tensor_mult(self.__permeability_inv,g) + util.grad(p)) ,
234                                      X=p*util.kronecker(self.__pde_k.getDomain()))
235                #self.__pde_k.setValue(Y=0.5*util.grad(p), 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.3500

  ViewVC Help
Powered by ViewVC 1.1.26