用于断裂和可变形多孔介质的Python模拟工具___下载.zip

""" Module contains two sets of tests: One on MVEM itself, one on MVEM with gravity. """ import unittest import numpy as np import scipy.sparse as sps import porepy as pp class BasicsTest(unittest.TestCase): def _matrix(self, g, perm, bc): solver = pp.MVEM(keyword="flow") data = pp.initialize_default_data( g, {}, "flow", {"second_order_tensor": perm, "bc": bc} ) solver.discretize(g, data) return solver.assemble_matrix(g, data).todense() # ------------------------------------------------------------------------------# def test_dual_vem_1d_iso(self): g = pp.CartGrid(3, 1) g.compute_geometry() kxx = np.ones(g.num_cells) perm = pp.SecondOrderTensor(kxx) bf = g.tags["domain_boundary_faces"].nonzero()[0] bc = pp.BoundaryCondition(g, bf, bf.size * ["dir"]) M = self._matrix(g, perm, bc) M_known = 1e-2 * np.array( [ [25, -(8 + 1 / 3.0), 0, 0, 1e2, 0, 0], [-(8 + 1 / 3.0), 50, -(8 + 1 / 3.0), 0, -1e2, 1e2, 0], [0, -(8 + 1 / 3.0), 50, -(8 + 1 / 3.0), 0, -1e2, 1e2], [0, 0, -(8 + 1 / 3.0), 25, 0, 0, -1e2], [1e2, -1e2, 0, 0, 0, 0, 0], [0, 1e2, -1e2, 0, 0, 0, 0], [0, 0, 1e2, -1e2, 0, 0, 0], ] ) rtol = 1e-15 atol = rtol self.assertTrue(np.allclose(M, M.T, rtol, atol)) self.assertTrue(np.allclose(M, M_known, rtol, atol)) # ------------------------------------------------------------------------------# def test_dual_vem_1d_ani(self): g = pp.CartGrid(3, 1) g.compute_geometry() kxx = np.sin(g.cell_centers[0, :]) + 1 perm = pp.SecondOrderTensor( kxx, kyy=np.ones(g.num_cells), kzz=np.ones(g.num_cells) ) bf = g.tags["domain_boundary_faces"].nonzero()[0] bc = pp.BoundaryCondition(g, bf, bf.size * ["dir"]) M = self._matrix(g, perm, bc) M_known = 1e-2 * np.array( [ [21.4427334468001192, -7.14757781560004, 0, 0, 1e2, 0, 0], [ -7.14757781560004, 38.3411844655562319, -5.6328170062520355, 0, -1e2, 1e2, 0, ], [ 0, -5.6328170062520355, 31.2648069176318977, -4.788785299625264, 0, -1e2, 1e2, ], [0, 0, -4.7887852996252649, 14.3663558988757991, 0, 0, -1e2], [1e2, -1e2, 0, 0, 0, 0, 0], [0, 1e2, -1e2, 0, 0, 0, 0], [0, 0, 1e2, -1e2, 0, 0, 0], ] ) rtol = 1e-15 atol = rtol self.assertTrue(np.allclose(M, M.T, rtol, atol)) self.assertTrue(np.allclose(M, M_known, rtol, atol)) # ------------------------------------------------------------------------------# def test_dual_vem_2d_iso_cart(self): g = pp.CartGrid([2, 1], [1, 1]) g.compute_geometry() kxx = np.ones(g.num_cells) perm = pp.SecondOrderTensor(kxx=kxx, kyy=kxx, kzz=kxx) bf = g.tags["domain_boundary_faces"].nonzero()[0] bc = pp.BoundaryCondition(g, bf, bf.size * ["dir"]) M = self._matrix(g, perm, bc) # Matrix computed with an already validated code (MRST) M_known = np.array( [ [0.625, -0.375, 0, 0, 0, 0, 0, 1, 0], [-0.375, 1.25, -0.375, 0, 0, 0, 0, -1, 1], [0, -0.375, 0.625, 0, 0, 0, 0, 0, -1], [0, 0, 0, 1, 0, 0, 0, 1, 0], [0, 0, 0, 0, 1, 0, 0, 0, 1], [0, 0, 0, 0, 0, 1, 0, -1, 0], [0, 0, 0, 0, 0, 0, 1, 0, -1], [1, -1, 0, 1, 0, -1, 0, 0, 0], [0, 1, -1, 0, 1, 0, -1, 0, 0], ] ) rtol = 1e-15 atol = rtol self.assertTrue(np.allclose(M, M.T, rtol, atol)) self.assertTrue(np.allclose(M, M_known, rtol, atol)) # ------------------------------------------------------------------------------# def test_dual_vem_2d_ani_cart(self): g = pp.CartGrid([2, 1], [1, 1]) g.compute_geometry() kxx = np.square(g.cell_centers[1, :]) + 1 kyy = np.square(g.cell_centers[0, :]) + 1 kxy = -np.multiply(g.cell_centers[0, :], g.cell_centers[1, :]) perm = pp.SecondOrderTensor(kxx=kxx, kyy=kyy, kxy=kxy) bf = g.tags["domain_boundary_faces"].nonzero()[0] bc = pp.BoundaryCondition(g, bf, bf.size * ["dir"]) M = self._matrix(g, perm, bc) # Matrix computed with an already validated code (MRST) M_known = np.array( [ [ 0.625000000000000, -0.422619047619048, 0, 0.023809523809524, 0, 0.023809523809524, 0, 1, 0, ], [ -0.422619047619048, 1.267241379310345, -0.426724137931035, 0.023809523809524, 0.051724137931034, 0.023809523809524, 0.051724137931034, -1, 1, ], [ 0, -0.426724137931035, 0.642241379310345, 0, 0.051724137931034, 0, 0.051724137931034, 0, -1, ], [ 0.023809523809524, 0.023809523809524, 0, 1, 0, -0.047619047619048, 0, 1, 0, ], [ 0, 0.051724137931034, 0.051724137931034, 0, 0.879310344827586, 0, -0.189655172413793, 0, 1, ], [ 0.023809523809524, 0.023809523809524, 0, -0.047619047619048, 0, 1, 0, -1, 0, ], [ 0, 0.051724137931034, 0.051724137931034, 0, -0.189655172413793, 0, 0.879310344827586, 0, -1, ], [1, -1, 0, 1, 0, -1, 0, 0, 0], [0, 1, -1, 0, 1, 0, -1, 0, 0], ] ) rtol = 1e-15 atol = rtol self.assertTrue(np.allclose(M, M.T, rtol, atol)) self.assertTrue(np.allclose(M, M_known, rtol, atol)) # ------------------------------------------------------------------------------# def test_dual_vem_2d_iso_simplex(self): g = pp.StructuredTriangleGrid([1, 1], [1, 1]) g.compute_geometry() kxx = np.ones(g.num_cells) perm = pp.SecondOrderTensor(kxx=kxx, kyy=kxx) bf = g.tags["domain_boundary_faces"].nonzero()[0] bc = pp.BoundaryCondition(g, bf, bf.size * ["dir"]) M = self._matrix(g, perm, bc) # Matrix computed with an already validated code (MRST) faces = np.arange(5) M_known = np.array