Merge branch 'develop' into make_examples_working

Conflicts:
    .gitignore
    src/porepy/fracs/meshing.py
    src/porepy/numerics/fv/mpfa.py
    src/porepy/numerics/fv/tpfa.py
    src/porepy/numerics/fv/transport/upwind_coupling.py
    src/porepy/numerics/mixed_dim/coupler.py
    test/integration/test_tpfaMultiDim.py
    test/unit/test_upwind_elimination.py

.gitignore

@@ -22,3 +22,6 @@
 *.o
 *.py_o
 *.so
+*.coverage
+.cache/*
+bin/*

examples/example8/aa_test_advection_diffusion_coupling.py

@@ -125,6 +125,7 @@ def add_data_advection_diffusion(gb, domain, tol):
 #------------------------------------------------------------------------------#
 
 
+do_save = False
 folder = os.path.dirname(os.path.realpath(__file__)) + "/"
 export_folder = folder + 'advection_diffusion_coupling'
 tol = 1e-3
@@ -164,7 +165,8 @@ def add_data_advection_diffusion(gb, domain, tol):
     d["p"] = darcy.discr.extract_p(g, d["up"])
     d["P0u"] = darcy.discr.project_u(g, discharge, d)
 
-exporter.export_vtk(gb, 'darcy', ["p", "P0u"], folder=export_folder)
+if do_save:
+    exporter.export_vtk(gb, 'darcy', ["p", "P0u"], folder=export_folder)
 
 #################################################################
 
@@ -184,6 +186,6 @@ def add_data_advection_diffusion(gb, domain, tol):
 theta = sps.linalg.spsolve(D + U, rhs_u + rhs_d)
 diffusion.split(gb, "temperature", theta)
 
-exporter.export_vtk(gb, 'advection_diffusion', [
-                    "temperature"], folder=export_folder)
-
+if do_save:
+    exporter.export_vtk(gb, 'advection_diffusion', [
+        "temperature"], folder=export_folder)

src/porepy/fracs/meshing.py

@@ -8,6 +8,7 @@
 """
 import numpy as np
 import scipy.sparse as sps
+import warnings
 
 from porepy.fracs import structured, simplex, split_grid
 from porepy.fracs.fractures import Intersection
@@ -23,6 +24,14 @@ def simplex_grid(fracs, domain, **kwargs):
     Main function for grid generation. Creates a fractured simiplex grid in 2
     or 3 dimensions.
 
+    NOTE: For some fracture networks, what appears to be a bug in Gmsh leads to
+    surface grids with cells that does not have a corresponding face in the 3d
+    grid. The problem may have been resolved (at least partly) by newer
+    versions of Gmsh, but can still be an issue for our purposes. If this
+    behavior is detected, an assertion error is raised. To avoid the issue,
+    and go on with a surface mesh that likely is problematic, kwargs should
+    contain a keyword ensure_matching_face_cell=False.
+
     Parameters
     ----------
     fracs (list of np.ndarray): One list item for each fracture. Each item
@@ -81,10 +90,11 @@ def simplex_grid(fracs, domain, **kwargs):
     tag_faces(grids)
 
     # Assemble grids in a bucket
-    gb = assemble_in_bucket(grids)
+    gb = assemble_in_bucket(grids, **kwargs)
     gb.compute_geometry()
     # Split the grids.
-    split_grid.split_fractures(gb)
+    split_grid.split_fractures(gb, **kwargs)
+    gb.assign_node_ordering()
     return gb
 
 #------------------------------------------------------------------------------#
@@ -211,6 +221,7 @@ def from_gmsh(file_name, dim, **kwargs):
 
 #------------------------------------------------------------------------------#
 
+
 def cart_grid(fracs, nx, **kwargs):
     """
     Creates a cartesian fractured GridBucket in 2- or 3-dimensions.
@@ -275,6 +286,7 @@ def cart_grid(fracs, nx, **kwargs):
 
     # Split grid.
     split_grid.split_fractures(gb, **kwargs)
+    gb.assign_node_ordering()
     return gb
 
 
@@ -356,7 +368,7 @@ def nodes_per_face(g):
     return n_per_face
 
 
-def assemble_in_bucket(grids):
+def assemble_in_bucket(grids, **kwargs):
     """
     Create a GridBucket from a list of grids.
     Parameters
@@ -392,7 +404,7 @@ def assemble_in_bucket(grids):
 
             for lg in grids[dim + 1]:
                 cell_2_face, cell = obtain_interdim_mappings(
-                    lg, fn, n_per_face)
+                    lg, fn, n_per_face, **kwargs)
                 face_cells = sps.csc_matrix(
                     (np.array([True] * cell.size), (cell, cell_2_face)),
                     (lg.num_cells, hg.num_faces))
@@ -404,10 +416,20 @@ def assemble_in_bucket(grids):
     return bucket
 
 
-def obtain_interdim_mappings(lg, fn, n_per_face):
+def obtain_interdim_mappings(lg, fn, n_per_face,
+                             ensure_matching_face_cell=True, **kwargs):
     """
     Find mappings between faces in higher dimension and cells in the lower
     dimension
+
+    Parameters:
+        lg: Lower dimensional grid.
+        fn: Higher dimensional face-node relation.
+        n_per_face: Number of nodes per face in the higher-dimensional grid.
+        ensure_matching_face_cell: Boolean, defaults to True. If True, an
+            assertion is made that all lower-dimensional cells corresponds to a
+            higher dimensional cell.
+
     """
     if lg.dim > 0:
         cn_loc = lg.cell_nodes().indices.reshape((n_per_face,
@@ -426,6 +448,14 @@ def obtain_interdim_mappings(lg, fn, n_per_face):
     # An element in cell_2_face gives, for all cells in the
     # lower-dimensional grid, the index of the corresponding face
     # in the higher-dimensional structure.
-
+    if not (np.all(is_mem) or np.all(~is_mem)):
+        if ensure_matching_face_cell:
+            raise ValueError(
+            '''Either all cells should have a corresponding face in a higher
+            dim grid or no cells should have a corresponding face in a higher
+            dim grid. This likely is related to gmsh behavior. ''')
+        else:
+            warnings.warn('''Found inconsistency between cells and higher
+                          dimensional faces. Continuing, faces crossed''')
     low_dim_cell = np.where(is_mem)[0]
     return cell_2_face, low_dim_cell

src/porepy/numerics/compressible/problems.py

@@ -2,88 +2,110 @@
 import scipy.sparse as sps
 
 from porepy.grids import structured
-from porepy.numerics.fv import tpfa
-from porepy.numerics.compressible import solvers
+from porepy.numerics.parabolic import *
+from porepy.numerics.fv import tpfa, mass_matrix, fvutils
+from porepy.numerics.mixed_dim.coupler import Coupler
 from porepy.params.data import Parameters
 from porepy.params import tensor
 from porepy.params import bc
 from porepy.viz.exporter import export_vtk, export_pvd
 
 
-class SlightlyCompressible():
-    """
-    Base-class for slightly compressible flow. Initialize all needed
-    attributes for a slightly compressible solver.
-    """
+class SlightlyCompressible(ParabolicProblem):
+    '''
+    Inherits from ParabolicProblem
+    This class solves equations of the type:
+    phi *c_p dp/dt  - \nabla K \nabla p = q
 
-    def __init__(self):
-        self.solver = solvers.Implicit(self)
-        self.solver.parameters['store_results'] = True
-        self.parameters = {'file_name': 'pressure'}
-        self.parameters['folder_name'] = 'results'
-        self.data = dict()
-    #---------Discretization---------------
+    Init:
+    - gb (Grid/GridBucket) Grid or grid bucket for the problem
+    - physics (string) Physics key word. See Parameters class for valid physics
 
-    def flux_disc(self):
-        """
-        Returns the flux discretization. 
-        """
-        return tpfa.Tpfa()
+    functions:
+    discharge(): computes the discharges and saves it in the grid bucket as 'p'
+    Also see functions from ParabolicProblem
+
+    Example:
+    # We create a problem with standard data
+
+    gb = meshing.cart_grid([], [10,10], physdims=[1,1])
+    for g, d in gb:
+        d['problem'] = SlightlyCompressibleData(g, d)
+    problem = SlightlyCompressible(gb)
+    problem.solve()
+   '''
+
+    def __init__(self, gb, physics='flow'):
+        ParabolicProblem.__init__(self, gb, physics)
 
-    def solve(self):
+    def space_disc(self):
+        return self.diffusive_disc(), self.source_disc()
+
+    def time_disc(self):
         """
-        Call the solver
+        Returns the time discretization.
         """
-        self.data = self.solver.solve()
-        return self.data
-
-    #-----Parameters------------
-    def porosity(self):
-        return np.ones(self.grid().num_cells)
+        class TimeDisc(mass_matrix.MassMatrix):
+            def __init__(self, deltaT):
+                self.deltaT = deltaT
+
+            def matrix_rhs(self, g, data):
+                lhs, rhs = mass_matrix.MassMatrix.matrix_rhs(self, g, data)
+                return lhs * data['compressibility'], rhs * data['compressibility']
+        single_dim_discr = TimeDisc(self.time_step())
+        multi_dim_discr = Coupler(single_dim_discr)
+        return multi_dim_discr
+
+    def discharge(self):
+        self.diffusive_disc().split(self.grid(), 'p', self._solver.p)
+        fvutils.compute_discharges(self.grid())
+
+
+class SlightlyCompressibleData(ParabolicData):
+    '''
+    Inherits from ParabolicData
+    Base class for assigning valid data for a slighly compressible problem.
+    Init:
+    - g    (Grid) Grid that data should correspond to
+    - d    (dictionary) data dictionary that data will be assigned to
+    - physics (string) Physics key word. See Parameters class for valid physics
+
+    Functions:
+        compressibility: (float) the compressibility of the fluid
+        permeability: (tensor.SecondOrder) The permeability tensor for the rock.
+                      Setting the permeability is equivalent to setting
+                      the ParabolicData.diffusivity() function. 
+    Example:
+    # We set an inflow and outflow boundary condition by overloading the
+    # bc_val term
+    class ExampleData(SlightlyCompressibleData):
+        def __init__(g, d):
+            SlightlyCompressibleData.__init__(self, g, d)
+        def bc_val(self):
+            left = self.grid().nodes[0] < 1e-6
+            right = self.grid().nodes[0] > 1 - 1e-6
+            val = np.zeros(g.num_faces)
+            val[left] = 1
+            val[right] = -1
+            return val
+    gb = meshing.cart_grid([], [10,10], physdims=[1,1])
+    for g, d in gb:
+        d['problem'] = ExampleData(g, d)
+    '''
+
+    def __init__(self, g, data, physics='flow'):
+        ParabolicData.__init__(self, g, data, physics)
+
+    def _set_data(self):
+        ParabolicData._set_data(self)
+        self.data()['compressibility'] = self.compressibility()
 
     def compressibility(self):
-        return 1
+        return 1.0
 
     def permeability(self):
         kxx = np.ones(self.grid().num_cells)
         return tensor.SecondOrder(self.grid().dim, kxx)
 
-    #--------Inn/outflow terms---------------
-
-    def initial_pressure(self):
-        return np.zeros(self.grid().num_cells)
-
-    def source(self, t):
-        return np.zeros(self.g.num_cells)
-
-    def bc(self):
-        dir_bound = np.array([])
-        return bc.BoundaryCondition(self.grid(), dir_bound,
-                                    ['dir'] * dir_bound.size)
-
-    def bc_val(self, t):
-        return np.zeros(self.grid().num_faces)
-
-    #---------Overloaded Functions-----------------
-    def grid(self):
-        raise NotImplementedError('subclass must overload function grid()')
-
-    #--------Time stepping------------
-    def time_step(self):
-        return 1.0
-
-    def end_time(self):
-        return 1.0
-
-    def save_results(self):
-        pressures = self.data['pressure']
-        times = np.array(self.data['times'])
-        folder = self.parameters['folder_name']
-        f_name = self.parameters['file_name']
-        for i, p in enumerate(pressures):
-            data_to_plot = {'pressure': p}
-            export_vtk(
-                self.grid(), f_name, data_to_plot, time_step=i, folder=folder)
-
-        export_pvd(
-            self.grid(), self.parameters['file_name'], times, folder=folder)
+    def diffusivity(self):
+        return self.permeability()