Python LinearPDE.getSolverOptions方法代码示例

# 需要导入模块: from esys.escript.linearPDEs import LinearPDE [as 别名]
# 或者: from esys.escript.linearPDEs.LinearPDE import getSolverOptions [as 别名]
    def setUpPDE(self):
        """
        Creates and returns the underlying PDE.

        :rtype: `LinearPDE`
        """
        if self.__pde is None:
            if not HAVE_DIRECT:
                raise ValueError("Either this build of escript or the current MPI configuration does not support direct solvers.")
            pde=LinearPDE(self.__domain, numEquations=2)
            D=pde.createCoefficient('D')
            A=pde.createCoefficient('A')
            A[0,:,0,:]=kronecker(self.__domain.getDim())
            A[1,:,1,:]=kronecker(self.__domain.getDim())
            pde.setValue(A=A, D=D)
            if self.__fixAtBottom:
                DIM=self.__domain.getDim()
                z = self.__domain.getX()[DIM-1]
                pde.setValue(q=whereZero(z-self.__BX[DIM-1][0])*[1,1])

            pde.getSolverOptions().setSolverMethod(SolverOptions.DIRECT)
            pde.getSolverOptions().setTolerance(self.__tol)
            pde.setSymmetryOff()
        else:
            pde=self.__pde
            pde.resetRightHandSideCoefficients()
        return pde

# 需要导入模块: from esys.escript.linearPDEs import LinearPDE [as 别名]
# 或者: from esys.escript.linearPDEs.LinearPDE import getSolverOptions [as 别名]
    def setUpPDE(self):
        """
        Return the underlying PDE.

        :rtype: `LinearPDE`
        """
        if self.__pde is None:
            dom=self.__domain
            x = dom.getX()
            DIM=dom.getDim()
            q=whereZero(x[DIM-1]-inf(x[DIM-1]))
            for i in range(DIM-1):
                xi=x[i]
                q+=whereZero(xi-inf(xi))+whereZero(xi-sup(xi))

            pde=LinearPDE(dom, numEquations=1)
            pde.getSolverOptions().setTolerance(self.__tol)
            pde.setSymmetryOn()
            A=pde.createCoefficient('A')
            X=pde.createCoefficient('X')
            pde.setValue(A=A, X=X, q=q)

        else:
            pde=self.__pde
            pde.resetRightHandSideCoefficients()
        return pde

# 需要导入模块: from esys.escript.linearPDEs import LinearPDE [as 别名]
# 或者: from esys.escript.linearPDEs.LinearPDE import getSolverOptions [as 别名]
    def getPotentialNumeric(self):
        """
        Returns 3 list each made up of a number of list containing primary, secondary and total
        potentials diferences. Each of the lists contain a list for each value of n.
        """
        primCon=self.primaryConductivity
        coords=self.domain.getX()
        tol=1e-8
        pde=LinearPDE(self.domain, numEquations=1)
        pde.getSolverOptions().setTolerance(tol)
        pde.setSymmetryOn()
        primaryPde=LinearPDE(self.domain, numEquations=1)
        primaryPde.getSolverOptions().setTolerance(tol)
        primaryPde.setSymmetryOn()
        DIM=self.domain.getDim()
        x=self.domain.getX()
        q=whereZero(x[DIM-1]-inf(x[DIM-1]))
        for i in xrange(DIM-1):
            xi=x[i]
            q+=whereZero(xi-inf(xi))+whereZero(xi-sup(xi))
        A = self.secondaryConductivity * kronecker(self.domain)
        APrimary = self.primaryConductivity * kronecker(self.domain)
        pde.setValue(A=A,q=q)
        primaryPde.setValue(A=APrimary,q=q)

        delPhiSecondaryList = []
        delPhiPrimaryList = []
        delPhiTotalList = []
        for i in range(1,self.n+1): # 1 to n
            maxR = self.numElectrodes - 1 - (2*i) #max amount of readings that will fit in the survey
            delPhiSecondary = []
            delPhiPrimary = []
            delPhiTotal = []
            for j in range(maxR):
                y_dirac=Scalar(0,DiracDeltaFunctions(self.domain))
                y_dirac.setTaggedValue(self.electrodeTags[j],self.current)
                y_dirac.setTaggedValue(self.electrodeTags[j + ((2*i) + 1)],-self.current)
                self.sources.append([self.electrodeTags[j], self.electrodeTags[j + ((2*i) + 1)]])
                primaryPde.setValue(y_dirac=y_dirac)
                numericPrimaryPot = primaryPde.getSolution()
                X=(primCon-self.secondaryConductivity) * grad(numericPrimaryPot)
                pde.setValue(X=X)
                u=pde.getSolution()
                loc=Locator(self.domain,[self.electrodes[j+i],self.electrodes[j+i+1]])
                self.samples.append([self.electrodeTags[j+i],self.electrodeTags[j+i+1]])
                valPrimary=loc.getValue(numericPrimaryPot)
                valSecondary=loc.getValue(u)
                delPhiPrimary.append(valPrimary[1]-valPrimary[0])
                delPhiSecondary.append(valSecondary[1]-valSecondary[0])
                delPhiTotal.append(delPhiPrimary[j]+delPhiSecondary[j])
            delPhiPrimaryList.append(delPhiPrimary)
            delPhiSecondaryList.append(delPhiSecondary)
            delPhiTotalList.append(delPhiTotal)

        self.delPhiPrimaryList=delPhiPrimaryList
        self.delPhiSecondaryList=delPhiSecondaryList
        self.delPhiTotalList = delPhiTotalList
        return [delPhiPrimaryList, delPhiSecondaryList, delPhiTotalList]

# 需要导入模块: from esys.escript.linearPDEs import LinearPDE [as 别名]
# 或者: from esys.escript.linearPDEs.LinearPDE import getSolverOptions [as 别名]
    def getPotential(self):
        """
        returns a list containing 3 lists one for each the primary, secondary
        and total potential.
        """


        primCon=self.primaryConductivity
        coords=self.domain.getX()
        pde=LinearPDE(self.domain, numEquations=1)
        tol=1e-8
        pde.getSolverOptions().setTolerance(tol)
        pde.setSymmetryOn()

        DIM=self.domain.getDim()
        x=self.domain.getX()
        q=whereZero(x[DIM-1]-inf(x[DIM-1]))
        for i in xrange(DIM-1):
            xi=x[i]
            q+=whereZero(xi-inf(xi))+whereZero(xi-sup(xi))
        A = self.secondaryConductivity * kronecker(self.domain)
        pde.setValue(A=A,q=q)

        delPhiSecondary = []
        delPhiPrimary = []
        delPhiTotal = []
        if(len(self.electrodes[0])==3):

            for i in range(self.numElectrodes-1):
                analyticRs=Data(0,(3,),ContinuousFunction(self.domain))
                analyticRs[0]=(coords[0]-self.electrodes[i][0])
                analyticRs[1]=(coords[1]-self.electrodes[i][1])
                analyticRs[2]=(coords[2])
                rsMag=(analyticRs[0]**2+analyticRs[1]**2+analyticRs[2]**2)**0.5
                analyticPrimaryPot=(self.current*(1./primCon))/(2*pi*(rsMag+(whereZero(rsMag)*0.0000001))) #the magic number 0.0000001 is to avoid devide by 0
                analyticRsPolePower=(analyticRs[0]**2+analyticRs[1]**2+analyticRs[2]**2)**1.5
                analyticRsPolePower = analyticRsPolePower+(whereZero(analyticRsPolePower)*0.0000001)
                gradUPrimary = Data(0,(3,),ContinuousFunction(self.domain))
                gradUPrimary[0] =(self.current/(2*pi*primCon)) * (analyticRs[0]/analyticRsPolePower)
                gradUPrimary[1] =(self.current/(2*pi*primCon)) * (analyticRs[1]/analyticRsPolePower)
                gradUPrimary[2] =(self.current/(2*pi*primCon)) * (analyticRs[2]/analyticRsPolePower)
                gradUPrimary=-gradUPrimary
                X=(primCon-self.secondaryConductivity) * gradUPrimary
                pde.setValue(X=X)
                u=pde.getSolution()
                loc=Locator(self.domain,self.electrodes[i+1])
                delPhiSecondary.append(loc.getValue(u))
                delPhiPrimary.append(loc.getValue(analyticPrimaryPot))
        else:
            raise NotImplementedError("2d forward model is not yet implemented")

        self.delPhiSecondary = delPhiSecondary
        self.delPhiPrimary = delPhiPrimary
        for i in range(len(delPhiPrimary)):
            delPhiTotal.append(delPhiPrimary[i] + delPhiSecondary[i])
        self.delPhiTotal=delPhiTotal
        return [delPhiPrimary, delPhiSecondary, delPhiTotal]

# 需要导入模块: from esys.escript.linearPDEs import LinearPDE [as 别名]
# 或者: from esys.escript.linearPDEs.LinearPDE import getSolverOptions [as 别名]
def runTest(dom, n=1, a=0, b=0):
   print("================== TEST : n= %s a=%s b=%s ================"%(n,a,b))
   DIM=dom.getDim()
   normal=dom.getNormal()
   mypde=LinearPDE(dom, numEquations=n, numSolutions=n)
   x=dom.getX()
   A=mypde.createCoefficient("A")
   D=mypde.createCoefficient("D")
   Y=mypde.createCoefficient("Y")
   y=mypde.createCoefficient("y")
   q=mypde.createCoefficient("q")
   if n==1:
      u_ref=Scalar(0.,Solution(dom))
      for j in range(DIM):
         q+=whereZero(x[j])
         A[j,j]=1
      y+=sin(sqrt(2.))*normal[0]
      Y+=b*x[0]*sin(sqrt(2.))
      D+=b
      u_ref=x[0]*sin(sqrt(2.))
   else:
      u_ref=Data(0.,(n,),Solution(dom))
      for i in range(n):
          for j in range(DIM):
             q[i]+=whereZero(x[j])
             A[i,j,i,j]=1
             if j == i%DIM: y[i]+=sin(i+sqrt(2.))*normal[j]
          for j in range(n):
                if i==j:
                   D[i,i]=b+(n-1)*a
                   Y[i]+=b*x[i%DIM]*sin(i+sqrt(2.))
                else:
                   D[i,j]=-a
                   Y[i]+=a*(x[i%DIM]*sin(i+sqrt(2.))-x[j%DIM]*sin(j+sqrt(2.)))
          u_ref[i]=x[i%DIM]*sin(i+sqrt(2.))
          
#    - u_{j,ii} + b*u_i+ a*sum_{k<>j}  (u_i-u_k) = F_j
          
   mypde.setValue(A=A, D=D, y=y, q=q, r=u_ref, Y=Y)
   mypde.getSolverOptions().setVerbosityOn()
   mypde.getSolverOptions().setTolerance(TOL)
   mypde.setSymmetryOn()
   u=mypde.getSolution()
   error=Lsup(u-u_ref)/Lsup(u_ref)
   print("error = ",error)
   if error > 10*TOL: print("XXXXXXXXXX Error to large ! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX")

# 需要导入模块: from esys.escript.linearPDEs import LinearPDE [as 别名]
# 或者: from esys.escript.linearPDEs.LinearPDE import getSolverOptions [as 别名]
    def getPDE(self, system):
        dim = self.domain.getDim()
        if system:
            pde=LinearPDE(self.domain, numEquations=dim)
        else:
            pde=LinearPDE(self.domain, numEquations=1)

        self._setCoefficients(pde, system)
        so = pde.getSolverOptions()
        so.setPackage(self.package)
        so.setSolverMethod(self.method)
        so.setPreconditioner(self.preconditioner)
        so.setTolerance(self.SOLVER_TOL)
        so.setVerbosity(self.SOLVER_VERBOSE)
        self._setSolverOptions(so)
        return pde, self._getSolution(system), self._getGrad(system)

# 需要导入模块: from esys.escript.linearPDEs import LinearPDE [as 别名]
# 或者: from esys.escript.linearPDEs.LinearPDE import getSolverOptions [as 别名]
   def test_StrongCoupled4(self):
        x=Solution(self.domain).getX()
        # --- set exact solution ----
        u_ex=Data(1.,(4,),Solution(self.domain))
        g_ex=Data(0.,(4,self.domain.getDim()), Solution(self.domain))
        A=Data(0.,(4,self.domain.getDim(),4,self.domain.getDim()), Function(self.domain))
        for i in range(self.domain.getDim()):
           u_ex[0]+= 1*(i+1)*x[i]
           g_ex[0,i]=1*(i+1)
           u_ex[1]+= 2*(i+1)*x[i]
           g_ex[1,i]=2*(i+1)
           u_ex[2]+= 3*(i+1)*x[i]
           g_ex[2,i]=3*(i+1)
           u_ex[3]+= 4*(i+1)*x[i]
           g_ex[3,i]=4*(i+1)
           A[0,i,0,i]=1
           A[1,i,1,i]=1
           A[2,i,2,i]=1
           A[3,i,3,i]=1

        a=-1./4.*0.99

        Y=Data(0.,(4,),Function(self.domain))
        Y[0]=     u_ex[0]+ a * u_ex[1]+ a * u_ex[2]+ a * u_ex[3]
        Y[1]= a * u_ex[0]+     u_ex[1]+ a * u_ex[2]+ a * u_ex[3]
        Y[2]= a * u_ex[0]+ a * u_ex[1]+     u_ex[2]+ a * u_ex[3]
        Y[3]= a * u_ex[0]+ a * u_ex[1]+ a * u_ex[2]+     u_ex[3]
        # create PDE:
        pde=LinearPDE(self.domain,numEquations=4)
        pde.setSymmetryOn()
        mask=whereZero(x[0])*[1,1,1,1]
        pde.setValue(r=u_ex,q=mask)
        pde.setValue(A=A,
                     D=kronecker(4)+a*(numpy.ones((4,4))-kronecker(4)),
                     Y=Y,
                     y=matrixmult(g_ex,self.domain.getNormal()))
        # -------- get the solution ---------------------------
        pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
        pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
        pde.getSolverOptions().setPreconditioner(SolverOptions.AMG)
        pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
        if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
        if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
        if MAX_LEVEL!=None: pde.getSolverOptions().setLevelMax(MAX_LEVEL)
        u=pde.getSolution()
        # -------- test the solution ---------------------------
        error=Lsup(u-u_ex)/Lsup(u_ex)
        self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)

# 需要导入模块: from esys.escript.linearPDEs import LinearPDE [as 别名]
# 或者: from esys.escript.linearPDEs.LinearPDE import getSolverOptions [as 别名]
class SplitRegularization(CostFunction):
    """
    The regularization term for the level set function ``m`` within the cost
    function J for an inversion:

    *J(m)=1/2 * sum_k integrate( mu[k] * ( w0[k] * m_k**2 * w1[k,i] * m_{k,i}**2) + sum_l<k mu_c[l,k] wc[l,k] * | curl(m_k) x curl(m_l) |^2*

    where w0[k], w1[k,i] and  wc[k,l] are non-negative weighting factors and
    mu[k] and mu_c[l,k] are trade-off factors which may be altered
    during the inversion. The weighting factors are normalized such that their
    integrals over the domain are constant:

    *integrate(w0[k] + inner(w1[k,:],1/L[:]**2))=scale[k]* volume(domain)*
    *integrate(wc[l,k]*1/L**4)=scale_c[k]* volume(domain) *

    """

    def __init__(
        self,
        domain,
        numLevelSets=1,
        w0=None,
        w1=None,
        wc=None,
        location_of_set_m=Data(),
        useDiagonalHessianApproximation=False,
        tol=1e-8,
        coordinates=None,
        scale=None,
        scale_c=None,
    ):
        """
        initialization.

        :param domain: domain
        :type domain: `Domain`
        :param numLevelSets: number of level sets
        :type numLevelSets: ``int``
        :param w0: weighting factor for the m**2 term. If not set zero is assumed.
        :type w0: ``Scalar`` if ``numLevelSets`` == 1 or `Data` object of shape
                  (``numLevelSets`` ,) if ``numLevelSets`` > 1
        :param w1: weighting factor for the grad(m_i) terms. If not set zero is assumed
        :type w1: ``Vector`` if ``numLevelSets`` == 1 or `Data` object of shape
                  (``numLevelSets`` , DIM) if ``numLevelSets`` > 1
        :param wc: weighting factor for the cross gradient terms. If not set
                   zero is assumed. Used for the case if ``numLevelSets`` > 1
                   only. Only values ``wc[l,k]`` in the lower triangle (l<k)
                   are used.
        :type wc: `Data` object of shape (``numLevelSets`` , ``numLevelSets``)
        :param location_of_set_m: marks location of zero values of the level
                                  set function ``m`` by a positive entry.
        :type location_of_set_m: ``Scalar`` if ``numLevelSets`` == 1 or `Data`
                object of shape (``numLevelSets`` ,) if ``numLevelSets`` > 1
        :param useDiagonalHessianApproximation: if True cross gradient terms
                    between level set components are ignored when calculating
                    approximations of the inverse of the Hessian Operator.
                    This can speed-up the calculation of the inverse but may
                    lead to an increase of the number of iteration steps in the
                    inversion.
        :type useDiagonalHessianApproximation: ``bool``
        :param tol: tolerance when solving the PDE for the inverse of the
                    Hessian Operator
        :type tol: positive ``float``

        :param coordinates: defines coordinate system to be used
        :type coordinates: ReferenceSystem` or `SpatialCoordinateTransformation`
        :param scale: weighting factor for level set function variation terms.
                      If not set one is used.
        :type scale: ``Scalar`` if ``numLevelSets`` == 1 or `Data` object of
                     shape (``numLevelSets`` ,) if ``numLevelSets`` > 1
        :param scale_c: scale for the cross gradient terms. If not set
                   one is assumed. Used for the case if ``numLevelSets`` > 1
                   only. Only values ``scale_c[l,k]`` in the lower triangle
                   (l<k) are used.
        :type scale_c: `Data` object of shape (``numLevelSets``,``numLevelSets``)

        """
        if w0 is None and w1 is None:
            raise ValueError("Values for w0 or for w1 must be given.")
        if wc is None and numLevelSets > 1:
            raise ValueError("Values for wc must be given.")

        self.__pre_input = None
        self.__pre_args = None
        self.logger = logging.getLogger("inv.%s" % self.__class__.__name__)
        self.__domain = domain
        DIM = self.__domain.getDim()
        self.__numLevelSets = numLevelSets
        self.__trafo = makeTransformation(domain, coordinates)
        self.__pde = LinearPDE(self.__domain, numEquations=self.__numLevelSets, numSolutions=self.__numLevelSets)
        self.__pde.getSolverOptions().setTolerance(tol)
        self.__pde.setSymmetryOn()
        self.__pde.setValue(A=self.__pde.createCoefficient("A"), D=self.__pde.createCoefficient("D"))
        try:
            self.__pde.setValue(q=location_of_set_m)
        except IllegalCoefficientValue:
            raise ValueError("Unable to set location of fixed level set function.")

        # =========== check the shape of the scales: ========================
        if scale is None:
#.........这里部分代码省略.........

# 需要导入模块: from esys.escript.linearPDEs import LinearPDE [as 别名]
# 或者: from esys.escript.linearPDEs.LinearPDE import getSolverOptions [as 别名]
    rtime_inc = tend / 20.0  # time increment to record
    # Check to make sure number of time steps is not too large.
    print("Time step size= ", h, "Expected number of outputs= ", tend / h)

    U0 = 0.005  # amplitude of point source
    # want a spherical source in the middle of area
    xc = [500, 500]  # with reference to mx,my this is the source location

    ####################################################DOMAIN CONSTRUCTION
    mydomain = Rectangle(l0=mx, l1=my, n0=ndx, n1=ndy)  # create the domain
    x = mydomain.getX()  # get the node locations of the domain

    ##########################################################ESTABLISH PDE
    mypde = LinearPDE(mydomain)  # create pde
    # turn lumping on for more efficient solving
    mypde.getSolverOptions().setSolverMethod(SolverOptions.HRZ_LUMPING)
    mypde.setSymmetryOn()  # turn symmetry on
    mypde.setValue(D=1.0)  # set the value of D in the general form to 1.

    ############################################FIRST TIME STEPS AND SOURCE
    # define small radius around point xc
    src_radius = 25.0
    print("src_radius = ", src_radius)
    # set initial values for first two time steps with source terms
    u = U0 * (cos(length(x - xc) * 3.1415 / src_radius) + 1) * whereNegative(length(x - xc) - src_radius)
    u_m1 = u
    # plot source shape
    cut_loc = []  # where the cross section of the source along x will be
    src_cut = []  # where the cross section of the source will be
    # create locations for source cross section
    for i in range(ndx // 2 - ndx // 10, ndx // 2 + ndx // 10):

# 需要导入模块: from esys.escript.linearPDEs import LinearPDE [as 别名]
# 或者: from esys.escript.linearPDEs.LinearPDE import getSolverOptions [as 别名]
class MultiScale(object):
   """
   problem description:
   -(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: old/current stress tensor
   :var Y: vector, body force
   :var y: vector, Neumann bc traction
   :var q: vector, Dirichlet bc mask
   :var r: vector, Dirichlet bc value
   """
   def __init__(self,domain,ng=1,useMPI=False,np=1,random=False,rtol=1.e-2,usePert=False,pert=-2.e-6,verbose=False):
      """
      initialization of the problem, i.e. model constructor
      :param domain: type Domain, domain of the problem
      :param ng: type integer, number of Gauss points
      :param useMPI: type boolean, use MPI or not
      :param np: type integer, number of processors
      :param random: type boolean, if or not use random density field
      :param rtol: type float, relevative tolerance for global convergence
      :param usePert: type boolean, if or not use perturbation method
      :param pert: type float, perturbated strain applied to DEM to obtain tangent operator
      :param verbose: type boolean, if or not print messages during calculation
      """
      self.__domain=domain
      self.__pde=LinearPDE(domain,numEquations=self.__domain.getDim(),numSolutions=self.__domain.getDim())
      self.__pde.getSolverOptions().setSolverMethod(SolverOptions.DIRECT)
      self.__pde.setSymmetryOn()
      #self.__pde.getSolverOptions().setTolerance(rtol**2)
      #self.__pde.getSolverOptions().setPackage(SolverOptions.UMFPACK)
      self.__numGaussPoints=ng
      self.__rtol=rtol
      self.__usepert=usePert
      self.__pert=pert
      self.__verbose=verbose
      self.__pool=get_pool(mpi=useMPI,threads=np)
      self.__scenes=self.__pool.map(initLoad,range(ng))
      self.__strain=escript.Tensor(0,escript.Function(self.__domain))
      self.__stress=escript.Tensor(0,escript.Function(self.__domain))
      self.__S=escript.Tensor4(0,escript.Function(self.__domain))
      
      if self.__usepert:
         s = self.__pool.map(getStressTensor,self.__scenes)
         t = self.__pool.map(getTangentOperator,zip(self.__scenes,repeat(pert)))
         for i in xrange(ng):
            self.__stress.setValueOfDataPoint(i,s[i])
            self.__S.setValueOfDataPoint(i,t[i])
      else:
         st = self.__pool.map(getStressAndTangent2D,self.__scenes)
         for i in xrange(ng):
            self.__stress.setValueOfDataPoint(i,st[i][0])
            self.__S.setValueOfDataPoint(i,st[i][1])
                     
   def initialize(self, b=escript.Data(), f=escript.Data(), specified_u_mask=escript.Data(), specified_u_val=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 specified_u_mask: type vector, mask of location for Dirichlet boundary
      :param specified_u_val: type vector, specified displacement for Dirichlet boundary
      """
      self.__pde.setValue(Y=b,y=f,q=specified_u_mask,r=specified_u_val)
      
   def getDomain(self):
      """
      return model domain
      """
      return self.__domain
      
   def getRelTolerance(self):
      """
      return relative tolerance for convergence
      type float
      """
      return self.__rtol
  
   def getCurrentPacking(self,pos=(),time=0,prefix=''):
      if len(pos) == 0:
         # output all Gauss points packings
         self.__pool.map(outputPack,zip(self.__scenes,repeat(time),repeat(prefix)))
      else:
         # output selected Gauss points packings
         scene = [self.__scenes[i] for i in pos]
         self.__pool.map(outputPack,zip(scene,repeat(time),repeat(prefix)))
   
   def getLocalVoidRatio(self):
      void=escript.Scalar(0,escript.Function(self.__domain))
      e = self.__pool.map(getVoidRatio2D,self.__scenes)
      for i in xrange(self.__numGaussPoints):
         void.setValueOfDataPoint(i,e[i])
      return void
   
   def getLocalAvgRotation(self):
      rot=escript.Scalar(0,escript.Function(self.__domain))
      r = self.__pool.map(avgRotation2D,self.__scenes)
      for i in xrange(self.__numGaussPoints):
         rot.setValueOfDataPoint(i,r[i])
#.........这里部分代码省略.........

# 需要导入模块: from esys.escript.linearPDEs import LinearPDE [as 别名]
# 或者: from esys.escript.linearPDEs.LinearPDE import getSolverOptions [as 别名]
#~ pl.colorbar()
#~ pl.savefig("tester2.png")
#~ 
#~ tester=fbottom*wherePositive(bottom)
#~ tester=tester.toListOfTuples()
#~ tester=np.reshape(tester,(ndx+1,ndy+1))
#~ pl.clf()
#~ pl.imshow(tester)
#~ pl.colorbar()
#~ pl.savefig("tester3.png")


# ... open new PDE ...
mypde=LinearPDE(domain)
print(mypde.isUsingLumping())
print(mypde.getSolverOptions())
#mypde.getSolverOptions().setSolverMethod(SolverOptions.LUMPING)
mypde.setSymmetryOn()
kmat = kronecker(domain)
mypde.setValue(D=kmat*rho)

# define small radius around point xc
# Lsup(x) returns the maximum value of the argument x
src_radius = 50#2*Lsup(domain.getSize())
print("src_radius = ",src_radius)

dunit=numpy.array([0.,1.]) # defines direction of point source
#~ dunit=(x-xc)
#~ absrc=length(dunit)
#~ dunit=dunit/maximum(absrc,1e-10)

# 需要导入模块: from esys.escript.linearPDEs import LinearPDE [as 别名]
# 或者: from esys.escript.linearPDEs.LinearPDE import getSolverOptions [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))))
#.........这里部分代码省略.........

# 需要导入模块: from esys.escript.linearPDEs import LinearPDE [as 别名]
# 或者: from esys.escript.linearPDEs.LinearPDE import getSolverOptions [as 别名]
   def test_Poisson4Classic(self):
        x=Solution(self.domain).getX()
        # --- set exact solution ----
        u_ex=Data(1.,(4,),Solution(self.domain))
        g_ex=Data(0.,(4,self.domain.getDim()), Solution(self.domain))
        A=Data(0.,(4,self.domain.getDim(),4,self.domain.getDim()), Function(self.domain))
        for i in range(self.domain.getDim()):
           u_ex[0]+= 1*(i+1) *x[i]
           g_ex[0,i]=1*(i+1)
           u_ex[1]+= 2*(i+1)*x[i]
           g_ex[1,i]=2*(i+1)
           u_ex[2]+= 3*(i+1)*x[i]
           g_ex[2,i]=3*(i+1)
           u_ex[3]+= 4*(i+1)*x[i]
           g_ex[3,i]=4*(i+1)
           A[0,i,0,i]=1
           A[1,i,1,i]=1
           A[2,i,2,i]=1
           A[3,i,3,i]=1

        # create PDE:
        pde=LinearPDE(self.domain,numEquations=4)
        pde.setSymmetryOn()
        mask=whereZero(x[0])*[1,1,1,1]
        pde.setValue(r=u_ex,q=mask)
        pde.setValue(A=A,y=matrixmult(g_ex,self.domain.getNormal()))
        # -------- get the solution ---------------------------
        pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
        pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
        pde.getSolverOptions().setPreconditioner(SolverOptions.AMG)
        pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
        if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
        if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
        if MAX_LEVEL!=None: pde.getSolverOptions().setLevelMax(MAX_LEVEL)
        pde.getSolverOptions().setAMGInterpolation(SolverOptions.CLASSIC_INTERPOLATION)

        u=pde.getSolution()
        # -------- test the solution ---------------------------
        error=Lsup(u-u_ex)/Lsup(u_ex)
        self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)

# 需要导入模块: from esys.escript.linearPDEs import LinearPDE [as 别名]
# 或者: from esys.escript.linearPDEs.LinearPDE import getSolverOptions [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:
#.........这里部分代码省略.........

# 需要导入模块: from esys.escript.linearPDEs import LinearPDE [as 别名]
# 或者: from esys.escript.linearPDEs.LinearPDE import getSolverOptions [as 别名]
class MultiScale(object):
   """
   problem description:
   D_ij u_j = -X_{ij,j} + Y_i
   Neumann boundary: d_ij u_j = n_j X_{ij} + y_i
   Dirichlet boundary: u_i = r_i where q_i > 0
   :var u: unknown vector, displacement
   :var X: old/current stress tensor
   :var Y: vector, body force
   :var y: vector, Neumann bc traction
   :var q: vector, Dirichlet bc mask
   :var r: vector, Dirichlet bc value
   """
   def __init__(self,domain,dim,ng=1,useMPI=False,np=1,rho=2.35e3,mIds=False,\
                FEDENodeMap=False,DE_ext=False,FEDEBoundMap=False,conf=False):
      """
      initialization of the problem, i.e. model constructor
      :param domain: type Domain, domain of the problem
      :param ng: type integer, number of Gauss points
      :param useMPI: type boolean, use MPI or not
      :param np: type integer, number of processors
      :param rho: type float, density of material
      :param mIds: a list contains membrane node IDs
      :param FEDENodeMap: a dictionary with FE and DE boundary node IDs in keys and values
      :param DE_ext: name of file which saves initial state of exterior DE domain
      :param FEDEBoundMap: a dictionary with FE and DE boundary element IDs in keys and values (deprecated)
      :param conf: type float, conf pressure on membrane
      """
      self.__domain=domain
      self.__pde=LinearPDE(self.__domain,numEquations=dim,numSolutions=dim)
      self.__pde.getSolverOptions().setSolverMethod(SolverOptions.HRZ_LUMPING)
      self.__pde.setSymmetryOn()
      self.__numGaussPoints=ng
      self.__rho=rho
      self.__mIds=mIds
      self.__FEDENodeMap=FEDENodeMap
      self.__FEDEBoundMap=FEDEBoundMap
      self.__conf=conf
      self.__pool=get_pool(mpi=useMPI,threads=np)
      self.__scenes=self.__pool.map(initLoad,range(ng))
      self.__strain=escript.Tensor(0,escript.Function(self.__domain))
      self.__stress=escript.Tensor(0,escript.Function(self.__domain))
      self.__u=escript.Vector(0,escript.Solution(self.__domain))
      self.__u_last=escript.Vector(0,escript.Solution(self.__domain))
      self.__u_t=escript.Vector(0,escript.Solution(self.__domain))
      self.__dt=0
      # ratio between timesteps in internal DE and FE domain
      self.__nsOfDE_int=1
      # if FEDENodeMap is given, employ exterior DE domain
      if self.__FEDENodeMap:
         self.__sceneExt=self.__pool.apply(initLoadExt,(DE_ext,))
         # ratio between timesteps in external DE and FE domain 
         self.__nsOfDE_ext=1
         # get interface nodal forces as boundary condition
         self.__FEf = self.__pool.apply(initNbc,(self.__sceneExt,self.__conf,mIds,FEDENodeMap))
         """
         # get interface nodal tractions as boundary condition (deprecated)
         self.__Nbc=escript.Vector(0,escript.Solution(self.__domain))
         FENbc = self.__pool.apply(initNbc,(self.__sceneExt,self.__conf,mIds,FEDENodeMap))
         for FEid in FENbc.keys():
            self.__Nbc.setValueOfDataPoint(FEid,FENbc[FEid])
         """
      # get stress tensor at material points
      s = self.__pool.map(getStress2D,self.__scenes)
      for i in xrange(ng):
         self.__stress.setValueOfDataPoint(i,s[i])
                     
   def initialize(self, b=escript.Data(), f=escript.Data(), specified_u_mask=escript.Data(), specified_u_val=escript.Data(), dt=0):
      """
      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 specified_u_mask: type vector, mask of location for Dirichlet boundary
      :param specified_u_val: type vector, specified displacement for Dirichlet boundary
      """
      self.__pde.setValue(Y=b,y=f,q=specified_u_mask,r=specified_u_val)
      # if FEDENodeMap is given
      if self.__FEDENodeMap:
			# assign dt_FE/dt_DE_ext to self.__nsOfDE_ext
         dt_ext = self.__pool.apply(getScenetDt,(self.__sceneExt,))
         self.__nsOfDE_ext = int(round(dt/dt_ext))
         print "Ratio between time step in FE and exterior DE domain: %1.1e"%self.__nsOfDE_ext
      dt_int = self.__pool.map(getScenetDt,self.__scenes)
      # assign a list of dt_FE/dt_DE_int to self.__nsOfDE_int
      self.__nsOfDE_int = numpy.round(numpy.array(dt)/dt_int).astype(int)
      print "Maximum ratio between time step in FE and interior DE domains: %1.1e"%max(self.__nsOfDE_int)
      if dt == 0:
         raise RuntimeError,"Time step in FE domain is not given"
      self.__dt = dt
      print "Time step in FE domain:%1.1e"%self.__dt

   def getDomain(self):
      """
      return model domain
      """
      return self.__domain

## below is Ning Guo's original script ##
    
   def getCurrentPacking(self,pos=(),time=0,prefix=''):
#.........这里部分代码省略.........