/[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 2959 by jfenwick, Thu Jan 28 02:03:15 2010 UTC revision 2960 by gross, Tue Mar 2 07:54:11 2010 UTC
# Line 53  class DarcyFlow(object): Line 53  class DarcyFlow(object):
53          initializes the Darcy flux problem          initializes the Darcy flux problem
54          :param domain: domain of the problem          :param domain: domain of the problem
55          :type domain: `Domain`          :type domain: `Domain`
56            :param weight: the weighting factor for the incempressiblity equation *u_{i,i} = f* in the variational principle.
57            If not present an apprppriate weight is chosen.
58            :type  weight: positive sacalar
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            """
64            self.domain=domain
65            if weight == None:
66               self.__l=None
67               self.__update_weight=True
68            else:
69               self.__update_weight=False
70               self.__l=weight
71            self.__pde_v=LinearPDESystem(domain)
72            if useReduced: self.__pde_v.setReducedOrderOn()
73            self.__pde_v.setSymmetryOn()
74            self.__pde_p=LinearSinglePDE(domain)
75            self.__pde_p.setSymmetryOn()
76            if useReduced: self.__pde_p.setReducedOrderOn()
77            self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
78            self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
79            self.setTolerance()
80            self.setAbsoluteTolerance()
81        self.__adaptSubTolerance=adaptSubTolerance
82        self.verbose=False
83        def getSolverOptionsFlux(self):
84        """
85        Returns the solver options used to solve the flux problems
86        
87        *(I+D^* weight D)u=F*
88        
89        :return: `SolverOptions`
90        """
91        return self.__pde_v.getSolverOptions()
92    
93        def setSolverOptionsFlux(self, options=None):
94        """
95        Sets the solver options used to solve the flux problems
96        
97        *(I+D^*  weight  D)u=F*
98        
99        If ``options`` is not present, the options are reset to default
100        :param options: `SolverOptions`
101        :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
102        """
103        return self.__pde_v.setSolverOptions(options)
104        def getSolverOptionsPressure(self):
105        """
106        Returns the solver options used to solve the pressure problems
107        
108        *(Q^* S Q)p=Q^*G*
109        
110             as a preconditioner where S=weight*k^2/dx^2
111    
112        :return: `SolverOptions`
113        """
114        return self.__pde_p.getSolverOptions()
115        def setSolverOptionsPressure(self, options=None):
116        """
117        Sets the solver options used to solve the pressure problems
118        
119        *(Q^* S Q)p=Q^*G*
120        
121            was a preconditioner here S=weight*k^2/dx^2
122        
123        If ``options`` is not present, the options are reset to default
124    
125        :param options: `SolverOptions`
126        :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
127        """
128        return self.__pde_p.setSolverOptions(options)
129    
130        def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
131            """
132            assigns values to model parameters
133    
134            :param f: volumetic sources/sinks
135            :type f: scalar value on the domain (e.g. `Data`)
136            :param g: flux sources/sinks
137            :type g: vector values on the domain (e.g. `Data`)
138            :param location_of_fixed_pressure: mask for locations where pressure is fixed
139            :type location_of_fixed_pressure: scalar value on the domain (e.g. `Data`)
140            :param location_of_fixed_flux:  mask for locations where flux is fixed.
141            :type location_of_fixed_flux: vector values on the domain (e.g. `Data`)
142            :param permeability: permeability tensor. If scalar ``s`` is given the tensor with
143                                 ``s`` on the main diagonal is used. If vector ``v`` is given the tensor with
144                                 ``v`` on the main diagonal is used.
145            :type permeability: scalar, vector or tensor values on the domain (e.g. `Data`)
146    
147            :note: the values of parameters which are not set by calling ``setValue`` are not altered.
148            :note: at any point on the boundary of the domain the pressure (``location_of_fixed_pressure`` >0)
149                   or the normal component of the flux (``location_of_fixed_flux[i]>0`` if direction of the normal
150                   is along the *x_i* axis.
151            """
152            if f !=None:
153               f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
154               if f.isEmpty():
155                   f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
156               else:
157                   if f.getRank()>0: raise ValueError,"illegal rank of f."
158               self.__f=f
159            if g !=None:
160               g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
161               if g.isEmpty():
162                 g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
163               else:
164                 if not g.getShape()==(self.domain.getDim(),):
165                   raise ValueError,"illegal shape of g"
166               self.__g=g
167    
168            if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)
169            if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux)
170    
171            if permeability!=None:
172               perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
173               if perm.getRank()==0:
174                   perm_inv=1./perm*util.kronecker(self.domain.getDim())
175                   perm=perm*util.kronecker(self.domain.getDim())
176               elif perm.getRank()==1:
177                   perm2=perm
178                   perm=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A"))
179                   perm_inv=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A"))
180                   for i in range(self.domain.getDim()):
181                          p=perm2[i]
182                          perm[i,i]=p
183                          perm_inv[i,i]=1/p
184               elif perm.getRank()==2:
185                   perm_inv=util.inverse(perm)
186               else:
187                  raise ValueError,"illegal rank of permeability."
188               self.__permeability=perm
189               self.__permeability_inv=perm_inv
190               if self.__update_weight: self.__l = util.longestEdge(self.domain)**2/3.14**2*util.length(self.__permeability_inv)
191               self.__pde_v.setValue(D=self.__permeability_inv, A=self.__l*util.outer(util.kronecker(self.domain),util.kronecker(self.domain)))
192               s=self.__pde_p.getFunctionSpaceForCoefficient("A").getSize()
193               self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability)*self.__l/s**2)
194    
195        def getFlux(self,p=None, fixed_flux=Data()):
196            """
197            returns the flux for a given pressure ``p`` where the flux is equal to ``fixed_flux``
198            on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
199            Note that ``g`` and ``f`` are used, see `setValue`.
200    
201            :param p: pressure.
202            :type p: scalar value on the domain (e.g. `Data`).
203            :param fixed_flux: flux on the locations of the domain marked be ``location_of_fixed_flux``.
204            :type fixed_flux: vector values on the domain (e.g. `Data`).
205            :return: flux
206            :rtype: `Data`
207            :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}*
208                   for the permeability *k_{ij}*
209            """
210        self.setSubProblemTolerance()
211            self.__pde_v.setValue(X=util.sqrt(self.__l)*self.__f*util.kronecker(self.domain), r=fixed_flux)
212            g2=util.tensor_mult(self.__permeability_inv,self.__g)
213            if p == None:
214               self.__pde_v.setValue(Y=g2)
215            else:
216               self.__pde_v.setValue(Y=g2-util.grad(p))
217            return self.__pde_v.getSolution()
218    
219        def __Aprod(self,dp):
220              if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
221              Gdp=util.grad(dp)
222              self.__pde_v.setValue(Y=-Gdp,X=Data(), r=Data())
223              du=self.__pde_v.getSolution()
224              # self.__pde_v.getOperator().saveMM("proj.mm")
225              return ArithmeticTuple(util.tensor_mult(self.__permeability,Gdp),-du)
226    
227        def __Msolve_PCG(self,r):
228          if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"
229              self.__pde_p.setValue(X=r[0]-r[1], Y=Data(), r=Data())
230              # self.__pde_p.getOperator().saveMM("prec.mm")
231              return self.__pde_p.getSolution()
232    
233        def __inner_PCG(self,p,r):
234             return util.integrate(util.inner(util.grad(p), r[0]-r[1]))
235    
236        def __L2(self,v):
237             return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))
238    
239        def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):
240             """
241             solves the problem.
242    
243             The iteration is terminated if the residual norm is less then self.getTolerance().
244    
245             :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.
246             :type u0: vector value on the domain (e.g. `Data`).
247             :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.
248             :type p0: scalar value on the domain (e.g. `Data`).
249             :param verbose: if set some information on iteration progress are printed
250             :type verbose: ``bool``
251             :return: flux and pressure
252             :rtype: ``tuple`` of `Data`.
253    
254             :note: The problem is solved as a least squares form
255    
256             *(K^[-1]+D^* weight D)u+G p=D^* sqrt(weight) f + K^[-1]g*
257             *G^*u+G^* K Gp=G^*g*
258    
259             where *D* is the *div* operator and *(Gp)_i=p_{,i}* for the permeability *K=k_{ij}*.
260             We eliminate the flux form the problem by setting
261    
262             *u=(K^[-1]+D^* weight D)^{-1}(D^* sqrt(weight) f + K^[-1]g* -G p ) with u=u0 on location_of_fixed_flux
263    
264             form the first equation. Inserted into the second equation we get
265    
266             *G^*(K-(K^[-1]+D^* weight D)^{-1})Gp= G^*(g-(K^[-1]+D^* weight D)^{-1}(D^* sqrt(weight) f + K^[-1]g))* with p=p0  on location_of_fixed_pressure
267    
268             which is solved using the PCG method (precondition is *G^*K^2*weight*l^2/dx^2 G*. In each iteration step
269             PDEs with operator *K^[-1]+D^* weight D* and with *G^*K^2*weight*l^2/dx^2 G* needs to be solved using a sub iteration scheme.
270             """
271             self.verbose=verbose
272             rtol=self.getTolerance()
273             atol=self.getAbsoluteTolerance()
274         self.setSubProblemTolerance()
275             num_corrections=0
276             converged=False
277             p=p0
278             norm_r=None
279             v=self.getFlux(p, fixed_flux=u0)
280             while not converged:
281                   Gp=util.grad(p)
282                   KGp=util.tensor_mult(self.__permeability,Gp)
283                   norm_v=util.integrate(util.inner(v,util.tensor_mult(self.__permeability_inv,v)))**0.5
284                   norm_Gp=util.integrate(util.inner(Gp,KGp))**0.5
285                   if self.verbose:
286                       print "DarcyFlux: L2: g-v-K*grad(p) = %e."%(util.integrate(util.length(self.__g-util.interpolate(v,Function(self.domain))-KGp)**2)**(0.5),)
287                       print "DarcyFlux: L2: f-div(v) = %e."%(util.integrate((self.__f-util.div(v))**2)**(0.5),)
288                       print "DarcyFlux: K^{-1}-norm of v = %e."%norm_v
289                       print "DarcyFlux: K-norm of grad(p) = %e."%norm_Gp
290                   if norm_v == 0.:
291                      if norm_Gp == 0.:
292                         return v,p
293                      else:
294                        fac=norm_Gp
295                   else:
296                      if norm_Gp == 0.:
297                        fac=norm_v
298                      else:
299                        fac=2./(1./norm_v+1./norm_Gp)
300                   #fac=max(norm_Gp, norm_v)
301                   ATOL=atol+rtol*fac
302                   if self.verbose: print "DarcyFlux: absolute tolerance ATOL = %e (fac= %e)."%(ATOL,fac)
303                   if norm_r == None or norm_r>ATOL:
304                       if num_corrections>max_num_corrections:
305                             raise ValueError,"maximum number of correction steps reached."
306                       # initial residual is r=G^*(self.__g-KGp - v)
307                       p,r, norm_r=PCG(ArithmeticTuple(self.__g-KGp,v),
308                                       self.__Aprod,
309                                       p,
310                                       self.__Msolve_PCG,
311                                       self.__inner_PCG,
312                                       atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
313                       if self.verbose: print "DarcyFlux: residual norm = %e."%norm_r
314                       v=r[1]
315                       num_corrections+=1
316                   else:
317                       if self.verbose: print "DarcyFlux: stopping criterium reached."
318                       converged=True
319             return v,p
320    
321              
322    
323        def setTolerance(self,rtol=1e-4):
324            """
325            sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if
326    
327            *|g-v-K gard(p)|_PCG <= atol + rtol * (1/|K^{-1/2}u|_0 +1/|K^{1/2}grad(p)|_0)^(-1}*
328    
329            where ``atol`` is an absolut tolerance (see `setAbsoluteTolerance`).
330    
331            :param rtol: relative tolerance for the pressure
332            :type rtol: non-negative ``float``
333            """
334            if rtol<0:
335                raise ValueError,"Relative tolerance needs to be non-negative."
336            self.__rtol=rtol
337        def getTolerance(self):
338            """
339            returns the relative tolerance
340    
341            :return: current relative tolerance
342            :rtype: ``float``
343            """
344            return self.__rtol
345    
346        def setAbsoluteTolerance(self,atol=0.):
347            """
348            sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if
349    
350             *|g-v-K gard(p)|_PCG <= atol + rtol * (1/|K^{-1/2}u|_0 +1/|K^{1/2}grad(p)|_0)^(-1}*
351    
352            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}*.
353    
354            :param atol: absolute tolerance for the pressure
355            :type atol: non-negative ``float``
356            """
357            if atol<0:
358                raise ValueError,"Absolute tolerance needs to be non-negative."
359            self.__atol=atol
360        def getAbsoluteTolerance(self):
361           """
362           returns the absolute tolerance
363          
364           :return: current absolute tolerance
365           :rtype: ``float``
366           """
367           return self.__atol
368        def getSubProblemTolerance(self):
369        """
370        Returns a suitable subtolerance
371        @type: ``float``
372        """
373        return max(util.EPSILON**(0.75),self.getTolerance()**2)
374    
375        def setSubProblemTolerance(self):
376             """
377             Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected.
378             """
379         if self.__adaptSubTolerance:
380             sub_tol=self.getSubProblemTolerance()
381                 self.getSolverOptionsFlux().setTolerance(sub_tol)
382             self.getSolverOptionsFlux().setAbsoluteTolerance(0.)
383             self.getSolverOptionsPressure().setTolerance(sub_tol)
384             self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
385             if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol
386    
387    
388    class DarcyFlowOld(object):
389        """
390        solves the problem
391    
392        *u_i+k_{ij}*p_{,j} = g_i*
393        *u_{i,i} = f*
394    
395        where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,
396    
397        :note: The problem is solved in a least squares formulation.
398        """
399    
400        def __init__(self, domain, weight=None, useReduced=False, adaptSubTolerance=True):
401            """
402            initializes the Darcy flux problem
403            :param domain: domain of the problem
404            :type domain: `Domain`
405      :param useReduced: uses reduced oreder on flux and pressure      :param useReduced: uses reduced oreder on flux and pressure
406      :type useReduced: ``bool``      :type useReduced: ``bool``
407      :param adaptSubTolerance: switches on automatic subtolerance selection      :param adaptSubTolerance: switches on automatic subtolerance selection
# Line 296  class DarcyFlow(object): Line 645  class DarcyFlow(object):
645                 ATOL=(atol+rtol*fac)                 ATOL=(atol+rtol*fac)
646                 if self.verbose:                 if self.verbose:
647                      print "DarcyFlux: L2 norm of v = %e."%norm_v                      print "DarcyFlux: L2 norm of v = %e."%norm_v
648                      print "DarcyFlux: L2 norm of k.grad(p) = %e."%norm_Qp                      print "DarcyFlux: L2 norm of k.util.grad(p) = %e."%norm_Qp
649                      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,Function(self.domain))-Qp)**2)**(0.5),)
650                      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),)
651                      print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL                      print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL

Legend:
Removed from v.2959  
changed lines
  Added in v.2960

  ViewVC Help
Powered by ViewVC 1.1.26