/[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 3502 - (hide annotations)
Thu Apr 28 05:06:24 2011 UTC (8 years, 3 months ago) by gross
Original Path: trunk/escript/py_src/flows.py
File MIME type: text/x-python
File size: 26407 byte(s)
Darcy flow solver is now supporting 4 solvers
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     y=-util.inner(self.domain.getNormal(),u0 * self.location_of_fixed_flux),
269     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     dp=self.__pde_p.getSolution()
347     p=GMRES(dp,
348     self.__STAB_Aprod,
349     p0,
350     self.__inner,
351     atol=self.__norm(p0+dp)*self.getTolerance() ,
352     rtol=0.,
353     iter_max=max_iter,
354     iter_restart=iter_restart,
355     verbose=self.verbose,P_R=None)
356    
357     u=self.getFlux(p, u0)
358     return u,p
359 gross 2264
360 gross 3502 def __solve_SYMSTAB(self, u0, p0, max_iter, iter_restart):
361     # p0 is used as an initial guess
362     u=self.getFlux(p0, u0)
363     self.__pde_p.setValue( Y= self.__f,
364     X= 0.5*(self.__g + u - util.tensor_mult(self.__permeability,util.grad(p0)) ),
365     y= - util.inner(self.domain.getNormal(), u),
366     r=escript.Data())
367     dp=self.__pde_p.getSolution()
368     p=GMRES(dp,
369     self.__SYMSTAB_Aprod,
370     p0,
371     self.__inner,
372     atol=self.__norm(p0+dp)*self.getTolerance() ,
373     rtol=0.,
374     iter_max=max_iter,
375     iter_restart=iter_restart,
376     verbose=self.verbose,P_R=None)
377    
378     u=self.getFlux(p, u0)
379     return u,p
380    
381     def __L2(self,v):
382     return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))
383    
384     def __norm(self,r):
385     return util.sqrt(self.__inner(r,r))
386    
387     def __inner(self,r,s):
388     return util.integrate(util.inner(r,s), escript.Function(self.domain))
389    
390     def __STAB_Aprod(self,p):
391     gp=util.grad(p)
392     self.__pde_v.setValue(Y=-0.5*gp,
393     X=-p*util.kronecker(self.__pde_v.getDomain()),
394     y= p * self.domain.getNormal(),
395     r=escript.Data())
396     u = -self.__pde_v.getSolution()
397     self.__pde_p.setValue(Y=util.div(u),
398     X=0.5*(u+util.tensor_mult(self.__permeability,gp)),
399     y=escript.Data(),
400     r=escript.Data())
401    
402     return self.__pde_p.getSolution()
403    
404     def __SYMSTAB_Aprod(self,p):
405     gp=util.grad(p)
406     self.__pde_v.setValue(Y=0.5*gp ,
407     X=escript.Data(),
408     y=escript.Data(),
409     r=escript.Data())
410     u = -self.__pde_v.getSolution()
411     self.__pde_p.setValue(Y=escript.Data(),
412     X=0.5*(-u+util.tensor_mult(self.__permeability,gp)),
413     y= util.inner(self.domain.getNormal(), u),
414     r=escript.Data())
415    
416     return self.__pde_p.getSolution()
417    
418    
419 gross 1414 class StokesProblemCartesian(HomogeneousSaddlePointProblem):
420 gross 2251 """
421 gross 2264 solves
422 gross 1414
423 gross 2208 -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
424     u_{i,i}=0
425 gross 1414
426 gross 2208 u=0 where fixed_u_mask>0
427     eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j
428 gross 1414
429 gross 2264 if surface_stress is not given 0 is assumed.
430 gross 1414
431 gross 2251 typical usage:
432 gross 1414
433 gross 2208 sp=StokesProblemCartesian(domain)
434     sp.setTolerance()
435     sp.initialize(...)
436     v,p=sp.solve(v0,p0)
437 gross 2251 """
438 gross 2719 def __init__(self,domain,**kwargs):
439 gross 2100 """
440 gross 2208 initialize the Stokes Problem
441 gross 2100
442 gross 2793 The approximation spaces used for velocity (=Solution(domain)) and pressure (=ReducedSolution(domain)) must be
443     LBB complient, for instance using quadratic and linear approximation on the same element or using linear approximation
444     with macro elements for the pressure.
445    
446     :param domain: domain of the problem.
447 jfenwick 2625 :type domain: `Domain`
448 gross 2100 """
449 gross 2719 HomogeneousSaddlePointProblem.__init__(self,**kwargs)
450 gross 1414 self.domain=domain
451 gross 3502 self.__pde_v=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
452     self.__pde_v.setSymmetryOn()
453 gross 2474
454 gross 1414 self.__pde_prec=LinearPDE(domain)
455     self.__pde_prec.setReducedOrderOn()
456     self.__pde_prec.setSymmetryOn()
457    
458 gross 2415 self.__pde_proj=LinearPDE(domain)
459     self.__pde_proj.setReducedOrderOn()
460     self.__pde_proj.setValue(D=1)
461     self.__pde_proj.setSymmetryOn()
462    
463 gross 2474 def getSolverOptionsVelocity(self):
464     """
465     returns the solver options used solve the equation for velocity.
466    
467 jfenwick 2625 :rtype: `SolverOptions`
468 gross 2474 """
469 gross 3502 return self.__pde_v.getSolverOptions()
470 gross 2474 def setSolverOptionsVelocity(self, options=None):
471     """
472     set the solver options for solving the equation for velocity.
473    
474 jfenwick 2625 :param options: new solver options
475     :type options: `SolverOptions`
476 gross 2474 """
477 gross 3502 self.__pde_v.setSolverOptions(options)
478 gross 2474 def getSolverOptionsPressure(self):
479     """
480     returns the solver options used solve the equation for pressure.
481 jfenwick 2625 :rtype: `SolverOptions`
482 gross 2474 """
483     return self.__pde_prec.getSolverOptions()
484     def setSolverOptionsPressure(self, options=None):
485     """
486     set the solver options for solving the equation for pressure.
487 jfenwick 2625 :param options: new solver options
488     :type options: `SolverOptions`
489 gross 2474 """
490     self.__pde_prec.setSolverOptions(options)
491 gross 2415
492 gross 2474 def setSolverOptionsDiv(self, options=None):
493     """
494     set the solver options for solving the equation to project the divergence of
495     the velocity onto the function space of presure.
496    
497 jfenwick 2625 :param options: new solver options
498     :type options: `SolverOptions`
499 gross 2474 """
500 artak 2689 self.__pde_proj.setSolverOptions(options)
501 gross 2474 def getSolverOptionsDiv(self):
502     """
503     returns the solver options for solving the equation to project the divergence of
504     the velocity onto the function space of presure.
505    
506 jfenwick 2625 :rtype: `SolverOptions`
507 gross 2474 """
508 artak 2689 return self.__pde_proj.getSolverOptions()
509 gross 2474
510 gross 2793 def updateStokesEquation(self, v, p):
511     """
512     updates the Stokes equation to consider dependencies from ``v`` and ``p``
513     :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values.
514     """
515     pass
516     def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None):
517     """
518     assigns new values to the model parameters.
519    
520     :param f: external force
521     :type f: `Vector` object in `FunctionSpace` `Function` or similar
522     :param fixed_u_mask: mask of locations with fixed velocity.
523     :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
524     :param eta: viscosity
525     :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
526     :param surface_stress: normal surface stress
527     :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
528     :param stress: initial stress
529     :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
530     """
531     if eta !=None:
532     k=util.kronecker(self.domain.getDim())
533     kk=util.outer(k,k)
534 jfenwick 3432 self.eta=util.interpolate(eta, escript.Function(self.domain))
535 gross 2793 self.__pde_prec.setValue(D=1/self.eta)
536 gross 3502 self.__pde_v.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))
537 gross 2793 if restoration_factor!=None:
538     n=self.domain.getNormal()
539 gross 3502 self.__pde_v.setValue(d=restoration_factor*util.outer(n,n))
540 gross 2793 if fixed_u_mask!=None:
541 gross 3502 self.__pde_v.setValue(q=fixed_u_mask)
542 gross 2793 if f!=None: self.__f=f
543     if surface_stress!=None: self.__surface_stress=surface_stress
544     if stress!=None: self.__stress=stress
545    
546 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):
547 gross 2208 """
548     assigns values to the model parameters
549 gross 2100
550 jfenwick 2625 :param f: external force
551     :type f: `Vector` object in `FunctionSpace` `Function` or similar
552     :param fixed_u_mask: mask of locations with fixed velocity.
553     :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
554     :param eta: viscosity
555     :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
556     :param surface_stress: normal surface stress
557     :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
558     :param stress: initial stress
559     :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
560 gross 2208 """
561 gross 2793 self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
562 gross 1414
563 gross 2719 def Bv(self,v,tol):
564 gross 2251 """
565     returns inner product of element p and div(v)
566 gross 1414
567 jfenwick 2625 :param v: a residual
568     :return: inner product of element p and div(v)
569     :rtype: ``float``
570 gross 2100 """
571 gross 2793 self.__pde_proj.setValue(Y=-util.div(v))
572 gross 2719 self.getSolverOptionsDiv().setTolerance(tol)
573     self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
574     out=self.__pde_proj.getSolution()
575     return out
576 gross 2208
577 gross 2445 def inner_pBv(self,p,Bv):
578     """
579     returns inner product of element p and Bv=-div(v)
580    
581 jfenwick 2625 :param p: a pressure increment
582     :param Bv: a residual
583     :return: inner product of element p and Bv=-div(v)
584     :rtype: ``float``
585 gross 2445 """
586 jfenwick 3432 return util.integrate(util.interpolate(p,escript.Function(self.domain))*util.interpolate(Bv, escript.Function(self.domain)))
587 gross 2445
588 gross 2251 def inner_p(self,p0,p1):
589 gross 2100 """
590 gross 2251 Returns inner product of p0 and p1
591 gross 1414
592 jfenwick 2625 :param p0: a pressure
593     :param p1: a pressure
594     :return: inner product of p0 and p1
595     :rtype: ``float``
596 gross 2100 """
597 jfenwick 3432 s0=util.interpolate(p0, escript.Function(self.domain))
598     s1=util.interpolate(p1, escript.Function(self.domain))
599 gross 2100 return util.integrate(s0*s1)
600 artak 1550
601 gross 2251 def norm_v(self,v):
602 gross 2100 """
603 gross 2251 returns the norm of v
604 gross 2208
605 jfenwick 2625 :param v: a velovity
606     :return: norm of v
607     :rtype: non-negative ``float``
608 gross 2100 """
609 gross 2719 return util.sqrt(util.integrate(util.length(util.grad(v))**2))
610 gross 2100
611 gross 2793
612 gross 2719 def getDV(self, p, v, tol):
613 gross 1414 """
614 gross 2251 return the value for v for a given p (overwrite)
615    
616 jfenwick 2625 :param p: a pressure
617 gross 2719 :param v: a initial guess for the value v to return.
618     :return: dv given as *Adv=(f-Av-B^*p)*
619 gross 1414 """
620 gross 2793 self.updateStokesEquation(v,p)
621 gross 3502 self.__pde_v.setValue(Y=self.__f, y=self.__surface_stress)
622 gross 2719 self.getSolverOptionsVelocity().setTolerance(tol)
623     self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
624 gross 2100 if self.__stress.isEmpty():
625 gross 3502 self.__pde_v.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
626 gross 2100 else:
627 gross 3502 self.__pde_v.setValue(X=self.__stress+p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
628     out=self.__pde_v.getSolution()
629 gross 2208 return out
630 gross 1414
631 gross 2445 def norm_Bv(self,Bv):
632 gross 2251 """
633     Returns Bv (overwrite).
634    
635 jfenwick 2625 :rtype: equal to the type of p
636     :note: boundary conditions on p should be zero!
637 gross 2251 """
638 jfenwick 3432 return util.sqrt(util.integrate(util.interpolate(Bv, escript.Function(self.domain))**2))
639 gross 2251
640 gross 2719 def solve_AinvBt(self,p, tol):
641 gross 2251 """
642 gross 2719 Solves *Av=B^*p* with accuracy `tol`
643 gross 2251
644 jfenwick 2625 :param p: a pressure increment
645     :return: the solution of *Av=B^*p*
646     :note: boundary conditions on v should be zero!
647 gross 2251 """
648 gross 3502 self.__pde_v.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))
649     out=self.__pde_v.getSolution()
650 gross 2251 return out
651    
652 gross 2719 def solve_prec(self,Bv, tol):
653 gross 2251 """
654 jfenwick 2625 applies preconditioner for for *BA^{-1}B^** to *Bv*
655     with accuracy `self.getSubProblemTolerance()`
656 gross 2251
657 jfenwick 2625 :param Bv: velocity increment
658     :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*
659     :note: boundary conditions on p are zero.
660 gross 2251 """
661 gross 2445 self.__pde_prec.setValue(Y=Bv)
662 gross 2719 self.getSolverOptionsPressure().setTolerance(tol)
663     self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
664     out=self.__pde_prec.getSolution()
665     return out

  ViewVC Help
Powered by ViewVC 1.1.26