Wrote some comments for mpsa

src/porepy/numerics/fv/mpfa.py

@@ -162,7 +162,9 @@ def discretize(self, g, data):
         a = param.aperture
         robin_weight = param.get_robin_weight(self)
 
-        trm, bound_flux, bp_cell, bp_face = mpfa(g, k, bnd, apertures=a, robin_weight=robin_weight)
+        trm, bound_flux, bp_cell, bp_face = mpfa(
+            g, k, bnd, apertures=a, robin_weight=robin_weight
+        )
         data["flux"] = trm
         data["bound_flux"] = bound_flux
         data["bound_pressure_cell"] = bp_cell

src/porepy/numerics/fv/mpsa.py

@@ -988,7 +988,7 @@ def mpsa_elasticity(
     ncsym_rob = subcell_topology.pair_over_subfaces_nd(ncsym_all + ncasym)
     del ncasym
     ncsym_rob = bound_exclusion.keep_robin_nd(ncsym_rob)
-    
+
     if robin_weight is None:
         if ncsym_rob.shape[0] != 0:
             raise ValueError(
@@ -1004,7 +1004,15 @@ def mpsa_elasticity(
         g, subcell_topology, eta, num_sub_cells, bound_exclusion, robin_weight
     )
 
-    # Pair the forces from each side
+    # Pair the forces from each side. For the Neumann faces this does in fact
+    # lead to an inconsistency because the asymetric parts does not cancell.
+    # We can add the asymmetric term for the Neumann faces, as is done above
+    # for the Robin faces, however, in some cases this might lead to sigular
+    # local matrices (when a subcell only have Neumann faces). We need to
+    # come back and adress this issue.
+
+    # ncsym * G is in fact (due to pair_over_subfaces)
+    # ncsym_L * G_L + ncsym_R * G_R for the left and right faces.
     ncsym = subcell_topology.pair_over_subfaces_nd(ncsym_all)
     del ncsym_all
     ncsym = bound_exclusion.exclude_robin_dirichlet_nd(ncsym)

test/unit/test_parameters.py

@@ -190,11 +190,11 @@ def test_porosity_assertion_vector(self):
     #####
 
     def _validate_robin_weight(self, p, physics, val=1):
-        if physics=='flow':
+        if physics == "flow":
             self.assertTrue(np.allclose(p.robin_weight_flow, val))
-        elif physics=='transport':
+        elif physics == "transport":
             self.assertTrue(np.allclose(p.robin_weight_transport, val))
-        elif physics=='mechanics':
+        elif physics == "mechanics":
             self.assertTrue(np.allclose(p.robin_weight_mechanics, val))
         else:
             self.assertTrue(False)
@@ -206,6 +206,7 @@ def test_robin_weight_default(self, val=1):
             self._validate_robin_weight(p, name)
 
         # Set by scalar
+
     def test_robin_weight_set(self):
         p = Parameters(self.g)
         for name in p.known_physics: