Python linearPDEs.LinearPDE类代码示例

    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

      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 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]

      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 test_PDE2D(self):
        dx_tests=0.1
        sigma0=1.
        electrodes=[(0.5-2*dx_tests,1.), (0.5-dx_tests,1.), (0.5+dx_tests,1.), (0.5+2*dx_tests,1.)]
        domain=finRectangle(20,20, d1=mpisize,  diracPoints=electrodes, diracTags=["sl0", "sl1", "sr0", "sr1"] )
        loc=Locator(domain,electrodes[2:])

        # this creates some reference Data:
        x=domain.getX()
        q=whereZero(x[0]-inf(x[0]))+whereZero(x[0]-sup(x[0]))+whereZero(x[1]-inf(x[1]))
        ppde=LinearPDE(domain, numEquations=1)
        s=Scalar(0.,DiracDeltaFunctions(domain))
        s.setTaggedValue("sl0" ,1.)
        s.setTaggedValue("sl1",-1.)
        ppde.setValue(A=kronecker(2)*sigma0, q=q, y_dirac=s)
        pp=ppde.getSolution()
        uu=loc(pp)

        # arguments for DcRes
        current = 10.
        sourceInfo = [ "sl0",  "sl1" ]
        sampleTags = [ ("sr0", "sr1") ]

        sigmaPrimary=7.
        phiPrimary=pp*current*sigma0/sigmaPrimary

        uuscale=1-current*sigma0/sigmaPrimary
        delphi_in = [ (uu[1]-uu[0]) * uuscale]

        acw=DcRes(domain, loc, delphi_in, sampleTags,  phiPrimary, sigmaPrimary)

        self.assertLess(Lsup(phiPrimary-acw.getPrimaryPotential()), 1.e-10 * Lsup(acw.getPrimaryPotential()))

        SIGMA=10. # matches current
        args0=acw.getArguments(SIGMA)
        p=args0[0]
        u=args0[1]

        # true secondary potential
        pps=pp-phiPrimary
        self.assertLess(Lsup(p-pps), 1.e-6*Lsup(pps))

        # test return values at electrodes:
        self.assertLess(abs(u[0]-uu[0]*uuscale), 1.e-6 * abs(uu[0]*uuscale))
        self.assertLess(abs(u[1]-uu[1]*uuscale), 1.e-6 * abs(uu[1]*uuscale))

        # this sould be zero
        dd=acw.getDefect(SIGMA, *args0)
        self.assertTrue( dd >= 0.)
        self.assertTrue( dd <= 1e-7 )

 def getMaxEigenvalue(self):
    """
    return the max eigenvalue of the model
    type: float
    """
    dom = self.getDomain()
    dim = dom.getDim()
    pdeKu_P = LinearPDE(dom,numEquations=dim,numSolutions=dim)
    T = self.getCurrentTangent()
    pdeKu_P.setValue(A=T)
    pdeKu_P.getOperator().saveMM('tempMTX.mtx')
    mtx = mmread('tempMTX.mtx')
    maxEigVal = max(eigs(mtx,k=1,return_eigenvectors=False,which='LR'))
    return maxEigVal.real

 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 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)

 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)

 def __init__(self,domain,ng=1,useMPI=False,np=1,random=False,rtol=1.e-2,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, relevant tolerance for global convergence
    :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.__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))
    
    st = self.__pool.map(getStressAndTangent,self.__scenes)
    for i in xrange(ng):
       self.__stress.setValueOfDataPoint(i,st[i][0])
       self.__S.setValueOfDataPoint(i,st[i][1])

    def __init__(self, domain, phi_f=None, phi=None, L_g=None, 
                       perm_f_0=None, perm_f_1=None, perm_f_2=None,
                       k_w =None, k_g=None, mu_w =None, mu_g =None,
                       rho_w =None, rho_g=None, sigma=0, A_mg=0, f_rg=1.,   
                         wells=[], g=9.81*U.m/U.sec**2):
       """
       set up
       
       :param domain: domain
       :type domain: `esys.escript.Domain`
       :param phi_f: porosity of the fractured rock as function of the gas pressure
       :type phi_f: `MaterialPropertyWithDifferential`
       :param phi: total porosity if None phi=1.
       :type phi: scalar or None
       :param L_g: gas adsorption as function of gas pressure
       :type L_g: `MaterialPropertyWithDifferential`
       :param S_fg: gas saturation in the fractured rock as function of the capillary pressure p_fc=S_fg- p_f
       :type S_fg: `MaterialPropertyWithDifferential`
       :param perm_f_0: permeability in the x_0 direction in the fractured media
       :type perm_f_0: scalar 
       :param perm_f_1: permeability in the x_1 direction in the fractured media
       :type perm_f_1: scalar
       :param perm_f_2: permeability in the x_2 direction in the fractured media
       :type perm_f_2: scalar
       :param k_w: relative permeability of water as function of water saturation
       :type k_w: `MaterialProperty`
       :param k_g: relative permeability of gas as function of gas saturation
       :type k_g: `MaterialProperty`
       :param mu_w: viscosity of water as function of water pressure
       :type mu_w: `MaterialProperty`
       :param mu_g: viscosity of gas as function of gas pressure
       :type mu_g: `MaterialProperty`
       :param rho_w: density of water as function of water pressure
       :type rho_w: `MaterialPropertyWithDifferential`
       :param rho_g: density of gas as function of gas pressure
       :type rho_g: `MaterialPropertyWithDifferential`
       :param wells : list of wells
       :type wells: list of `Well`
       :param sigma: shape factor for gas matrix diffusion 
       :param A_mg: diffusion constant for gas adsorption
       :param f_rg: gas re-adsorption factor
       """
 
       DualPorosity.__init__(self, domain,
                            phi_f=phi_f, phi_m=None, phi=phi,
                            S_fg=None, S_mg=None, 
                            perm_f_0=perm_f_0, perm_f_1=perm_f_1, perm_f_2=perm_f_2,
                            perm_m_0=None , perm_m_1=None, perm_m_2=None, 
                            k_w =k_w, k_g=k_g, mu_w =mu_w, mu_g =mu_g,
                            rho_w =rho_w, rho_g=rho_g, 
                            wells=wells, g=g)
       self.L_g=L_g
       self.sigma = sigma 
       self.A_mg = A_mg
       self.f_rg  = f_rg
       self.__pde=LinearPDE(self.domain, numEquations=3, numSolutions =3)

 def __update(self,tag,tag_value,value):
     if self.__pde==None:
        self.__pde=LinearPDE(self.domain,numSolutions=1)
     mask=Scalar(0.,Function(self.domain))
     mask.setTaggedValue(tag,1.)
     self.__pde.setValue(Y=mask)
     mask=wherePositive(abs(self.__pde.getRightHandSide()))
     value*=(1.-mask)
     value+=tag_value*mask
     return 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 __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())
    try:
          self.__pde.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.__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,list(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,list(zip(self.__scenes,repeat(pert))))
       for i in range(ng):
          self.__stress.setValueOfDataPoint(i,s[i])
          self.__S.setValueOfDataPoint(i,t[i])
    else:
       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])

    def gzpot(p, y, x, *args):
        #rho, rhox, rhoy, R = p
        rhox=args[0]/2.; rhoy=args[1]/2.
        rho, R, z =p
        #Domain related.
        mx = args[0]; my = args[1];
        #PDE related
        G=6.67300*10E-11

        #DOMAIN CONSTRUCTION
        domain=ReadGmsh('data/example10m/example10m.msh',2)
        domx=Solution(domain).getX()
        mask=wherePositive(R-length(domx-rholoc))
        rhoe=rho*mask
        kro=kronecker(domain)

        q=whereZero(domx[1]-my)+whereZero(domx[1])+whereZero(domx[0])+whereZero(domx[0]-mx)
        #ESCRIPT PDE CONSTRUCTION
        mypde=LinearPDE(domain)
        mypde.setValue(A=kro,Y=4.*np.pi*G*rhoe,q=q,r=0.0)
        mypde.setSymmetryOn()
        sol=mypde.getSolution()

        g_field=grad(sol) #The graviational accelleration g.
        g_fieldz=g_field*[0,1] #The vertical component of the g field.
        gz=length(g_fieldz) #The magnitude of the vertical component.

        #MODEL SIZE SAMPLING
        sol_escgz=[]
        sol_escx=[]
        for i in range(0,len(x)):
            sol_escgz.append([x[i],rhoy+z])

        sample=[] # array to hold values
        rec=Locator(gz.getFunctionSpace(),sol_escgz) #location to record
        psol=rec.getValue(gz)

        err = np.sum((np.array(y) - np.array(psol))**2.)
        print("Lsup= ",Lsup(np.array(psol)-np.array(sol_angz))/Lsup(np.array(psol)))
        return err