Python LinearPDE.setSolverMethod方法代码示例

# 需要导入模块: from esys.escript.linearPDEs import LinearPDE [as 别名]
# 或者: from esys.escript.linearPDEs.LinearPDE import setSolverMethod [as 别名]
def wavePropagation(dom,rho,mu,lmbd,eta):
   x=Function(dom).getX()
   # ... open new PDE ...
   mypde=LinearPDE(dom)
   mypde.setSolverMethod(LinearPDE.LUMPING)
   k=kronecker(Function(dom))
   mypde.setValue(D=k*rho)

   dt=(1./5.)*inf(dom.getSize()/sqrt((2*mu+lmbd)/rho))
   if output: print("time step size = ",dt)
   # ... set initial values ....
   n=0
   t=0
   t_write=0.
   n_write=0
   # initial value of displacement at point source is constant (U0=0.01)
   # for first two time steps
   u=Vector(0.,Solution(dom))
   v=Vector(0.,Solution(dom))
   a=Vector(0.,Solution(dom))
   a2=Vector(0.,Solution(dom))
   v=Vector(0.,Solution(dom))

   if not os.path.isdir(WORKDIR): os.mkdir(WORKDIR)

   starttime = time.clock()
   while t<t_end and n<n_end:
     if output: print(n+1,"-th time step t ",t+dt," max u and F: ",Lsup(u), end=' ')
     # prediction:
     u_pr=u+dt*v+(dt**2/2)*a+(dt**3/6)*a2
     v_pr=v+dt*a+(dt**2/2)*a2
     a_pr=a+dt*a2
     # ... get current stress ....
     eps=symmetric(grad(u_pr))
     stress=lmbd*trace(eps)*k+2*mu*eps
     # ... force due to event:
     if abs(t-tc)<5*tc_length:
        F=exp(-((t-tc)/tc_length)**2)*exp(-(length(x-xc)/src_radius)**2)*event
        if output: print(Lsup(F))
     else:
        if output: print(0.)
     # ... get new acceleration ....
     mypde.setValue(X=-stress,Y=F-eta*v_pr)
     a=mypde.getSolution()
     # ... get new displacement ...
     da=a-a_pr
     u=u_pr+(dt**2/12.)*da
     v=v_pr+(5*dt/12.)*da
     a2+=da/dt
     # ... save current acceleration in units of gravity and displacements 
     if output:
          if t>=t_write: 
             saveVTK(os.path.join(WORKDIR,"disp.%i.vtu"%n_write),displacement=u, amplitude=length(u))
             t_write+=dt_write
             n_write+=1
     t+=dt
     n+=1

   endtime = time.clock()
   totaltime = endtime-starttime
   global netotal
   print(">>number of elements: %s, total time: %s, per time step: %s <<"%(netotal,totaltime,totaltime/n))

# 需要导入模块: from esys.escript.linearPDEs import LinearPDE [as 别名]
# 或者: from esys.escript.linearPDEs.LinearPDE import setSolverMethod [as 别名]
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)
right = whereZero(xx-l0)
front = whereZero(yy)

# 需要导入模块: from esys.escript.linearPDEs import LinearPDE [as 别名]
# 或者: from esys.escript.linearPDEs.LinearPDE import setSolverMethod [as 别名]
class Mechanics(Model):
      """
      base class for mechanics models in updated lagrangean framework

      :ivar domain: domain (in)
      :ivar internal_force: =Data()
      :ivar external_force: =Data()
      :ivar prescribed_velocity: =Data()
      :ivar location_prescribed_velocity: =Data()
      :ivar temperature:  = None
      :ivar expansion_coefficient:  = 0.
      :ivar bulk_modulus: =1.
      :ivar shear_modulus: =1.
      :ivar rel_tol: =1.e-3
      :ivar abs_tol: =1.e-15
      :ivar max_iter: =10
      :ivar displacement: =None
      :ivar stress: =None
      """
      SAFTY_FACTOR_ITERATION=1./100.
      def __init__(self,**kwargs):
         """
         set up the model
         
         :keyword debug: debug flag
         :type debug: ``bool``
         """
         super(Mechanics, self).__init__(self,**kwargs)
         self.declareParameter(domain=None, \
                               displacement=None, \
                               stress=None, \
                               velocity=None, \
                               internal_force=None, \
                               external_force=None, \
                               prescribed_velocity=None, \
                               location_prescribed_velocity=None, \
                               temperature = None, \
                               expansion_coefficient = 0., \
                               bulk_modulus=2., \
                               shear_modulus=1., \
                               rel_tol=1.e-3,abs_tol=1.e-15,max_iter=10)
         self.__iter=0

      def doInitialization(self):
           """
           initialize model
           """
           if not self.displacement: self.displacement=Vector(0.,ContinuousFunction(self.domain))
           if not self.velocity: self.velocity=Vector(0.,ContinuousFunction(self.domain))
           if not self.stress: self.stress=Tensor(0.,ContinuousFunction(self.domain))
           if not self.internal_force: self.internal_force = Data()
           if not self.external_force: self.external_force = Data()
           if not self.prescribed_velocity: self.prescribed_velocity = Data()
           if not self.location_prescribed_velocity: self.location_prescribed_velocity =Data()
           # save the old values:
           self.__stress_safe=self.stress
           self.__temperature_safe=self.temperature
           self.__displacement_safe=self.displacement
           self.__velocity_safe=self.velocity
           self.__velocity_old=None
           self.__old_dt=None
           self.__very_old_dt=None
           # get node cooridnates and apply initial displacement
           self.__x=self.domain.getX()
           self.domain.setX(self.__x+self.displacement)
           # open PDE:
           self.__pde=LinearPDE(self.domain)
           self.__pde.setSolverMethod(self.__pde.DIRECT)
           self.__solver_options=self.__pde.getSolverOptions()
           self.__solver_options.setSolverMethod(self.__solver_options.DIRECT)
           self.__solver_options.setVerbosity(self.debug)

           # self.__pde.setSymmetryOn()

      def doStepPreprocessing(self,dt):
            """
            step up pressure iteration

            if run within a time dependend problem extrapolation of pressure from previous time steps is used to
            get an initial guess (that needs some work!!!!!!!)
            """
            # reset iteration counters:
            self.__iter=0
            self.__diff=self.UNDEF_DT
            # set initial guesses for the iteration:
            self.displacement=self.__displacement_safe
            self.stress=self.__stress_safe
            self.velocity=self.__velocity_safe
            # update geometry
            self.domain.setX(self.__x+self.displacement)

      def doStep(self,dt):
          """
          """
          self.__iter+=1
          k3=kronecker(self.domain)
          # set new thermal stress increment
          if self.temperature == None:
             self.deps_th=0.
          else:
#.........这里部分代码省略.........