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

revision 1481 by artak, Wed Apr 9 00:45:47 2008 UTC revision 1639 by gross, Mon Jul 14 08:55:25 2008 UTC
# Line 36  __date__="\$Date:\$" Line 36  __date__="\$Date:\$"
36  from escript import *  from escript import *
37  import util  import util
38  from linearPDEs import LinearPDE  from linearPDEs import LinearPDE
40
42        """        """
# Line 63  class StokesProblemCartesian(Homogeneous Line 63  class StokesProblemCartesian(Homogeneous
63           self.vol=util.integrate(1.,Function(self.domain))           self.vol=util.integrate(1.,Function(self.domain))
64           self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())           self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
65           self.__pde_u.setSymmetryOn()           self.__pde_u.setSymmetryOn()
66           self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)           # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)
67
68           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
69           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
# Line 82  class StokesProblemCartesian(Homogeneous Line 82  class StokesProblemCartesian(Homogeneous
82          for j in range(self.domain.getDim()):          for j in range(self.domain.getDim()):
83              A[i,j,j,i] += 1.              A[i,j,j,i] += 1.
84              A[i,j,i,j] += 1.              A[i,j,i,j] += 1.
85      self.__pde_prec.setValue(D=1./self.eta)      self.__pde_prec.setValue(D=1/self.eta)
87
88        def B(self,arg):        def B(self,arg):
# Line 96  class StokesProblemCartesian(Homogeneous Line 96  class StokesProblemCartesian(Homogeneous
96           s1=util.interpolate(p1,Function(self.domain))           s1=util.interpolate(p1,Function(self.domain))
97           return util.integrate(s0*s1)           return util.integrate(s0*s1)
98
99          def inner_a(self,a0,a1):
100             p0=util.interpolate(a0[1],Function(self.domain))
101             p1=util.interpolate(a1[1],Function(self.domain))
102             alfa=(1/self.vol)*util.integrate(p0)
103             beta=(1/self.vol)*util.integrate(p1)
106         #print "NORM",alfa,beta,util.integrate((p0-alfa)*(p1-beta))+util.integrate(util.inner(v0,v1))
107             return util.integrate((p0-alfa)*(p1-beta)+((1/self.eta)**2)*util.inner(v0,v1))
108
109
110        def getStress(self,u):        def getStress(self,u):
112           return 2.*self.eta*util.symmetric(mg)           return 2.*self.eta*util.symmetric(mg)
113          def getEtaEffective(self):
114             return self.eta
115
116        def solve_A(self,u,p):        def solve_A(self,u,p):
117           """           """
# Line 108  class StokesProblemCartesian(Homogeneous Line 121  class StokesProblemCartesian(Homogeneous
121           self.__pde_u.setValue(X=-self.getStress(u)-p*util.kronecker(self.domain))           self.__pde_u.setValue(X=-self.getStress(u)-p*util.kronecker(self.domain))
122           return  self.__pde_u.getSolution(verbose=self.show_details)           return  self.__pde_u.getSolution(verbose=self.show_details)
123
124
125        def solve_prec(self,p):        def solve_prec(self,p):
126         #proj=Projector(domain=self.domain, reduce = True, fast=False)
127             self.__pde_prec.setTolerance(self.getSubProblemTolerance())
128             self.__pde_prec.setValue(Y=p)
129             q=self.__pde_prec.getSolution(verbose=self.show_details)
130             return q
131
132          def solve_prec1(self,p):
133         #proj=Projector(domain=self.domain, reduce = True, fast=False)
134           self.__pde_prec.setTolerance(self.getSubProblemTolerance())           self.__pde_prec.setTolerance(self.getSubProblemTolerance())
135           self.__pde_prec.setValue(Y=p)           self.__pde_prec.setValue(Y=p)
136           q=self.__pde_prec.getSolution(verbose=self.show_details)           q=self.__pde_prec.getSolution(verbose=self.show_details)
137         q0=util.interpolate(q,Function(self.domain))
138             print util.inf(q*q0),util.sup(q*q0)
139             q-=(1/self.vol)*util.integrate(q0)
140             print util.inf(q*q0),util.sup(q*q0)
141           return q           return q
142
143        def stoppingcriterium(self,Bv,v,p):        def stoppingcriterium(self,Bv,v,p):
144            n_r=util.sqrt(self.inner(Bv,Bv))            n_r=util.sqrt(self.inner(Bv,Bv))
146            if self.verbose: print "PCG step %s: L2(div(v)) = %s, Lsup(v)=%s"%(self.iter,n_r,n_v)            if self.verbose: print "PCG step %s: L2(div(v)) = %s, L2(grad(v))=%s"%(self.iter,n_r,n_v)
147              if self.iter == 0: self.__n_v=n_v;
148              self.__n_v, n_v_old =n_v, self.__n_v
149            self.iter+=1            self.iter+=1
150            if n_r <= self.vol**(1./2.-1./self.domain.getDim())*n_v*self.getTolerance():            print abs(n_v_old-self.__n_v), n_v, self.getTolerance()
151              if self.iter>1 and n_r <= n_v*self.getTolerance() and abs(n_v_old-self.__n_v) <= n_v * self.getTolerance():
152                if self.verbose: print "PCG terminated after %s steps."%self.iter                if self.verbose: print "PCG terminated after %s steps."%self.iter
153                return True                return True
154            else:            else:
155                return False                return False
156        def stoppingcriterium_GMRES(self,norm_r,norm_b):        def stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None):
157            if self.verbose: print "GMRES step %s: L2(r) = %s, L2(b)*TOL=%s"%(self.iter,norm_r,norm_b*self.getTolerance())        if TOL==None:
158                 TOL=self.getTolerance()
159              if self.verbose: print "%s step %s: L2(r) = %s, L2(b)*TOL=%s"%(solver,self.iter,norm_r,norm_b*TOL)
160            self.iter+=1            self.iter+=1
161            if norm_r <= norm_b*self.getTolerance():
162                if self.verbose: print "GMRES terminated after %s steps."%self.iter            if norm_r <= norm_b*TOL:
163                  if self.verbose: print "%s terminated after %s steps."%(solver,self.iter)
164                return True                return True
165            else:            else:
166                return False                return False
167
168        def stoppingcriterium_MINRES(self,norm_r,norm_Ax):
if self.verbose: print "MINRES step %s: L2(r) = %s, L2(b)*TOL=%s"%(self.iter,norm_r,norm_Ax*self.getTolerance())
self.iter+=1
if norm_r <= norm_Ax*self.getTolerance():
if self.verbose: print "MINRES terminated after %s steps."%self.iter
return True
else:
return False

Legend:
 Removed from v.1481 changed lines Added in v.1639