Python LinearPDE.setReducedOrderOn方法代码示例

# 需要导入模块: from esys.escript.linearPDEs import LinearPDE [as 别名]
# 或者: from esys.escript.linearPDEs.LinearPDE import setReducedOrderOn [as 别名]
class SimpleStokesProblem(SaddlePointProblem):
      """
      simple example of saddle point problem
      """
      def __init__(self,domain):
         super(SimpleStokesProblem, self).__init__(self)

         self.__pde_u=LinearPDE(domain)
         self.__pde_u.setSymmetryOn()
         self.__pde_u.setValue(A=identityTensor4(dom))

         self.__pde_p=LinearPDE(domain)
         self.__pde_p.setReducedOrderOn()
         self.__pde_p.setSymmetryOn()
         self.__pde_p.setValue(D=1.)

      def initialize(self,f=Data(),fixed_u_mask=Data()):
         self.__pde_u.setValue(q=fixed_u_mask,Y=f)
      def inner(self,p0,p1):
         return integrate(p0*p1,Function(self.__pde_p.getDomain()))

      def solve_f(self,u,p,tol=1.e-8):
         self.__pde_u.setTolerance(tol)
         self.__pde_u.setValue(X=grad(u)+p*kronecker(self.__pde_u.getDomain()))
         return  self.__pde_u.getSolution()
      def solve_g(self,u,tol=1.e-8):
         self.__pde_p.setTolerance(tol)
         self.__pde_p.setValue(X=-u) 
         dp=self.__pde_p.getSolution()
         return  dp

# 需要导入模块: from esys.escript.linearPDEs import LinearPDE [as 别名]
# 或者: from esys.escript.linearPDEs.LinearPDE import setReducedOrderOn [as 别名]
class StokesProblem(SaddlePointProblem):
      """
      simple example of saddle point problem
      """
      def __init__(self,domain,debug=False):
         super(StokesProblem, self).__init__(self,debug)
         self.domain=domain
         self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
         self.__pde_u.setSymmetryOn()

         self.__pde_p=LinearPDE(domain)
         self.__pde_p.setReducedOrderOn()
         self.__pde_p.setSymmetryOn()

      def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1):
         self.eta=eta
         A =self.__pde_u.createCoefficientOfGeneralPDE("A")
         for i in range(self.domain.getDim()):
           for j in range(self.domain.getDim()):
             A[i,j,j,i] += self.eta
             A[i,j,i,j] += self.eta
         self.__pde_p.setValue(D=1/self.eta)
         self.__pde_u.setValue(A=A,q=fixed_u_mask,Y=f)

      def inner(self,p0,p1):
         return integrate(p0*p1,Function(self.__pde_p.getDomain()))

      def solve_f(self,u,p,tol=1.e-8):
         self.__pde_u.setTolerance(tol)
         g=grad(u)
         self.__pde_u.setValue(X=self.eta*symmetric(g)+p*kronecker(self.__pde_u.getDomain()))
         return  self.__pde_u.getSolution()

      def solve_g(self,u,tol=1.e-8):
         self.__pde_p.setTolerance(tol)
         self.__pde_p.setValue(X=-u) 
         dp=self.__pde_p.getSolution()
         return  dp

# 需要导入模块: from esys.escript.linearPDEs import LinearPDE [as 别名]
# 或者: from esys.escript.linearPDEs.LinearPDE import setReducedOrderOn [as 别名]
class TemperatureAdvection(Model):
       """

       The conservation of internal heat energy is given by

       *rho c_p ( dT/dt+v[j] * grad(T)[j])-grad(\kappa grad(T)_{,i}=Q*

       *n_i \kappa T_{,i}=0*

       it is assummed that *\rho c_p* is constant in time.

       solved by Taylor Galerkin method

       """
       def __init__(self,**kwargs):
           super(TemperatureAdvection, self).__init__(**kwargs)
           self.declareParameter(domain=None, \
                                 temperature=1., \
                                 velocity=numpy.zeros([3]),
                                 density=1., \
                                 heat_capacity=1., \
                                 thermal_permabilty=1., \
                                 # reference_temperature=0., \
                                 # radiation_coefficient=0., \
                                 thermal_source=0., \
                                 fixed_temperature=0.,
                                 location_fixed_temperature=Data(),
                                 safety_factor=0.1)

       def doInitialization(self):
           self.__pde=LinearPDE(self.domain)
           self.__pde.setSymmetryOn()
           self.__pde.setReducedOrderOn()
           self.__pde.getSolverOptions().setSolverMethod(SolverOptions.LUMPING)
           self.__pde.setValue(D=self.heat_capacity*self.density)

       def getSafeTimeStepSize(self,dt):
           """
           returns new step size
           """
           h=self.domain.getSize()
           return self.safety_factor*inf(h**2/(h*abs(self.heat_capacity*self.density)*length(self.velocity)+self.thermal_permabilty))

       def G(self,T,alpha):
           """
           tangential operator for taylor galerikin
           """
           g=grad(T)
           self.__pde.setValue(X=-self.thermal_permabilty*g, \
                               Y=self.thermal_source-self.__rhocp*inner(self.velocity,g), \
                               r=(self.__fixed_T-self.temperature)*alpha,\
                               q=self.location_fixed_temperature)
           return self.__pde.getSolution()
           

       def doStepPostprocessing(self,dt):
           """
           perform taylor galerkin step
           """
           T=self.temperature
           self.__rhocp=self.heat_capacity*self.density
           self.__fixed_T=self.fixed_temperature
           self.temperature=dt*self.G(dt/2*self.G(T,1./dt)+T,1./dt)+T
           self.trace("Temperature range is %e %e"%(inf(self.temperature),sup(self.temperature)))

# 需要导入模块: from esys.escript.linearPDEs import LinearPDE [as 别名]
# 或者: from esys.escript.linearPDEs.LinearPDE import setReducedOrderOn [as 别名]
eta2 = 1.0e2         # VISCOSITY OF THE FLUID ON TOP
penalty = 1.0e3      # PENALTY FACTOR FOT THE PENALTY METHOD
g=10.                # GRAVITY
t_step = 0
t_step_end = 2000
reinit_max = 30      # NUMBER OF ITERATIONS DURING THE REINITIALISATION PROCEDURE
reinit_each = 3      # NUMBER OF TIME STEPS BETWEEN TWO REINITIALISATIONS
h = Lsup(mesh.getSize())
numDim = mesh.getDim()
smooth = h*2.0       # SMOOTHING PARAMETER FOR THE TRANSITION ACROSS THE INTERFACE

### DEFINITION OF THE PDE ###
velocityPDE = LinearPDE(mesh, numEquations=numDim)

advectPDE = LinearPDE(mesh)
advectPDE.setReducedOrderOn()
advectPDE.setValue(D=1.0)
advectPDE.setSolverMethod(solver=LinearPDE.DIRECT)

reinitPDE = LinearPDE(mesh, numEquations=1)
reinitPDE.setReducedOrderOn()
reinitPDE.setSolverMethod(solver=LinearPDE.LUMPING)
my_proj=Projector(mesh)

### BOUNDARY CONDITIONS ###
xx = mesh.getX()[0]
yy = mesh.getX()[1]
zz = mesh.getX()[2]
top = whereZero(zz-l1)
bottom = whereZero(zz)
left = whereZero(xx)

# 需要导入模块: from esys.escript.linearPDEs import LinearPDE [as 别名]
# 或者: from esys.escript.linearPDEs.LinearPDE import setReducedOrderOn [as 别名]
class MultiScale(object):
   """
   problem description
   1. displacement:
   -(A_{ijkl} u_{k,l})_{,j} = -X_{ij,j} + Y_i
   Neumann boundary: n_j A_{ijkl} u_{k,l} = n_j X_{ij} + y_i
   Dirichlet boundary: u_i = r_i where q_i > 0
   :var u: unknown vector, displacement
   :var A: elastic tensor / tangent operator
   :var X: tensor, minus old stress
   :var Y: vector, gamma - grad(p)
   :var y: vector, Neumann bc traction
   :var q: vector, Dirichlet bc mask
   :var r: vector, Dirichlet bc value
   2. pore pressure:
   -(A_{ij} p_{,j})_{,i} + D p = Y
   Neumann boundary: n_j A_{jl} p_{,l} = y
   Dirichlet boundary: p = r where q > 0
   :var p: unknown scalar, pore pressure
   :var A: permeability tensor
   :var D: scalar, n / (K_f dt)
   :var Y: scalar, -dot(u)_{i,i} + n p_pre / (K_f dt)
   :var y: scalar, Neumann bc flux
   :var q: scalar, Dirichlet bc mask
   :var r: scalar, Dirichlet bc value
   """
   def __init__(self,domain,pore0=0.,perm=1.e-5,kf=2.2e9,dt=0.001,ng=1,useMPI=False,np=1,rtol=1.e-2):
      """
      initialization of the problem, i.e. model constructor
      :param domain: type Domain, domain of the problem
      :param pore0: type float, initial pore pressure
      :param perm: type float, d^2/(150 mu_f) in KC equation
      :param kf: type float, bulk modulus of the fluid
      :param dt: type float, time step for calculation
      :param ng: type integer, number of Gauss points
      :param useMPI: type boolean, use MPI or not
      :param np: type integer, number of processors
      :param rtol: type float, relevative tolerance for global convergence
      """
      self.__domain=domain
      self.__upde=LinearPDE(domain,numEquations=domain.getDim(),numSolutions=domain.getDim())
      self.__ppde=LinearPDE(domain,numEquations=1,numSolutions=1)
      # use reduced interpolation for pore pressure
      self.__ppde.setReducedOrderOn()
      try:
            self.__upde.getSolverOptions().setSolverMethod(SolverOptions.DIRECT)
            self.__ppde.getSolverOptions().setSolverMethod(SolverOptions.DIRECT)
      except:
            #import time
            print("=======================================================================")
            print("For better performance compile python-escript with direct solver method")
            print("=======================================================================")
            input("Press Enter to continue...")
            #time.sleep(5)
      self.__upde.setSymmetryOn()
      self.__ppde.setSymmetryOn()
      
      self.__dt=dt
      self.__bulkFluid=kf
      self.__numGaussPoints=ng
      self.__rtol=rtol
      self.__stress=escript.Tensor(0,escript.Function(domain))
      self.__S=escript.Tensor4(0,escript.Function(domain))
      self.__pool=get_pool(mpi=useMPI,threads=np)
      self.__scenes=self.__pool.map(initLoad,list(range(ng)))
      st = self.__pool.map(getStressAndTangent2D,self.__scenes)
      for i in range(ng):
         self.__stress.setValueOfDataPoint(i,st[i][0])
         self.__S.setValueOfDataPoint(i,st[i][1])
      self.__strain=escript.Tensor(0,escript.Function(domain))
      self.__pore=escript.Scalar(pore0,escript.ReducedSolution(domain))
      self.__pgauss=util.interpolate(pore0,escript.Function(domain))
      self.__permeability=perm
      self.__meanStressRate=escript.Scalar(0,escript.Function(domain))
      self.__r=escript.Vector(0,escript.Solution(domain)) #Dirichlet BC for u

   def initialize(self, b=escript.Data(), f=escript.Data(), umsk=escript.Data(), uvalue=escript.Data(), flux=escript.Data(), pmsk=escript.Data(), pvalue=escript.Data()):
      """
      initialize the model for each time step, e.g. assign parameters
      :param b: type vector, body force on FunctionSpace, e.g. gravity
      :param f: type vector, boundary traction on FunctionSpace (FunctionOnBoundary)
      :param umsk: type vector, mask of location for Dirichlet boundary
      :param uvalue: type vector, specified displacement for Dirichlet boundary
      """
      self.__upde.setValue(Y=b,y=f,q=umsk,r=uvalue)
      self.__ppde.setValue(y=flux,q=pmsk,r=pvalue)
      self.__r=uvalue                  
   
   def getDomain(self):
      """
      return model domain
      """
      return self.__domain
      
   def setTimeStep(self,dt=1.0):
      self.__dt = dt

   def getCurrentPacking(self,pos=(),time=0,prefix=''):
      if len(pos) == 0: # output all Gauss points packings
         self.__pool.map(outputPack,list(zip(self.__scenes,repeat(time),repeat(prefix))))
#.........这里部分代码省略.........