Merge remote-tracking branch 'origin/primal_contact' into primal_contact

examples/papers/Readme.md

@@ -4,17 +4,18 @@ This folder contains paper examples of PorePy. The intention is that the results
 We try to keep the examples updated as the code changes, but may not always succeed due to lack of time, or simply because the code has changed so much that compatability is no longer prioritized. The papers should at least work with the commit that merged them into the main branch.
 
 # Overview
-* [arXiv_1610_02905](./arXiv_1610_02905/) paper "*Dual virtual element method for discrete fractures networks*" by Alessio Fumagalli, Eirik Keilegavlen. See [arXiv pre-print](https://arxiv.org/abs/1610.02905) or [journal article](https://epubs.siam.org/doi/abs/10.1137/16M1098231).
-* [arXiv_1709_03519](./arXiv_1709_03519/) paper "*Dual virtual element method in presence of an inclusion*" by Alessio Fumagalli. See [arXiv pre-print](https://arxiv.org/abs/1709.03519) or [journal article](https://www.sciencedirect.com/science/article/pii/S0893965918301812).
-* [arXiv_1711_01818](./arXiv_1711_01818/) paper "*Dual Virtual Element Methods for Discrete Fracture Matrix Models*" by Alessio Fumagalli, Eirik Keilegavlen. See [arXiv pre-print](https://arxiv.org/abs/1711.01818) or [journal article](https://ogst.ifpenergiesnouvelles.fr/articles/ogst/full_html/2019/01/ogst170210/ogst170210.html).
-* [arXiv_1712_00460](./arXiv_1712_00460) paper "*PorePy: An Open-Source Simulation Tool for Flow and Transport in Deformable Fractured Rocks*" by Eirik Keilegavlen, Alessio Fumagalli, Runar Berge, Ivar Stefansson, Inga Berre. See [arXiv pre-print](https://arxiv.org/abs/1712.00460).
-* [arXiv_1712_07392](./arXiv_1712_07392) paper "*Implementation of mixed-dimensional models for flow in fractured porous media*" by Eirik Keilegavlen, Alessio Fumagalli, Runar Berge, Ivar Stefansson. See [arXiv pre-print](https://arxiv.org/abs/1712.07392) or [journal article](https://www.springerprofessional.de/en/implementation-of-mixed-dimensional-models-for-flow-in-fractured/16377424).
-* [arXiv_1712_08479](./arXiv_1712_08479) paper "*Hybrid-Dimensional Finite Volume Discretizations for Fractured Porous Media*" by Ivar Stefansson, Inga Berre, Eirik Keilegavlen.
-* [arXiv_1803_01732](./arXiv_1803_01732) paper "*Conforming, non-conforming and non-matching discretization couplings in discrete fracture network simulations*" by Alessio Fumagalli, Eirik Keilegavlen, and Stefano Scialo'. See [arXiv pre-print](https://arxiv.org/abs/1803.01732) or [journal article](https://www.sciencedirect.com/science/article/pii/S0021999118306508).
-* [arXiv_1802_05961](./arXiv_1802_05961) paper "*Unified approach to discretization of flow in fractured porous media*" by Jan M. Nordbotten, Wietse M. Boon, Alessio Fumagalli, and  Eirik Keilegavlen. See [arXiv pre-print](https://arxiv.org/abs/1802.05961) or [journal article](https://link.springer.com/article/10.1007/s10596-018-9778-9).
-* [arXiv_1809_06926](./arXiv_1809_06926) paper "*Call for participation: Verification benchmarks for single-phase flow in three-dimensional fractured porous media*" by Inga Berre, Wietse Boon, Bernd Flemisch, Alessio Fumagalli, Dennis Glaser, Eirik Keilegavlen, Anna Scotti, Ivar Stefansson, Alexandru Tatomir. See [arXiv pre-print](https://arxiv.org/abs/1809.06926).
-* [arXiv_1810_12761](./arXiv_1810_12761) paper "*A multiscale flux basis for mortar mixed discretizations of reduced Darcy-Forchheimer fracture models*" by Elyes Ahmed, Alessio Fumagalli, Ana Budiša. See [arXiv pre-print](https://arxiv.org/abs/1810.12761) or [journal article](https://www.sciencedirect.com/science/article/pii/S0045782519303044).
-* [arXiv_1903_01117](./arXiv_1903_01117) paper "*A multi-layer reduced model for flow in porous media with a fault and surrounding damage zones*" by Alessio Fumagalli, Anna Scotti. See [arXiv pre-print](https://arxiv.org/abs/1903.01117).
+* [arXiv:1610.02905](./arXiv_1610_02905/) paper "*Dual virtual element method for discrete fractures networks*" by Alessio Fumagalli, Eirik Keilegavlen. See [arXiv pre-print](https://arxiv.org/abs/1610.02905) or [journal article](https://epubs.siam.org/doi/abs/10.1137/16M1098231).
+* [arXiv:1709.03519](./arXiv_1709_03519/) paper "*Dual virtual element method in presence of an inclusion*" by Alessio Fumagalli. See [arXiv pre-print](https://arxiv.org/abs/1709.03519) or [journal article](https://www.sciencedirect.com/science/article/pii/S0893965918301812).
+* [arXiv:1711.01818](./arXiv_1711_01818/) paper "*Dual Virtual Element Methods for Discrete Fracture Matrix Models*" by Alessio Fumagalli, Eirik Keilegavlen. See [arXiv pre-print](https://arxiv.org/abs/1711.01818) or [journal article](https://ogst.ifpenergiesnouvelles.fr/articles/ogst/full_html/2019/01/ogst170210/ogst170210.html).
+* [arXiv:1712.00460](./arXiv_1712_00460) paper "*PorePy: An Open-Source Simulation Tool for Flow and Transport in Deformable Fractured Rocks*" by Eirik Keilegavlen, Alessio Fumagalli, Runar Berge, Ivar Stefansson, Inga Berre. See [arXiv pre-print](https://arxiv.org/abs/1712.00460).
+* [arXiv:1712.07392](./arXiv_1712_07392) paper "*Implementation of mixed-dimensional models for flow in fractured porous media*" by Eirik Keilegavlen, Alessio Fumagalli, Runar Berge, Ivar Stefansson. See [arXiv pre-print](https://arxiv.org/abs/1712.07392) or [journal article](https://www.springerprofessional.de/en/implementation-of-mixed-dimensional-models-for-flow-in-fractured/16377424).
+* [arXiv:1712.08479](./arXiv_1712_08479) paper "*Hybrid-Dimensional Finite Volume Discretizations for Fractured Porous Media*" by Ivar Stefansson, Inga Berre, Eirik Keilegavlen.
+* [arXiv:1803.01732](./arXiv_1803_01732) paper "*Conforming, non-conforming and non-matching discretization couplings in discrete fracture network simulations*" by Alessio Fumagalli, Eirik Keilegavlen, and Stefano Scialo'. See [arXiv pre-print](https://arxiv.org/abs/1803.01732) or [journal article](https://www.sciencedirect.com/science/article/pii/S0021999118306508).
+* [arXiv:1802.05961](./arXiv_1802_05961) paper "*Unified approach to discretization of flow in fractured porous media*" by Jan M. Nordbotten, Wietse M. Boon, Alessio Fumagalli, and  Eirik Keilegavlen. See [arXiv pre-print](https://arxiv.org/abs/1802.05961) or [journal article](https://link.springer.com/article/10.1007/s10596-018-9778-9).
+* [arXiv:1809.06926](./arXiv_1809_06926) paper "*Call for participation: Verification benchmarks for single-phase flow in three-dimensional fractured porous media*" by Inga Berre, Wietse Boon, Bernd Flemisch, Alessio Fumagalli, Dennis Glaser, Eirik Keilegavlen, Anna Scotti, Ivar Stefansson, Alexandru Tatomir. See [arXiv pre-print](https://arxiv.org/abs/1809.06926).
+* [arXiv:1810.12761](./arXiv_1810_12761) paper "*A multiscale flux basis for mortar mixed discretizations of reduced Darcy-Forchheimer fracture models*" by Elyes Ahmed, Alessio Fumagalli, Ana Budiša. See [arXiv pre-print](https://arxiv.org/abs/1810.12761) or [journal article](https://www.sciencedirect.com/science/article/pii/S0045782519303044).
+* [arXiv:1903.01117](./arXiv_1903_01117) paper "*A multi-layer reduced model for flow in porous media with a fault and surrounding damage zones*" by Alessio Fumagalli, Anna Scotti. See [arXiv pre-print](https://arxiv.org/abs/1903.01117).
 
 # Separate repository
-* [repo](https://github.com/alessiofumagalli/multiscale_timedependent) paper "*Robust linear domain decomposition schemes for reduced non-linear fracture flow models*" by Elyes Ahmed, Alessio Fumagalli, Ana Budiša, Eirik Keilegavlen, Jan M. Nordbotten, A. Radu Florin. See [arXiv pre-print](https://arxiv.org/abs/1906.05831).
+* [arXiv:1906.05831 [math.NA]](https://github.com/alessiofumagalli/multiscale_timedependent) paper "*Robust linear domain decomposition schemes for reduced non-linear fracture flow models*" by Elyes Ahmed, Alessio Fumagalli, Ana Budiša, Eirik Keilegavlen, Jan M. Nordbotten, A. Radu Florin. See [arXiv pre-print](https://arxiv.org/abs/1906.05831).
+* [arXiv:1907.12514 [math.NA]](https://github.com/alessiofumagalli/dfn_transport) paper "*Analysis of conforming, non-matching, and polygonal methods for Darcy and advection-diffusion-reaction simulations in discrete fracture networks*" by Andrea Borio, Alessio Fumagalli, Stefano Scialò. See [arXiv pre-print](https://arxiv.org/abs/1907.12514).

src/porepy/models/contact_mechanics_biot_model.py

@@ -350,7 +350,7 @@ def assign_variables(self):
     def discretize_biot(self, gb):
         """
         To save computational time, the full Biot equation (without contact mechanics)
-        is discretized once. This computing the same terms multiple times.
+        is discretized once. This is to avoid computing the same terms multiple times.
         """
         g = gb.grids_of_dimension(gb.dim_max())[0]
         d = gb.node_props(g)
@@ -459,4 +459,4 @@ def run_biot(setup, atol=1e-10):
         assembler.distribute_variable(sol)
         setup.export_step()
         errors.append(newton_errors)
-    setup.export_pvd()
+    setup.export_pvd()

src/porepy/numerics/contact_mechanics/contact_conditions.py

@@ -233,8 +233,11 @@ def discretize(self, g_h, g_l, data_h, data_l, data_edge):
             else:  # should never happen
                 raise AssertionError("Should not get here")
 
-            # Scale equations (helps iterative solver)
-            # TODO: Find out what happens here
+            # Depending on the state of the system, the weights in the tangential direction may
+            # become huge or tiny compared to the other equations. This will
+            # impede convergence of an iterative solver for the linearized
+            # system. As a partial remedy, rescale the condition to become
+            # closer to unity.
             w_diag = np.diag(loc_displacement_weight) + np.diag(loc_traction_weight)
             W_inv = np.diag(1 / w_diag)
             loc_displacement_weight = W_inv.dot(loc_displacement_weight)
@@ -403,9 +406,6 @@ def set_projections(gb):
     It is assumed that the surface is planar.
 
     """
-
-    ambient_dim = gb.dim_max()
-
     # Information on the vector normal to the surface is not available directly
     # from the surface grid (it could be constructed from the surface geometry,
     # which spans the tangential plane). We instead get the normal vector from
@@ -415,7 +415,7 @@ def set_projections(gb):
 
         mg = d_m["mortar_grid"]
         # Only consider edges where the lower-dimensional neighbor is of co-dimension 1
-        if not mg.dim == ambient_dim - 1:
+        if not mg.dim == (gb.dim_max() - 1):
             continue
 
         # Neigboring grids
@@ -436,15 +436,16 @@ def set_projections(gb):
         # Ensure all normal vectors on the relevant surface points outwards
         unit_normal[:, faces_on_surface] *= sgn
 
-        # Now we need to pick out a normal vector of the higher dimensional grid
+        # Now we need to pick out *one*  normal vector of the higher dimensional grid
         # which coincides with this mortar grid. This could probably have been
         # done with face tags, but we instead project the normal vectors onto the
-        # mortar grid to kill off all irrelevant grids.
+        # mortar grid to kill off all irrelevant faces. Restriction to a single
+        # normal vector is done in the construction of the projection object
+        # (below).
+        # NOTE: Use a single normal vector to span the tangential and normal space,
+        # thus assuming the surface is planar.
         outwards_unit_vector_mortar = mg.master_to_mortar_int().dot(unit_normal.T).T
 
-        # Use a single normal vector to span the tangential and normal space,
-        # assuming the surface is planar.
-
         # NOTE: The normal vector is based on the first cell in the mortar grid,
         # and will be pointing from that cell towards the other side of the
         # mortar grid. This defines the positive direction in the normal direction.

src/porepy/numerics/interface_laws/contact_mechanics_interface_laws.py

@@ -3,9 +3,6 @@
 """
 Implementation of contact conditions for fracture mechanics, using a primal formulation.
 
-
-The primal formulation is conceptually similar, but mathematically different from,
-the dual formulation, currently located in elliptic_interface_laws.
 """
 
 import numpy as np
@@ -164,7 +161,8 @@ def discretize(self, g_h, g_l, data_h, data_l, data_edge):
             # The final equations should relate to continuity of the normal froces
             matrix_dictionary_edge[self.SURFACE_DISCRETIZATION_KEY] = A
 
-        # Discretization of the contact mechanics is done externally.
+        # Discretization of the contact mechanics is done by a ColumbContact
+        # object. 
         # The resulting equations are located at the lower-dimensional grid,
         # however, the discretization is inherently linked to the mortar grid.
         # It is therefore constructed here.
@@ -306,7 +304,9 @@ def assemble_matrix_rhs(
 
         # This is first a stress balance: stress from the higher dimensional
         # domain (both interior and bound_stress) should match with the contact stress:
-        # -\lambda_slave + \lambda_mortar = 0.
+        #
+        #     traction_slave + traction_master = 0
+        #
         # Optionally, a diffusion term can be added in the tangential direction
         # of the stresses, this is currently under implementation.
 
@@ -323,8 +323,8 @@ def assemble_matrix_rhs(
         ## First, we obtain \lambda_mortar = stress * u_master + bound_stress * u_mortar
         # Stress contribution from the higher dimensional domain, projected onto
         # the mortar grid
-        # Switch the direction of the vectors, so that for all faces, a positive
-        # force points into the fracture surface.
+        # Switch the direction of the vectors to obtain the traction as defined
+        # by the outwards pointing normal vector.
         stress_from_master = (
             mg.master_to_mortar_int(nd=ambient_dimension)
             * sign_switcher
@@ -396,7 +396,7 @@ class MatrixScalarToForceBalance:
     This class adds the matrix scalar (pressure) contribution to the force balance posed
     on the mortar grid by PrimalContactCoupling.
 
-    We account for the grad P contribution to the forces on the higher-dimensional
+    We account for the pressure contribution to the forces on the higher-dimensional
     internal boundary, i.e. the last term of:
 
         boundary_traction_hat = stress * u_hat + bound_stress * u_mortar + gradP * p_hat

src/porepy/numerics/mixed_dim/assembler.py

@@ -1044,7 +1044,7 @@ def merge_variable(self, var):
             else:
                 data = self.gb.node_props(g)
             if var_name == var:
-                loc_value = data[var_name]
+                loc_value = data[pp.STATE][var_name]
             else:
                 loc_value = 0
             values[dof[bi] : dof[bi + 1]] = loc_value

src/porepy/numerics/vem/dual_elliptic.py

@@ -33,7 +33,7 @@ def project_flux(gb, discr, flux, P0_flux, mortar_key="mortar_solution"):
 
         # we need to recover the flux from the mortar variable before
         # the projection, only lower dimensional edges need to be considered.
-        edge_flux = np.zeros(d[flux].size)
+        edge_flux = np.zeros(d[pp.STATE][flux].size)
         faces = g.tags["fracture_faces"]
         if np.any(faces):
             # recover the sign of the flux, since the mortar is assumed
@@ -47,10 +47,10 @@ def project_flux(gb, discr, flux, P0_flux, mortar_key="mortar_solution"):
                     continue
                 # project the mortar variable back to the higher dimensional
                 # problem
-                # edge_flux += sign * g_m.mortar_to_master_int() * d_e[mortar_key]
-                edge_flux += sign * g_m.master_to_mortar_avg().T * d_e[mortar_key]
+                # edge_flux += sign * g_m.mortar_to_master_int() * d_e[pp.STATE][mortar_key]
+                edge_flux += sign * g_m.master_to_mortar_avg().T * d_e[pp.STATE][mortar_key]
 
-        d[P0_flux] = discr.project_flux(g, edge_flux + d[flux], d)
+        d[pp.STATE][P0_flux] = discr.project_flux(g, edge_flux + d[pp.STATE][flux], d)
 
 
 # ------------------------------------------------------------------------------#