/[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 3510 - (hide annotations)
Fri May 13 06:09:49 2011 UTC (8 years, 4 months ago) by gross
File MIME type: text/x-python
File size: 26432 byte(s)
some fixes for the compilation without BOOMERAMG
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 3502 class DarcyFlow(object):
41 gross 3499 """
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 gross 3502 :cvar SIMPLE: simple solver
50     :cvar POST: solver using global postprocessing of flux
51     :cvar STAB: solver uses (non-symmetric) stabilization
52     :cvar SYMSTAB: solver uses symmetric stabilization
53 gross 3499 """
54 gross 3502 SIMPLE="SIMPLE"
55     POST="POST"
56     STAB="STAB"
57     SYMSTAB="SYMSTAB"
58     def __init__(self, domain, useReduced=False, solver="SYMSTAB", verbose=False, w=1.):
59 gross 3499 """
60     initializes the Darcy flux problem
61     :param domain: domain of the problem
62     :type domain: `Domain`
63     :param useReduced: uses reduced oreder on flux and pressure
64     :type useReduced: ``bool``
65 gross 3502 :param solver: solver method
66     :type solver: in [`DarcyFlow.SIMPLE`, `DarcyFlow.POST', `DarcyFlow.STAB`, `DarcyFlow.SYMSTAB` ]
67     :param verbose: if ``True`` some information on the iteration progress are printed.
68     :type verbose: ``bool``
69     :param w: weighting factor for `DarcyFlow.POST` solver
70     :type w: ``float``
71    
72 gross 3499 """
73     self.domain=domain
74 gross 3502 self.solver=solver
75 gross 3499 self.useReduced=useReduced
76 gross 3502 self.verbose=verbose
77 gross 3499 self.scale=1.
78 gross 3500
79 gross 3499
80 gross 3502 self.__pde_v=LinearPDESystem(domain)
81     self.__pde_v.setSymmetryOn()
82     if self.useReduced: self.__pde_v.setReducedOrderOn()
83     self.__pde_p=LinearSinglePDE(domain)
84     self.__pde_p.setSymmetryOn()
85     if self.useReduced: self.__pde_p.setReducedOrderOn()
86 gross 3499
87 gross 3502 if self.solver == self.SIMPLE:
88     if self.verbose: print "DarcyFlow: simple solver is used."
89     self.__pde_v.setValue(D=util.kronecker(self.domain.getDim()))
90     elif self.solver == self.POST:
91     self.w=w
92     if util.inf(w)<0.:
93     raise ValueError,"Weighting factor must be non-negative."
94     if self.verbose: print "DarcyFlow: global postprocessing of flux is used."
95     elif self.solver == self.STAB:
96     if self.verbose: print "DarcyFlow: (non-symmetric) stabilization is used."
97     elif self.solver == self.SYMSTAB:
98     if self.verbose: print "DarcyFlow: symmetric stabilization is used."
99     else:
100     raise ValueError,"unknown solver %s"%self.solver
101     self.__f=escript.Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))
102     self.__g=escript.Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
103     self.location_of_fixed_pressure = escript.Scalar(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
104     self.location_of_fixed_flux = escript.Vector(0, self.__pde_v.getFunctionSpaceForCoefficient("q"))
105     self.setTolerance()
106    
107    
108 gross 3499 def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
109     """
110     assigns values to model parameters
111    
112     :param f: volumetic sources/sinks
113     :type f: scalar value on the domain (e.g. `escript.Data`)
114     :param g: flux sources/sinks
115     :type g: vector values on the domain (e.g. `escript.Data`)
116     :param location_of_fixed_pressure: mask for locations where pressure is fixed
117     :type location_of_fixed_pressure: scalar value on the domain (e.g. `escript.Data`)
118     :param location_of_fixed_flux: mask for locations where flux is fixed.
119     :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)
120     :param permeability: permeability tensor. If scalar ``s`` is given the tensor with ``s`` on the main diagonal is used.
121 gross 3502 :type permeability: scalar or symmetric tensor values on the domain (e.g. `escript.Data`)
122 gross 3499
123     :note: the values of parameters which are not set by calling ``setValue`` are not altered.
124     :note: at any point on the boundary of the domain the pressure
125     (``location_of_fixed_pressure`` >0) or the normal component of the
126     flux (``location_of_fixed_flux[i]>0``) if direction of the normal
127     is along the *x_i* axis.
128    
129     """
130 gross 3502 if location_of_fixed_pressure!=None:
131     self.location_of_fixed_pressure=util.wherePositive(location_of_fixed_pressure)
132     self.__pde_p.setValue(q=self.location_of_fixed_pressure)
133     if location_of_fixed_flux!=None:
134     self.location_of_fixed_flux=util.wherePositive(location_of_fixed_flux)
135     self.__pde_v.setValue(q=self.location_of_fixed_flux)
136 gross 3499
137    
138 gross 3502 # pressure is rescaled by the factor 1/self.scale
139 gross 3499 if permeability!=None:
140 gross 3502
141     perm=util.interpolate(permeability,self.__pde_v.getFunctionSpaceForCoefficient("A"))
142 gross 3499 V=util.vol(self.domain)
143 gross 3502 l=V**(1./self.domain.getDim())
144    
145 gross 3499 if perm.getRank()==0:
146     perm_inv=(1./perm)
147 gross 3502 self.scale=util.integrate(perm_inv)/V*l
148 gross 3499 perm_inv=perm_inv*((1./self.scale)*util.kronecker(self.domain.getDim()))
149     perm=perm*(self.scale*util.kronecker(self.domain.getDim()))
150 gross 3502
151    
152 gross 3499 elif perm.getRank()==2:
153     perm_inv=util.inverse(perm)
154 gross 3502 self.scale=util.sqrt(util.integrate(util.length(perm_inv)**2)/V)*l
155 gross 3499 perm_inv*=(1./self.scale)
156     perm=perm*self.scale
157     else:
158     raise ValueError,"illegal rank of permeability."
159    
160     self.__permeability=perm
161     self.__permeability_inv=perm_inv
162 gross 3502 if self.verbose: print "DarcyFlow: scaling factor for pressure is %e."%self.scale
163    
164     if self.solver == self.SIMPLE:
165     self.__pde_p.setValue(A=self.__permeability)
166     elif self.solver == self.POST:
167     self.__pde_p.setValue(A=self.__permeability)
168     k=util.kronecker(self.domain.getDim())
169     self.lamb = self.w*util.length(perm_inv)*l
170     self.__pde_v.setValue(D=self.__permeability_inv, A=self.lamb*self.domain.getSize()*util.outer(k,k))
171     elif self.solver == self.STAB:
172 gross 3499 self.__pde_p.setValue(A=0.5*self.__permeability)
173 gross 3502 self.__pde_v.setValue(D=0.5*self.__permeability_inv)
174     elif self.solver == self.SYMSTAB:
175     self.__pde_p.setValue(A=0.5*self.__permeability)
176     self.__pde_v.setValue(D=0.5*self.__permeability_inv)
177    
178 gross 3499 if g != None:
179 gross 3502 g=util.interpolate(g, self.__pde_v.getFunctionSpaceForCoefficient("Y"))
180 gross 3499 if g.isEmpty():
181 gross 3502 g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
182 gross 3499 else:
183     if not g.getShape()==(self.domain.getDim(),): raise ValueError,"illegal shape of g"
184 gross 3502 self.__g=g
185 gross 3499 if f !=None:
186 gross 3502 f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
187     if f.isEmpty():
188     f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
189 gross 3499 else:
190     if f.getRank()>0: raise ValueError,"illegal rank of f."
191 gross 3502 self.__f=f
192 gross 3071 def getSolverOptionsFlux(self):
193     """
194     Returns the solver options used to solve the flux problems
195     :return: `SolverOptions`
196     """
197 gross 3502 return self.__pde_v.getSolverOptions()
198 gross 3071
199     def setSolverOptionsFlux(self, options=None):
200     """
201     Sets the solver options used to solve the flux problems
202     If ``options`` is not present, the options are reset to default
203     :param options: `SolverOptions`
204     """
205     return self.__pde_v.setSolverOptions(options)
206    
207     def getSolverOptionsPressure(self):
208     """
209     Returns the solver options used to solve the pressure problems
210     :return: `SolverOptions`
211     """
212     return self.__pde_p.getSolverOptions()
213    
214     def setSolverOptionsPressure(self, options=None):
215     """
216     Sets the solver options used to solve the pressure problems
217     If ``options`` is not present, the options are reset to default
218    
219     :param options: `SolverOptions`
220     :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
221     """
222     return self.__pde_p.setSolverOptions(options)
223    
224     def setTolerance(self,rtol=1e-4):
225     """
226 gross 3502 sets the relative tolerance ``rtol`` for the pressure for the stabelized solvers.
227 gross 3071
228     :param rtol: relative tolerance for the pressure
229     :type rtol: non-negative ``float``
230     """
231     if rtol<0:
232     raise ValueError,"Relative tolerance needs to be non-negative."
233     self.__rtol=rtol
234 gross 3502
235 gross 3071 def getTolerance(self):
236     """
237     returns the relative tolerance
238     :return: current relative tolerance
239     :rtype: ``float``
240     """
241     return self.__rtol
242 gross 3502
243     def solve(self,u0,p0, max_iter=100, iter_restart=20):
244 gross 3071 """
245 gross 3502 solves the problem.
246 gross 3071
247 gross 3502 The iteration is terminated if the residual norm is less then self.getTolerance().
248 gross 2960
249 gross 3502 :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.
250     :type u0: vector value on the domain (e.g. `escript.Data`).
251     :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.
252     :type p0: scalar value on the domain (e.g. `escript.Data`).
253     :param max_iter: maximum number of (outer) iteration steps for the stabilization solvers,
254     :type max_iter: ``int``
255     :param iter_restart: number of steps after which the iteration is restarted. The larger ``iter_restart`` the larger the required memory.
256     A small value for ``iter_restart`` may require a large number of iteration steps or may even lead to a failure
257     of the iteration. ``iter_restart`` is relevant for the stabilization solvers only.
258     :type iter_restart: ``int``
259     :return: flux and pressure
260     :rtype: ``tuple`` of `escript.Data`.
261 gross 2960
262 gross 3071 """
263 gross 3502 # rescale initial guess:
264     p0=p0/self.scale
265     if self.solver == self.SIMPLE or self.solver == self.POST :
266     self.__pde_p.setValue(X=self.__g ,
267     Y=self.__f,
268 gross 3510 y= - util.inner(self.domain.getNormal(),u0 * self.location_of_fixed_flux),
269 gross 3502 r=p0)
270     p=self.__pde_p.getSolution()
271     u = self.getFlux(p, u0)
272     elif self.solver == self.STAB:
273     u,p = self.__solve_STAB(u0,p0, max_iter, iter_restart)
274     elif self.solver == self.SYMSTAB:
275     u,p = self.__solve_SYMSTAB(u0,p0, max_iter, iter_restart)
276 gross 2474
277 gross 3502 if self.verbose:
278     KGp=util.tensor_mult(self.__permeability,util.grad(p))
279     def_p=self.__g-(u+KGp)
280     def_v=self.__f-util.div(u, self.__pde_v.getFunctionSpaceForCoefficient("X"))
281     print "DarcyFlux: |g-u-K*grad(p)|_2 = %e (|u|_2 = %e)."%(self.__L2(def_p),self.__L2(u))
282     print "DarcyFlux: |f-div(u)|_2 = %e (|grad(u)|_2 = %e)."%(self.__L2(def_v),self.__L2(util.grad(u)))
283     #rescale result
284     p=p*self.scale
285     return u,p
286    
287     def getFlux(self,p, u0=None):
288 gross 2100 """
289 gross 3502 returns the flux for a given pressure ``p`` where the flux is equal to ``u0``
290 jfenwick 2625 on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
291 gross 3502 Notice that ``g`` and ``f`` are used, see `setValue`.
292 gross 2264
293 jfenwick 2625 :param p: pressure.
294 gross 3440 :type p: scalar value on the domain (e.g. `escript.Data`).
295 gross 3502 :param u0: flux on the locations of the domain marked be ``location_of_fixed_flux``.
296     :type u0: vector values on the domain (e.g. `escript.Data`) or ``None``
297 jfenwick 2625 :return: flux
298 gross 3440 :rtype: `escript.Data`
299 gross 2264 """
300 gross 3502 if self.solver == self.SIMPLE or self.solver == self.POST :
301     KGp=util.tensor_mult(self.__permeability,util.grad(p))
302     self.__pde_v.setValue(Y=self.__g-KGp, X=escript.Data())
303     if u0 == None:
304     self.__pde_v.setValue(r=escript.Data())
305     else:
306     self.__pde_v.setValue(r=u0)
307     u= self.__pde_v.getSolution()
308     elif self.solver == self.POST:
309     self.__pde_v.setValue(Y=util.tensor_mult(self.__permeability_inv,self.__g)-util.grad(p),
310     X=self.lamb * self.__f * util.kronecker(self.domain.getDim()))
311     if u0 == None:
312     self.__pde_v.setValue(r=escript.Data())
313     else:
314     self.__pde_v.setValue(r=u0)
315     u= self.__pde_v.getSolution()
316     elif self.solver == self.STAB:
317     gp=util.grad(p)
318     self.__pde_v.setValue(Y=0.5*(util.tensor_mult(self.__permeability_inv,self.__g)+gp),
319     X= p * util.kronecker(self.domain.getDim()),
320     y= - p * self.domain.getNormal())
321     if u0 == None:
322     self.__pde_v.setValue(r=escript.Data())
323     else:
324     self.__pde_v.setValue(r=u0)
325     u= self.__pde_v.getSolution()
326     elif self.solver == self.SYMSTAB:
327     gp=util.grad(p)
328     self.__pde_v.setValue(Y=0.5*(util.tensor_mult(self.__permeability_inv,self.__g)-gp),
329     X= escript.Data() ,
330     y= escript.Data() )
331     if u0 == None:
332     self.__pde_v.setValue(r=escript.Data())
333     else:
334     self.__pde_v.setValue(r=u0)
335     u= self.__pde_v.getSolution()
336     return u
337    
338    
339     def __solve_STAB(self, u0, p0, max_iter, iter_restart):
340     # p0 is used as an initial guess
341     u=self.getFlux(p0, u0)
342     self.__pde_p.setValue( Y=self.__f-util.div(u),
343     X=0.5*(self.__g - u - util.tensor_mult(self.__permeability,util.grad(p0)) ),
344     y= escript.Data(),
345     r=escript.Data())
346 gross 3510
347 gross 3502 dp=self.__pde_p.getSolution()
348     p=GMRES(dp,
349     self.__STAB_Aprod,
350     p0,
351     self.__inner,
352     atol=self.__norm(p0+dp)*self.getTolerance() ,
353     rtol=0.,
354     iter_max=max_iter,
355     iter_restart=iter_restart,
356     verbose=self.verbose,P_R=None)
357    
358     u=self.getFlux(p, u0)
359     return u,p
360 gross 2264
361 gross 3502 def __solve_SYMSTAB(self, u0, p0, max_iter, iter_restart):
362     # p0 is used as an initial guess
363     u=self.getFlux(p0, u0)
364     self.__pde_p.setValue( Y= self.__f,
365     X= 0.5*(self.__g + u - util.tensor_mult(self.__permeability,util.grad(p0)) ),
366 gross 3510 y= - util.inner(self.domain.getNormal(), u),
367 gross 3502 r=escript.Data())
368     dp=self.__pde_p.getSolution()
369 gross 3510
370     print dp
371     print p0+dp
372    
373 gross 3502 p=GMRES(dp,
374     self.__SYMSTAB_Aprod,
375     p0,
376     self.__inner,
377     atol=self.__norm(p0+dp)*self.getTolerance() ,
378     rtol=0.,
379     iter_max=max_iter,
380     iter_restart=iter_restart,
381     verbose=self.verbose,P_R=None)
382    
383     u=self.getFlux(p, u0)
384     return u,p
385    
386     def __L2(self,v):
387     return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))
388    
389     def __norm(self,r):
390     return util.sqrt(self.__inner(r,r))
391    
392     def __inner(self,r,s):
393     return util.integrate(util.inner(r,s), escript.Function(self.domain))
394    
395     def __STAB_Aprod(self,p):
396     gp=util.grad(p)
397     self.__pde_v.setValue(Y=-0.5*gp,
398     X=-p*util.kronecker(self.__pde_v.getDomain()),
399     y= p * self.domain.getNormal(),
400     r=escript.Data())
401     u = -self.__pde_v.getSolution()
402     self.__pde_p.setValue(Y=util.div(u),
403     X=0.5*(u+util.tensor_mult(self.__permeability,gp)),
404     y=escript.Data(),
405     r=escript.Data())
406    
407     return self.__pde_p.getSolution()
408    
409     def __SYMSTAB_Aprod(self,p):
410     gp=util.grad(p)
411     self.__pde_v.setValue(Y=0.5*gp ,
412     X=escript.Data(),
413     y=escript.Data(),
414     r=escript.Data())
415     u = -self.__pde_v.getSolution()
416     self.__pde_p.setValue(Y=escript.Data(),
417     X=0.5*(-u+util.tensor_mult(self.__permeability,gp)),
418 gross 3510 y=escript.Data(),
419 gross 3502 r=escript.Data())
420    
421     return self.__pde_p.getSolution()
422    
423    
424 gross 1414 class StokesProblemCartesian(HomogeneousSaddlePointProblem):
425 gross 2251 """
426 gross 2264 solves
427 gross 1414
428 gross 2208 -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
429     u_{i,i}=0
430 gross 1414
431 gross 2208 u=0 where fixed_u_mask>0
432     eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j
433 gross 1414
434 gross 2264 if surface_stress is not given 0 is assumed.
435 gross 1414
436 gross 2251 typical usage:
437 gross 1414
438 gross 2208 sp=StokesProblemCartesian(domain)
439     sp.setTolerance()
440     sp.initialize(...)
441     v,p=sp.solve(v0,p0)
442 gross 2251 """
443 gross 2719 def __init__(self,domain,**kwargs):
444 gross 2100 """
445 gross 2208 initialize the Stokes Problem
446 gross 2100
447 gross 2793 The approximation spaces used for velocity (=Solution(domain)) and pressure (=ReducedSolution(domain)) must be
448     LBB complient, for instance using quadratic and linear approximation on the same element or using linear approximation
449     with macro elements for the pressure.
450    
451     :param domain: domain of the problem.
452 jfenwick 2625 :type domain: `Domain`
453 gross 2100 """
454 gross 2719 HomogeneousSaddlePointProblem.__init__(self,**kwargs)
455 gross 1414 self.domain=domain
456 gross 3502 self.__pde_v=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
457     self.__pde_v.setSymmetryOn()
458 gross 2474
459 gross 1414 self.__pde_prec=LinearPDE(domain)
460     self.__pde_prec.setReducedOrderOn()
461     self.__pde_prec.setSymmetryOn()
462    
463 gross 2415 self.__pde_proj=LinearPDE(domain)
464     self.__pde_proj.setReducedOrderOn()
465     self.__pde_proj.setValue(D=1)
466     self.__pde_proj.setSymmetryOn()
467    
468 gross 2474 def getSolverOptionsVelocity(self):
469     """
470     returns the solver options used solve the equation for velocity.
471    
472 jfenwick 2625 :rtype: `SolverOptions`
473 gross 2474 """
474 gross 3502 return self.__pde_v.getSolverOptions()
475 gross 2474 def setSolverOptionsVelocity(self, options=None):
476     """
477     set the solver options for solving the equation for velocity.
478    
479 jfenwick 2625 :param options: new solver options
480     :type options: `SolverOptions`
481 gross 2474 """
482 gross 3502 self.__pde_v.setSolverOptions(options)
483 gross 2474 def getSolverOptionsPressure(self):
484     """
485     returns the solver options used solve the equation for pressure.
486 jfenwick 2625 :rtype: `SolverOptions`
487 gross 2474 """
488     return self.__pde_prec.getSolverOptions()
489     def setSolverOptionsPressure(self, options=None):
490     """
491     set the solver options for solving the equation for pressure.
492 jfenwick 2625 :param options: new solver options
493     :type options: `SolverOptions`
494 gross 2474 """
495     self.__pde_prec.setSolverOptions(options)
496 gross 2415
497 gross 2474 def setSolverOptionsDiv(self, options=None):
498     """
499     set the solver options for solving the equation to project the divergence of
500     the velocity onto the function space of presure.
501    
502 jfenwick 2625 :param options: new solver options
503     :type options: `SolverOptions`
504 gross 2474 """
505 artak 2689 self.__pde_proj.setSolverOptions(options)
506 gross 2474 def getSolverOptionsDiv(self):
507     """
508     returns the solver options for solving the equation to project the divergence of
509     the velocity onto the function space of presure.
510    
511 jfenwick 2625 :rtype: `SolverOptions`
512 gross 2474 """
513 artak 2689 return self.__pde_proj.getSolverOptions()
514 gross 2474
515 gross 2793 def updateStokesEquation(self, v, p):
516     """
517     updates the Stokes equation to consider dependencies from ``v`` and ``p``
518     :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values.
519     """
520     pass
521     def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None):
522     """
523     assigns new values to the model parameters.
524    
525     :param f: external force
526     :type f: `Vector` object in `FunctionSpace` `Function` or similar
527     :param fixed_u_mask: mask of locations with fixed velocity.
528     :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
529     :param eta: viscosity
530     :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
531     :param surface_stress: normal surface stress
532     :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
533     :param stress: initial stress
534     :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
535     """
536     if eta !=None:
537     k=util.kronecker(self.domain.getDim())
538     kk=util.outer(k,k)
539 jfenwick 3432 self.eta=util.interpolate(eta, escript.Function(self.domain))
540 gross 2793 self.__pde_prec.setValue(D=1/self.eta)
541 gross 3502 self.__pde_v.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))
542 gross 2793 if restoration_factor!=None:
543     n=self.domain.getNormal()
544 gross 3502 self.__pde_v.setValue(d=restoration_factor*util.outer(n,n))
545 gross 2793 if fixed_u_mask!=None:
546 gross 3502 self.__pde_v.setValue(q=fixed_u_mask)
547 gross 2793 if f!=None: self.__f=f
548     if surface_stress!=None: self.__surface_stress=surface_stress
549     if stress!=None: self.__stress=stress
550    
551 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):
552 gross 2208 """
553     assigns values to the model parameters
554 gross 2100
555 jfenwick 2625 :param f: external force
556     :type f: `Vector` object in `FunctionSpace` `Function` or similar
557     :param fixed_u_mask: mask of locations with fixed velocity.
558     :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
559     :param eta: viscosity
560     :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
561     :param surface_stress: normal surface stress
562     :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
563     :param stress: initial stress
564     :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
565 gross 2208 """
566 gross 2793 self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
567 gross 1414
568 gross 2719 def Bv(self,v,tol):
569 gross 2251 """
570     returns inner product of element p and div(v)
571 gross 1414
572 jfenwick 2625 :param v: a residual
573     :return: inner product of element p and div(v)
574     :rtype: ``float``
575 gross 2100 """
576 gross 2793 self.__pde_proj.setValue(Y=-util.div(v))
577 gross 2719 self.getSolverOptionsDiv().setTolerance(tol)
578     self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
579     out=self.__pde_proj.getSolution()
580     return out
581 gross 2208
582 gross 2445 def inner_pBv(self,p,Bv):
583     """
584     returns inner product of element p and Bv=-div(v)
585    
586 jfenwick 2625 :param p: a pressure increment
587     :param Bv: a residual
588     :return: inner product of element p and Bv=-div(v)
589     :rtype: ``float``
590 gross 2445 """
591 jfenwick 3432 return util.integrate(util.interpolate(p,escript.Function(self.domain))*util.interpolate(Bv, escript.Function(self.domain)))
592 gross 2445
593 gross 2251 def inner_p(self,p0,p1):
594 gross 2100 """
595 gross 2251 Returns inner product of p0 and p1
596 gross 1414
597 jfenwick 2625 :param p0: a pressure
598     :param p1: a pressure
599     :return: inner product of p0 and p1
600     :rtype: ``float``
601 gross 2100 """
602 jfenwick 3432 s0=util.interpolate(p0, escript.Function(self.domain))
603     s1=util.interpolate(p1, escript.Function(self.domain))
604 gross 2100 return util.integrate(s0*s1)
605 artak 1550
606 gross 2251 def norm_v(self,v):
607 gross 2100 """
608 gross 2251 returns the norm of v
609 gross 2208
610 jfenwick 2625 :param v: a velovity
611     :return: norm of v
612     :rtype: non-negative ``float``
613 gross 2100 """
614 gross 2719 return util.sqrt(util.integrate(util.length(util.grad(v))**2))
615 gross 2100
616 gross 2793
617 gross 2719 def getDV(self, p, v, tol):
618 gross 1414 """
619 gross 2251 return the value for v for a given p (overwrite)
620    
621 jfenwick 2625 :param p: a pressure
622 gross 2719 :param v: a initial guess for the value v to return.
623     :return: dv given as *Adv=(f-Av-B^*p)*
624 gross 1414 """
625 gross 2793 self.updateStokesEquation(v,p)
626 gross 3502 self.__pde_v.setValue(Y=self.__f, y=self.__surface_stress)
627 gross 2719 self.getSolverOptionsVelocity().setTolerance(tol)
628     self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
629 gross 2100 if self.__stress.isEmpty():
630 gross 3502 self.__pde_v.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
631 gross 2100 else:
632 gross 3502 self.__pde_v.setValue(X=self.__stress+p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
633     out=self.__pde_v.getSolution()
634 gross 2208 return out
635 gross 1414
636 gross 2445 def norm_Bv(self,Bv):
637 gross 2251 """
638     Returns Bv (overwrite).
639    
640 jfenwick 2625 :rtype: equal to the type of p
641     :note: boundary conditions on p should be zero!
642 gross 2251 """
643 jfenwick 3432 return util.sqrt(util.integrate(util.interpolate(Bv, escript.Function(self.domain))**2))
644 gross 2251
645 gross 2719 def solve_AinvBt(self,p, tol):
646 gross 2251 """
647 gross 2719 Solves *Av=B^*p* with accuracy `tol`
648 gross 2251
649 jfenwick 2625 :param p: a pressure increment
650     :return: the solution of *Av=B^*p*
651     :note: boundary conditions on v should be zero!
652 gross 2251 """
653 gross 3502 self.__pde_v.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))
654     out=self.__pde_v.getSolution()
655 gross 2251 return out
656    
657 gross 2719 def solve_prec(self,Bv, tol):
658 gross 2251 """
659 jfenwick 2625 applies preconditioner for for *BA^{-1}B^** to *Bv*
660     with accuracy `self.getSubProblemTolerance()`
661 gross 2251
662 jfenwick 2625 :param Bv: velocity increment
663     :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*
664     :note: boundary conditions on p are zero.
665 gross 2251 """
666 gross 2445 self.__pde_prec.setValue(Y=Bv)
667 gross 2719 self.getSolverOptionsPressure().setTolerance(tol)
668     self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
669     out=self.__pde_prec.getSolution()
670     return out

  ViewVC Help
Powered by ViewVC 1.1.26