update example 4

examples/papers/dfn_comparison/example_4_data.py

@@ -0,0 +1,71 @@
+import numpy as np
+
+from porepy.params import tensor
+from porepy.params.bc import BoundaryCondition
+from porepy.params.data import Parameters
+
+#------------------------------------------------------------------------------#
+
+def add_data(gb, domain, direction, tol):
+    """
+    Define the permeability, apertures, boundary conditions
+    """
+    gb.add_node_props(['param'])
+
+    # Aavatsmark_transmissibilities only for tpfa intra-dimensional coupling
+
+    for g, d in gb:
+        d['Aavatsmark_transmissibilities'] = True
+
+        if g.dim < 2:
+            continue
+
+        param = Parameters(g)
+
+        # Permeability
+        kxx = np.ones(g.num_cells)
+        param.set_tensor("flow", tensor.SecondOrder(3, kxx))
+
+        # Source term
+        param.set_source("flow", np.zeros(g.num_cells))
+
+        # Boundaries
+        bound_faces = g.get_boundary_faces()
+        if bound_faces.size != 0:
+            bound_face_centers = g.face_centers[:, bound_faces]
+
+            if direction == 'left_right':
+                left = bound_face_centers[0, :] < domain['xmin'] + tol
+                right = bound_face_centers[0, :] > domain['xmax'] - tol
+                bc_dir = np.logical_or(left, right)
+                bc_one = right
+            elif direction == 'bottom_top':
+                bottom = bound_face_centers[2, :] < domain['zmin'] + tol
+                top = bound_face_centers[2, :] > domain['zmax'] - tol
+                bc_dir = np.logical_or(top, bottom)
+                bc_one = top
+            elif direction == 'back_front':
+                back = bound_face_centers[1, :] < domain['ymin'] + tol
+                front = bound_face_centers[1, :] > domain['ymax'] - tol
+                bc_dir = np.logical_or(back, front)
+                bc_one = front
+
+            labels = np.array(['neu'] * bound_faces.size)
+            labels[bc_dir] = 'dir'
+
+            bc_val = np.zeros(g.num_faces)
+            bc_val[bound_faces[bc_one]] = 1
+
+            param.set_bc("flow", BoundaryCondition(g, bound_faces, labels))
+            param.set_bc_val("flow", bc_val)
+        else:
+            param.set_bc("flow", BoundaryCondition(
+                g, np.empty(0), np.empty(0)))
+
+        d['param'] = param
+
+    gb.add_edge_prop('Aavatsmark_transmissibilities')
+    for _, d in gb.edges_props():
+        d['Aavatsmark_transmissibilities'] = True
+
+#------------------------------------------------------------------------------#

examples/papers/dfn_comparison/example_4_vem.py

@@ -1,17 +1,66 @@
-from porepy.viz.exporter import Exporter
 import scipy.sparse as sps
+import pickle
+
+from porepy.viz.exporter import Exporter
+
+from porepy.grids.grid import FaceTag
+#from porepy.grids import coarsening as co
+
+from porepy.numerics.vem import vem_dual, vem_source
 
-import example_4_create_grid
+#import example_4_create_grid
+import example_4_data
 
 #------------------------------------------------------------------------------#
 
-folder_export = 'viz/'
-file_export = 'vem'
+def main(conforming, direction):
+    file_export = 'solution'
+    tol = 1e-4
+
+    if conforming:
+        gb = pickle.load(open('dfn_conf.grid' ,'rb'))
+        folder_export = './example_4_vem_conf_'+direction+"/"
+    else:
+        gb = pickle.load(open('dfn_non-conf.grid' ,'rb'))
+        folder_export = './example_4_vem_non-conf_'+direction+"/"
+
+    domain = {'xmin': -800, 'xmax': 600,
+              'ymin': 100, 'ymax': 1500,
+              'zmin': -100, 'zmax': 1000}
+
+    internal_flag = FaceTag.FRACTURE
+    [g.remove_face_tag_if_tag(FaceTag.BOUNDARY, internal_flag) for g, _ in gb]
+
+    example_4_data.add_data(gb, domain, direction, tol)
 
-# consider gmsh 2.11
-conforming = True
-gb = example_4_create_grid.create(conforming)
+    # Choose and define the solvers and coupler
+    solver_flow = vem_dual.DualVEMDFN(gb.dim_max(), 'flow')
+    A_flow, b_flow = solver_flow.matrix_rhs(gb)
 
-save = Exporter(gb, file_export, folder_export)
-save.write_vtk()
+    solver_source = vem_source.IntegralDFN(gb.dim_max(), 'flow')
+    A_source, b_source = solver_source.matrix_rhs(gb)
 
+    up = sps.linalg.spsolve(A_flow+A_source, b_flow+b_source)
+    solver_flow.split(gb, "up", up)
+
+    gb.add_node_props(["discharge", "p", "P0u"])
+    solver_flow.extract_u(gb, "up", "discharge")
+    solver_flow.extract_p(gb, "up", "p")
+    solver_flow.project_u(gb, "discharge", "P0u")
+
+    save = Exporter(gb, file_export, folder_export)
+    save.write_vtk(["p", "P0u"])
+
+#------------------------------------------------------------------------------#
+
+#conforming = True
+#main(conforming, 'left_right')
+#main(conforming, 'bottom_top')
+#main(conforming, 'back_front')
+
+conforming = False
+main(conforming, 'left_right')
+main(conforming, 'bottom_top')
+main(conforming, 'back_front')
+
+#------------------------------------------------------------------------------#