Diff of /trunk/escript/py_src/flows.py

revision 2208 by gross, Mon Jan 12 06:37:07 2009 UTC revision 2251 by gross, Fri Feb 6 06:50:39 2009 UTC
# Line 324  class DarcyFlow(object): Line 324  class DarcyFlow(object):
324            return self.__pde_p.getSolution(verbose=self.show_details)            return self.__pde_p.getSolution(verbose=self.show_details)
325
327        """       """
328        solves       solves
329
330            -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}            -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
331                  u_{i,i}=0                  u_{i,i}=0
# Line 333  class StokesProblemCartesian(Homogeneous Line 333  class StokesProblemCartesian(Homogeneous
334            eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j            eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j
335
336        if surface_stress is not given 0 is assumed.       if surface_stress is not given 0 is assumed.
337
338        typical usage:       typical usage:
339
340              sp=StokesProblemCartesian(domain)              sp=StokesProblemCartesian(domain)
341              sp.setTolerance()              sp.setTolerance()
342              sp.initialize(...)              sp.initialize(...)
343              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
344        """       """
345        def __init__(self,domain,**kwargs):       def __init__(self,domain,**kwargs):
346           """           """
347           initialize the Stokes Problem           initialize the Stokes Problem
348
# Line 363  class StokesProblemCartesian(Homogeneous Line 363  class StokesProblemCartesian(Homogeneous
363           # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)           # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)
364           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
365
self.__pde_proj.setReducedOrderOn()
self.__pde_proj.setSymmetryOn()
self.__pde_proj.setValue(D=1.)

367          """          """
368          assigns values to the model parameters          assigns values to the model parameters
369
# Line 393  class StokesProblemCartesian(Homogeneous Line 388  class StokesProblemCartesian(Homogeneous
388              A[i,j,j,i] += 1.              A[i,j,j,i] += 1.
389              A[i,j,i,j] += 1.              A[i,j,i,j] += 1.
390      self.__pde_prec.setValue(D=1/self.eta)      self.__pde_prec.setValue(D=1/self.eta)
392            self.__f=f
393            self.__surface_stress=surface_stress
394          self.__stress=stress          self.__stress=stress
395
396        def B(self,v):  #===============================================================================================================
397          """       def inner_pBv(self,p,v):
398          returns div(v)           """
399          @rtype: equal to the type of p           returns inner product of element p and div(v)
400
401          @note: boundary conditions on p should be zero!           @param p: a pressure increment
402          """           @param v: a residual
403          if self.show_details: print "apply divergence:"           @return: inner product of element p and div(v)
self.__pde_proj.setValue(Y=-util.div(v))
self.__pde_proj.setTolerance(self.getSubProblemTolerance())
return self.__pde_proj.getSolution(verbose=self.show_details)

def inner_pBv(self,p,Bv):
"""
returns inner product of element p and Bv  (overwrite)

@type p: equal to the type of p
@type Bv: equal to the type of result of operator B
404           @rtype: C{float}           @rtype: C{float}

@rtype: equal to the type of p
405           """           """
406           s0=util.interpolate(p,Function(self.domain))           return util.integrate(-p*util.div(v))
s1=util.interpolate(Bv,Function(self.domain))
return util.integrate(s0*s1)
407
408        def inner_p(self,p0,p1):       def inner_p(self,p0,p1):
409           """           """
410           returns inner product of element p0 and p1  (overwrite)           Returns inner product of p0 and p1

@type p0: equal to the type of p
@type p1: equal to the type of p
@rtype: C{float}
411
412           @rtype: equal to the type of p           @param p0: a pressure
413             @param p1: a pressure
414             @return: inner product of p0 and p1
415             @rtype: C{float}
416           """           """
417           s0=util.interpolate(p0/self.eta,Function(self.domain))           s0=util.interpolate(p0/self.eta,Function(self.domain))
418           s1=util.interpolate(p1/self.eta,Function(self.domain))           s1=util.interpolate(p1/self.eta,Function(self.domain))
419           return util.integrate(s0*s1)           return util.integrate(s0*s1)
420
421        def inner_v(self,v0,v1):       def norm_v(self,v):
422           """           """
423           returns inner product of two element v0 and v1  (overwrite)           returns the norm of v

@type v0: equal to the type of v
@type v1: equal to the type of v
@rtype: C{float}
424
425           @rtype: equal to the type of v           @param v: a velovity
426             @return: norm of v
427             @rtype: non-negative C{float}
428           """           """
return util.integrate(util.inner(gv0,gv1))
430
431        def solve_A(self,u,p):       def getV(self, p, v0):
432           """           """
433           solves Av=f-Au-B^*p (v=0 on fixed_u_mask)           return the value for v for a given p (overwrite)
434
435             @param p: a pressure
436             @param v0: a initial guess for the value v to return.
437             @return: v given as M{v= A^{-1} (f-B^*p)}
438           """           """
if self.show_details: print "solve for velocity:"
439           self.__pde_u.setTolerance(self.getSubProblemTolerance())           self.__pde_u.setTolerance(self.getSubProblemTolerance())
440             self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress, r=v0)
441           if self.__stress.isEmpty():           if self.__stress.isEmpty():
443           else:           else:
445           out=self.__pde_u.getSolution(verbose=self.show_details)           out=self.__pde_u.getSolution(verbose=self.show_details)
446           return  out           return  out
447
448        def solve_prec(self,p):
449           if self.show_details: print "apply preconditioner:"           raise NotImplementedError,"no v calculation implemented."
450
451
452         def norm_Bv(self,v):
453            """
454            Returns Bv (overwrite).
455
456            @rtype: equal to the type of p
457            @note: boundary conditions on p should be zero!
458            """
459            return util.sqrt(util.integrate(util.div(v)**2))
460
461         def solve_AinvBt(self,p):
462             """
463             Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}
464
465             @param p: a pressure increment
466             @return: the solution of M{Av=B^*p}
467             @note: boundary conditions on v should be zero!
468             """
469             self.__pde_u.setTolerance(self.getSubProblemTolerance())
470             self.__pde_u.setValue(Y=Data(), y=Data(), r=Data(),X=-p*util.kronecker(self.domain))
471             out=self.__pde_u.getSolution(verbose=self.show_details)
472             return  out
473
474         def solve_precB(self,v):
475             """
476             applies preconditioner for for M{BA^{-1}B^*} to M{Bv}
477             with accuracy L{self.getSubProblemTolerance()} (overwrite).
478
479             @param v: velocity increment
480             @return: M{p=P(Bv)} where M{P^{-1}} is an approximation of M{BA^{-1}B^*}
481             @note: boundary conditions on p are zero.
482             """
483             self.__pde_prec.setValue(Y=-util.div(v))
484           self.__pde_prec.setTolerance(self.getSubProblemTolerance())           self.__pde_prec.setTolerance(self.getSubProblemTolerance())
485           self.__pde_prec.setValue(Y=p)           return self.__pde_prec.getSolution(verbose=self.show_details)
q=self.__pde_prec.getSolution(verbose=self.show_details)
return q

Legend:
 Removed from v.2208 changed lines Added in v.2251