Python LinearPDE.getSolution方法代码示例

# 需要导入模块: from esys.escript.linearPDEs import LinearPDE [as 别名]
# 或者: from esys.escript.linearPDEs.LinearPDE import getSolution [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 getSolution [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 getSolution [as 别名]
   def test_Poisson2Classic(self):
        x=Solution(self.domain).getX()
        # --- set exact solution ----
        u_ex=Data(1.,(2,),Solution(self.domain))
        g_ex=Data(0.,(2,self.domain.getDim()), Solution(self.domain))
        A=Data(0.,(2,self.domain.getDim(),2,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)
           A[0,i,0,i]=1
           A[1,i,1,i]=1

        # create PDE:
        pde=LinearPDE(self.domain,numEquations=2)
        pde.setSymmetryOn()
        mask=whereZero(x[0])*[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 getSolution [as 别名]
   def test_PoissonWithDirectInterpolation(self):
        x=Solution(self.domain).getX()
        # --- set exact solution ----
        u_ex=Scalar(1,Solution(self.domain))
        g_ex=Vector(0.,Solution(self.domain))
        for i in range(self.domain.getDim()):
           u_ex+=(i+1)*x[i]
           g_ex[i]=(i+1)

        # create PDE:
        pde=LinearPDE(self.domain,numEquations=1)
        pde.setSymmetryOn()
        mask=whereZero(x[0])
        pde.setValue(r=u_ex,q=mask)
        pde.setValue(A=kronecker(self.domain),y=inner(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().setNumPreSweeps(3)
        pde.getSolverOptions().setNumPostSweeps(3)
        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.DIRECT_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 getSolution [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 getSolution [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 getSolution [as 别名]
    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 )

# 需要导入模块: from esys.escript.linearPDEs import LinearPDE [as 别名]
# 或者: from esys.escript.linearPDEs.LinearPDE import getSolution [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 getSolution [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 getSolution [as 别名]
class SlippingFault(SaddlePointProblem):
      """
      simple example of saddle point problem
      """
      def __init__(self,domain):
         super(SlippingFault, self).__init__(self)
         self.domain=domain
         self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
         self.__pde_u.setSymmetryOn()

      def initialize(self,density=1.,lmbd=1., mu=1., traction=Data(),fixed_u_mask=Data(), slip=0.):
         d=self.domain.getDim()
         self.slip=slip
         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] += mu
             A[i,j,i,j] += mu
             A[i,i,j,j] += lmbd
         self.__pde_u.setValue(A=A,q=fixed_u_mask,Y=-kronecker(Function(self.domain))[d-1]*g*density,y=traction)

      def inner(self,p0,p1):
         return integrate(inner(p0,p1),FunctionOnContactZero(self.domain))

      def solve_f(self,u,p,tol=1.e-8):
         self.__pde_u.setTolerance(tol)
         self.__pde_u.setValue(y_contact=-p)
         # print "p:",inf(p),sup(p)
         # print "u:",inf(u),sup(u)
         self.__pde_u.setValue(y_contact=-p)
         return  self.__pde_u.getSolution()

      def solve_g(self,u,tol=1.e-8):
         dp=Vector(0.,FunctionOnContactZero(self.domain))
         h=FunctionOnContactZero(self.domain).getSize()
         # print jump(u)-self.slip
         dp[0]=(self.slip[0]-jump(u[0]))*lam_mu/h
         dp[1]=(self.slip[1]-jump(u[1]))*lam_mu/h
         dp[2]=(self.slip[2]-jump(u[2]))*lam_mu/h
         return  dp

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

# 需要导入模块: from esys.escript.linearPDEs import LinearPDE [as 别名]
# 或者: from esys.escript.linearPDEs.LinearPDE import getSolution [as 别名]

#.........这里部分代码省略.........
      """
      return current strain
      type: Tensor on FunctionSpace
      """
      return self.__strain

   def getCurrentFlux(self):
      """
      return current Darcy flux
      type: Vector on FunctionSpace
      """
      n=self.getEquivalentPorosity()
      perm=self.__permeability*n**3/(1.-n)**2
      flux=-perm*util.grad(self.__pore)
      return flux
   
   def exitSimulation(self):
      """finish the whole simulation, exit"""
      self.__pool.close()
   
   def solveSolid(self, p_iter_gauss=escript.Data(), iter_max=50):
      """
      solve the pde for displacement using Newton-Ralphson scheme
      """
      k=util.kronecker(self.__domain)
      p=p_iter_gauss*k
      iterate=0
      rtol=self.__rtol
      stress_safe=self.__stress
      s_safe=self.__S
      x_safe=self.__domain.getX()
      self.__upde.setValue(A=s_safe, X=-stress_safe+p, r=self.__r)
      #residual0=util.L2(self.__pde.getRightHandSide()) # using force error
      u=self.__upde.getSolution()  # trial solution, displacement
      D=util.grad(u)               # trial strain tensor
      # !!!!!! obtain stress and tangent operator from DEM part
      update_stress,update_s,update_scenes=self.applyStrain_getStressTangentDEM(st=D)
      err=util.Lsup(u) # initial error before iteration
      converged=(err<1.e-12)
      while (not converged) and (iterate<iter_max):
         #if iterate>iter_max:
         #   raise RuntimeError,"Convergence for Newton-Raphson failed after %s steps."%(iter_max)
         iterate+=1
         self.__domain.setX(x_safe+u)
         self.__upde.setValue(A=update_s,X=-update_stress+p,r=escript.Data())
         #residual=util.L2(self.__pde.getRightHandSide())
         du=self.__upde.getSolution()
         u+=du
         l,d=util.L2(u),util.L2(du)
         err=d/l # displacement error, alternatively using force error 'residual'
         converged=(err<rtol)
         if err>rtol**3: # only update DEM parts when error is large enough
            self.__domain.setX(x_safe)
            D=util.grad(u)
            update_stress,update_s,update_scenes=self.applyStrain_getStressTangentDEM(st=D)
      """reset domain geometry to original until global convergence"""
      self.__domain.setX(x_safe)
      return u,D,update_stress,update_s,update_scenes
   
   def solve(self, globalIter=10, solidIter=50):
      """
      solve the coupled PDE using fixed-stress split method,
      call solveSolid to get displacement
      """
      rtol=self.__rtol
      x_safe=self.__domain.getX()

# 需要导入模块: from esys.escript.linearPDEs import LinearPDE [as 别名]
# 或者: from esys.escript.linearPDEs.LinearPDE import getSolution [as 别名]
    ###########################SAVING THE VALUE AT A LOC FOR EACH TIME STEP
    u_rec0 = []  # array to hold values
    rec = Locator(mydomain, [250.0, 250.0])  # location to record
    u_rec = rec.getValue(u)
    u_rec0.append(u_rec)  # get the first two time steps

    ####################################################ITERATION VARIABLES
    n = 0  # iteration counter
    t = 0  # time counter
    ##############################################################ITERATION
    while t < tend:
        g = grad(u)
        pres = csq * g  # get current pressure
        mypde.setValue(X=-pres)  # set values in pde
        accel = mypde.getSolution()  # get new acceleration
        u_p1 = (2.0 * u - u_m1) + h * h * accel  # calculate the displacement for the next time step
        u_m1 = u
        u = u_p1  # shift values back one time step for next iteration
        # save current displacement, acceleration and pressure
        if t >= rtime:
            saveVTK(
                os.path.join(savepath, "ex07b.%i.vtu" % n),
                displacement=length(u),
                acceleration=length(accel),
                tensor=pres,
            )
            rtime = rtime + rtime_inc  # increment data save time
        u_rec0.append(rec.getValue(u))  # location specific recording
        # increment loop values
        t = t + h

# 需要导入模块: from esys.escript.linearPDEs import LinearPDE [as 别名]
# 或者: from esys.escript.linearPDEs.LinearPDE import getSolution [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 getSolution [as 别名]
    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):
        cut_loc.append(xstep*i)
        src_cut.append([xstep*i,xc[1]])
    # locate the nearest nodes to the points in src_cut
    src=Locator(mydomain,src_cut)
    src_cut=src.getValue(u) #retrieve the values from the nodes
    # plot the x locations vs value and save the figure
    pl.plot(cut_loc,src_cut)
    pl.axis([xc[0]-src_radius*3,xc[0]+src_radius*3,0.,2.*U0])
    pl.savefig(os.path.join(savepath,"source_line.png"))

    ####################################################ITERATION VARIABLES
    n=0 # iteration counter
    t=0 # time counter
    ##############################################################ITERATION
    while t<tend:
        g=grad(u); pres=csq*h*h*g # get current pressure
        mypde.setValue(X=-pres,Y=(2.*u-u_m1)) # set values in pde
        u_p1 = mypde.getSolution() # get the new displacement
        u_m1=u; u=u_p1 # shift values back one time step for next iteration
            # save current displacement, acceleration and pressure
        if (t >= rtime):
            saveVTK(os.path.join(savepath,"ex07a.%i.vtu"%n),displacement=length(u),tensor=pres)
            rtime=rtime+rtime_inc #increment data save time
        # increment loop values
        t=t+h; n=n+1
        print("time step %d, t=%s"%(n,t))