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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3462 - (hide annotations)
Mon Feb 7 02:01:50 2011 UTC (8 years, 6 months ago) by caltinay
Original Path: trunk/escript/py_src/flows.py
File MIME type: text/x-python
File size: 44083 byte(s)
Removed debug output.

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

  ViewVC Help
Powered by ViewVC 1.1.26