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

revision 1658 by gross, Mon Jul 14 08:55:25 2008 UTC revision 1659 by gross, Fri Jul 18 02:28:13 2008 UTC
# Line 38  import util Line 38  import util
38  from linearPDEs import LinearPDE  from linearPDEs import LinearPDE
39  from pdetools import HomogeneousSaddlePointProblem,Projector  from pdetools import HomogeneousSaddlePointProblem,Projector
40
42          """
43          solves
44
45              -(eta*(u_{i,j}+u_{j,i}))_j - p_i = f_i
46                    u_{i,i}=0
47
48              u=0 where  fixed_u_mask>0
49              eta*(u_{i,j}+u_{j,i})*n_j=surface_stress
50
51          if surface_stress is not give 0 is assumed.
52
53          typical usage:
54
55                sp=StokesProblemCartesian(domain)
56                sp.setTolerance()
57                sp.initialize(...)
58                v,p=sp.solve(v0,p0)
59          """
60          def __init__(self,domain,**kwargs):
62             self.domain=domain
63             self.vol=util.integrate(1.,Function(self.domain))
64             self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
65             self.__pde_u.setSymmetryOn()
66             self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)
67
68             self.__pde_proj=LinearPDE(domain,numEquations=1,numSolutions=1)
69             self.__pde_proj.setReducedOrderOn()
70             self.__pde_proj.setSymmetryOn()
71             # self.__pde_proj.setSolverMethod(LinearPDE.LUMPING)
72
74            self.eta=eta
75            A =self.__pde_u.createCoefficientOfGeneralPDE("A")
76        self.__pde_u.setValue(A=Data())
77            for i in range(self.domain.getDim()):
78            for j in range(self.domain.getDim()):
79                A[i,j,j,i] += 1.
80                A[i,j,i,j] += 1.
81            # self.__inv_eta=util.interpolate(self.eta,ReducedFunction(self.domain))
83
84            self.__pde_proj.setValue(D=1/eta)
85            self.__pde_proj.setValue(Y=1.)
86            self.__inv_eta=util.interpolate(self.__pde_proj.getSolution(),ReducedFunction(self.domain))
87
88          def B(self,arg):
89             a=util.div(arg, ReducedFunction(self.domain))
90             return a-util.integrate(a)/self.vol
91
92          def inner(self,p0,p1):
93             return util.integrate(p0*p1)
94
95          def getStress(self,u):
97             return 2.*self.eta*util.symmetric(mg)
98          def getEtaEffective(self):
99             return self.eta
100
101          def solve_A(self,u,p):
102             """
103             solves Av=f-Au-B^*p (v=0 on fixed_u_mask)
104             """
105             self.__pde_u.setTolerance(self.getSubProblemTolerance())
106             self.__pde_u.setValue(X=-self.getStress(u),X_reduced=-p*util.kronecker(self.domain))
107             return  self.__pde_u.getSolution(verbose=self.show_details)
108
109
110          def solve_prec(self,p):
111            a=self.__inv_eta*p
112            return a-util.integrate(a)/self.vol
113
114          def stoppingcriterium(self,Bv,v,p):
115              n_r=util.sqrt(self.inner(Bv,Bv))
117              if self.verbose: print "PCG step %s: L2(div(v)) = %s, L2(grad(v))=%s"%(self.iter,n_r,n_v) , util.Lsup(v)
118              if self.iter == 0: self.__n_v=n_v;
119              self.__n_v, n_v_old =n_v, self.__n_v
120              self.iter+=1
121              if self.iter>1 and n_r <= n_v*self.getTolerance() and abs(n_v_old-self.__n_v) <= n_v * self.getTolerance():
122                  if self.verbose: print "PCG terminated after %s steps."%self.iter
123                  return True
124              else:
125                  return False
126
127
129        """        """
130        solves        solves
# Line 103  class StokesProblemCartesian(Homogeneous Line 190  class StokesProblemCartesian(Homogeneous
190           beta=(1/self.vol)*util.integrate(p1)           beta=(1/self.vol)*util.integrate(p1)
#print "NORM",alfa,beta,util.integrate((p0-alfa)*(p1-beta))+util.integrate(util.inner(v0,v1))
193           return util.integrate((p0-alfa)*(p1-beta)+((1/self.eta)**2)*util.inner(v0,v1))           return util.integrate((p0-alfa)*(p1-beta)+((1/self.eta)**2)*util.inner(v0,v1))
194
195
# Line 147  class StokesProblemCartesian(Homogeneous Line 233  class StokesProblemCartesian(Homogeneous
233            if self.iter == 0: self.__n_v=n_v;            if self.iter == 0: self.__n_v=n_v;
234            self.__n_v, n_v_old =n_v, self.__n_v            self.__n_v, n_v_old =n_v, self.__n_v
235            self.iter+=1            self.iter+=1
print abs(n_v_old-self.__n_v), n_v, self.getTolerance()
236            if self.iter>1 and n_r <= n_v*self.getTolerance() and abs(n_v_old-self.__n_v) <= n_v * self.getTolerance():            if self.iter>1 and n_r <= n_v*self.getTolerance() and abs(n_v_old-self.__n_v) <= n_v * self.getTolerance():
237                if self.verbose: print "PCG terminated after %s steps."%self.iter                if self.verbose: print "PCG terminated after %s steps."%self.iter
238                return True                return True

Legend:
 Removed from v.1658 changed lines Added in v.1659