/[escript]/branches/symbolic_from_3470/escript/py_src/flows.py
ViewVC logotype

Diff of /branches/symbolic_from_3470/escript/py_src/flows.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3867 by caltinay, Tue Jan 31 04:55:05 2012 UTC revision 3868 by caltinay, Thu Mar 15 06:07:08 2012 UTC
# Line 86  class DarcyFlow(object): Line 86  class DarcyFlow(object):
86    
87        if self.solver  == self.EVAL:        if self.solver  == self.EVAL:
88           self.__pde_v=None           self.__pde_v=None
89       if self.verbose: print("DarcyFlow: simple solver is used.")           if self.verbose: print("DarcyFlow: simple solver is used.")
90    
91        elif self.solver  == self.POST:        elif self.solver  == self.POST:
92       if util.inf(w)<0.:           if util.inf(w)<0.:
93          raise ValueError("Weighting factor must be non-negative.")              raise ValueError("Weighting factor must be non-negative.")
94       if self.verbose: print("DarcyFlow: global postprocessing of flux is used.")           if self.verbose: print("DarcyFlow: global postprocessing of flux is used.")
95           self.__pde_v=LinearPDESystem(domain)           self.__pde_v=LinearPDESystem(domain)
96           self.__pde_v.setSymmetryOn()           self.__pde_v.setSymmetryOn()
97           if self.useReduced: self.__pde_v.setReducedOrderOn()           if self.useReduced: self.__pde_v.setReducedOrderOn()
98       self.w=w           self.w=w
99           self.l=util.vol(self.domain)**(1./self.domain.getDim()) # length scale           self.l=util.vol(self.domain)**(1./self.domain.getDim()) # length scale
100    
101        elif self.solver  == self.SMOOTH:        elif self.solver  == self.SMOOTH:
102           self.__pde_v=LinearPDESystem(domain)           self.__pde_v=LinearPDESystem(domain)
103           self.__pde_v.setSymmetryOn()           self.__pde_v.setSymmetryOn()
104           if self.useReduced: self.__pde_v.setReducedOrderOn()           if self.useReduced: self.__pde_v.setReducedOrderOn()
105       if self.verbose: print("DarcyFlow: flux smoothing is used.")           if self.verbose: print("DarcyFlow: flux smoothing is used.")
106       self.w=0           self.w=0
107    
108        self.__f=escript.Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))        self.__f=escript.Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))
109        self.__g=escript.Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))        self.__g=escript.Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
# Line 143  class DarcyFlow(object): Line 143  class DarcyFlow(object):
143                            
144        if permeability!=None:        if permeability!=None:
145            
146       perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))           perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
147           self.perm_scale=util.Lsup(util.length(perm))           self.perm_scale=util.Lsup(util.length(perm))
148       if self.verbose: print(("DarcyFlow: permeability scaling factor = %e."%self.perm_scale))           if self.verbose: print(("DarcyFlow: permeability scaling factor = %e."%self.perm_scale))
149           perm=perm*(1./self.perm_scale)           perm=perm*(1./self.perm_scale)
150                    
151       if perm.getRank()==0:           if perm.getRank()==0:
152    
153          perm_inv=(1./perm)              perm_inv=(1./perm)
154          perm_inv=perm_inv*util.kronecker(self.domain.getDim())              perm_inv=perm_inv*util.kronecker(self.domain.getDim())
155          perm=perm*util.kronecker(self.domain.getDim())              perm=perm*util.kronecker(self.domain.getDim())
156                    
157                    
158       elif perm.getRank()==2:           elif perm.getRank()==2:
159          perm_inv=util.inverse(perm)              perm_inv=util.inverse(perm)
160       else:           else:
161          raise ValueError("illegal rank of permeability.")              raise ValueError("illegal rank of permeability.")
162                    
163       self.__permeability=perm           self.__permeability=perm
164       self.__permeability_inv=perm_inv           self.__permeability_inv=perm_inv
165            
166           #====================           #====================
167       self.__pde_p.setValue(A=self.__permeability)           self.__pde_p.setValue(A=self.__permeability)
168           if self.solver  == self.EVAL:           if self.solver  == self.EVAL:
169                pass # no extra work required                pass # no extra work required
170           elif self.solver  == self.POST:           elif self.solver  == self.POST:
171          k=util.kronecker(self.domain.getDim())                k=util.kronecker(self.domain.getDim())
172          self.omega = self.w*util.length(perm_inv)*self.l*self.domain.getSize()                self.omega = self.w*util.length(perm_inv)*self.l*self.domain.getSize()
173          self.__pde_v.setValue(D=self.__permeability_inv, A=self.omega*util.outer(k,k))                self.__pde_v.setValue(D=self.__permeability_inv, A=self.omega*util.outer(k,k))
174           elif self.solver  == self.SMOOTH:           elif self.solver  == self.SMOOTH:
175          self.__pde_v.setValue(D=self.__permeability_inv)              self.__pde_v.setValue(D=self.__permeability_inv)
176    
177        if g != None:        if g != None:
178      g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))          g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
179      if g.isEmpty():          if g.isEmpty():
180            g=Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))               g=Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
181      else:          else:
182          if not g.getShape()==(self.domain.getDim(),): raise ValueError("illegal shape of g")               if not g.getShape()==(self.domain.getDim(),): raise ValueError("illegal shape of g")
183      self.__g=g          self.__g=g
184        if f !=None:        if f !=None:
185       f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("Y"))           f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
186       if f.isEmpty():                 if f.isEmpty():      
187            f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))               f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
188       else:           else:
189           if f.getRank()>0: raise ValueError("illegal rank of f.")               if f.getRank()>0: raise ValueError("illegal rank of f.")
190       self.__f=f           self.__f=f
191    
192     def getSolverOptionsFlux(self):     def getSolverOptionsFlux(self):
193        """        """
# Line 263  class DarcyFlow(object): Line 263  class DarcyFlow(object):
263          elif self.solver  == self.POST or self.solver  == self.SMOOTH:          elif self.solver  == self.POST or self.solver  == self.SMOOTH:
264              self.__pde_v.setValue(Y=util.tensor_mult(self.__permeability_inv,self.__g * 1./self.perm_scale)-util.grad(p))              self.__pde_v.setValue(Y=util.tensor_mult(self.__permeability_inv,self.__g * 1./self.perm_scale)-util.grad(p))
265              if u0 == None:              if u0 == None:
266             self.__pde_v.setValue(r=escript.Data())                 self.__pde_v.setValue(r=escript.Data())
267          else:              else:
268                 if not isinstance(u0, escript.Data) : u0 = escript.Vector(u0, escript.Solution(self.domain))                 if not isinstance(u0, escript.Data) : u0 = escript.Vector(u0, escript.Solution(self.domain))
269             self.__pde_v.setValue(r=1./self.perm_scale * u0)                 self.__pde_v.setValue(r=1./self.perm_scale * u0)
270              u= self.__pde_v.getSolution() * self.perm_scale                 u= self.__pde_v.getSolution() * self.perm_scale
271      return u          return u
272                
273  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
274       """       """
# Line 313  class StokesProblemCartesian(Homogeneous Line 313  class StokesProblemCartesian(Homogeneous
313    
314           self.__pde_proj=LinearPDE(domain)           self.__pde_proj=LinearPDE(domain)
315           self.__pde_proj.setReducedOrderOn()           self.__pde_proj.setReducedOrderOn()
316       self.__pde_proj.setValue(D=1)           self.__pde_proj.setValue(D=1)
317           self.__pde_proj.setSymmetryOn()           self.__pde_proj.setSymmetryOn()
318    
319       def getSolverOptionsVelocity(self):       def getSolverOptionsVelocity(self):
# Line 322  class StokesProblemCartesian(Homogeneous Line 322  class StokesProblemCartesian(Homogeneous
322            
323       :rtype: `SolverOptions`       :rtype: `SolverOptions`
324       """       """
325       return self.__pde_v.getSolverOptions()           return self.__pde_v.getSolverOptions()
326       def setSolverOptionsVelocity(self, options=None):       def setSolverOptionsVelocity(self, options=None):
327           """           """
328       set the solver options for solving the equation for velocity.       set the solver options for solving the equation for velocity.
# Line 336  class StokesProblemCartesian(Homogeneous Line 336  class StokesProblemCartesian(Homogeneous
336       returns the solver options used  solve the equation for pressure.       returns the solver options used  solve the equation for pressure.
337       :rtype: `SolverOptions`       :rtype: `SolverOptions`
338       """       """
339       return self.__pde_prec.getSolverOptions()           return self.__pde_prec.getSolverOptions()
340       def setSolverOptionsPressure(self, options=None):       def setSolverOptionsPressure(self, options=None):
341           """           """
342       set the solver options for solving the equation for pressure.       set the solver options for solving the equation for pressure.
343       :param options: new solver  options       :param options: new solver  options
344       :type options: `SolverOptions`       :type options: `SolverOptions`
345       """       """
346       self.__pde_prec.setSolverOptions(options)           self.__pde_prec.setSolverOptions(options)
347    
348       def setSolverOptionsDiv(self, options=None):       def setSolverOptionsDiv(self, options=None):
349           """           """
# Line 353  class StokesProblemCartesian(Homogeneous Line 353  class StokesProblemCartesian(Homogeneous
353       :param options: new solver options       :param options: new solver options
354       :type options: `SolverOptions`       :type options: `SolverOptions`
355       """       """
356       self.__pde_proj.setSolverOptions(options)           self.__pde_proj.setSolverOptions(options)
357       def getSolverOptionsDiv(self):       def getSolverOptionsDiv(self):
358           """           """
359       returns the solver options for solving the equation to project the divergence of       returns the solver options for solving the equation to project the divergence of
# Line 361  class StokesProblemCartesian(Homogeneous Line 361  class StokesProblemCartesian(Homogeneous
361            
362       :rtype: `SolverOptions`       :rtype: `SolverOptions`
363       """       """
364       return self.__pde_proj.getSolverOptions()           return self.__pde_proj.getSolverOptions()
365    
366       def updateStokesEquation(self, v, p):       def updateStokesEquation(self, v, p):
367           """           """
# Line 388  class StokesProblemCartesian(Homogeneous Line 388  class StokesProblemCartesian(Homogeneous
388              k=util.kronecker(self.domain.getDim())              k=util.kronecker(self.domain.getDim())
389              kk=util.outer(k,k)              kk=util.outer(k,k)
390              self.eta=util.interpolate(eta, escript.Function(self.domain))              self.eta=util.interpolate(eta, escript.Function(self.domain))
391          self.__pde_prec.setValue(D=1/self.eta)              self.__pde_prec.setValue(D=1/self.eta)
392              self.__pde_v.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))              self.__pde_v.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))
393          if restoration_factor!=None:          if restoration_factor!=None:
394              n=self.domain.getNormal()              n=self.domain.getNormal()
# Line 412  class StokesProblemCartesian(Homogeneous Line 412  class StokesProblemCartesian(Homogeneous
412          :param surface_stress: normal surface stress          :param surface_stress: normal surface stress
413          :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar          :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
414          :param stress: initial stress          :param stress: initial stress
415      :type stress: `Tensor` object on `FunctionSpace` `Function` or similar          :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
416          """          """
417          self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)          self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
418    
# Line 425  class StokesProblemCartesian(Homogeneous Line 425  class StokesProblemCartesian(Homogeneous
425           :rtype: ``float``           :rtype: ``float``
426           """           """
427           self.__pde_proj.setValue(Y=-util.div(v))           self.__pde_proj.setValue(Y=-util.div(v))
428       self.getSolverOptionsDiv().setTolerance(tol)           self.getSolverOptionsDiv().setTolerance(tol)
429       self.getSolverOptionsDiv().setAbsoluteTolerance(0.)           self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
430           out=self.__pde_proj.getSolution()           out=self.__pde_proj.getSolution()
431           return out           return out
432    
# Line 475  class StokesProblemCartesian(Homogeneous Line 475  class StokesProblemCartesian(Homogeneous
475           """           """
476           self.updateStokesEquation(v,p)           self.updateStokesEquation(v,p)
477           self.__pde_v.setValue(Y=self.__f, y=self.__surface_stress)           self.__pde_v.setValue(Y=self.__f, y=self.__surface_stress)
478       self.getSolverOptionsVelocity().setTolerance(tol)           self.getSolverOptionsVelocity().setTolerance(tol)
479       self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)           self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
480           if self.__stress.isEmpty():           if self.__stress.isEmpty():
481              self.__pde_v.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))              self.__pde_v.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
482           else:           else:
# Line 515  class StokesProblemCartesian(Homogeneous Line 515  class StokesProblemCartesian(Homogeneous
515           :note: boundary conditions on p are zero.           :note: boundary conditions on p are zero.
516           """           """
517           self.__pde_prec.setValue(Y=Bv)           self.__pde_prec.setValue(Y=Bv)
518       self.getSolverOptionsPressure().setTolerance(tol)           self.getSolverOptionsPressure().setTolerance(tol)
519       self.getSolverOptionsPressure().setAbsoluteTolerance(0.)           self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
520           out=self.__pde_prec.getSolution()           out=self.__pde_prec.getSolution()
521           return out           return out

Legend:
Removed from v.3867  
changed lines
  Added in v.3868

  ViewVC Help
Powered by ViewVC 1.1.26