Python LinearPDE.createCoefficient方法代码示例

# 需要导入模块: from esys.escript.linearPDEs import LinearPDE [as 别名]
# 或者: from esys.escript.linearPDEs.LinearPDE import createCoefficient [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 createCoefficient [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 createCoefficient [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 createCoefficient [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 createCoefficient [as 别名]

#.........这里部分代码省略.........
         E_fps=  -  phi_f * rho_fw 
         E_fsp= S_fg * ( rho_fg * dphi_fdp + phi_f * drho_fgdp )
         E_fss= phi_f * rho_fg 
        
        
         
         F_fw=0.
         F_fg=0.
         D_fw=0.
         D_fg=0.

         prod_index=Scalar(0., DiracDeltaFunctions(self.domain))
         geo_fac=Scalar(0., DiracDeltaFunctions(self.domain))
         for I in self.wells:
             prod_index.setTaggedValue(I.name, I.getProductivityIndex() )
             geo_fac.setTaggedValue(I.name, I.geo_fac)

         F_fp_Y = A_fw * prod_index * BHP  * geo_fac
         F_fs_Y = A_fg * prod_index * BHP * geo_fac
         D_fpp =  A_fw * prod_index * geo_fac
         D_fsp =  A_fg * prod_index * geo_fac

         
         if self.domain.getDim() == 3:
            F_fp_X = ( A_fw * self.g * rho_fw * self.perm_f[2] ) * kronecker(self.domain)[2]
            F_fs_X = ( A_fg * self.g * rho_fg * self.perm_f[2] ) * kronecker(self.domain)[2]
         else:
            F_fp_X = 0 * kronecker(self.domain)[1]
            F_fs_X = 0 * kronecker(self.domain)[1]
            
         F_mg_Y = H * L_g


         D=self.__pde.createCoefficient("D")
         A=self.__pde.createCoefficient("A")
         Y=self.__pde.createCoefficient("Y")
         X=self.__pde.createCoefficient("X") 
         
         d_dirac=self.__pde.createCoefficient("d_dirac")
         y_dirac=self.__pde.createCoefficient("y_dirac")


        
         D[0,0]=E_fpp
         D[0,1]=E_fps
         d_dirac[0,0]=dt * D_fpp
         
         D[1,0]=E_fsp
         D[1,1]=E_fss 
         D[1,2]= rho_g_surf
         d_dirac[1,0]=dt * D_fsp
         
         D[0,2]= -dt * H * dL_gdp * 0
         D[2,2]= 1 + dt * H
         
         
         for i in range(self.domain.getDim()):
            A[0,i,0,i] = dt * A_fw * self.perm_f[i]
            A[1,i,1,i] = dt * A_fg * self.perm_f[i]

         X[0,:]=  dt * F_fp_X
         X[1,:]=  dt * F_fs_X

         Y[0] = E_fpp *  p_f_old + E_fps * S_fg_old 
         Y[1] = E_fsp *  p_f_old + E_fss * S_fg_old + rho_g_surf * c_mg_old 
         Y[2] = c_mg_old                                                    + dt * ( F_mg_Y -  H * dL_gdp * p_f *0 )