/[escript]/trunk/escriptcore/py_src/flows.py
ViewVC logotype

Annotation of /trunk/escriptcore/py_src/flows.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2960 - (hide annotations)
Tue Mar 2 07:54:11 2010 UTC (9 years, 5 months ago) by gross
Original Path: trunk/escript/py_src/flows.py
File MIME type: text/x-python
File size: 39649 byte(s)
new DarcyFlux solver
1 ksteube 1809 ########################################################
2 gross 1414 #
3 jfenwick 2881 # Copyright (c) 2003-2010 by University of Queensland
4 ksteube 1809 # Earth Systems Science Computational Center (ESSCC)
5     # http://www.uq.edu.au/esscc
6 gross 1414 #
7 ksteube 1809 # Primary Business: Queensland, Australia
8     # Licensed under the Open Software License version 3.0
9     # http://www.opensource.org/licenses/osl-3.0.php
10 gross 1414 #
11 ksteube 1809 ########################################################
12 gross 1414
13 jfenwick 2881 __copyright__="""Copyright (c) 2003-2010 by University of Queensland
14 ksteube 1809 Earth Systems Science Computational Center (ESSCC)
15     http://www.uq.edu.au/esscc
16     Primary Business: Queensland, Australia"""
17     __license__="""Licensed under the Open Software License version 3.0
18     http://www.opensource.org/licenses/osl-3.0.php"""
19 jfenwick 2344 __url__="https://launchpad.net/escript-finley"
20 ksteube 1809
21 gross 1414 """
22     Some models for flow
23    
24 jfenwick 2625 :var __author__: name of author
25     :var __copyright__: copyrights
26     :var __license__: licence agreement
27     :var __url__: url entry point on documentation
28     :var __version__: version
29     :var __date__: date of the version
30 gross 1414 """
31    
32     __author__="Lutz Gross, l.gross@uq.edu.au"
33    
34     from escript import *
35     import util
36 gross 2474 from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions
37 gross 2264 from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
38 gross 1414
39 gross 2100 class DarcyFlow(object):
40     """
41 gross 2264 solves the problem
42 gross 1659
43 jfenwick 2625 *u_i+k_{ij}*p_{,j} = g_i*
44     *u_{i,i} = f*
45 gross 1659
46 jfenwick 2625 where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,
47 gross 1659
48 jfenwick 2625 :note: The problem is solved in a least squares formulation.
49 gross 2100 """
50 gross 1659
51 gross 2474 def __init__(self, domain, weight=None, useReduced=False, adaptSubTolerance=True):
52 gross 2100 """
53 gross 2208 initializes the Darcy flux problem
54 jfenwick 2625 :param domain: domain of the problem
55     :type domain: `Domain`
56 gross 2960 :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 jfenwick 2625 :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 gross 2100 """
64     self.domain=domain
65 gross 2351 if weight == None:
66 gross 2960 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
406     :type useReduced: ``bool``
407     :param adaptSubTolerance: switches on automatic subtolerance selection
408     :type adaptSubTolerance: ``bool``
409     """
410     self.domain=domain
411     if weight == None:
412 gross 2354 s=self.domain.getSize()
413     self.__l=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2
414 gross 2486 # self.__l=(3.*util.longestEdge(self.domain))**2
415 artak 2689 #self.__l=(0.1*util.longestEdge(self.domain)*s/util.sup(s))**2
416 gross 2351 else:
417     self.__l=weight
418 gross 2100 self.__pde_v=LinearPDESystem(domain)
419 gross 2208 if useReduced: self.__pde_v.setReducedOrderOn()
420     self.__pde_v.setSymmetryOn()
421 gross 2267 self.__pde_v.setValue(D=util.kronecker(domain), A=self.__l*util.outer(util.kronecker(domain),util.kronecker(domain)))
422 gross 2100 self.__pde_p=LinearSinglePDE(domain)
423     self.__pde_p.setSymmetryOn()
424 gross 2208 if useReduced: self.__pde_p.setReducedOrderOn()
425 gross 2100 self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
426     self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
427 gross 2264 self.setTolerance()
428     self.setAbsoluteTolerance()
429 gross 2474 self.__adaptSubTolerance=adaptSubTolerance
430     self.verbose=False
431     def getSolverOptionsFlux(self):
432     """
433     Returns the solver options used to solve the flux problems
434    
435 jfenwick 2625 *(I+D^*D)u=F*
436 gross 2474
437 jfenwick 2625 :return: `SolverOptions`
438 gross 2474 """
439     return self.__pde_v.getSolverOptions()
440     def setSolverOptionsFlux(self, options=None):
441     """
442     Sets the solver options used to solve the flux problems
443    
444 jfenwick 2625 *(I+D^*D)u=F*
445 gross 2474
446 jfenwick 2625 If ``options`` is not present, the options are reset to default
447     :param options: `SolverOptions`
448     :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
449 gross 2474 """
450     return self.__pde_v.setSolverOptions(options)
451     def getSolverOptionsPressure(self):
452     """
453     Returns the solver options used to solve the pressure problems
454    
455 jfenwick 2625 *(Q^*Q)p=Q^*G*
456 gross 2474
457 jfenwick 2625 :return: `SolverOptions`
458 gross 2474 """
459     return self.__pde_p.getSolverOptions()
460     def setSolverOptionsPressure(self, options=None):
461     """
462     Sets the solver options used to solve the pressure problems
463    
464 jfenwick 2625 *(Q^*Q)p=Q^*G*
465 gross 2474
466 jfenwick 2625 If ``options`` is not present, the options are reset to default
467     :param options: `SolverOptions`
468     :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
469 gross 2474 """
470     return self.__pde_p.setSolverOptions(options)
471 gross 1659
472 gross 2100 def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
473     """
474 gross 2208 assigns values to model parameters
475 gross 1659
476 jfenwick 2625 :param f: volumetic sources/sinks
477     :type f: scalar value on the domain (e.g. `Data`)
478     :param g: flux sources/sinks
479     :type g: vector values on the domain (e.g. `Data`)
480     :param location_of_fixed_pressure: mask for locations where pressure is fixed
481     :type location_of_fixed_pressure: scalar value on the domain (e.g. `Data`)
482     :param location_of_fixed_flux: mask for locations where flux is fixed.
483     :type location_of_fixed_flux: vector values on the domain (e.g. `Data`)
484     :param permeability: permeability tensor. If scalar ``s`` is given the tensor with
485     ``s`` on the main diagonal is used. If vector ``v`` is given the tensor with
486     ``v`` on the main diagonal is used.
487     :type permeability: scalar, vector or tensor values on the domain (e.g. `Data`)
488 gross 1659
489 jfenwick 2625 :note: the values of parameters which are not set by calling ``setValue`` are not altered.
490     :note: at any point on the boundary of the domain the pressure (``location_of_fixed_pressure`` >0)
491     or the normal component of the flux (``location_of_fixed_flux[i]>0`` if direction of the normal
492     is along the *x_i* axis.
493 gross 2100 """
494 gross 2264 if f !=None:
495 gross 2100 f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
496     if f.isEmpty():
497     f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
498     else:
499     if f.getRank()>0: raise ValueError,"illegal rank of f."
500 gross 2267 self.__f=f
501 gross 2264 if g !=None:
502 gross 2100 g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
503     if g.isEmpty():
504     g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
505     else:
506     if not g.getShape()==(self.domain.getDim(),):
507     raise ValueError,"illegal shape of g"
508     self.__g=g
509 gross 1659
510 gross 2100 if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)
511     if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux)
512 gross 1659
513 gross 2100 if permeability!=None:
514     perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
515     if perm.getRank()==0:
516     perm=perm*util.kronecker(self.domain.getDim())
517     elif perm.getRank()==1:
518     perm, perm2=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A")), perm
519     for i in range(self.domain.getDim()): perm[i,i]=perm2[i]
520     elif perm.getRank()==2:
521     pass
522     else:
523     raise ValueError,"illegal rank of permeability."
524     self.__permeability=perm
525     self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))
526 gross 1659
527 gross 2264 def setTolerance(self,rtol=1e-4):
528     """
529 jfenwick 2625 sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if
530 gross 1659
531 jfenwick 2625 *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
532 gross 2264
533 jfenwick 2625 where ``atol`` is an absolut tolerance (see `setAbsoluteTolerance`), *|f|^2 = integrate(length(f)^2)* and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.
534 gross 2264
535 jfenwick 2625 :param rtol: relative tolerance for the pressure
536     :type rtol: non-negative ``float``
537 gross 2100 """
538 gross 2264 if rtol<0:
539     raise ValueError,"Relative tolerance needs to be non-negative."
540     self.__rtol=rtol
541     def getTolerance(self):
542 gross 2100 """
543 gross 2264 returns the relative tolerance
544 gross 1659
545 jfenwick 2625 :return: current relative tolerance
546     :rtype: ``float``
547 caltinay 2169 """
548 gross 2264 return self.__rtol
549 gross 1659
550 gross 2264 def setAbsoluteTolerance(self,atol=0.):
551 gross 2208 """
552 jfenwick 2625 sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if
553 gross 1659
554 jfenwick 2625 *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
555 gross 2208
556 jfenwick 2625 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}*.
557 gross 2208
558 jfenwick 2625 :param atol: absolute tolerance for the pressure
559     :type atol: non-negative ``float``
560 gross 2264 """
561     if atol<0:
562     raise ValueError,"Absolute tolerance needs to be non-negative."
563     self.__atol=atol
564     def getAbsoluteTolerance(self):
565     """
566     returns the absolute tolerance
567    
568 jfenwick 2625 :return: current absolute tolerance
569     :rtype: ``float``
570 gross 2264 """
571     return self.__atol
572     def getSubProblemTolerance(self):
573 gross 2474 """
574     Returns a suitable subtolerance
575 jfenwick 2625 @type: ``float``
576 gross 2474 """
577     return max(util.EPSILON**(0.75),self.getTolerance()**2)
578     def setSubProblemTolerance(self):
579 gross 2264 """
580 gross 2474 Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected.
581 gross 2264 """
582 gross 2474 if self.__adaptSubTolerance:
583     sub_tol=self.getSubProblemTolerance()
584     self.getSolverOptionsFlux().setTolerance(sub_tol)
585     self.getSolverOptionsFlux().setAbsoluteTolerance(0.)
586     self.getSolverOptionsPressure().setTolerance(sub_tol)
587     self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
588     if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol
589 gross 2208
590 gross 2474 def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):
591 gross 2264 """
592     solves the problem.
593 gross 1659
594 gross 2208 The iteration is terminated if the residual norm is less then self.getTolerance().
595 gross 1659
596 jfenwick 2625 :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.
597     :type u0: vector value on the domain (e.g. `Data`).
598     :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.
599     :type p0: scalar value on the domain (e.g. `Data`).
600     :param verbose: if set some information on iteration progress are printed
601     :type verbose: ``bool``
602     :return: flux and pressure
603     :rtype: ``tuple`` of `Data`.
604 gross 2264
605 jfenwick 2625 :note: The problem is solved as a least squares form
606 gross 2100
607 jfenwick 2625 *(I+D^*D)u+Qp=D^*f+g*
608     *Q^*u+Q^*Qp=Q^*g*
609 gross 2100
610 jfenwick 2625 where *D* is the *div* operator and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.
611 gross 2264 We eliminate the flux form the problem by setting
612 caltinay 2169
613 jfenwick 2625 *u=(I+D^*D)^{-1}(D^*f-g-Qp)* with u=u0 on location_of_fixed_flux
614 caltinay 2169
615 gross 2208 form the first equation. Inserted into the second equation we get
616 caltinay 2169
617 jfenwick 2625 *Q^*(I-(I+D^*D)^{-1})Qp= Q^*(g-(I+D^*D)^{-1}(D^*f+g))* with p=p0 on location_of_fixed_pressure
618 gross 2264
619 jfenwick 2625 which is solved using the PCG method (precondition is *Q^*Q*). In each iteration step
620     PDEs with operator *I+D^*D* and with *Q^*Q* needs to be solved using a sub iteration scheme.
621 gross 2208 """
622 gross 2386 self.verbose=verbose
623 gross 2264 rtol=self.getTolerance()
624     atol=self.getAbsoluteTolerance()
625 gross 2474 self.setSubProblemTolerance()
626 gross 2264 num_corrections=0
627     converged=False
628     p=p0
629     norm_r=None
630     while not converged:
631 gross 2474 v=self.getFlux(p, fixed_flux=u0)
632 gross 2264 Qp=self.__Q(p)
633     norm_v=self.__L2(v)
634     norm_Qp=self.__L2(Qp)
635     if norm_v == 0.:
636     if norm_Qp == 0.:
637     return v,p
638     else:
639     fac=norm_Qp
640     else:
641     if norm_Qp == 0.:
642     fac=norm_v
643     else:
644     fac=2./(1./norm_v+1./norm_Qp)
645     ATOL=(atol+rtol*fac)
646     if self.verbose:
647     print "DarcyFlux: L2 norm of v = %e."%norm_v
648 gross 2960 print "DarcyFlux: L2 norm of k.util.grad(p) = %e."%norm_Qp
649 gross 2486 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),)
651 gross 2264 print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
652     if norm_r == None or norm_r>ATOL:
653     if num_corrections>max_num_corrections:
654     raise ValueError,"maximum number of correction steps reached."
655 gross 2354 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)
656 gross 2264 num_corrections+=1
657     else:
658     converged=True
659     return v,p
660     def __L2(self,v):
661     return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))
662    
663     def __Q(self,p):
664     return util.tensor_mult(self.__permeability,util.grad(p))
665    
666     def __Aprod(self,dp):
667 gross 2474 if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
668 gross 2264 Qdp=self.__Q(dp)
669     self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())
670 gross 2474 du=self.__pde_v.getSolution()
671 gross 2370 # self.__pde_v.getOperator().saveMM("proj.mm")
672 gross 2264 return Qdp+du
673     def __inner_GMRES(self,r,s):
674     return util.integrate(util.inner(r,s))
675    
676 gross 2100 def __inner_PCG(self,p,r):
677 gross 2264 return util.integrate(util.inner(self.__Q(p), r))
678 gross 2100
679     def __Msolve_PCG(self,r):
680 gross 2474 if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"
681 gross 2264 self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=Data(), r=Data())
682 gross 2370 # self.__pde_p.getOperator().saveMM("prec.mm")
683 gross 2474 return self.__pde_p.getSolution()
684 gross 2100
685 gross 2474 def getFlux(self,p=None, fixed_flux=Data()):
686 gross 2264 """
687 jfenwick 2625 returns the flux for a given pressure ``p`` where the flux is equal to ``fixed_flux``
688     on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
689     Note that ``g`` and ``f`` are used, see `setValue`.
690 gross 2264
691 jfenwick 2625 :param p: pressure.
692     :type p: scalar value on the domain (e.g. `Data`).
693     :param fixed_flux: flux on the locations of the domain marked be ``location_of_fixed_flux``.
694     :type fixed_flux: vector values on the domain (e.g. `Data`).
695     :return: flux
696     :rtype: `Data`
697     :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}*
698     for the permeability *k_{ij}*
699 gross 2264 """
700 gross 2474 self.setSubProblemTolerance()
701 gross 2264 g=self.__g
702     f=self.__f
703 gross 2267 self.__pde_v.setValue(X=self.__l*f*util.kronecker(self.domain), r=fixed_flux)
704 gross 2264 if p == None:
705     self.__pde_v.setValue(Y=g)
706     else:
707     self.__pde_v.setValue(Y=g-self.__Q(p))
708 gross 2474 return self.__pde_v.getSolution()
709 gross 2264
710 gross 1414 class StokesProblemCartesian(HomogeneousSaddlePointProblem):
711 gross 2251 """
712 gross 2264 solves
713 gross 1414
714 gross 2208 -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
715     u_{i,i}=0
716 gross 1414
717 gross 2208 u=0 where fixed_u_mask>0
718     eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j
719 gross 1414
720 gross 2264 if surface_stress is not given 0 is assumed.
721 gross 1414
722 gross 2251 typical usage:
723 gross 1414
724 gross 2208 sp=StokesProblemCartesian(domain)
725     sp.setTolerance()
726     sp.initialize(...)
727     v,p=sp.solve(v0,p0)
728 gross 2251 """
729 gross 2719 def __init__(self,domain,**kwargs):
730 gross 2100 """
731 gross 2208 initialize the Stokes Problem
732 gross 2100
733 gross 2793 The approximation spaces used for velocity (=Solution(domain)) and pressure (=ReducedSolution(domain)) must be
734     LBB complient, for instance using quadratic and linear approximation on the same element or using linear approximation
735     with macro elements for the pressure.
736    
737     :param domain: domain of the problem.
738 jfenwick 2625 :type domain: `Domain`
739 gross 2100 """
740 gross 2719 HomogeneousSaddlePointProblem.__init__(self,**kwargs)
741 gross 1414 self.domain=domain
742     self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
743     self.__pde_u.setSymmetryOn()
744 gross 2474
745 gross 1414 self.__pde_prec=LinearPDE(domain)
746     self.__pde_prec.setReducedOrderOn()
747     self.__pde_prec.setSymmetryOn()
748    
749 gross 2415 self.__pde_proj=LinearPDE(domain)
750     self.__pde_proj.setReducedOrderOn()
751     self.__pde_proj.setValue(D=1)
752     self.__pde_proj.setSymmetryOn()
753    
754 gross 2474 def getSolverOptionsVelocity(self):
755     """
756     returns the solver options used solve the equation for velocity.
757    
758 jfenwick 2625 :rtype: `SolverOptions`
759 gross 2474 """
760     return self.__pde_u.getSolverOptions()
761     def setSolverOptionsVelocity(self, options=None):
762     """
763     set the solver options for solving the equation for velocity.
764    
765 jfenwick 2625 :param options: new solver options
766     :type options: `SolverOptions`
767 gross 2474 """
768     self.__pde_u.setSolverOptions(options)
769     def getSolverOptionsPressure(self):
770     """
771     returns the solver options used solve the equation for pressure.
772 jfenwick 2625 :rtype: `SolverOptions`
773 gross 2474 """
774     return self.__pde_prec.getSolverOptions()
775     def setSolverOptionsPressure(self, options=None):
776     """
777     set the solver options for solving the equation for pressure.
778 jfenwick 2625 :param options: new solver options
779     :type options: `SolverOptions`
780 gross 2474 """
781     self.__pde_prec.setSolverOptions(options)
782 gross 2415
783 gross 2474 def setSolverOptionsDiv(self, options=None):
784     """
785     set the solver options for solving the equation to project the divergence of
786     the velocity onto the function space of presure.
787    
788 jfenwick 2625 :param options: new solver options
789     :type options: `SolverOptions`
790 gross 2474 """
791 artak 2689 self.__pde_proj.setSolverOptions(options)
792 gross 2474 def getSolverOptionsDiv(self):
793     """
794     returns the solver options for solving the equation to project the divergence of
795     the velocity onto the function space of presure.
796    
797 jfenwick 2625 :rtype: `SolverOptions`
798 gross 2474 """
799 artak 2689 return self.__pde_proj.getSolverOptions()
800 gross 2474
801 gross 2793 def updateStokesEquation(self, v, p):
802     """
803     updates the Stokes equation to consider dependencies from ``v`` and ``p``
804     :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values.
805     """
806     pass
807     def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None):
808     """
809     assigns new values to the model parameters.
810    
811     :param f: external force
812     :type f: `Vector` object in `FunctionSpace` `Function` or similar
813     :param fixed_u_mask: mask of locations with fixed velocity.
814     :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
815     :param eta: viscosity
816     :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
817     :param surface_stress: normal surface stress
818     :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
819     :param stress: initial stress
820     :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
821     """
822     if eta !=None:
823     k=util.kronecker(self.domain.getDim())
824     kk=util.outer(k,k)
825     self.eta=util.interpolate(eta, Function(self.domain))
826     self.__pde_prec.setValue(D=1/self.eta)
827     self.__pde_u.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))
828     if restoration_factor!=None:
829     n=self.domain.getNormal()
830     self.__pde_u.setValue(d=restoration_factor*util.outer(n,n))
831     if fixed_u_mask!=None:
832     self.__pde_u.setValue(q=fixed_u_mask)
833     if f!=None: self.__f=f
834     if surface_stress!=None: self.__surface_stress=surface_stress
835     if stress!=None: self.__stress=stress
836    
837 gross 2620 def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data(), restoration_factor=0):
838 gross 2208 """
839     assigns values to the model parameters
840 gross 2100
841 jfenwick 2625 :param f: external force
842     :type f: `Vector` object in `FunctionSpace` `Function` or similar
843     :param fixed_u_mask: mask of locations with fixed velocity.
844     :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
845     :param eta: viscosity
846     :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
847     :param surface_stress: normal surface stress
848     :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
849     :param stress: initial stress
850     :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
851 gross 2208 """
852 gross 2793 self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
853 gross 1414
854 gross 2719 def Bv(self,v,tol):
855 gross 2251 """
856     returns inner product of element p and div(v)
857 gross 1414
858 jfenwick 2625 :param v: a residual
859     :return: inner product of element p and div(v)
860     :rtype: ``float``
861 gross 2100 """
862 gross 2793 self.__pde_proj.setValue(Y=-util.div(v))
863 gross 2719 self.getSolverOptionsDiv().setTolerance(tol)
864     self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
865     out=self.__pde_proj.getSolution()
866     return out
867 gross 2208
868 gross 2445 def inner_pBv(self,p,Bv):
869     """
870     returns inner product of element p and Bv=-div(v)
871    
872 jfenwick 2625 :param p: a pressure increment
873     :param Bv: a residual
874     :return: inner product of element p and Bv=-div(v)
875     :rtype: ``float``
876 gross 2445 """
877     return util.integrate(util.interpolate(p,Function(self.domain))*util.interpolate(Bv,Function(self.domain)))
878    
879 gross 2251 def inner_p(self,p0,p1):
880 gross 2100 """
881 gross 2251 Returns inner product of p0 and p1
882 gross 1414
883 jfenwick 2625 :param p0: a pressure
884     :param p1: a pressure
885     :return: inner product of p0 and p1
886     :rtype: ``float``
887 gross 2100 """
888 gross 2719 s0=util.interpolate(p0,Function(self.domain))
889     s1=util.interpolate(p1,Function(self.domain))
890 gross 2100 return util.integrate(s0*s1)
891 artak 1550
892 gross 2251 def norm_v(self,v):
893 gross 2100 """
894 gross 2251 returns the norm of v
895 gross 2208
896 jfenwick 2625 :param v: a velovity
897     :return: norm of v
898     :rtype: non-negative ``float``
899 gross 2100 """
900 gross 2719 return util.sqrt(util.integrate(util.length(util.grad(v))**2))
901 gross 2100
902 gross 2793
903 gross 2719 def getDV(self, p, v, tol):
904 gross 1414 """
905 gross 2251 return the value for v for a given p (overwrite)
906    
907 jfenwick 2625 :param p: a pressure
908 gross 2719 :param v: a initial guess for the value v to return.
909     :return: dv given as *Adv=(f-Av-B^*p)*
910 gross 1414 """
911 gross 2793 self.updateStokesEquation(v,p)
912 gross 2719 self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress)
913     self.getSolverOptionsVelocity().setTolerance(tol)
914     self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
915 gross 2100 if self.__stress.isEmpty():
916 gross 2719 self.__pde_u.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
917 gross 2100 else:
918 gross 2719 self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
919 gross 2474 out=self.__pde_u.getSolution()
920 gross 2208 return out
921 gross 1414
922 gross 2445 def norm_Bv(self,Bv):
923 gross 2251 """
924     Returns Bv (overwrite).
925    
926 jfenwick 2625 :rtype: equal to the type of p
927     :note: boundary conditions on p should be zero!
928 gross 2251 """
929 gross 2445 return util.sqrt(util.integrate(util.interpolate(Bv,Function(self.domain))**2))
930 gross 2251
931 gross 2719 def solve_AinvBt(self,p, tol):
932 gross 2251 """
933 gross 2719 Solves *Av=B^*p* with accuracy `tol`
934 gross 2251
935 jfenwick 2625 :param p: a pressure increment
936     :return: the solution of *Av=B^*p*
937     :note: boundary conditions on v should be zero!
938 gross 2251 """
939 gross 2719 self.__pde_u.setValue(Y=Data(), y=Data(), X=-p*util.kronecker(self.domain))
940 gross 2474 out=self.__pde_u.getSolution()
941 gross 2251 return out
942    
943 gross 2719 def solve_prec(self,Bv, tol):
944 gross 2251 """
945 jfenwick 2625 applies preconditioner for for *BA^{-1}B^** to *Bv*
946     with accuracy `self.getSubProblemTolerance()`
947 gross 2251
948 jfenwick 2625 :param Bv: velocity increment
949     :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*
950     :note: boundary conditions on p are zero.
951 gross 2251 """
952 gross 2445 self.__pde_prec.setValue(Y=Bv)
953 gross 2719 self.getSolverOptionsPressure().setTolerance(tol)
954     self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
955     out=self.__pde_prec.getSolution()
956     return out

  ViewVC Help
Powered by ViewVC 1.1.26