Merge branch 'develop' into frac_mesh_3d

.gitignore

@@ -11,3 +11,6 @@
 *.vrml
 *~
 *.swo
+*.coverage
+.cache/*
+bin/*

examples/example8/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
@@ -166,7 +167,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)
 
 #################################################################
 
@@ -192,6 +194,6 @@ def add_data_advection_diffusion(gb, domain, tol):
 theta = sps.linalg.spsolve(D + U, rhs_u + rhs_d)
 diffusion_coupler.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

@@ -85,11 +85,13 @@ def simplex_grid(fracs=None, domain=None, network=None, **kwargs):
     gb = assemble_in_bucket(grids)
     gb.compute_geometry()
     # Split the grids.
-    split_grid.split_fractures(gb)
+    split_grid.split_fractures(gb, **kwargs)
+    gb.assign_node_ordering()
     return gb
 
 #------------------------------------------------------------------------------#
 
+
 def from_gmsh(file_name, dim, **kwargs):
     """
     Import an already generated grid from gmsh.
@@ -140,6 +142,7 @@ def from_gmsh(file_name, dim, **kwargs):
 
 #------------------------------------------------------------------------------#
 
+
 def cart_grid(fracs, nx, **kwargs):
     """
     Creates a cartesian fractured GridBucket in 2- or 3-dimensions.
@@ -204,6 +207,7 @@ def cart_grid(fracs, nx, **kwargs):
 
     # Split grid.
     split_grid.split_fractures(gb, **kwargs)
+    gb.assign_node_ordering()
     return gb
 
 
@@ -342,6 +346,9 @@ 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.
-
+    assert np.all(is_mem) or np.all(~is_mem),\
+        '''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
+        might be a problem with gmsh or how we read the gmsh output, not sure.. '''
     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()