/[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 2474 - (hide annotations)
Tue Jun 16 06:32:15 2009 UTC (10 years, 3 months ago) by gross
File MIME type: text/x-python
File size: 22893 byte(s)
linearPDEs has is now using the SolverOptions class to talk to PASO
1 ksteube 1809 ########################################################
2 gross 1414 #
3 ksteube 1809 # Copyright (c) 2003-2008 by University of Queensland
4     # Earth Systems Science Computational Center (ESSCC)
5     # http://www.uq.edu.au/esscc
6 gross 1414 #
7 ksteube 1809 # Primary Business: Queensland, Australia
8     # Licensed under the Open Software License version 3.0
9     # http://www.opensource.org/licenses/osl-3.0.php
10 gross 1414 #
11 ksteube 1809 ########################################################
12 gross 1414
13 ksteube 1809 __copyright__="""Copyright (c) 2003-2008 by University of Queensland
14     Earth Systems Science Computational Center (ESSCC)
15     http://www.uq.edu.au/esscc
16     Primary Business: Queensland, Australia"""
17     __license__="""Licensed under the Open Software License version 3.0
18     http://www.opensource.org/licenses/osl-3.0.php"""
19 jfenwick 2344 __url__="https://launchpad.net/escript-finley"
20 ksteube 1809
21 gross 1414 """
22     Some models for flow
23    
24     @var __author__: name of author
25     @var __copyright__: copyrights
26     @var __license__: licence agreement
27     @var __url__: url entry point on documentation
28     @var __version__: version
29     @var __date__: date of the version
30     """
31    
32     __author__="Lutz Gross, l.gross@uq.edu.au"
33    
34     from escript import *
35     import util
36 gross 2474 from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions
37 gross 2264 from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
38 gross 1414
39 gross 2100 class DarcyFlow(object):
40     """
41 gross 2264 solves the problem
42 gross 1659
43 gross 2100 M{u_i+k_{ij}*p_{,j} = g_i}
44     M{u_{i,i} = f}
45 gross 1659
46 gross 2264 where M{p} represents the pressure and M{u} the Darcy flux. M{k} represents the permeability,
47 gross 1659
48 gross 2100 @note: The problem is solved in a least squares formulation.
49     """
50 gross 1659
51 gross 2474 def __init__(self, domain, weight=None, useReduced=False, adaptSubTolerance=True):
52 gross 2100 """
53 gross 2208 initializes the Darcy flux problem
54 gross 2100 @param domain: domain of the problem
55     @type domain: L{Domain}
56 gross 2474 @param useReduced: uses reduced oreder on flux and pressure
57     @type useReduced: C{bool}
58     @param adaptSubTolerance: switches on automatic subtolerance selection
59     @type adaptSubTolerance: C{bool}
60 gross 2100 """
61     self.domain=domain
62 gross 2351 if weight == None:
63 gross 2354 s=self.domain.getSize()
64     self.__l=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2
65 gross 2351 else:
66     self.__l=weight
67 gross 2100 self.__pde_v=LinearPDESystem(domain)
68 gross 2208 if useReduced: self.__pde_v.setReducedOrderOn()
69     self.__pde_v.setSymmetryOn()
70 gross 2267 self.__pde_v.setValue(D=util.kronecker(domain), A=self.__l*util.outer(util.kronecker(domain),util.kronecker(domain)))
71 gross 2100 self.__pde_p=LinearSinglePDE(domain)
72     self.__pde_p.setSymmetryOn()
73 gross 2208 if useReduced: self.__pde_p.setReducedOrderOn()
74 gross 2100 self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
75     self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
76 gross 2264 self.setTolerance()
77     self.setAbsoluteTolerance()
78 gross 2474 self.__adaptSubTolerance=adaptSubTolerance
79     self.verbose=False
80     def getSolverOptionsFlux(self):
81     """
82     Returns the solver options used to solve the flux problems
83    
84     M{(I+D^*D)u=F}
85    
86     @return: L{SolverOptions}
87     """
88     return self.__pde_v.getSolverOptions()
89     def setSolverOptionsFlux(self, options=None):
90     """
91     Sets the solver options used to solve the flux problems
92    
93     M{(I+D^*D)u=F}
94    
95     If C{options} is not present, the options are reset to default
96     @param options: L{SolverOptions}
97     @note: if the adaption of subtolerance is choosen, the tolerance set by C{options} will be overwritten before the solver is called.
98     """
99     return self.__pde_v.setSolverOptions(options)
100     def getSolverOptionsPressure(self):
101     """
102     Returns the solver options used to solve the pressure problems
103    
104     M{(Q^*Q)p=Q^*G}
105    
106     @return: L{SolverOptions}
107     """
108     return self.__pde_p.getSolverOptions()
109     def setSolverOptionsPressure(self, options=None):
110     """
111     Sets the solver options used to solve the pressure problems
112    
113     M{(Q^*Q)p=Q^*G}
114    
115     If C{options} is not present, the options are reset to default
116     @param options: L{SolverOptions}
117     @note: if the adaption of subtolerance is choosen, the tolerance set by C{options} will be overwritten before the solver is called.
118     """
119     return self.__pde_p.setSolverOptions(options)
120 gross 1659
121 gross 2100 def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
122     """
123 gross 2208 assigns values to model parameters
124 gross 1659
125 gross 2208 @param f: volumetic sources/sinks
126     @type f: scalar value on the domain (e.g. L{Data})
127 gross 2100 @param g: flux sources/sinks
128 gross 2208 @type g: vector values on the domain (e.g. L{Data})
129 gross 2100 @param location_of_fixed_pressure: mask for locations where pressure is fixed
130 gross 2208 @type location_of_fixed_pressure: scalar value on the domain (e.g. L{Data})
131     @param location_of_fixed_flux: mask for locations where flux is fixed.
132     @type location_of_fixed_flux: vector values on the domain (e.g. L{Data})
133 gross 2264 @param permeability: permeability tensor. If scalar C{s} is given the tensor with
134     C{s} on the main diagonal is used. If vector C{v} is given the tensor with
135 gross 2208 C{v} on the main diagonal is used.
136     @type permeability: scalar, vector or tensor values on the domain (e.g. L{Data})
137 gross 1659
138 gross 2208 @note: the values of parameters which are not set by calling C{setValue} are not altered.
139     @note: at any point on the boundary of the domain the pressure (C{location_of_fixed_pressure} >0)
140     or the normal component of the flux (C{location_of_fixed_flux[i]>0} if direction of the normal
141     is along the M{x_i} axis.
142 gross 2100 """
143 gross 2264 if f !=None:
144 gross 2100 f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
145     if f.isEmpty():
146     f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
147     else:
148     if f.getRank()>0: raise ValueError,"illegal rank of f."
149 gross 2267 self.__f=f
150 gross 2264 if g !=None:
151 gross 2100 g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
152     if g.isEmpty():
153     g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
154     else:
155     if not g.getShape()==(self.domain.getDim(),):
156     raise ValueError,"illegal shape of g"
157     self.__g=g
158 gross 1659
159 gross 2100 if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)
160     if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux)
161 gross 1659
162 gross 2100 if permeability!=None:
163     perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
164     if perm.getRank()==0:
165     perm=perm*util.kronecker(self.domain.getDim())
166     elif perm.getRank()==1:
167     perm, perm2=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A")), perm
168     for i in range(self.domain.getDim()): perm[i,i]=perm2[i]
169     elif perm.getRank()==2:
170     pass
171     else:
172     raise ValueError,"illegal rank of permeability."
173     self.__permeability=perm
174     self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))
175 gross 1659
176 gross 2264 def setTolerance(self,rtol=1e-4):
177     """
178     sets the relative tolerance C{rtol} used to terminate the solution process. The iteration is terminated if
179 gross 1659
180 gross 2264 M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }
181    
182     where C{atol} is an absolut tolerance (see L{setAbsoluteTolerance}), M{|f|^2 = integrate(length(f)^2)} and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.
183    
184     @param rtol: relative tolerance for the pressure
185     @type rtol: non-negative C{float}
186 gross 2100 """
187 gross 2264 if rtol<0:
188     raise ValueError,"Relative tolerance needs to be non-negative."
189     self.__rtol=rtol
190     def getTolerance(self):
191 gross 2100 """
192 gross 2264 returns the relative tolerance
193 gross 1659
194 gross 2264 @return: current relative tolerance
195     @rtype: C{float}
196 caltinay 2169 """
197 gross 2264 return self.__rtol
198 gross 1659
199 gross 2264 def setAbsoluteTolerance(self,atol=0.):
200 gross 2208 """
201 gross 2264 sets the absolute tolerance C{atol} used to terminate the solution process. The iteration is terminated if
202 gross 1659
203 gross 2264 M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }
204 gross 2208
205 gross 2264 where C{rtol} is an absolut tolerance (see L{setTolerance}), M{|f|^2 = integrate(length(f)^2)} and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.
206 gross 2208
207 gross 2264 @param atol: absolute tolerance for the pressure
208     @type atol: non-negative C{float}
209     """
210     if atol<0:
211     raise ValueError,"Absolute tolerance needs to be non-negative."
212     self.__atol=atol
213     def getAbsoluteTolerance(self):
214     """
215     returns the absolute tolerance
216    
217     @return: current absolute tolerance
218     @rtype: C{float}
219     """
220     return self.__atol
221     def getSubProblemTolerance(self):
222 gross 2474 """
223     Returns a suitable subtolerance
224     @type: C{float}
225     """
226     return max(util.EPSILON**(0.75),self.getTolerance()**2)
227     def setSubProblemTolerance(self):
228 gross 2264 """
229 gross 2474 Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected.
230 gross 2264 """
231 gross 2474 if self.__adaptSubTolerance:
232     sub_tol=self.getSubProblemTolerance()
233     self.getSolverOptionsFlux().setTolerance(sub_tol)
234     self.getSolverOptionsFlux().setAbsoluteTolerance(0.)
235     self.getSolverOptionsPressure().setTolerance(sub_tol)
236     self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
237     if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol
238 gross 2208
239 gross 2474 def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):
240 gross 2264 """
241     solves the problem.
242 gross 1659
243 gross 2208 The iteration is terminated if the residual norm is less then self.getTolerance().
244 gross 1659
245 gross 2208 @param u0: initial guess for the flux. At locations in the domain marked by C{location_of_fixed_flux} the value of C{u0} is kept unchanged.
246     @type u0: vector value on the domain (e.g. L{Data}).
247     @param p0: initial guess for the pressure. At locations in the domain marked by C{location_of_fixed_pressure} the value of C{p0} is kept unchanged.
248     @type p0: scalar value on the domain (e.g. L{Data}).
249     @param verbose: if set some information on iteration progress are printed
250     @type verbose: C{bool}
251     @return: flux and pressure
252     @rtype: C{tuple} of L{Data}.
253 gross 2264
254 gross 2208 @note: The problem is solved as a least squares form
255 gross 2100
256 gross 2208 M{(I+D^*D)u+Qp=D^*f+g}
257     M{Q^*u+Q^*Qp=Q^*g}
258 gross 2100
259 gross 2264 where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.
260     We eliminate the flux form the problem by setting
261 caltinay 2169
262 gross 2208 M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} with u=u0 on location_of_fixed_flux
263 caltinay 2169
264 gross 2208 form the first equation. Inserted into the second equation we get
265 caltinay 2169
266 gross 2208 M{Q^*(I-(I+D^*D)^{-1})Qp= Q^*(g-(I+D^*D)^{-1}(D^*f+g))} with p=p0 on location_of_fixed_pressure
267 gross 2264
268     which is solved using the PCG method (precondition is M{Q^*Q}). In each iteration step
269 gross 2208 PDEs with operator M{I+D^*D} and with M{Q^*Q} needs to be solved using a sub iteration scheme.
270     """
271 gross 2386 self.verbose=verbose
272 gross 2264 rtol=self.getTolerance()
273     atol=self.getAbsoluteTolerance()
274 gross 2474 self.setSubProblemTolerance()
275    
276 gross 2264 num_corrections=0
277     converged=False
278     p=p0
279     norm_r=None
280     while not converged:
281 gross 2474 v=self.getFlux(p, fixed_flux=u0)
282 gross 2264 Qp=self.__Q(p)
283     norm_v=self.__L2(v)
284     norm_Qp=self.__L2(Qp)
285     if norm_v == 0.:
286     if norm_Qp == 0.:
287     return v,p
288     else:
289     fac=norm_Qp
290     else:
291     if norm_Qp == 0.:
292     fac=norm_v
293     else:
294     fac=2./(1./norm_v+1./norm_Qp)
295     ATOL=(atol+rtol*fac)
296     if self.verbose:
297     print "DarcyFlux: L2 norm of v = %e."%norm_v
298     print "DarcyFlux: L2 norm of k.grad(p) = %e."%norm_Qp
299     print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
300     if norm_r == None or norm_r>ATOL:
301     if num_corrections>max_num_corrections:
302     raise ValueError,"maximum number of correction steps reached."
303 gross 2354 p,r, norm_r=PCG(self.__g-util.interpolate(v,Function(self.domain))-Qp,self.__Aprod,p,self.__Msolve_PCG,self.__inner_PCG,atol=0.5*ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
304 gross 2264 num_corrections+=1
305     else:
306     converged=True
307     return v,p
308     def __L2(self,v):
309     return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))
310    
311     def __Q(self,p):
312     return util.tensor_mult(self.__permeability,util.grad(p))
313    
314     def __Aprod(self,dp):
315 gross 2474 if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
316 gross 2264 Qdp=self.__Q(dp)
317     self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())
318 gross 2474 du=self.__pde_v.getSolution()
319 gross 2370 # self.__pde_v.getOperator().saveMM("proj.mm")
320 gross 2264 return Qdp+du
321     def __inner_GMRES(self,r,s):
322     return util.integrate(util.inner(r,s))
323    
324 gross 2100 def __inner_PCG(self,p,r):
325 gross 2264 return util.integrate(util.inner(self.__Q(p), r))
326 gross 2100
327     def __Msolve_PCG(self,r):
328 gross 2474 if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"
329 gross 2264 self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=Data(), r=Data())
330 gross 2370 # self.__pde_p.getOperator().saveMM("prec.mm")
331 gross 2474 return self.__pde_p.getSolution()
332 gross 2100
333 gross 2474 def getFlux(self,p=None, fixed_flux=Data()):
334 gross 2264 """
335     returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux}
336     on locations where C{location_of_fixed_flux} is positive (see L{setValue}).
337     Note that C{g} and C{f} are used, see L{setValue}.
338    
339     @param p: pressure.
340     @type p: scalar value on the domain (e.g. L{Data}).
341     @param fixed_flux: flux on the locations of the domain marked be C{location_of_fixed_flux}.
342     @type fixed_flux: vector values on the domain (e.g. L{Data}).
343     @param tol: relative tolerance to be used.
344     @type tol: positive C{float}.
345     @return: flux
346     @rtype: L{Data}
347     @note: the method uses the least squares solution M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}}
348     for the permeability M{k_{ij}}
349     """
350 gross 2474 self.setSubProblemTolerance()
351 gross 2264 g=self.__g
352     f=self.__f
353 gross 2267 self.__pde_v.setValue(X=self.__l*f*util.kronecker(self.domain), r=fixed_flux)
354 gross 2264 if p == None:
355     self.__pde_v.setValue(Y=g)
356     else:
357     self.__pde_v.setValue(Y=g-self.__Q(p))
358 gross 2474 return self.__pde_v.getSolution()
359 gross 2264
360 gross 1414 class StokesProblemCartesian(HomogeneousSaddlePointProblem):
361 gross 2251 """
362 gross 2264 solves
363 gross 1414
364 gross 2208 -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
365     u_{i,i}=0
366 gross 1414
367 gross 2208 u=0 where fixed_u_mask>0
368     eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j
369 gross 1414
370 gross 2264 if surface_stress is not given 0 is assumed.
371 gross 1414
372 gross 2251 typical usage:
373 gross 1414
374 gross 2208 sp=StokesProblemCartesian(domain)
375     sp.setTolerance()
376     sp.initialize(...)
377     v,p=sp.solve(v0,p0)
378 gross 2251 """
379 gross 2474 def __init__(self,domain,adaptSubTolerance=True, **kwargs):
380 gross 2100 """
381 gross 2208 initialize the Stokes Problem
382 gross 2100
383 gross 2208 @param domain: domain of the problem. The approximation order needs to be two.
384 gross 2100 @type domain: L{Domain}
385 gross 2474 @param adaptSubTolerance: If True the tolerance for subproblem is set automatically.
386     @type adaptSubTolerance: C{bool}
387 gross 2208 @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.
388 gross 2100 """
389 gross 2474 HomogeneousSaddlePointProblem.__init__(self,adaptSubTolerance=adaptSubTolerance,**kwargs)
390 gross 1414 self.domain=domain
391     self.vol=util.integrate(1.,Function(self.domain))
392     self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
393     self.__pde_u.setSymmetryOn()
394 gross 2474
395 gross 1414 self.__pde_prec=LinearPDE(domain)
396     self.__pde_prec.setReducedOrderOn()
397     self.__pde_prec.setSymmetryOn()
398    
399 gross 2415 self.__pde_proj=LinearPDE(domain)
400     self.__pde_proj.setReducedOrderOn()
401     self.__pde_proj.setValue(D=1)
402     self.__pde_proj.setSymmetryOn()
403    
404 gross 2474 def getSolverOptionsVelocity(self):
405     """
406     returns the solver options used solve the equation for velocity.
407    
408     @rtype: L{SolverOptions}
409     """
410     return self.__pde_u.getSolverOptions()
411     def setSolverOptionsVelocity(self, options=None):
412     """
413     set the solver options for solving the equation for velocity.
414    
415     @param options: new solver options
416     @type options: L{SolverOptions}
417     """
418     self.__pde_u.setSolverOptions(options)
419     def getSolverOptionsPressure(self):
420     """
421     returns the solver options used solve the equation for pressure.
422     @rtype: L{SolverOptions}
423     """
424     return self.__pde_prec.getSolverOptions()
425     def setSolverOptionsPressure(self, options=None):
426     """
427     set the solver options for solving the equation for pressure.
428     @param options: new solver options
429     @type options: L{SolverOptions}
430     """
431     self.__pde_prec.setSolverOptions(options)
432 gross 2415
433 gross 2474 def setSolverOptionsDiv(self, options=None):
434     """
435     set the solver options for solving the equation to project the divergence of
436     the velocity onto the function space of presure.
437    
438     @param options: new solver options
439     @type options: L{SolverOptions}
440     """
441     self.__pde_prec.setSolverOptions(options)
442     def getSolverOptionsDiv(self):
443     """
444     returns the solver options for solving the equation to project the divergence of
445     the velocity onto the function space of presure.
446    
447     @rtype: L{SolverOptions}
448     """
449     return self.__pde_prec.getSolverOptions()
450     def setSubProblemTolerance(self):
451     """
452     Updates the tolerance for subproblems
453     """
454     if self.adaptSubTolerance():
455     sub_tol=self.getSubProblemTolerance()
456     self.getSolverOptionsDiv().setTolerance(sub_tol)
457     self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
458     self.getSolverOptionsPressure().setTolerance(sub_tol)
459     self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
460     self.getSolverOptionsVelocity().setTolerance(sub_tol)
461     self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
462    
463    
464 gross 2251 def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):
465 gross 2208 """
466     assigns values to the model parameters
467 gross 2100
468 gross 2208 @param f: external force
469     @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar
470     @param fixed_u_mask: mask of locations with fixed velocity.
471     @type fixed_u_mask: L{Vector} object on L{FunctionSpace} L{Solution} or similar
472     @param eta: viscosity
473     @type eta: L{Scalar} object on L{FunctionSpace} L{Function} or similar
474     @param surface_stress: normal surface stress
475     @type eta: L{Vector} object on L{FunctionSpace} L{FunctionOnBoundary} or similar
476     @param stress: initial stress
477     @type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar
478     @note: All values needs to be set.
479     """
480     self.eta=eta
481     A =self.__pde_u.createCoefficient("A")
482     self.__pde_u.setValue(A=Data())
483     for i in range(self.domain.getDim()):
484     for j in range(self.domain.getDim()):
485 gross 2264 A[i,j,j,i] += 1.
486 gross 2208 A[i,j,i,j] += 1.
487 gross 2264 self.__pde_prec.setValue(D=1/self.eta)
488 gross 2251 self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask)
489     self.__f=f
490     self.__surface_stress=surface_stress
491 gross 2208 self.__stress=stress
492 gross 1414
493 gross 2445 def Bv(self,v):
494 gross 2251 """
495     returns inner product of element p and div(v)
496 gross 1414
497 gross 2251 @param p: a pressure increment
498     @param v: a residual
499     @return: inner product of element p and div(v)
500     @rtype: C{float}
501 gross 2100 """
502 gross 2445 self.__pde_proj.setValue(Y=-util.div(v))
503     return self.__pde_proj.getSolution()
504 gross 2208
505 gross 2445 def inner_pBv(self,p,Bv):
506     """
507     returns inner product of element p and Bv=-div(v)
508    
509     @param p: a pressure increment
510     @param v: a residual
511     @return: inner product of element p and Bv=-div(v)
512     @rtype: C{float}
513     """
514     return util.integrate(util.interpolate(p,Function(self.domain))*util.interpolate(Bv,Function(self.domain)))
515    
516 gross 2251 def inner_p(self,p0,p1):
517 gross 2100 """
518 gross 2251 Returns inner product of p0 and p1
519 gross 1414
520 gross 2251 @param p0: a pressure
521     @param p1: a pressure
522     @return: inner product of p0 and p1
523 gross 2208 @rtype: C{float}
524 gross 2100 """
525     s0=util.interpolate(p0/self.eta,Function(self.domain))
526     s1=util.interpolate(p1/self.eta,Function(self.domain))
527     return util.integrate(s0*s1)
528 artak 1550
529 gross 2251 def norm_v(self,v):
530 gross 2100 """
531 gross 2251 returns the norm of v
532 gross 2208
533 gross 2251 @param v: a velovity
534     @return: norm of v
535     @rtype: non-negative C{float}
536 gross 2100 """
537 gross 2251 return util.sqrt(util.integrate(util.length(util.grad(v))))
538 gross 2100
539 gross 2251 def getV(self, p, v0):
540 gross 1414 """
541 gross 2251 return the value for v for a given p (overwrite)
542    
543     @param p: a pressure
544 gross 2264 @param v0: a initial guess for the value v to return.
545 gross 2251 @return: v given as M{v= A^{-1} (f-B^*p)}
546 gross 1414 """
547 gross 2251 self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress, r=v0)
548 gross 2100 if self.__stress.isEmpty():
549 gross 2251 self.__pde_u.setValue(X=p*util.kronecker(self.domain))
550 gross 2100 else:
551 gross 2251 self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain))
552 gross 2474 out=self.__pde_u.getSolution()
553 gross 2208 return out
554 gross 1414
555 gross 2445 def norm_Bv(self,Bv):
556 gross 2251 """
557     Returns Bv (overwrite).
558    
559     @rtype: equal to the type of p
560     @note: boundary conditions on p should be zero!
561     """
562 gross 2445 return util.sqrt(util.integrate(util.interpolate(Bv,Function(self.domain))**2))
563 gross 2251
564     def solve_AinvBt(self,p):
565     """
566     Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}
567    
568     @param p: a pressure increment
569 gross 2264 @return: the solution of M{Av=B^*p}
570 gross 2251 @note: boundary conditions on v should be zero!
571     """
572     self.__pde_u.setValue(Y=Data(), y=Data(), r=Data(),X=-p*util.kronecker(self.domain))
573 gross 2474 out=self.__pde_u.getSolution()
574 gross 2251 return out
575    
576 gross 2445 def solve_prec(self,Bv):
577 gross 2251 """
578     applies preconditioner for for M{BA^{-1}B^*} to M{Bv}
579 gross 2445 with accuracy L{self.getSubProblemTolerance()}
580 gross 2251
581     @param v: velocity increment
582     @return: M{p=P(Bv)} where M{P^{-1}} is an approximation of M{BA^{-1}B^*}
583     @note: boundary conditions on p are zero.
584     """
585 gross 2445 self.__pde_prec.setValue(Y=Bv)
586 gross 2474 return self.__pde_prec.getSolution()

  ViewVC Help
Powered by ViewVC 1.1.26