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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3501 - (hide annotations)
Wed Apr 13 04:07:53 2011 UTC (8 years, 5 months ago) by gross
File MIME type: text/x-python
File size: 55236 byte(s)


1 gross 3071 # -*- coding: utf-8 -*-
2 ksteube 1809 ########################################################
3 gross 1414 #
4 jfenwick 2881 # Copyright (c) 2003-2010 by University of Queensland
5 ksteube 1809 # Earth Systems Science Computational Center (ESSCC)
6     # http://www.uq.edu.au/esscc
7 gross 1414 #
8 ksteube 1809 # Primary Business: Queensland, Australia
9     # Licensed under the Open Software License version 3.0
10     # http://www.opensource.org/licenses/osl-3.0.php
11 gross 1414 #
12 ksteube 1809 ########################################################
13 gross 1414
14 jfenwick 2881 __copyright__="""Copyright (c) 2003-2010 by University of Queensland
15 ksteube 1809 Earth Systems Science Computational Center (ESSCC)
16     http://www.uq.edu.au/esscc
17     Primary Business: Queensland, Australia"""
18     __license__="""Licensed under the Open Software License version 3.0
19     http://www.opensource.org/licenses/osl-3.0.php"""
20 jfenwick 2344 __url__="https://launchpad.net/escript-finley"
21 ksteube 1809
22 gross 1414 """
23     Some models for flow
24    
25 jfenwick 2625 :var __author__: name of author
26     :var __copyright__: copyrights
27     :var __license__: licence agreement
28     :var __url__: url entry point on documentation
29     :var __version__: version
30     :var __date__: date of the version
31 gross 1414 """
32    
33     __author__="Lutz Gross, l.gross@uq.edu.au"
34    
35 jfenwick 3432 import escript
36 gross 1414 import util
37 gross 2474 from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions
38 gross 2264 from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
39 gross 1414
40 gross 3499
41 gross 3501
42     class DarcyFlowNew(object):
43 gross 3499 """
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 gross 3500 def __init__(self, domain, useReduced=False, useVPIteration=True, *args, **kargs):
54 gross 3499 """
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 gross 3500 self.useVPIteration=useVPIteration
66 gross 3499 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 gross 3500 self.setTolerance()
97    
98    
99 gross 3499 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 gross 3500 self.__pde_k.setValue(Y=-0.5*util.grad(p), X=-p*util.kronecker(self.__pde_k.getDomain()) )
106 gross 3499 du=self.__pde_k.getSolution()
107 gross 3500 self.__pde_p.setValue(Y=-util.div(du), X=0.5*(-du+util.tensor_mult(self.__permeability,util.grad(p))))
108 gross 3499 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 gross 3501 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 gross 3499
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 gross 3500 def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10, iter_restart=20):
200 gross 3499 """
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 gross 3500 rtol=self.getTolerance()
223 gross 3501 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 gross 3499 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 gross 3501 p_init,
231 gross 3499 self.__inner_GMRES,
232     atol=0,
233 gross 3500 rtol=rtol,
234     iter_max=max_iter,
235     iter_restart=iter_restart, verbose=self.verbose,P_R=None)
236 gross 3499 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 gross 3500 #self.__pde_k.setValue(Y=0.5*util.grad(p), X=p*util.kronecker(self.__pde_k.getDomain()) )
239 gross 3499 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 gross 3501
285 gross 2100 class DarcyFlow(object):
286 gross 3071 """
287     solves the problem
288    
289     *u_i+k_{ij}*p_{,j} = g_i*
290     *u_{i,i} = f*
291    
292     where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,
293    
294     :note: The problem is solved in a least squares formulation.
295     """
296    
297 gross 3501 def __init__(self, domain, useReduced=False, adaptSubTolerance=True, solveForFlux=False, useVPIteration=True, weighting_scale=1.):
298 gross 3071 """
299     initializes the Darcy flux problem
300     :param domain: domain of the problem
301     :type domain: `Domain`
302     :param useReduced: uses reduced oreder on flux and pressure
303     :type useReduced: ``bool``
304     :param adaptSubTolerance: switches on automatic subtolerance selection
305     :type adaptSubTolerance: ``bool``
306     :param solveForFlux: if True the solver solves for the flux (do not use!)
307     :type solveForFlux: ``bool``
308 gross 3440 :param useVPIteration: if True altenative iteration over v and p is performed. Otherwise V and P are calculated in a single PDE.
309     :type useVPIteration: ``bool``
310 gross 3071 """
311     self.domain=domain
312 gross 3440 self.useVPIteration=useVPIteration
313 gross 3071 self.useReduced=useReduced
314 gross 3440 self.weighting_scale=weighting_scale
315     if self.useVPIteration:
316     self.solveForFlux=solveForFlux
317     self.__adaptSubTolerance=adaptSubTolerance
318     self.verbose=False
319    
320     self.__pde_k=LinearPDESystem(domain)
321     self.__pde_k.setSymmetryOn()
322     if self.useReduced: self.__pde_k.setReducedOrderOn()
323    
324     self.__pde_p=LinearSinglePDE(domain)
325     self.__pde_p.setSymmetryOn()
326     if self.useReduced: self.__pde_p.setReducedOrderOn()
327     self.setTolerance()
328     self.setAbsoluteTolerance()
329     else:
330     self.__pde_k=LinearPDE(self.domain, numEquations=self.domain.getDim()+1)
331     self.__pde_k.setSymmetryOn()
332     if self.useReduced: self.__pde_k.setReducedOrderOn()
333     C=self.__pde_k.createCoefficient("C")
334     B=self.__pde_k.createCoefficient("B")
335     for i in range(self.domain.getDim()):
336     C[i,self.domain.getDim(),i]=1
337     B[self.domain.getDim(),i,i]=1
338     self.__pde_k.setValue(C=C, B=B)
339     self.__f=escript.Scalar(0,self.__pde_k.getFunctionSpaceForCoefficient("X"))
340     self.__g=escript.Vector(0,self.__pde_k.getFunctionSpaceForCoefficient("Y"))
341 gross 3071
342     def getSolverOptionsFlux(self):
343     """
344     Returns the solver options used to solve the flux problems
345    
346     *K^{-1} u=F*
347    
348     :return: `SolverOptions`
349     """
350     return self.__pde_k.getSolverOptions()
351    
352     def setSolverOptionsFlux(self, options=None):
353     """
354     Sets the solver options used to solve the flux problems
355    
356     *K^{-1}u=F*
357    
358     If ``options`` is not present, the options are reset to default
359    
360     :param options: `SolverOptions`
361     :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
362     """
363     return self.__pde_v.setSolverOptions(options)
364    
365     def getSolverOptionsPressure(self):
366     """
367     Returns the solver options used to solve the pressure problems
368    
369     *(Q^* K Q)p=-Q^*G*
370    
371     :return: `SolverOptions`
372     """
373     return self.__pde_p.getSolverOptions()
374    
375     def setSolverOptionsPressure(self, options=None):
376     """
377     Sets the solver options used to solve the pressure problems
378    
379     *(Q^* K Q)p=-Q^*G*
380    
381     If ``options`` is not present, the options are reset to default
382    
383     :param options: `SolverOptions`
384     :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
385     """
386     return self.__pde_p.setSolverOptions(options)
387 gross 1659
388    
389 gross 3071 def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
390     """
391     assigns values to model parameters
392 jfenwick 3360
393 gross 3071 :param f: volumetic sources/sinks
394 gross 3440 :type f: scalar value on the domain (e.g. `escript.Data`)
395 gross 3071 :param g: flux sources/sinks
396 gross 3440 :type g: vector values on the domain (e.g. `escript.Data`)
397 gross 3071 :param location_of_fixed_pressure: mask for locations where pressure is fixed
398 gross 3440 :type location_of_fixed_pressure: scalar value on the domain (e.g. `escript.Data`)
399 gross 3071 :param location_of_fixed_flux: mask for locations where flux is fixed.
400 gross 3440 :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)
401 jfenwick 3360 :param permeability: permeability tensor. If scalar ``s`` is given the tensor with ``s`` on the main diagonal is used.
402 gross 3440 :type permeability: scalar or tensor values on the domain (e.g. `escript.Data`)
403 jfenwick 3360
404 gross 3071 :note: the values of parameters which are not set by calling ``setValue`` are not altered.
405 caltinay 3452 :note: at any point on the boundary of the domain the pressure
406     (``location_of_fixed_pressure`` >0) or the normal component of the
407     flux (``location_of_fixed_flux[i]>0``) if direction of the normal
408     is along the *x_i* axis.
409 jfenwick 3360
410 gross 3071 """
411 gross 3440 if self.useVPIteration:
412 gross 3501 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=util.wherePositive(location_of_fixed_flux))
414 gross 3440 else:
415     q=self.__pde_k.getCoefficient("q")
416     if q.isEmpty(): q=self.__pde_k.createCoefficient("q")
417 gross 3501 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()]=util.wherePositive(location_of_fixed_flux)
419 gross 3440 self.__pde_k.setValue(q=q)
420 gross 3071
421 gross 3440 # flux is rescaled by the factor mean value(perm_inv)*length where length**self.domain.getDim()=vol(self.domain)
422 gross 3501 V=util.vol(self.domain)
423 gross 3071 if permeability!=None:
424 gross 3440 perm=util.interpolate(permeability,self.__pde_k.getFunctionSpaceForCoefficient("A"))
425 gross 3501
426 gross 3071 if perm.getRank()==0:
427 gross 3440 perm_inv=(1./perm)
428     if self.useVPIteration:
429     self.scale=1.
430     else:
431     self.scale=util.integrate(perm_inv)*V**(1./self.domain.getDim()-1.)
432    
433     perm_inv=perm_inv*((1./self.scale)*util.kronecker(self.domain.getDim()))
434     perm=perm*(self.scale*util.kronecker(self.domain.getDim()))
435 gross 3071 elif perm.getRank()==2:
436     perm_inv=util.inverse(perm)
437 gross 3440 if self.useVPIteration:
438     self.scale=1.
439     else:
440     self.scale=util.sqrt(util.integrate(util.length(perm_inv)**2)*V**(2./self.domain.getDim()-1.)/self.domain.getDim())
441     perm_inv*=(1./self.scale)
442     perm=perm*self.scale
443 gross 3071 else:
444     raise ValueError,"illegal rank of permeability."
445 gross 2960
446 gross 3071 self.__permeability=perm
447     self.__permeability_inv=perm_inv
448 gross 3440
449 gross 3501 self.__l2 = V**(1./self.domain.getDim())*util.length(self.__permeability_inv)*self.domain.getSize()*self.weighting_scale
450 gross 3440 if self.useVPIteration:
451     if self.solveForFlux:
452     self.__pde_k.setValue(D=self.__permeability_inv)
453     else:
454     self.__pde_k.setValue(D=self.__permeability_inv, A=self.__l2*util.outer(util.kronecker(self.domain),util.kronecker(self.domain)))
455     self.__pde_p.setValue(A=self.__permeability)
456     else:
457     D=self.__pde_k.createCoefficient("D")
458     A=self.__pde_k.createCoefficient("A")
459     D[:self.domain.getDim(),:self.domain.getDim()]=self.__permeability_inv
460     for i in range(self.domain.getDim()):
461     for j in range(self.domain.getDim()):
462     A[i,i,j,j]=self.__l2
463     A[self.domain.getDim(),:,self.domain.getDim(),:]=self.__permeability
464     self.__pde_k.setValue(A=A, D=D)
465     if g !=None:
466     g=util.interpolate(g, self.__pde_k.getFunctionSpaceForCoefficient("Y"))
467     if g.isEmpty():
468     g=Vector(0,self.__pde_k.getFunctionSpaceForCoefficient("Y"))
469 gross 3071 else:
470 gross 3440 if not g.getShape()==(self.domain.getDim(),):
471     raise ValueError,"illegal shape of g"
472     self.__g=g
473     elif permeability!=None:
474     X
475     if f !=None:
476     f=util.interpolate(f, self.__pde_k.getFunctionSpaceForCoefficient("X"))
477     if f.isEmpty():
478     f=Scalar(0,self.__pde_k.getFunctionSpaceForCoefficient("X"))
479     else:
480     if f.getRank()>0: raise ValueError,"illegal rank of f."
481     self.__f=f
482 gross 2960
483 gross 3071
484 gross 3440 def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):
485     """
486     solves the problem.
487    
488     The iteration is terminated if the residual norm is less then self.getTolerance().
489 gross 2960
490 gross 3440 :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.
491     :type u0: vector value on the domain (e.g. `escript.Data`).
492     :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.
493     :type p0: scalar value on the domain (e.g. `escript.Data`).
494     :param verbose: if set some information on iteration progress are printed
495     :type verbose: ``bool``
496     :return: flux and pressure
497     :rtype: ``tuple`` of `escript.Data`.
498 gross 2960
499 gross 3440 :note: The problem is solved as a least squares form
500 caltinay 3452 *(K^[-1]+D^* l2 D)u+G p=D^* l2 * f + K^[-1]g*
501     *G^*u+*G^* K Gp=G^*g*
502     where *D* is the *div* operator and *(Gp)_i=p_{,i}* for the permeability *K=k_{ij}*.
503 gross 3440 """
504     self.verbose=verbose
505     if self.useVPIteration:
506     return self.__solveVP(u0,p0,max_iter,max_num_corrections)
507 gross 3071 else:
508 gross 3440 X=self.__pde_k.createCoefficient("X")
509     Y=self.__pde_k.createCoefficient("Y")
510     Y[:self.domain.getDim()]=self.scale*util.tensor_mult(self.__permeability_inv,self.__g)
511     rtmp=self.__f * self.__l2 * self.scale
512     for i in range(self.domain.getDim()): X[i,i]=rtmp
513     X[self.domain.getDim(),:]=self.__g*self.scale
514     r=self.__pde_k.createCoefficient("r")
515     r[:self.domain.getDim()]=u0*self.scale
516     r[self.domain.getDim()]=p0
517     self.__pde_k.setValue(X=X, Y=Y, r=r)
518     self.__pde_k.getSolverOptions().setVerbosity(self.verbose)
519     #self.__pde_k.getSolverOptions().setPreconditioner(self.__pde_k.getSolverOptions().AMG)
520     self.__pde_k.getSolverOptions().setSolverMethod(self.__pde_k.getSolverOptions().DIRECT)
521     U=self.__pde_k.getSolution()
522     # self.__pde_k.getOperator().saveMM("k.mm")
523     u=U[:self.domain.getDim()]*(1./self.scale)
524     p=U[self.domain.getDim()]
525     if self.verbose:
526     KGp=util.tensor_mult(self.__permeability,util.grad(p)/self.scale)
527     def_p=self.__g-(u+KGp)
528     def_v=self.__f-util.div(u, self.__pde_k.getFunctionSpaceForCoefficient("X"))
529     print "DarcyFlux: L2: g-v-K*grad(p) = %e (v = %e)."%(self.__L2(def_p),self.__L2(u))
530     print "DarcyFlux: L2: f-div(v) = %e (grad(v) = %e)."%(self.__L2(def_v),self.__L2(util.grad(u)))
531     return u,p
532 gross 2960
533 gross 3440 def __solveVP(self,u0,p0, max_iter=100, max_num_corrections=10):
534 gross 3071 """
535     solves the problem.
536    
537 caltinay 3452 The iteration is terminated if the residual norm is less than self.getTolerance().
538 gross 2960
539 gross 3071 :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.
540 gross 3440 :type u0: vector value on the domain (e.g. `escript.Data`).
541 gross 3071 :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.
542 gross 3440 :type p0: scalar value on the domain (e.g. `escript.Data`).
543 gross 3071 :return: flux and pressure
544 gross 3440 :rtype: ``tuple`` of `escript.Data`.
545 gross 2960
546 gross 3071 :note: The problem is solved as a least squares form
547 caltinay 3452 *(K^[-1]+D^* (DKD^*)^[-1] D)u+G p=D^* (DKD^*)^[-1] f + K^[-1]g*
548     *G^*u+*G^* K Gp=G^*g*
549     where *D* is the *div* operator and *(Gp)_i=p_{,i}* for the permeability *K=k_{ij}*.
550 gross 3071 """
551     rtol=self.getTolerance()
552     atol=self.getAbsoluteTolerance()
553     self.setSubProblemTolerance()
554     num_corrections=0
555     converged=False
556     norm_r=None
557    
558     # Eliminate the hydrostatic pressure:
559     if self.verbose: print "DarcyFlux: calculate hydrostatic pressure component."
560     self.__pde_p.setValue(X=self.__g, r=p0, y=-util.inner(self.domain.getNormal(),u0))
561     p0=self.__pde_p.getSolution()
562     g2=self.__g - util.tensor_mult(self.__permeability, util.grad(p0))
563     norm_g2=util.integrate(util.inner(g2,util.tensor_mult(self.__permeability_inv,g2)))**0.5
564 gross 2960
565 gross 3071 p=p0*0
566     if self.solveForFlux:
567     v=u0.copy()
568     else:
569     v=self.__getFlux(p, u0, f=self.__f, g=g2)
570 gross 2960
571 gross 3071 while not converged and norm_g2 > 0:
572     Gp=util.grad(p)
573     KGp=util.tensor_mult(self.__permeability,Gp)
574     if self.verbose:
575     def_p=g2-(v+KGp)
576     def_v=self.__f-util.div(v)
577     print "DarcyFlux: L2: g-v-K*grad(p) = %e (v = %e)."%(self.__L2(def_p),self.__L2(v))
578     print "DarcyFlux: L2: f-div(v) = %e (grad(v) = %e)."%(self.__L2(def_v),self.__L2(util.grad(v)))
579     print "DarcyFlux: K^{-1}-norm of v = %e."%util.integrate(util.inner(v,util.tensor_mult(self.__permeability_inv,v)))**0.5
580     print "DarcyFlux: K^{-1}-norm of g2 = %e."%norm_g2
581     print "DarcyFlux: K-norm of grad(dp) = %e."%util.integrate(util.inner(Gp,KGp))**0.5
582     ATOL=atol+rtol*norm_g2
583     if self.verbose: print "DarcyFlux: absolute tolerance ATOL = %e."%(ATOL,)
584     if norm_r == None or norm_r>ATOL:
585     if num_corrections>max_num_corrections:
586     raise ValueError,"maximum number of correction steps reached."
587    
588     if self.solveForFlux:
589     # initial residual is r=K^{-1}*(g-v-K*Gp)+D^*L^{-1}(f-Du)
590     v,r, norm_r=PCG(ArithmeticTuple(util.tensor_mult(self.__permeability_inv,g2-v)-Gp,self.__applWeight(v,self.__f),p),
591     self.__Aprod_v,
592     v,
593     self.__Msolve_PCG_v,
594     self.__inner_PCG_v,
595     atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
596     p=r[2]
597     else:
598     # initial residual is r=G^*(g2-KGp - v)
599     p,r, norm_r=PCG(ArithmeticTuple(g2-KGp,v),
600     self.__Aprod_p,
601     p,
602     self.__Msolve_PCG_p,
603     self.__inner_PCG_p,
604     atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
605     v=r[1]
606     if self.verbose: print "DarcyFlux: residual norm = %e."%norm_r
607     num_corrections+=1
608     else:
609     if self.verbose: print "DarcyFlux: stopping criterium reached."
610     converged=True
611     return v,p+p0
612 gross 3440
613     def __applWeight(self, v, f=None):
614     # solves L p = f-Dv with p = 0
615     if self.verbose: print "DarcyFlux: Applying weighting operator"
616     if f == None:
617     return -util.div(v)*self.__l2
618     else:
619     return (f-util.div(v))*self.__l2
620     def __getPressure(self, v, p0, g=None):
621     # solves (G*KG)p = G^(g-v) with p = p0 where location_of_fixed_pressure>0
622     if self.getSolverOptionsPressure().isVerbose() or self.verbose: print "DarcyFlux: Pressure update"
623     if g == None:
624     self.__pde_p.setValue(X=-v, r=p0)
625     else:
626     self.__pde_p.setValue(X=g-v, r=p0)
627     p=self.__pde_p.getSolution()
628     return p
629    
630     def __Aprod_v(self,dv):
631     # calculates: (a,b,c) = (K^{-1}(dv + KG * dp), L^{-1}Ddv, dp) with (G*KG)dp = - G^*dv
632     dp=self.__getPressure(dv, p0=escript.Data()) # dp = (G*KG)^{-1} (0-G^*dv)
633     a=util.tensor_mult(self.__permeability_inv,dv)+util.grad(dp) # a= K^{-1}u+G*dp
634     b= - self.__applWeight(dv) # b = - (D K D^*)^{-1} (0-Dv)
635     return ArithmeticTuple(a,b,-dp)
636    
637     def __Msolve_PCG_v(self,r):
638     # K^{-1} u = r[0] + D^*r[1] = K^{-1}(dv + KG * dp) + D^*L^{-1}Ddv
639     if self.getSolverOptionsFlux().isVerbose() or self.verbose: print "DarcyFlux: Applying preconditioner"
640     self.__pde_k.setValue(X=r[1]*util.kronecker(self.domain), Y=r[0], r=escript.Data())
641     return self.__pde_k.getSolution()
642    
643     def __inner_PCG_v(self,v,r):
644     return util.integrate(util.inner(v,r[0])+util.div(v)*r[1])
645    
646     def __Aprod_p(self,dp):
647     if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
648     Gdp=util.grad(dp)
649     self.__pde_k.setValue(Y=-Gdp,X=escript.Data(), r=escript.Data())
650     du=self.__pde_k.getSolution()
651     return ArithmeticTuple(util.tensor_mult(self.__permeability,Gdp),-du)
652    
653     def __getFlux(self,p, v0, f=None, g=None):
654     # solves (K^{-1}+D^*L^{-1} D) v = D^*L^{-1}f + K^{-1}g - Gp
655     if f!=None:
656     self.__pde_k.setValue(X=self.__applWeight(v0*0,self.__f)*util.kronecker(self.domain))
657     self.__pde_k.setValue(r=v0)
658     g2=util.tensor_mult(self.__permeability_inv,g)
659     if p == None:
660     self.__pde_k.setValue(Y=g2)
661     else:
662     self.__pde_k.setValue(Y=g2-util.grad(p))
663     return self.__pde_k.getSolution()
664    
665     #v=self.__getFlux(p, u0, f=self.__f, g=g2)
666     def __Msolve_PCG_p(self,r):
667     if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"
668     self.__pde_p.setValue(X=r[0]-r[1], Y=escript.Data(), r=escript.Data(), y=escript.Data())
669     return self.__pde_p.getSolution()
670    
671     def __inner_PCG_p(self,p,r):
672     return util.integrate(util.inner(util.grad(p), r[0]-r[1]))
673    
674     def __L2(self,v):
675     return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))
676    
677     def __L2_r(self,v):
678     return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.ReducedFunction(self.domain)))**2))
679    
680 gross 3071 def setTolerance(self,rtol=1e-4):
681     """
682     sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if
683 gross 2960
684 gross 3071 *|g-v-K gard(p)|_PCG <= atol + rtol * |K^{1/2}g2|_0*
685    
686     where ``atol`` is an absolut tolerance (see `setAbsoluteTolerance`).
687    
688     :param rtol: relative tolerance for the pressure
689     :type rtol: non-negative ``float``
690     """
691     if rtol<0:
692     raise ValueError,"Relative tolerance needs to be non-negative."
693     self.__rtol=rtol
694     def getTolerance(self):
695     """
696     returns the relative tolerance
697     :return: current relative tolerance
698     :rtype: ``float``
699     """
700     return self.__rtol
701 gross 2960
702 gross 3071 def setAbsoluteTolerance(self,atol=0.):
703     """
704     sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if
705    
706     *|g-v-K gard(p)|_PCG <= atol + rtol * |K^{1/2}g2|_0*
707 gross 2960
708    
709 gross 3071 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}*.
710 gross 2960
711 gross 3071 :param atol: absolute tolerance for the pressure
712     :type atol: non-negative ``float``
713     """
714     if atol<0:
715     raise ValueError,"Absolute tolerance needs to be non-negative."
716     self.__atol=atol
717     def getAbsoluteTolerance(self):
718     """
719     returns the absolute tolerance
720     :return: current absolute tolerance
721     :rtype: ``float``
722     """
723     return self.__atol
724     def getSubProblemTolerance(self):
725     """
726     Returns a suitable subtolerance
727     :type: ``float``
728     """
729 gross 3074 return max(util.EPSILON**(0.5),self.getTolerance()**2)
730 gross 2960
731 gross 3071 def setSubProblemTolerance(self):
732     """
733     Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected.
734     """
735     if self.__adaptSubTolerance:
736     sub_tol=self.getSubProblemTolerance()
737     self.getSolverOptionsFlux().setTolerance(sub_tol)
738     self.getSolverOptionsFlux().setAbsoluteTolerance(0.)
739     self.getSolverOptionsPressure().setTolerance(sub_tol)
740     self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
741     if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol
742 gross 2960
743    
744     class DarcyFlowOld(object):
745     """
746     solves the problem
747    
748     *u_i+k_{ij}*p_{,j} = g_i*
749     *u_{i,i} = f*
750    
751     where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,
752    
753     :note: The problem is solved in a least squares formulation.
754     """
755    
756     def __init__(self, domain, weight=None, useReduced=False, adaptSubTolerance=True):
757     """
758     initializes the Darcy flux problem
759     :param domain: domain of the problem
760     :type domain: `Domain`
761     :param useReduced: uses reduced oreder on flux and pressure
762     :type useReduced: ``bool``
763     :param adaptSubTolerance: switches on automatic subtolerance selection
764     :type adaptSubTolerance: ``bool``
765     """
766     self.domain=domain
767     if weight == None:
768 gross 2354 s=self.domain.getSize()
769 gross 3440 self.__l2=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2
770     # self.__l2=(3.*util.longestEdge(self.domain))**2
771     #self.__l2=(0.1*util.longestEdge(self.domain)*s/util.sup(s))**2
772 gross 2351 else:
773 gross 3440 self.__l2=weight
774 gross 2100 self.__pde_v=LinearPDESystem(domain)
775 gross 2208 if useReduced: self.__pde_v.setReducedOrderOn()
776     self.__pde_v.setSymmetryOn()
777 gross 3440 self.__pde_v.setValue(D=util.kronecker(domain), A=self.__l2*util.outer(util.kronecker(domain),util.kronecker(domain)))
778 gross 2100 self.__pde_p=LinearSinglePDE(domain)
779     self.__pde_p.setSymmetryOn()
780 gross 2208 if useReduced: self.__pde_p.setReducedOrderOn()
781 jfenwick 3432 self.__f=escript.Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
782     self.__g=escript.Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
783 gross 2264 self.setTolerance()
784     self.setAbsoluteTolerance()
785 gross 2474 self.__adaptSubTolerance=adaptSubTolerance
786     self.verbose=False
787     def getSolverOptionsFlux(self):
788     """
789     Returns the solver options used to solve the flux problems
790    
791 jfenwick 2625 *(I+D^*D)u=F*
792 gross 2474
793 jfenwick 2625 :return: `SolverOptions`
794 gross 2474 """
795     return self.__pde_v.getSolverOptions()
796     def setSolverOptionsFlux(self, options=None):
797     """
798     Sets the solver options used to solve the flux problems
799    
800 jfenwick 2625 *(I+D^*D)u=F*
801 gross 2474
802 jfenwick 2625 If ``options`` is not present, the options are reset to default
803     :param options: `SolverOptions`
804     :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
805 gross 2474 """
806     return self.__pde_v.setSolverOptions(options)
807     def getSolverOptionsPressure(self):
808     """
809     Returns the solver options used to solve the pressure problems
810    
811 jfenwick 2625 *(Q^*Q)p=Q^*G*
812 gross 2474
813 jfenwick 2625 :return: `SolverOptions`
814 gross 2474 """
815     return self.__pde_p.getSolverOptions()
816     def setSolverOptionsPressure(self, options=None):
817     """
818     Sets the solver options used to solve the pressure problems
819    
820 jfenwick 2625 *(Q^*Q)p=Q^*G*
821 gross 2474
822 jfenwick 2625 If ``options`` is not present, the options are reset to default
823     :param options: `SolverOptions`
824     :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
825 gross 2474 """
826     return self.__pde_p.setSolverOptions(options)
827 gross 1659
828 gross 2100 def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
829     """
830 gross 2208 assigns values to model parameters
831 gross 1659
832 jfenwick 2625 :param f: volumetic sources/sinks
833 gross 3440 :type f: scalar value on the domain (e.g. `escript.Data`)
834 jfenwick 2625 :param g: flux sources/sinks
835 gross 3440 :type g: vector values on the domain (e.g. `escript.Data`)
836 jfenwick 2625 :param location_of_fixed_pressure: mask for locations where pressure is fixed
837 gross 3440 :type location_of_fixed_pressure: scalar value on the domain (e.g. `escript.Data`)
838 jfenwick 2625 :param location_of_fixed_flux: mask for locations where flux is fixed.
839 gross 3440 :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)
840 jfenwick 2625 :param permeability: permeability tensor. If scalar ``s`` is given the tensor with
841     ``s`` on the main diagonal is used. If vector ``v`` is given the tensor with
842     ``v`` on the main diagonal is used.
843 gross 3440 :type permeability: scalar, vector or tensor values on the domain (e.g. `escript.Data`)
844 gross 1659
845 jfenwick 2625 :note: the values of parameters which are not set by calling ``setValue`` are not altered.
846     :note: at any point on the boundary of the domain the pressure (``location_of_fixed_pressure`` >0)
847     or the normal component of the flux (``location_of_fixed_flux[i]>0`` if direction of the normal
848     is along the *x_i* axis.
849 gross 2100 """
850 gross 2264 if f !=None:
851 gross 2100 f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
852     if f.isEmpty():
853 jfenwick 3432 f=escript.Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
854 gross 2100 else:
855     if f.getRank()>0: raise ValueError,"illegal rank of f."
856 gross 2267 self.__f=f
857 gross 2264 if g !=None:
858 gross 2100 g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
859     if g.isEmpty():
860 jfenwick 3432 g=escript.Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
861 gross 2100 else:
862     if not g.getShape()==(self.domain.getDim(),):
863     raise ValueError,"illegal shape of g"
864     self.__g=g
865 gross 1659
866 gross 2100 if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)
867     if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux)
868 gross 1659
869 gross 2100 if permeability!=None:
870     perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
871     if perm.getRank()==0:
872     perm=perm*util.kronecker(self.domain.getDim())
873     elif perm.getRank()==1:
874     perm, perm2=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A")), perm
875     for i in range(self.domain.getDim()): perm[i,i]=perm2[i]
876     elif perm.getRank()==2:
877     pass
878     else:
879     raise ValueError,"illegal rank of permeability."
880     self.__permeability=perm
881     self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))
882 gross 1659
883 gross 2264 def setTolerance(self,rtol=1e-4):
884     """
885 jfenwick 2625 sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if
886 gross 1659
887 jfenwick 2625 *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
888 gross 2264
889 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}*.
890 gross 2264
891 jfenwick 2625 :param rtol: relative tolerance for the pressure
892     :type rtol: non-negative ``float``
893 gross 2100 """
894 gross 2264 if rtol<0:
895     raise ValueError,"Relative tolerance needs to be non-negative."
896     self.__rtol=rtol
897     def getTolerance(self):
898 gross 2100 """
899 gross 2264 returns the relative tolerance
900 gross 1659
901 jfenwick 2625 :return: current relative tolerance
902     :rtype: ``float``
903 caltinay 2169 """
904 gross 2264 return self.__rtol
905 gross 1659
906 gross 2264 def setAbsoluteTolerance(self,atol=0.):
907 gross 2208 """
908 jfenwick 2625 sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if
909 gross 1659
910 jfenwick 2625 *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
911 gross 2208
912 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}*.
913 gross 2208
914 jfenwick 2625 :param atol: absolute tolerance for the pressure
915     :type atol: non-negative ``float``
916 gross 2264 """
917     if atol<0:
918     raise ValueError,"Absolute tolerance needs to be non-negative."
919     self.__atol=atol
920     def getAbsoluteTolerance(self):
921     """
922     returns the absolute tolerance
923    
924 jfenwick 2625 :return: current absolute tolerance
925     :rtype: ``float``
926 gross 2264 """
927     return self.__atol
928     def getSubProblemTolerance(self):
929 gross 2474 """
930     Returns a suitable subtolerance
931 jfenwick 2625 @type: ``float``
932 gross 2474 """
933     return max(util.EPSILON**(0.75),self.getTolerance()**2)
934     def setSubProblemTolerance(self):
935 gross 2264 """
936 gross 2474 Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected.
937 gross 2264 """
938 gross 2474 if self.__adaptSubTolerance:
939     sub_tol=self.getSubProblemTolerance()
940     self.getSolverOptionsFlux().setTolerance(sub_tol)
941     self.getSolverOptionsFlux().setAbsoluteTolerance(0.)
942     self.getSolverOptionsPressure().setTolerance(sub_tol)
943     self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
944     if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol
945 gross 2208
946 gross 2474 def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):
947 gross 2264 """
948     solves the problem.
949 gross 1659
950 gross 2208 The iteration is terminated if the residual norm is less then self.getTolerance().
951 gross 1659
952 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.
953 gross 3440 :type u0: vector value on the domain (e.g. `escript.Data`).
954 jfenwick 2625 :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.
955 gross 3440 :type p0: scalar value on the domain (e.g. `escript.Data`).
956 jfenwick 2625 :param verbose: if set some information on iteration progress are printed
957     :type verbose: ``bool``
958     :return: flux and pressure
959 gross 3440 :rtype: ``tuple`` of `escript.Data`.
960 gross 2264
961 jfenwick 2625 :note: The problem is solved as a least squares form
962 gross 2100
963 jfenwick 2625 *(I+D^*D)u+Qp=D^*f+g*
964     *Q^*u+Q^*Qp=Q^*g*
965 gross 2100
966 jfenwick 2625 where *D* is the *div* operator and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.
967 gross 2264 We eliminate the flux form the problem by setting
968 caltinay 2169
969 jfenwick 2625 *u=(I+D^*D)^{-1}(D^*f-g-Qp)* with u=u0 on location_of_fixed_flux
970 caltinay 2169
971 gross 2208 form the first equation. Inserted into the second equation we get
972 caltinay 2169
973 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
974 gross 2264
975 jfenwick 2625 which is solved using the PCG method (precondition is *Q^*Q*). In each iteration step
976     PDEs with operator *I+D^*D* and with *Q^*Q* needs to be solved using a sub iteration scheme.
977 gross 2208 """
978 gross 2386 self.verbose=verbose
979 gross 2264 rtol=self.getTolerance()
980     atol=self.getAbsoluteTolerance()
981 gross 2474 self.setSubProblemTolerance()
982 gross 2264 num_corrections=0
983     converged=False
984     p=p0
985     norm_r=None
986     while not converged:
987 gross 2474 v=self.getFlux(p, fixed_flux=u0)
988 gross 2264 Qp=self.__Q(p)
989     norm_v=self.__L2(v)
990     norm_Qp=self.__L2(Qp)
991     if norm_v == 0.:
992     if norm_Qp == 0.:
993     return v,p
994     else:
995     fac=norm_Qp
996     else:
997     if norm_Qp == 0.:
998     fac=norm_v
999     else:
1000     fac=2./(1./norm_v+1./norm_Qp)
1001     ATOL=(atol+rtol*fac)
1002     if self.verbose:
1003     print "DarcyFlux: L2 norm of v = %e."%norm_v
1004 gross 2960 print "DarcyFlux: L2 norm of k.util.grad(p) = %e."%norm_Qp
1005 jfenwick 3432 print "DarcyFlux: L2 defect u = %e."%(util.integrate(util.length(self.__g-util.interpolate(v,escript.Function(self.domain))-Qp)**2)**(0.5),)
1006 gross 2486 print "DarcyFlux: L2 defect div(v) = %e."%(util.integrate((self.__f-util.div(v))**2)**(0.5),)
1007 gross 2264 print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
1008     if norm_r == None or norm_r>ATOL:
1009     if num_corrections>max_num_corrections:
1010     raise ValueError,"maximum number of correction steps reached."
1011 jfenwick 3432 p,r, norm_r=PCG(self.__g-util.interpolate(v,escript.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)
1012 gross 2264 num_corrections+=1
1013     else:
1014     converged=True
1015     return v,p
1016     def __L2(self,v):
1017 jfenwick 3432 return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))
1018 gross 2264
1019     def __Q(self,p):
1020     return util.tensor_mult(self.__permeability,util.grad(p))
1021    
1022     def __Aprod(self,dp):
1023 gross 2474 if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
1024 gross 2264 Qdp=self.__Q(dp)
1025 jfenwick 3432 self.__pde_v.setValue(Y=-Qdp,X=escript.Data(), r=escript.Data())
1026 gross 2474 du=self.__pde_v.getSolution()
1027 gross 2264 return Qdp+du
1028     def __inner_GMRES(self,r,s):
1029     return util.integrate(util.inner(r,s))
1030    
1031 gross 2100 def __inner_PCG(self,p,r):
1032 gross 2264 return util.integrate(util.inner(self.__Q(p), r))
1033 gross 2100
1034     def __Msolve_PCG(self,r):
1035 gross 2474 if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"
1036 jfenwick 3432 self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=escript.Data(), r=escript.Data())
1037 gross 2474 return self.__pde_p.getSolution()
1038 gross 2100
1039 jfenwick 3432 def getFlux(self,p=None, fixed_flux=escript.Data()):
1040 gross 2264 """
1041 jfenwick 2625 returns the flux for a given pressure ``p`` where the flux is equal to ``fixed_flux``
1042     on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
1043     Note that ``g`` and ``f`` are used, see `setValue`.
1044 gross 2264
1045 jfenwick 2625 :param p: pressure.
1046 gross 3440 :type p: scalar value on the domain (e.g. `escript.Data`).
1047 jfenwick 2625 :param fixed_flux: flux on the locations of the domain marked be ``location_of_fixed_flux``.
1048 gross 3440 :type fixed_flux: vector values on the domain (e.g. `escript.Data`).
1049 jfenwick 2625 :return: flux
1050 gross 3440 :rtype: `escript.Data`
1051 jfenwick 2625 :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}*
1052     for the permeability *k_{ij}*
1053 gross 2264 """
1054 gross 2474 self.setSubProblemTolerance()
1055 gross 2264 g=self.__g
1056     f=self.__f
1057 gross 3440 self.__pde_v.setValue(X=self.__l2*f*util.kronecker(self.domain), r=fixed_flux)
1058 gross 2264 if p == None:
1059     self.__pde_v.setValue(Y=g)
1060     else:
1061     self.__pde_v.setValue(Y=g-self.__Q(p))
1062 gross 2474 return self.__pde_v.getSolution()
1063 gross 2264
1064 gross 1414 class StokesProblemCartesian(HomogeneousSaddlePointProblem):
1065 gross 2251 """
1066 gross 2264 solves
1067 gross 1414
1068 gross 2208 -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
1069     u_{i,i}=0
1070 gross 1414
1071 gross 2208 u=0 where fixed_u_mask>0
1072     eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j
1073 gross 1414
1074 gross 2264 if surface_stress is not given 0 is assumed.
1075 gross 1414
1076 gross 2251 typical usage:
1077 gross 1414
1078 gross 2208 sp=StokesProblemCartesian(domain)
1079     sp.setTolerance()
1080     sp.initialize(...)
1081     v,p=sp.solve(v0,p0)
1082 gross 2251 """
1083 gross 2719 def __init__(self,domain,**kwargs):
1084 gross 2100 """
1085 gross 2208 initialize the Stokes Problem
1086 gross 2100
1087 gross 2793 The approximation spaces used for velocity (=Solution(domain)) and pressure (=ReducedSolution(domain)) must be
1088     LBB complient, for instance using quadratic and linear approximation on the same element or using linear approximation
1089     with macro elements for the pressure.
1090    
1091     :param domain: domain of the problem.
1092 jfenwick 2625 :type domain: `Domain`
1093 gross 2100 """
1094 gross 2719 HomogeneousSaddlePointProblem.__init__(self,**kwargs)
1095 gross 1414 self.domain=domain
1096     self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
1097     self.__pde_u.setSymmetryOn()
1098 gross 2474
1099 gross 1414 self.__pde_prec=LinearPDE(domain)
1100     self.__pde_prec.setReducedOrderOn()
1101     self.__pde_prec.setSymmetryOn()
1102    
1103 gross 2415 self.__pde_proj=LinearPDE(domain)
1104     self.__pde_proj.setReducedOrderOn()
1105     self.__pde_proj.setValue(D=1)
1106     self.__pde_proj.setSymmetryOn()
1107    
1108 gross 2474 def getSolverOptionsVelocity(self):
1109     """
1110     returns the solver options used solve the equation for velocity.
1111    
1112 jfenwick 2625 :rtype: `SolverOptions`
1113 gross 2474 """
1114     return self.__pde_u.getSolverOptions()
1115     def setSolverOptionsVelocity(self, options=None):
1116     """
1117     set the solver options for solving the equation for velocity.
1118    
1119 jfenwick 2625 :param options: new solver options
1120     :type options: `SolverOptions`
1121 gross 2474 """
1122     self.__pde_u.setSolverOptions(options)
1123     def getSolverOptionsPressure(self):
1124     """
1125     returns the solver options used solve the equation for pressure.
1126 jfenwick 2625 :rtype: `SolverOptions`
1127 gross 2474 """
1128     return self.__pde_prec.getSolverOptions()
1129     def setSolverOptionsPressure(self, options=None):
1130     """
1131     set the solver options for solving the equation for pressure.
1132 jfenwick 2625 :param options: new solver options
1133     :type options: `SolverOptions`
1134 gross 2474 """
1135     self.__pde_prec.setSolverOptions(options)
1136 gross 2415
1137 gross 2474 def setSolverOptionsDiv(self, options=None):
1138     """
1139     set the solver options for solving the equation to project the divergence of
1140     the velocity onto the function space of presure.
1141    
1142 jfenwick 2625 :param options: new solver options
1143     :type options: `SolverOptions`
1144 gross 2474 """
1145 artak 2689 self.__pde_proj.setSolverOptions(options)
1146 gross 2474 def getSolverOptionsDiv(self):
1147     """
1148     returns the solver options for solving the equation to project the divergence of
1149     the velocity onto the function space of presure.
1150    
1151 jfenwick 2625 :rtype: `SolverOptions`
1152 gross 2474 """
1153 artak 2689 return self.__pde_proj.getSolverOptions()
1154 gross 2474
1155 gross 2793 def updateStokesEquation(self, v, p):
1156     """
1157     updates the Stokes equation to consider dependencies from ``v`` and ``p``
1158     :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values.
1159     """
1160     pass
1161     def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None):
1162     """
1163     assigns new values to the model parameters.
1164    
1165     :param f: external force
1166     :type f: `Vector` object in `FunctionSpace` `Function` or similar
1167     :param fixed_u_mask: mask of locations with fixed velocity.
1168     :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
1169     :param eta: viscosity
1170     :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
1171     :param surface_stress: normal surface stress
1172     :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
1173     :param stress: initial stress
1174     :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
1175     """
1176     if eta !=None:
1177     k=util.kronecker(self.domain.getDim())
1178     kk=util.outer(k,k)
1179 jfenwick 3432 self.eta=util.interpolate(eta, escript.Function(self.domain))
1180 gross 2793 self.__pde_prec.setValue(D=1/self.eta)
1181     self.__pde_u.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))
1182     if restoration_factor!=None:
1183     n=self.domain.getNormal()
1184     self.__pde_u.setValue(d=restoration_factor*util.outer(n,n))
1185     if fixed_u_mask!=None:
1186     self.__pde_u.setValue(q=fixed_u_mask)
1187     if f!=None: self.__f=f
1188     if surface_stress!=None: self.__surface_stress=surface_stress
1189     if stress!=None: self.__stress=stress
1190    
1191 jfenwick 3432 def initialize(self,f=escript.Data(),fixed_u_mask=escript.Data(),eta=1,surface_stress=escript.Data(),stress=escript.Data(), restoration_factor=0):
1192 gross 2208 """
1193     assigns values to the model parameters
1194 gross 2100
1195 jfenwick 2625 :param f: external force
1196     :type f: `Vector` object in `FunctionSpace` `Function` or similar
1197     :param fixed_u_mask: mask of locations with fixed velocity.
1198     :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
1199     :param eta: viscosity
1200     :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
1201     :param surface_stress: normal surface stress
1202     :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
1203     :param stress: initial stress
1204     :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
1205 gross 2208 """
1206 gross 2793 self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
1207 gross 1414
1208 gross 2719 def Bv(self,v,tol):
1209 gross 2251 """
1210     returns inner product of element p and div(v)
1211 gross 1414
1212 jfenwick 2625 :param v: a residual
1213     :return: inner product of element p and div(v)
1214     :rtype: ``float``
1215 gross 2100 """
1216 gross 2793 self.__pde_proj.setValue(Y=-util.div(v))
1217 gross 2719 self.getSolverOptionsDiv().setTolerance(tol)
1218     self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
1219     out=self.__pde_proj.getSolution()
1220     return out
1221 gross 2208
1222 gross 2445 def inner_pBv(self,p,Bv):
1223     """
1224     returns inner product of element p and Bv=-div(v)
1225    
1226 jfenwick 2625 :param p: a pressure increment
1227     :param Bv: a residual
1228     :return: inner product of element p and Bv=-div(v)
1229     :rtype: ``float``
1230 gross 2445 """
1231 jfenwick 3432 return util.integrate(util.interpolate(p,escript.Function(self.domain))*util.interpolate(Bv, escript.Function(self.domain)))
1232 gross 2445
1233 gross 2251 def inner_p(self,p0,p1):
1234 gross 2100 """
1235 gross 2251 Returns inner product of p0 and p1
1236 gross 1414
1237 jfenwick 2625 :param p0: a pressure
1238     :param p1: a pressure
1239     :return: inner product of p0 and p1
1240     :rtype: ``float``
1241 gross 2100 """
1242 jfenwick 3432 s0=util.interpolate(p0, escript.Function(self.domain))
1243     s1=util.interpolate(p1, escript.Function(self.domain))
1244 gross 2100 return util.integrate(s0*s1)
1245 artak 1550
1246 gross 2251 def norm_v(self,v):
1247 gross 2100 """
1248 gross 2251 returns the norm of v
1249 gross 2208
1250 jfenwick 2625 :param v: a velovity
1251     :return: norm of v
1252     :rtype: non-negative ``float``
1253 gross 2100 """
1254 gross 2719 return util.sqrt(util.integrate(util.length(util.grad(v))**2))
1255 gross 2100
1256 gross 2793
1257 gross 2719 def getDV(self, p, v, tol):
1258 gross 1414 """
1259 gross 2251 return the value for v for a given p (overwrite)
1260    
1261 jfenwick 2625 :param p: a pressure
1262 gross 2719 :param v: a initial guess for the value v to return.
1263     :return: dv given as *Adv=(f-Av-B^*p)*
1264 gross 1414 """
1265 gross 2793 self.updateStokesEquation(v,p)
1266 gross 2719 self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress)
1267     self.getSolverOptionsVelocity().setTolerance(tol)
1268     self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
1269 gross 2100 if self.__stress.isEmpty():
1270 gross 2719 self.__pde_u.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
1271 gross 2100 else:
1272 gross 2719 self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
1273 gross 2474 out=self.__pde_u.getSolution()
1274 gross 2208 return out
1275 gross 1414
1276 gross 2445 def norm_Bv(self,Bv):
1277 gross 2251 """
1278     Returns Bv (overwrite).
1279    
1280 jfenwick 2625 :rtype: equal to the type of p
1281     :note: boundary conditions on p should be zero!
1282 gross 2251 """
1283 jfenwick 3432 return util.sqrt(util.integrate(util.interpolate(Bv, escript.Function(self.domain))**2))
1284 gross 2251
1285 gross 2719 def solve_AinvBt(self,p, tol):
1286 gross 2251 """
1287 gross 2719 Solves *Av=B^*p* with accuracy `tol`
1288 gross 2251
1289 jfenwick 2625 :param p: a pressure increment
1290     :return: the solution of *Av=B^*p*
1291     :note: boundary conditions on v should be zero!
1292 gross 2251 """
1293 jfenwick 3432 self.__pde_u.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))
1294 gross 2474 out=self.__pde_u.getSolution()
1295 gross 2251 return out
1296    
1297 gross 2719 def solve_prec(self,Bv, tol):
1298 gross 2251 """
1299 jfenwick 2625 applies preconditioner for for *BA^{-1}B^** to *Bv*
1300     with accuracy `self.getSubProblemTolerance()`
1301 gross 2251
1302 jfenwick 2625 :param Bv: velocity increment
1303     :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*
1304     :note: boundary conditions on p are zero.
1305 gross 2251 """
1306 gross 2445 self.__pde_prec.setValue(Y=Bv)
1307 gross 2719 self.getSolverOptionsPressure().setTolerance(tol)
1308     self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
1309     out=self.__pde_prec.getSolution()
1310     return out

  ViewVC Help
Powered by ViewVC 1.1.26