/[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 3440 - (hide annotations)
Fri Jan 14 00:04:53 2011 UTC (8 years, 8 months ago) by gross
File MIME type: text/x-python
File size: 44138 byte(s)
a second AMG interpolation implemented
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 3440 print dir(escript)
41    
42 gross 2100 class DarcyFlow(object):
43 gross 3071 """
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 least squares formulation.
52     """
53    
54 gross 3440 def __init__(self, domain, useReduced=False, adaptSubTolerance=True, solveForFlux=False, useVPIteration=True, weighting_scale=0.1):
55 gross 3071 """
56     initializes the Darcy flux problem
57     :param domain: domain of the problem
58     :type domain: `Domain`
59     :param useReduced: uses reduced oreder on flux and pressure
60     :type useReduced: ``bool``
61     :param adaptSubTolerance: switches on automatic subtolerance selection
62     :type adaptSubTolerance: ``bool``
63     :param solveForFlux: if True the solver solves for the flux (do not use!)
64     :type solveForFlux: ``bool``
65 gross 3440 :param useVPIteration: if True altenative iteration over v and p is performed. Otherwise V and P are calculated in a single PDE.
66     :type useVPIteration: ``bool``
67 gross 3071 """
68     self.domain=domain
69 gross 3440 self.useVPIteration=useVPIteration
70 gross 3071 self.useReduced=useReduced
71 gross 3440 self.weighting_scale=weighting_scale
72     if self.useVPIteration:
73     self.solveForFlux=solveForFlux
74     self.__adaptSubTolerance=adaptSubTolerance
75     self.verbose=False
76    
77     self.__pde_k=LinearPDESystem(domain)
78     self.__pde_k.setSymmetryOn()
79     if self.useReduced: self.__pde_k.setReducedOrderOn()
80    
81     self.__pde_p=LinearSinglePDE(domain)
82     self.__pde_p.setSymmetryOn()
83     if self.useReduced: self.__pde_p.setReducedOrderOn()
84     self.setTolerance()
85     self.setAbsoluteTolerance()
86     else:
87     self.__pde_k=LinearPDE(self.domain, numEquations=self.domain.getDim()+1)
88     self.__pde_k.setSymmetryOn()
89     if self.useReduced: self.__pde_k.setReducedOrderOn()
90     C=self.__pde_k.createCoefficient("C")
91     B=self.__pde_k.createCoefficient("B")
92     for i in range(self.domain.getDim()):
93     C[i,self.domain.getDim(),i]=1
94     B[self.domain.getDim(),i,i]=1
95     self.__pde_k.setValue(C=C, B=B)
96     self.__f=escript.Scalar(0,self.__pde_k.getFunctionSpaceForCoefficient("X"))
97     self.__g=escript.Vector(0,self.__pde_k.getFunctionSpaceForCoefficient("Y"))
98 gross 3071
99     def getSolverOptionsFlux(self):
100     """
101     Returns the solver options used to solve the flux problems
102    
103     *K^{-1} u=F*
104    
105     :return: `SolverOptions`
106     """
107     return self.__pde_k.getSolverOptions()
108    
109     def setSolverOptionsFlux(self, options=None):
110     """
111     Sets the solver options used to solve the flux problems
112    
113     *K^{-1}u=F*
114    
115     If ``options`` is not present, the options are reset to default
116    
117     :param options: `SolverOptions`
118     :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
119     """
120     return self.__pde_v.setSolverOptions(options)
121    
122     def getSolverOptionsPressure(self):
123     """
124     Returns the solver options used to solve the pressure problems
125    
126     *(Q^* K Q)p=-Q^*G*
127    
128     :return: `SolverOptions`
129     """
130     return self.__pde_p.getSolverOptions()
131    
132     def setSolverOptionsPressure(self, options=None):
133     """
134     Sets the solver options used to solve the pressure problems
135    
136     *(Q^* K Q)p=-Q^*G*
137    
138     If ``options`` is not present, the options are reset to default
139    
140     :param options: `SolverOptions`
141     :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
142     """
143     return self.__pde_p.setSolverOptions(options)
144 gross 1659
145    
146 gross 3071 def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
147     """
148     assigns values to model parameters
149 jfenwick 3360
150 gross 3071 :param f: volumetic sources/sinks
151 gross 3440 :type f: scalar value on the domain (e.g. `escript.Data`)
152 gross 3071 :param g: flux sources/sinks
153 gross 3440 :type g: vector values on the domain (e.g. `escript.Data`)
154 gross 3071 :param location_of_fixed_pressure: mask for locations where pressure is fixed
155 gross 3440 :type location_of_fixed_pressure: scalar value on the domain (e.g. `escript.Data`)
156 gross 3071 :param location_of_fixed_flux: mask for locations where flux is fixed.
157 gross 3440 :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)
158 jfenwick 3360 :param permeability: permeability tensor. If scalar ``s`` is given the tensor with ``s`` on the main diagonal is used.
159 gross 3440 :type permeability: scalar or tensor values on the domain (e.g. `escript.Data`)
160 jfenwick 3360
161 gross 3071 :note: the values of parameters which are not set by calling ``setValue`` are not altered.
162     :note: at any point on the boundary of the domain the pressure (``location_of_fixed_pressure`` >0)
163 jfenwick 3360 or the normal component of the flux (``location_of_fixed_flux[i]>0``) if direction of the normal is along the *x_i* axis.
164    
165 gross 3071 """
166 gross 3440 if self.useVPIteration:
167     if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)
168     if location_of_fixed_flux!=None: self.__pde_k.setValue(q=location_of_fixed_flux)
169     else:
170     q=self.__pde_k.getCoefficient("q")
171     if q.isEmpty(): q=self.__pde_k.createCoefficient("q")
172     if location_of_fixed_pressure!=None: q[self.domain.getDim()]=location_of_fixed_pressure
173     if location_of_fixed_flux!=None: q[:self.domain.getDim()]=location_of_fixed_flux
174     self.__pde_k.setValue(q=q)
175 gross 3071
176 gross 3440 # flux is rescaled by the factor mean value(perm_inv)*length where length**self.domain.getDim()=vol(self.domain)
177 gross 3071 if permeability!=None:
178 gross 3440 perm=util.interpolate(permeability,self.__pde_k.getFunctionSpaceForCoefficient("A"))
179     V=util.vol(self.domain)
180 gross 3071 if perm.getRank()==0:
181 gross 3440 perm_inv=(1./perm)
182     if self.useVPIteration:
183     self.scale=1.
184     else:
185     self.scale=util.integrate(perm_inv)*V**(1./self.domain.getDim()-1.)
186    
187     perm_inv=perm_inv*((1./self.scale)*util.kronecker(self.domain.getDim()))
188     perm=perm*(self.scale*util.kronecker(self.domain.getDim()))
189 gross 3071 elif perm.getRank()==2:
190     perm_inv=util.inverse(perm)
191 gross 3440 if self.useVPIteration:
192     self.scale=1.
193     else:
194     self.scale=util.sqrt(util.integrate(util.length(perm_inv)**2)*V**(2./self.domain.getDim()-1.)/self.domain.getDim())
195     perm_inv*=(1./self.scale)
196     perm=perm*self.scale
197 gross 3071 else:
198     raise ValueError,"illegal rank of permeability."
199 gross 2960
200 gross 3071 self.__permeability=perm
201     self.__permeability_inv=perm_inv
202 gross 3440
203     self.__l2 =(util.longestEdge(self.domain)**2*util.length(self.__permeability_inv))*self.weighting_scale
204     if self.useVPIteration:
205     if self.solveForFlux:
206     self.__pde_k.setValue(D=self.__permeability_inv)
207     else:
208     self.__pde_k.setValue(D=self.__permeability_inv, A=self.__l2*util.outer(util.kronecker(self.domain),util.kronecker(self.domain)))
209     self.__pde_p.setValue(A=self.__permeability)
210     else:
211     D=self.__pde_k.createCoefficient("D")
212     A=self.__pde_k.createCoefficient("A")
213     D[:self.domain.getDim(),:self.domain.getDim()]=self.__permeability_inv
214     for i in range(self.domain.getDim()):
215     for j in range(self.domain.getDim()):
216     A[i,i,j,j]=self.__l2
217     A[self.domain.getDim(),:,self.domain.getDim(),:]=self.__permeability
218     self.__pde_k.setValue(A=A, D=D)
219     if g !=None:
220     g=util.interpolate(g, self.__pde_k.getFunctionSpaceForCoefficient("Y"))
221     if g.isEmpty():
222     g=Vector(0,self.__pde_k.getFunctionSpaceForCoefficient("Y"))
223 gross 3071 else:
224 gross 3440 if not g.getShape()==(self.domain.getDim(),):
225     raise ValueError,"illegal shape of g"
226     self.__g=g
227     elif permeability!=None:
228     X
229     if f !=None:
230     f=util.interpolate(f, self.__pde_k.getFunctionSpaceForCoefficient("X"))
231     if f.isEmpty():
232     f=Scalar(0,self.__pde_k.getFunctionSpaceForCoefficient("X"))
233     else:
234     if f.getRank()>0: raise ValueError,"illegal rank of f."
235     self.__f=f
236 gross 2960
237 gross 3071
238 gross 3440 def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):
239     """
240     solves the problem.
241    
242     The iteration is terminated if the residual norm is less then self.getTolerance().
243 gross 2960
244 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.
245     :type u0: vector value on the domain (e.g. `escript.Data`).
246     :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.
247     :type p0: scalar value on the domain (e.g. `escript.Data`).
248     :param verbose: if set some information on iteration progress are printed
249     :type verbose: ``bool``
250     :return: flux and pressure
251     :rtype: ``tuple`` of `escript.Data`.
252 gross 2960
253 gross 3440 :note: The problem is solved as a least squares form
254     *(K^[-1]+D^* l2 D)u+G p=D^* l2 * f + K^[-1]g*
255     *G^*u+*G^* K Gp=G^*g*
256     where *D* is the *div* operator and *(Gp)_i=p_{,i}* for the permeability *K=k_{ij}*.
257     """
258     self.verbose=verbose
259     if self.useVPIteration:
260     return self.__solveVP(u0,p0,max_iter,max_num_corrections)
261 gross 3071 else:
262 gross 3440 X=self.__pde_k.createCoefficient("X")
263     Y=self.__pde_k.createCoefficient("Y")
264     Y[:self.domain.getDim()]=self.scale*util.tensor_mult(self.__permeability_inv,self.__g)
265     rtmp=self.__f * self.__l2 * self.scale
266     for i in range(self.domain.getDim()): X[i,i]=rtmp
267     X[self.domain.getDim(),:]=self.__g*self.scale
268     r=self.__pde_k.createCoefficient("r")
269     r[:self.domain.getDim()]=u0*self.scale
270     r[self.domain.getDim()]=p0
271     self.__pde_k.setValue(X=X, Y=Y, r=r)
272     self.__pde_k.getSolverOptions().setVerbosity(self.verbose)
273     #self.__pde_k.getSolverOptions().setPreconditioner(self.__pde_k.getSolverOptions().AMG)
274     self.__pde_k.getSolverOptions().setSolverMethod(self.__pde_k.getSolverOptions().DIRECT)
275     U=self.__pde_k.getSolution()
276     # self.__pde_k.getOperator().saveMM("k.mm")
277     u=U[:self.domain.getDim()]*(1./self.scale)
278     p=U[self.domain.getDim()]
279     if self.verbose:
280     KGp=util.tensor_mult(self.__permeability,util.grad(p)/self.scale)
281     def_p=self.__g-(u+KGp)
282     def_v=self.__f-util.div(u, self.__pde_k.getFunctionSpaceForCoefficient("X"))
283     print "DarcyFlux: L2: g-v-K*grad(p) = %e (v = %e)."%(self.__L2(def_p),self.__L2(u))
284     print "DarcyFlux: L2: f-div(v) = %e (grad(v) = %e)."%(self.__L2(def_v),self.__L2(util.grad(u)))
285     return u,p
286 gross 2960
287 gross 3440 def __solveVP(self,u0,p0, max_iter=100, max_num_corrections=10):
288 gross 3071 """
289     solves the problem.
290    
291     The iteration is terminated if the residual norm is less then self.getTolerance().
292 gross 2960
293 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.
294 gross 3440 :type u0: vector value on the domain (e.g. `escript.Data`).
295 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.
296 gross 3440 :type p0: scalar value on the domain (e.g. `escript.Data`).
297 gross 3071 :param verbose: if set some information on iteration progress are printed
298     :type verbose: ``bool``
299     :return: flux and pressure
300 gross 3440 :rtype: ``tuple`` of `escript.Data`.
301 gross 2960
302 gross 3071 :note: The problem is solved as a least squares form
303     *(K^[-1]+D^* (DKD^*)^[-1] D)u+G p=D^* (DKD^*)^[-1] f + K^[-1]g*
304 jfenwick 3360 *G^*u+*G^* K Gp=G^*g*
305 gross 3071 where *D* is the *div* operator and *(Gp)_i=p_{,i}* for the permeability *K=k_{ij}*.
306     """
307     rtol=self.getTolerance()
308     atol=self.getAbsoluteTolerance()
309     self.setSubProblemTolerance()
310     num_corrections=0
311     converged=False
312     norm_r=None
313    
314     # Eliminate the hydrostatic pressure:
315     if self.verbose: print "DarcyFlux: calculate hydrostatic pressure component."
316     self.__pde_p.setValue(X=self.__g, r=p0, y=-util.inner(self.domain.getNormal(),u0))
317     p0=self.__pde_p.getSolution()
318     g2=self.__g - util.tensor_mult(self.__permeability, util.grad(p0))
319     norm_g2=util.integrate(util.inner(g2,util.tensor_mult(self.__permeability_inv,g2)))**0.5
320 gross 2960
321 gross 3071 p=p0*0
322     if self.solveForFlux:
323     v=u0.copy()
324     else:
325     v=self.__getFlux(p, u0, f=self.__f, g=g2)
326 gross 2960
327 gross 3071 while not converged and norm_g2 > 0:
328     Gp=util.grad(p)
329     KGp=util.tensor_mult(self.__permeability,Gp)
330     if self.verbose:
331     def_p=g2-(v+KGp)
332     def_v=self.__f-util.div(v)
333     print "DarcyFlux: L2: g-v-K*grad(p) = %e (v = %e)."%(self.__L2(def_p),self.__L2(v))
334     print "DarcyFlux: L2: f-div(v) = %e (grad(v) = %e)."%(self.__L2(def_v),self.__L2(util.grad(v)))
335     print "DarcyFlux: K^{-1}-norm of v = %e."%util.integrate(util.inner(v,util.tensor_mult(self.__permeability_inv,v)))**0.5
336     print "DarcyFlux: K^{-1}-norm of g2 = %e."%norm_g2
337     print "DarcyFlux: K-norm of grad(dp) = %e."%util.integrate(util.inner(Gp,KGp))**0.5
338     ATOL=atol+rtol*norm_g2
339     if self.verbose: print "DarcyFlux: absolute tolerance ATOL = %e."%(ATOL,)
340     if norm_r == None or norm_r>ATOL:
341     if num_corrections>max_num_corrections:
342     raise ValueError,"maximum number of correction steps reached."
343    
344     if self.solveForFlux:
345     # initial residual is r=K^{-1}*(g-v-K*Gp)+D^*L^{-1}(f-Du)
346     v,r, norm_r=PCG(ArithmeticTuple(util.tensor_mult(self.__permeability_inv,g2-v)-Gp,self.__applWeight(v,self.__f),p),
347     self.__Aprod_v,
348     v,
349     self.__Msolve_PCG_v,
350     self.__inner_PCG_v,
351     atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
352     p=r[2]
353     else:
354     # initial residual is r=G^*(g2-KGp - v)
355     p,r, norm_r=PCG(ArithmeticTuple(g2-KGp,v),
356     self.__Aprod_p,
357     p,
358     self.__Msolve_PCG_p,
359     self.__inner_PCG_p,
360     atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
361     v=r[1]
362     if self.verbose: print "DarcyFlux: residual norm = %e."%norm_r
363     num_corrections+=1
364     else:
365     if self.verbose: print "DarcyFlux: stopping criterium reached."
366     converged=True
367     return v,p+p0
368 gross 3440
369     def __applWeight(self, v, f=None):
370     # solves L p = f-Dv with p = 0
371     if self.verbose: print "DarcyFlux: Applying weighting operator"
372     if f == None:
373     return -util.div(v)*self.__l2
374     else:
375     return (f-util.div(v))*self.__l2
376     def __getPressure(self, v, p0, g=None):
377     # solves (G*KG)p = G^(g-v) with p = p0 where location_of_fixed_pressure>0
378     if self.getSolverOptionsPressure().isVerbose() or self.verbose: print "DarcyFlux: Pressure update"
379     if g == None:
380     self.__pde_p.setValue(X=-v, r=p0)
381     else:
382     self.__pde_p.setValue(X=g-v, r=p0)
383     p=self.__pde_p.getSolution()
384     return p
385    
386     def __Aprod_v(self,dv):
387     # calculates: (a,b,c) = (K^{-1}(dv + KG * dp), L^{-1}Ddv, dp) with (G*KG)dp = - G^*dv
388     dp=self.__getPressure(dv, p0=escript.Data()) # dp = (G*KG)^{-1} (0-G^*dv)
389     a=util.tensor_mult(self.__permeability_inv,dv)+util.grad(dp) # a= K^{-1}u+G*dp
390     b= - self.__applWeight(dv) # b = - (D K D^*)^{-1} (0-Dv)
391     return ArithmeticTuple(a,b,-dp)
392    
393     def __Msolve_PCG_v(self,r):
394     # K^{-1} u = r[0] + D^*r[1] = K^{-1}(dv + KG * dp) + D^*L^{-1}Ddv
395     if self.getSolverOptionsFlux().isVerbose() or self.verbose: print "DarcyFlux: Applying preconditioner"
396     self.__pde_k.setValue(X=r[1]*util.kronecker(self.domain), Y=r[0], r=escript.Data())
397     return self.__pde_k.getSolution()
398    
399     def __inner_PCG_v(self,v,r):
400     return util.integrate(util.inner(v,r[0])+util.div(v)*r[1])
401    
402     def __Aprod_p(self,dp):
403     if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
404     Gdp=util.grad(dp)
405     self.__pde_k.setValue(Y=-Gdp,X=escript.Data(), r=escript.Data())
406     du=self.__pde_k.getSolution()
407     return ArithmeticTuple(util.tensor_mult(self.__permeability,Gdp),-du)
408    
409     def __getFlux(self,p, v0, f=None, g=None):
410     # solves (K^{-1}+D^*L^{-1} D) v = D^*L^{-1}f + K^{-1}g - Gp
411     if f!=None:
412     self.__pde_k.setValue(X=self.__applWeight(v0*0,self.__f)*util.kronecker(self.domain))
413     self.__pde_k.setValue(r=v0)
414     g2=util.tensor_mult(self.__permeability_inv,g)
415     if p == None:
416     self.__pde_k.setValue(Y=g2)
417     else:
418     self.__pde_k.setValue(Y=g2-util.grad(p))
419     return self.__pde_k.getSolution()
420    
421     #v=self.__getFlux(p, u0, f=self.__f, g=g2)
422     def __Msolve_PCG_p(self,r):
423     if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"
424     self.__pde_p.setValue(X=r[0]-r[1], Y=escript.Data(), r=escript.Data(), y=escript.Data())
425     return self.__pde_p.getSolution()
426    
427     def __inner_PCG_p(self,p,r):
428     return util.integrate(util.inner(util.grad(p), r[0]-r[1]))
429    
430     def __L2(self,v):
431     return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))
432    
433     def __L2_r(self,v):
434     return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.ReducedFunction(self.domain)))**2))
435    
436 gross 3071 def setTolerance(self,rtol=1e-4):
437     """
438     sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if
439 gross 2960
440 gross 3071 *|g-v-K gard(p)|_PCG <= atol + rtol * |K^{1/2}g2|_0*
441    
442     where ``atol`` is an absolut tolerance (see `setAbsoluteTolerance`).
443    
444     :param rtol: relative tolerance for the pressure
445     :type rtol: non-negative ``float``
446     """
447     if rtol<0:
448     raise ValueError,"Relative tolerance needs to be non-negative."
449     self.__rtol=rtol
450     def getTolerance(self):
451     """
452     returns the relative tolerance
453     :return: current relative tolerance
454     :rtype: ``float``
455     """
456     return self.__rtol
457 gross 2960
458 gross 3071 def setAbsoluteTolerance(self,atol=0.):
459     """
460     sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if
461    
462     *|g-v-K gard(p)|_PCG <= atol + rtol * |K^{1/2}g2|_0*
463 gross 2960
464    
465 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}*.
466 gross 2960
467 gross 3071 :param atol: absolute tolerance for the pressure
468     :type atol: non-negative ``float``
469     """
470     if atol<0:
471     raise ValueError,"Absolute tolerance needs to be non-negative."
472     self.__atol=atol
473     def getAbsoluteTolerance(self):
474     """
475     returns the absolute tolerance
476     :return: current absolute tolerance
477     :rtype: ``float``
478     """
479     return self.__atol
480     def getSubProblemTolerance(self):
481     """
482     Returns a suitable subtolerance
483     :type: ``float``
484     """
485 gross 3074 return max(util.EPSILON**(0.5),self.getTolerance()**2)
486 gross 2960
487 gross 3071 def setSubProblemTolerance(self):
488     """
489     Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected.
490     """
491     if self.__adaptSubTolerance:
492     sub_tol=self.getSubProblemTolerance()
493     self.getSolverOptionsFlux().setTolerance(sub_tol)
494     self.getSolverOptionsFlux().setAbsoluteTolerance(0.)
495     self.getSolverOptionsPressure().setTolerance(sub_tol)
496     self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
497     if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol
498 gross 2960
499    
500     class DarcyFlowOld(object):
501     """
502     solves the problem
503    
504     *u_i+k_{ij}*p_{,j} = g_i*
505     *u_{i,i} = f*
506    
507     where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,
508    
509     :note: The problem is solved in a least squares formulation.
510     """
511    
512     def __init__(self, domain, weight=None, useReduced=False, adaptSubTolerance=True):
513     """
514     initializes the Darcy flux problem
515     :param domain: domain of the problem
516     :type domain: `Domain`
517     :param useReduced: uses reduced oreder on flux and pressure
518     :type useReduced: ``bool``
519     :param adaptSubTolerance: switches on automatic subtolerance selection
520     :type adaptSubTolerance: ``bool``
521     """
522     self.domain=domain
523     if weight == None:
524 gross 2354 s=self.domain.getSize()
525 gross 3440 self.__l2=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2
526     # self.__l2=(3.*util.longestEdge(self.domain))**2
527     #self.__l2=(0.1*util.longestEdge(self.domain)*s/util.sup(s))**2
528 gross 2351 else:
529 gross 3440 self.__l2=weight
530 gross 2100 self.__pde_v=LinearPDESystem(domain)
531 gross 2208 if useReduced: self.__pde_v.setReducedOrderOn()
532     self.__pde_v.setSymmetryOn()
533 gross 3440 self.__pde_v.setValue(D=util.kronecker(domain), A=self.__l2*util.outer(util.kronecker(domain),util.kronecker(domain)))
534 gross 2100 self.__pde_p=LinearSinglePDE(domain)
535     self.__pde_p.setSymmetryOn()
536 gross 2208 if useReduced: self.__pde_p.setReducedOrderOn()
537 jfenwick 3432 self.__f=escript.Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
538     self.__g=escript.Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
539 gross 2264 self.setTolerance()
540     self.setAbsoluteTolerance()
541 gross 2474 self.__adaptSubTolerance=adaptSubTolerance
542     self.verbose=False
543     def getSolverOptionsFlux(self):
544     """
545     Returns the solver options used to solve the flux problems
546    
547 jfenwick 2625 *(I+D^*D)u=F*
548 gross 2474
549 jfenwick 2625 :return: `SolverOptions`
550 gross 2474 """
551     return self.__pde_v.getSolverOptions()
552     def setSolverOptionsFlux(self, options=None):
553     """
554     Sets the solver options used to solve the flux problems
555    
556 jfenwick 2625 *(I+D^*D)u=F*
557 gross 2474
558 jfenwick 2625 If ``options`` is not present, the options are reset to default
559     :param options: `SolverOptions`
560     :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
561 gross 2474 """
562     return self.__pde_v.setSolverOptions(options)
563     def getSolverOptionsPressure(self):
564     """
565     Returns the solver options used to solve the pressure problems
566    
567 jfenwick 2625 *(Q^*Q)p=Q^*G*
568 gross 2474
569 jfenwick 2625 :return: `SolverOptions`
570 gross 2474 """
571     return self.__pde_p.getSolverOptions()
572     def setSolverOptionsPressure(self, options=None):
573     """
574     Sets the solver options used to solve the pressure problems
575    
576 jfenwick 2625 *(Q^*Q)p=Q^*G*
577 gross 2474
578 jfenwick 2625 If ``options`` is not present, the options are reset to default
579     :param options: `SolverOptions`
580     :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
581 gross 2474 """
582     return self.__pde_p.setSolverOptions(options)
583 gross 1659
584 gross 2100 def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
585     """
586 gross 2208 assigns values to model parameters
587 gross 1659
588 jfenwick 2625 :param f: volumetic sources/sinks
589 gross 3440 :type f: scalar value on the domain (e.g. `escript.Data`)
590 jfenwick 2625 :param g: flux sources/sinks
591 gross 3440 :type g: vector values on the domain (e.g. `escript.Data`)
592 jfenwick 2625 :param location_of_fixed_pressure: mask for locations where pressure is fixed
593 gross 3440 :type location_of_fixed_pressure: scalar value on the domain (e.g. `escript.Data`)
594 jfenwick 2625 :param location_of_fixed_flux: mask for locations where flux is fixed.
595 gross 3440 :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)
596 jfenwick 2625 :param permeability: permeability tensor. If scalar ``s`` is given the tensor with
597     ``s`` on the main diagonal is used. If vector ``v`` is given the tensor with
598     ``v`` on the main diagonal is used.
599 gross 3440 :type permeability: scalar, vector or tensor values on the domain (e.g. `escript.Data`)
600 gross 1659
601 jfenwick 2625 :note: the values of parameters which are not set by calling ``setValue`` are not altered.
602     :note: at any point on the boundary of the domain the pressure (``location_of_fixed_pressure`` >0)
603     or the normal component of the flux (``location_of_fixed_flux[i]>0`` if direction of the normal
604     is along the *x_i* axis.
605 gross 2100 """
606 gross 2264 if f !=None:
607 gross 2100 f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
608     if f.isEmpty():
609 jfenwick 3432 f=escript.Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
610 gross 2100 else:
611     if f.getRank()>0: raise ValueError,"illegal rank of f."
612 gross 2267 self.__f=f
613 gross 2264 if g !=None:
614 gross 2100 g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
615     if g.isEmpty():
616 jfenwick 3432 g=escript.Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
617 gross 2100 else:
618     if not g.getShape()==(self.domain.getDim(),):
619     raise ValueError,"illegal shape of g"
620     self.__g=g
621 gross 1659
622 gross 2100 if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)
623     if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux)
624 gross 1659
625 gross 2100 if permeability!=None:
626     perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
627     if perm.getRank()==0:
628     perm=perm*util.kronecker(self.domain.getDim())
629     elif perm.getRank()==1:
630     perm, perm2=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A")), perm
631     for i in range(self.domain.getDim()): perm[i,i]=perm2[i]
632     elif perm.getRank()==2:
633     pass
634     else:
635     raise ValueError,"illegal rank of permeability."
636     self.__permeability=perm
637     self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))
638 gross 1659
639 gross 2264 def setTolerance(self,rtol=1e-4):
640     """
641 jfenwick 2625 sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if
642 gross 1659
643 jfenwick 2625 *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
644 gross 2264
645 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}*.
646 gross 2264
647 jfenwick 2625 :param rtol: relative tolerance for the pressure
648     :type rtol: non-negative ``float``
649 gross 2100 """
650 gross 2264 if rtol<0:
651     raise ValueError,"Relative tolerance needs to be non-negative."
652     self.__rtol=rtol
653     def getTolerance(self):
654 gross 2100 """
655 gross 2264 returns the relative tolerance
656 gross 1659
657 jfenwick 2625 :return: current relative tolerance
658     :rtype: ``float``
659 caltinay 2169 """
660 gross 2264 return self.__rtol
661 gross 1659
662 gross 2264 def setAbsoluteTolerance(self,atol=0.):
663 gross 2208 """
664 jfenwick 2625 sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if
665 gross 1659
666 jfenwick 2625 *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
667 gross 2208
668 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}*.
669 gross 2208
670 jfenwick 2625 :param atol: absolute tolerance for the pressure
671     :type atol: non-negative ``float``
672 gross 2264 """
673     if atol<0:
674     raise ValueError,"Absolute tolerance needs to be non-negative."
675     self.__atol=atol
676     def getAbsoluteTolerance(self):
677     """
678     returns the absolute tolerance
679    
680 jfenwick 2625 :return: current absolute tolerance
681     :rtype: ``float``
682 gross 2264 """
683     return self.__atol
684     def getSubProblemTolerance(self):
685 gross 2474 """
686     Returns a suitable subtolerance
687 jfenwick 2625 @type: ``float``
688 gross 2474 """
689     return max(util.EPSILON**(0.75),self.getTolerance()**2)
690     def setSubProblemTolerance(self):
691 gross 2264 """
692 gross 2474 Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected.
693 gross 2264 """
694 gross 2474 if self.__adaptSubTolerance:
695     sub_tol=self.getSubProblemTolerance()
696     self.getSolverOptionsFlux().setTolerance(sub_tol)
697     self.getSolverOptionsFlux().setAbsoluteTolerance(0.)
698     self.getSolverOptionsPressure().setTolerance(sub_tol)
699     self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
700     if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol
701 gross 2208
702 gross 2474 def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):
703 gross 2264 """
704     solves the problem.
705 gross 1659
706 gross 2208 The iteration is terminated if the residual norm is less then self.getTolerance().
707 gross 1659
708 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.
709 gross 3440 :type u0: vector value on the domain (e.g. `escript.Data`).
710 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.
711 gross 3440 :type p0: scalar value on the domain (e.g. `escript.Data`).
712 jfenwick 2625 :param verbose: if set some information on iteration progress are printed
713     :type verbose: ``bool``
714     :return: flux and pressure
715 gross 3440 :rtype: ``tuple`` of `escript.Data`.
716 gross 2264
717 jfenwick 2625 :note: The problem is solved as a least squares form
718 gross 2100
719 jfenwick 2625 *(I+D^*D)u+Qp=D^*f+g*
720     *Q^*u+Q^*Qp=Q^*g*
721 gross 2100
722 jfenwick 2625 where *D* is the *div* operator and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.
723 gross 2264 We eliminate the flux form the problem by setting
724 caltinay 2169
725 jfenwick 2625 *u=(I+D^*D)^{-1}(D^*f-g-Qp)* with u=u0 on location_of_fixed_flux
726 caltinay 2169
727 gross 2208 form the first equation. Inserted into the second equation we get
728 caltinay 2169
729 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
730 gross 2264
731 jfenwick 2625 which is solved using the PCG method (precondition is *Q^*Q*). In each iteration step
732     PDEs with operator *I+D^*D* and with *Q^*Q* needs to be solved using a sub iteration scheme.
733 gross 2208 """
734 gross 2386 self.verbose=verbose
735 gross 2264 rtol=self.getTolerance()
736     atol=self.getAbsoluteTolerance()
737 gross 2474 self.setSubProblemTolerance()
738 gross 2264 num_corrections=0
739     converged=False
740     p=p0
741     norm_r=None
742     while not converged:
743 gross 2474 v=self.getFlux(p, fixed_flux=u0)
744 gross 2264 Qp=self.__Q(p)
745     norm_v=self.__L2(v)
746     norm_Qp=self.__L2(Qp)
747     if norm_v == 0.:
748     if norm_Qp == 0.:
749     return v,p
750     else:
751     fac=norm_Qp
752     else:
753     if norm_Qp == 0.:
754     fac=norm_v
755     else:
756     fac=2./(1./norm_v+1./norm_Qp)
757     ATOL=(atol+rtol*fac)
758     if self.verbose:
759     print "DarcyFlux: L2 norm of v = %e."%norm_v
760 gross 2960 print "DarcyFlux: L2 norm of k.util.grad(p) = %e."%norm_Qp
761 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),)
762 gross 2486 print "DarcyFlux: L2 defect div(v) = %e."%(util.integrate((self.__f-util.div(v))**2)**(0.5),)
763 gross 2264 print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
764     if norm_r == None or norm_r>ATOL:
765     if num_corrections>max_num_corrections:
766     raise ValueError,"maximum number of correction steps reached."
767 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)
768 gross 2264 num_corrections+=1
769     else:
770     converged=True
771     return v,p
772     def __L2(self,v):
773 jfenwick 3432 return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))
774 gross 2264
775     def __Q(self,p):
776     return util.tensor_mult(self.__permeability,util.grad(p))
777    
778     def __Aprod(self,dp):
779 gross 2474 if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
780 gross 2264 Qdp=self.__Q(dp)
781 jfenwick 3432 self.__pde_v.setValue(Y=-Qdp,X=escript.Data(), r=escript.Data())
782 gross 2474 du=self.__pde_v.getSolution()
783 gross 2264 return Qdp+du
784     def __inner_GMRES(self,r,s):
785     return util.integrate(util.inner(r,s))
786    
787 gross 2100 def __inner_PCG(self,p,r):
788 gross 2264 return util.integrate(util.inner(self.__Q(p), r))
789 gross 2100
790     def __Msolve_PCG(self,r):
791 gross 2474 if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"
792 jfenwick 3432 self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=escript.Data(), r=escript.Data())
793 gross 2474 return self.__pde_p.getSolution()
794 gross 2100
795 jfenwick 3432 def getFlux(self,p=None, fixed_flux=escript.Data()):
796 gross 2264 """
797 jfenwick 2625 returns the flux for a given pressure ``p`` where the flux is equal to ``fixed_flux``
798     on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
799     Note that ``g`` and ``f`` are used, see `setValue`.
800 gross 2264
801 jfenwick 2625 :param p: pressure.
802 gross 3440 :type p: scalar value on the domain (e.g. `escript.Data`).
803 jfenwick 2625 :param fixed_flux: flux on the locations of the domain marked be ``location_of_fixed_flux``.
804 gross 3440 :type fixed_flux: vector values on the domain (e.g. `escript.Data`).
805 jfenwick 2625 :return: flux
806 gross 3440 :rtype: `escript.Data`
807 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}*
808     for the permeability *k_{ij}*
809 gross 2264 """
810 gross 2474 self.setSubProblemTolerance()
811 gross 2264 g=self.__g
812     f=self.__f
813 gross 3440 self.__pde_v.setValue(X=self.__l2*f*util.kronecker(self.domain), r=fixed_flux)
814 gross 2264 if p == None:
815     self.__pde_v.setValue(Y=g)
816     else:
817     self.__pde_v.setValue(Y=g-self.__Q(p))
818 gross 2474 return self.__pde_v.getSolution()
819 gross 2264
820 gross 1414 class StokesProblemCartesian(HomogeneousSaddlePointProblem):
821 gross 2251 """
822 gross 2264 solves
823 gross 1414
824 gross 2208 -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
825     u_{i,i}=0
826 gross 1414
827 gross 2208 u=0 where fixed_u_mask>0
828     eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j
829 gross 1414
830 gross 2264 if surface_stress is not given 0 is assumed.
831 gross 1414
832 gross 2251 typical usage:
833 gross 1414
834 gross 2208 sp=StokesProblemCartesian(domain)
835     sp.setTolerance()
836     sp.initialize(...)
837     v,p=sp.solve(v0,p0)
838 gross 2251 """
839 gross 2719 def __init__(self,domain,**kwargs):
840 gross 2100 """
841 gross 2208 initialize the Stokes Problem
842 gross 2100
843 gross 2793 The approximation spaces used for velocity (=Solution(domain)) and pressure (=ReducedSolution(domain)) must be
844     LBB complient, for instance using quadratic and linear approximation on the same element or using linear approximation
845     with macro elements for the pressure.
846    
847     :param domain: domain of the problem.
848 jfenwick 2625 :type domain: `Domain`
849 gross 2100 """
850 gross 2719 HomogeneousSaddlePointProblem.__init__(self,**kwargs)
851 gross 1414 self.domain=domain
852     self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
853     self.__pde_u.setSymmetryOn()
854 gross 2474
855 gross 1414 self.__pde_prec=LinearPDE(domain)
856     self.__pde_prec.setReducedOrderOn()
857     self.__pde_prec.setSymmetryOn()
858    
859 gross 2415 self.__pde_proj=LinearPDE(domain)
860     self.__pde_proj.setReducedOrderOn()
861     self.__pde_proj.setValue(D=1)
862     self.__pde_proj.setSymmetryOn()
863    
864 gross 2474 def getSolverOptionsVelocity(self):
865     """
866     returns the solver options used solve the equation for velocity.
867    
868 jfenwick 2625 :rtype: `SolverOptions`
869 gross 2474 """
870     return self.__pde_u.getSolverOptions()
871     def setSolverOptionsVelocity(self, options=None):
872     """
873     set the solver options for solving the equation for velocity.
874    
875 jfenwick 2625 :param options: new solver options
876     :type options: `SolverOptions`
877 gross 2474 """
878     self.__pde_u.setSolverOptions(options)
879     def getSolverOptionsPressure(self):
880     """
881     returns the solver options used solve the equation for pressure.
882 jfenwick 2625 :rtype: `SolverOptions`
883 gross 2474 """
884     return self.__pde_prec.getSolverOptions()
885     def setSolverOptionsPressure(self, options=None):
886     """
887     set the solver options for solving the equation for pressure.
888 jfenwick 2625 :param options: new solver options
889     :type options: `SolverOptions`
890 gross 2474 """
891     self.__pde_prec.setSolverOptions(options)
892 gross 2415
893 gross 2474 def setSolverOptionsDiv(self, options=None):
894     """
895     set the solver options for solving the equation to project the divergence of
896     the velocity onto the function space of presure.
897    
898 jfenwick 2625 :param options: new solver options
899     :type options: `SolverOptions`
900 gross 2474 """
901 artak 2689 self.__pde_proj.setSolverOptions(options)
902 gross 2474 def getSolverOptionsDiv(self):
903     """
904     returns the solver options for solving the equation to project the divergence of
905     the velocity onto the function space of presure.
906    
907 jfenwick 2625 :rtype: `SolverOptions`
908 gross 2474 """
909 artak 2689 return self.__pde_proj.getSolverOptions()
910 gross 2474
911 gross 2793 def updateStokesEquation(self, v, p):
912     """
913     updates the Stokes equation to consider dependencies from ``v`` and ``p``
914     :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values.
915     """
916     pass
917     def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None):
918     """
919     assigns new values to the model parameters.
920    
921     :param f: external force
922     :type f: `Vector` object in `FunctionSpace` `Function` or similar
923     :param fixed_u_mask: mask of locations with fixed velocity.
924     :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
925     :param eta: viscosity
926     :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
927     :param surface_stress: normal surface stress
928     :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
929     :param stress: initial stress
930     :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
931     """
932     if eta !=None:
933     k=util.kronecker(self.domain.getDim())
934     kk=util.outer(k,k)
935 jfenwick 3432 self.eta=util.interpolate(eta, escript.Function(self.domain))
936 gross 2793 self.__pde_prec.setValue(D=1/self.eta)
937     self.__pde_u.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))
938     if restoration_factor!=None:
939     n=self.domain.getNormal()
940     self.__pde_u.setValue(d=restoration_factor*util.outer(n,n))
941     if fixed_u_mask!=None:
942     self.__pde_u.setValue(q=fixed_u_mask)
943     if f!=None: self.__f=f
944     if surface_stress!=None: self.__surface_stress=surface_stress
945     if stress!=None: self.__stress=stress
946    
947 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):
948 gross 2208 """
949     assigns values to the model parameters
950 gross 2100
951 jfenwick 2625 :param f: external force
952     :type f: `Vector` object in `FunctionSpace` `Function` or similar
953     :param fixed_u_mask: mask of locations with fixed velocity.
954     :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
955     :param eta: viscosity
956     :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
957     :param surface_stress: normal surface stress
958     :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
959     :param stress: initial stress
960     :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
961 gross 2208 """
962 gross 2793 self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
963 gross 1414
964 gross 2719 def Bv(self,v,tol):
965 gross 2251 """
966     returns inner product of element p and div(v)
967 gross 1414
968 jfenwick 2625 :param v: a residual
969     :return: inner product of element p and div(v)
970     :rtype: ``float``
971 gross 2100 """
972 gross 2793 self.__pde_proj.setValue(Y=-util.div(v))
973 gross 2719 self.getSolverOptionsDiv().setTolerance(tol)
974     self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
975     out=self.__pde_proj.getSolution()
976     return out
977 gross 2208
978 gross 2445 def inner_pBv(self,p,Bv):
979     """
980     returns inner product of element p and Bv=-div(v)
981    
982 jfenwick 2625 :param p: a pressure increment
983     :param Bv: a residual
984     :return: inner product of element p and Bv=-div(v)
985     :rtype: ``float``
986 gross 2445 """
987 jfenwick 3432 return util.integrate(util.interpolate(p,escript.Function(self.domain))*util.interpolate(Bv, escript.Function(self.domain)))
988 gross 2445
989 gross 2251 def inner_p(self,p0,p1):
990 gross 2100 """
991 gross 2251 Returns inner product of p0 and p1
992 gross 1414
993 jfenwick 2625 :param p0: a pressure
994     :param p1: a pressure
995     :return: inner product of p0 and p1
996     :rtype: ``float``
997 gross 2100 """
998 jfenwick 3432 s0=util.interpolate(p0, escript.Function(self.domain))
999     s1=util.interpolate(p1, escript.Function(self.domain))
1000 gross 2100 return util.integrate(s0*s1)
1001 artak 1550
1002 gross 2251 def norm_v(self,v):
1003 gross 2100 """
1004 gross 2251 returns the norm of v
1005 gross 2208
1006 jfenwick 2625 :param v: a velovity
1007     :return: norm of v
1008     :rtype: non-negative ``float``
1009 gross 2100 """
1010 gross 2719 return util.sqrt(util.integrate(util.length(util.grad(v))**2))
1011 gross 2100
1012 gross 2793
1013 gross 2719 def getDV(self, p, v, tol):
1014 gross 1414 """
1015 gross 2251 return the value for v for a given p (overwrite)
1016    
1017 jfenwick 2625 :param p: a pressure
1018 gross 2719 :param v: a initial guess for the value v to return.
1019     :return: dv given as *Adv=(f-Av-B^*p)*
1020 gross 1414 """
1021 gross 2793 self.updateStokesEquation(v,p)
1022 gross 2719 self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress)
1023     self.getSolverOptionsVelocity().setTolerance(tol)
1024     self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
1025 gross 2100 if self.__stress.isEmpty():
1026 gross 2719 self.__pde_u.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
1027 gross 2100 else:
1028 gross 2719 self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
1029 gross 2474 out=self.__pde_u.getSolution()
1030 gross 2208 return out
1031 gross 1414
1032 gross 2445 def norm_Bv(self,Bv):
1033 gross 2251 """
1034     Returns Bv (overwrite).
1035    
1036 jfenwick 2625 :rtype: equal to the type of p
1037     :note: boundary conditions on p should be zero!
1038 gross 2251 """
1039 jfenwick 3432 return util.sqrt(util.integrate(util.interpolate(Bv, escript.Function(self.domain))**2))
1040 gross 2251
1041 gross 2719 def solve_AinvBt(self,p, tol):
1042 gross 2251 """
1043 gross 2719 Solves *Av=B^*p* with accuracy `tol`
1044 gross 2251
1045 jfenwick 2625 :param p: a pressure increment
1046     :return: the solution of *Av=B^*p*
1047     :note: boundary conditions on v should be zero!
1048 gross 2251 """
1049 jfenwick 3432 self.__pde_u.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))
1050 gross 2474 out=self.__pde_u.getSolution()
1051 gross 2251 return out
1052    
1053 gross 2719 def solve_prec(self,Bv, tol):
1054 gross 2251 """
1055 jfenwick 2625 applies preconditioner for for *BA^{-1}B^** to *Bv*
1056     with accuracy `self.getSubProblemTolerance()`
1057 gross 2251
1058 jfenwick 2625 :param Bv: velocity increment
1059     :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*
1060     :note: boundary conditions on p are zero.
1061 gross 2251 """
1062 gross 2445 self.__pde_prec.setValue(Y=Bv)
1063 gross 2719 self.getSolverOptionsPressure().setTolerance(tol)
1064     self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
1065     out=self.__pde_prec.getSolution()
1066     return out

  ViewVC Help
Powered by ViewVC 1.1.26