From 5a45ad0a1ea71b0a6a2f86e9cb6debbb893629d9 Mon Sep 17 00:00:00 2001 From: Runar Date: Thu, 6 Sep 2018 09:17:54 +0200 Subject: [PATCH 01/26] Added robin-discretization for MPSA with added dofs --- src/porepy/__init__.py | 2 +- src/porepy/numerics/fv/mpsa.py | 175 +++++++++++++++++++++++++++++++++ src/porepy/params/data.py | 28 ++++++ 3 files changed, 204 insertions(+), 1 deletion(-) diff --git a/src/porepy/__init__.py b/src/porepy/__init__.py index 09f0a54599..ecf99a92e7 100644 --- a/src/porepy/__init__.py +++ b/src/porepy/__init__.py @@ -27,7 +27,7 @@ # Numerics # Control volume, elliptic -from porepy.numerics.fv.mpsa import Mpsa, FracturedMpsa +from porepy.numerics.fv.mpsa import Mpsa, FracturedMpsa, RobinMpsa from porepy.numerics.fv.tpfa import Tpfa, TpfaMixedDim from porepy.numerics.fv.mpfa import Mpfa, MpfaMixedDim from porepy.numerics.fv.biot import Biot diff --git a/src/porepy/numerics/fv/mpsa.py b/src/porepy/numerics/fv/mpsa.py index 43d7e91f12..f9d11eaed3 100644 --- a/src/porepy/numerics/fv/mpsa.py +++ b/src/porepy/numerics/fv/mpsa.py @@ -138,6 +138,181 @@ def rhs(self, g, bound_stress, bc_val, f): # ------------------------------------------------------------------------------# +class RobinMpsa(Mpsa): + """ + Subclass of MPSA for discretizing a domain with robin boundary conditinos. + Adds DOFs on each robin face. + """ + + def __init__(self, robin_faces, **kwargs): + """ + robin_faces is a function that takes a grid and return an array + of length grid.num_faces that are true on the robin faces and false otherwise + """ + Mpsa.__init__(self, **kwargs) + assert hasattr(self, "physics"), "Mpsa must assign physics" + self.robin_faces = robin_faces + + def ndof(self, g): + """ + Return the number of degrees of freedom associated to the method. + In this case number of cells pluss the number of robin face times + dimension (stress dof). + + Parameter + --------- + g: grid, or a subclass. + + Return + ------ + dof: the number of degrees of freedom. + + """ + return g.dim * (g.num_cells + np.sum(self.robin_faces(g))) + + def matrix_rhs(self, g, data, discretize=True): + """ + Return the matrix and right-hand side for a discretization of a second + order elliptic equation using a FV method with a multi-point stress + approximation with dofs added on the robin faces + + Parameters + ---------- + g : grid, or a subclass, with geometry fields computed. + data: dictionary to store the data. For details on necessary keywords, + see method discretize() + discretize (boolean, optional): default True. Whether to discetize + prior to matrix assembly. If False, data should already contain + discretization. + + Return + ------ + matrix: sparse csr (g.dim * (g_num_cells + {#of robin faces}, + g.dim * (g_num_cells + {#of robin faces})) + Discretization matrix. + rhs: array (g.dim * g_num_cells + {#of robin faces}) + Right-hand side which contains the boundary conditions and the scalar + source term. + """ + if discretize: + self.discretize_robin(g, data) + + b_e = data["b_e"] + A_e = data["A_e"] + + bc_val = data["param"].get_bc_val(self) + + rhs = b_e * bc_val + + return A_e, rhs + + def discretize_robin(self, g, data, faces=None, **kwargs): + """ + Discretize the vector elliptic equation by the multi-point stress and added + degrees of freedom on the robin faces + + The method computes fluxes over faces in terms of displacements in + adjacent cells (defined as the two cells sharing the face). + + The name of data in the input dictionary (data) are: + param : Parameter(Class). Contains the following parameters: + tensor : fourth_order_tensor + Permeability defined cell-wise. If not given a identity permeability + is assumed and a warning arised. + bc : boundary conditions (optional) + bc_val : dictionary (optional) + Values of the boundary conditions. The dictionary has at most the + following keys: 'dir' and 'neu', for Dirichlet and Neumann boundary + conditions, respectively. + apertures : (np.ndarray) (optional) apertures of the cells for scaling of + the face normals. + + Parameters + ---------- + g : grid, or a subclass, with geometry fields computed. + data: dictionary to store the data. + """ + + # dir_bound = g.get_all_boundary_faces() + # bound = bc.BoundaryCondition(g, dir_bound, ['dir'] * dir_bound.size) + + rob_faces = self.robin_faces(g) + + bound = data["param"].get_bc(self) + + if bound.bc_type == "scalar": + bound.is_dir[rob_faces] = True + bound.is_neu[rob_faces] = False + elif bound.bc_type == "vectorial": + bound.is_dir[:, rob_faces] = True + bound.is_neu[:, rob_faces] = False + else: + raise ValueError("Unknow boundary condition type: " + bound.bc_type) + if np.sum(bound.is_dir * bound.is_neu) !=0: + raise AssertionError('Found faces that are both dirichlet and neuman') + + # Discretize with normal mpsa + self.discretize(g, data, **kwargs) + stress, bound_stress = data["stress"], data["bound_stress"] + + # Create A and rhs + div = fvutils.vector_divergence(g) + a = div * stress + b = div * bound_stress + + # we find the matrix indices of the robin faces + rob_ind = mcolon.mcolon( + g.dim * rob_faces, g.dim * rob_faces + g.dim + ) + # We find the sign of the robin faces. + sgn_rob = _sign_matrix(g, rob_faces) + + # The displacement on the robin boundary face are considered unknowns, + # so we move them over to the lhs. The rhs now only consists of the + # external boundary faces + b_rob = b[:, rob_ind] + b_external = b.copy() + sparse_mat.zero_columns(b_external, rob_ind) + + bound_stress_external = bound_stress.copy().tocsc() + sparse_mat.zero_columns(bound_stress_external, rob_ind) + # We calculate the stress on the robin faces due to the internal cells + # and robin face displacements (the other boundaries are on the rhs) + roben_stress = sps.hstack( + ( + sgn_rob * stress[rob_ind, :], + (sgn_rob * bound_stress[rob_ind, :])[:, rob_ind] + ) + ) + # we add the displacement on the robin boundaries + alpha = data["param"].get_robin_factor()[rob_faces] + alpha = np.ravel([alpha]*g.dim, 'F') + roben_disp = sps.hstack( + ( + sps.csr_matrix((rob_ind.size, g.dim * g.num_cells)), + sps.diags(alpha, 0), + ) + ) + roben_condition = roben_stress + roben_disp + A = sps.vstack((sps.hstack((a, b_rob)), roben_condition), format="csr") + + # negative sign since we have moved b_external from lhs to rhs + d_b = -b_external + # sps.csr_matrix((int_b_left.size, g.num_faces * g.dim)) + d_t = ( + sps.csr_matrix( + (np.ones(rob_ind.size), (np.arange(rob_ind.size), rob_ind)), + (rob_ind.size, g.num_faces * g.dim), + ) + -sgn_rob * bound_stress_external[rob_ind] + ) + + b_matrix = sps.vstack((d_b, d_t), format="csr") + + data["b_e"] = b_matrix + data["A_e"] = A + + class FracturedMpsa(Mpsa): """ Subclass of MPSA for discretizing a fractured domain. Adds DOFs on each diff --git a/src/porepy/params/data.py b/src/porepy/params/data.py index 4b85d6a9f4..5505268867 100644 --- a/src/porepy/params/data.py +++ b/src/porepy/params/data.py @@ -765,3 +765,31 @@ def get_bc_val_mechanics(self): return np.zeros(self._num_faces * self.dim) bc_val_mechanics = property(get_bc_val_mechanics) + + + # ------------------------------------ + def get_robin_factor(self, default=1): + """ double or array-like + face-wise representation of robin_factor. Set as either a np.ndarary, or a + scalar (uniform) value. Always returned as np.ndarray. + """ + if not hasattr(self, "_robin_factor"): + return default * np.ones(self._num_faces * self.dim) + + if isinstance(self._robin_factor, np.ndarray): + # Hope that the user did not initialize as array with wrong size + return self._robin_factor + else: + return self._robin_factor * np.ones(self._num_faces) + + def set_robin_factor(self, val): + """ Set physics-specific robin factor. The robin factor is given by + T + robin_factor * u = g + where T is the traction, and g is given as a boundary value + + Parameters: + val : robin_factor, representing the robin factor + """ + self._robin_factor = val + + robin_factor = property(get_robin_factor, set_robin_factor) From 220b686ccc9786379c4211a3e47135a7232b9eb1 Mon Sep 17 00:00:00 2001 From: Runar Date: Thu, 6 Sep 2018 15:24:57 +0200 Subject: [PATCH 02/26] Added Robin boundary condition for mpfa discretization. Added this in the discretization itself, that is in the local systems --- src/porepy/numerics/fv/fvutils.py | 138 +++++++++++++++++++++++++++++- src/porepy/numerics/fv/mpfa.py | 59 ++++++++----- src/porepy/params/bc.py | 8 +- 3 files changed, 182 insertions(+), 23 deletions(-) diff --git a/src/porepy/numerics/fv/fvutils.py b/src/porepy/numerics/fv/fvutils.py index f4f023db71..1e9d12be5b 100644 --- a/src/porepy/numerics/fv/fvutils.py +++ b/src/porepy/numerics/fv/fvutils.py @@ -689,8 +689,42 @@ def __init__(self, subcell_topology, bound, nd): shape=(row_dir.size, num_subfno), ).tocsr() - elif self.bc_type == "vectorial": + col_dir = np.argwhere([not it for it in bound.is_rob[fno]]) + row_dir = np.arange(col_dir.size) + self.exclude_rob = sps.coo_matrix( + (np.ones(row_dir.size), (row_dir, col_dir.ravel("C"))), + shape=(row_dir.size, num_subfno), + ).tocsr() + + col_neu_dir = np.argwhere([not it for it in bound.is_neu[fno] | bound.is_dir[fno]]) + row_neu_dir = np.arange(col_neu_dir.size) + self._exclude_neu_dir = sps.coo_matrix( + (np.ones(row_neu_dir.size), (row_neu_dir, col_neu_dir.ravel("C"))), + shape=(row_neu_dir.size, num_subfno), + ).tocsr() + + col_neu_rob = np.argwhere([not it for it in bound.is_neu[fno] | bound.is_rob[fno]]) + row_neu_rob = np.arange(col_neu_rob.size) + self._exclude_neu_rob = sps.coo_matrix( + (np.ones(row_neu_rob.size), (row_neu_rob, col_neu_rob.ravel("C"))), + shape=(row_neu_rob.size, num_subfno), + ).tocsr() + + col_rob_dir = np.argwhere([not it for it in bound.is_dir[fno] | bound.is_rob[fno]]) + row_rob_dir = np.arange(col_rob_dir.size) + self._exclude_rob_dir = sps.coo_matrix( + (np.ones(row_rob_dir.size), (row_rob_dir, col_rob_dir.ravel("C"))), + shape=(row_rob_dir.size, num_subfno), + ).tocsr() + col_rob = np.argwhere([it for it in bound.is_rob[fno]]) + row_rob = np.arange(col_rob.size) + self._keep_robin = sps.coo_matrix( + (np.ones(row_rob.size), (row_rob, col_rob.ravel("C"))), + shape=(row_rob.size, num_subfno), + ).tocsr() + + elif self.bc_type == "vectorial": # Neumann col_neu_x = np.argwhere([not it for it in bound.is_neu[0, fno]]) row_neu_x = np.arange(col_neu_x.size) @@ -862,6 +896,27 @@ def exclude_neumann(self, other): return exclude_neu + def exclude_robin(self, other): + """ Mapping to exclude faces/components with robin boundary conditions from + local systems. + + Parameters: + other (scipy.sparse matrix): Matrix of local equations for + continuity of flux and pressure. + + Returns: + scipy.sparse matrix, with rows corresponding to faces/components with + robin conditions eliminated. + + """ + if self.bc_type == "scalar": + exclude_rob = self.exclude_rob * other + + elif self.bc_type == "vectorial": + raise NotImplementedError('can not exclude robin for vectorial bc') + + return exclude_rob + def exclude_neumann_x(self, other): """ Mapping to exclude components in the x-direction with Neumann boundary conditions from local systems. @@ -940,6 +995,87 @@ def exclude_neumann_nd(self, other): return exclude_neumann_nd * other + def exclude_neu_rob(self, other): + """ Mapping to exclude faces/components with Neumann and Robin boundary + conditions from local systems. + + Parameters: + other (scipy.sparse matrix): Matrix of local equations for + continuity of flux and pressure. + + Returns: + scipy.sparse matrix, with rows corresponding to faces/components with + Neumann conditions eliminated. + + """ + if self.bc_type == "scalar": + exclude_neu = self._exclude_neu_rob* other + + elif self.bc_type == "vectorial": + raise NotImplementedError() + return exclude_neu + + def exclude_neu_dir(self, other): + """ Mapping to exclude faces/components with Neumann and Dirichlet boundary + conditions from local systems. + + Parameters: + other (scipy.sparse matrix): Matrix of local equations for + continuity of flux and pressure. + + Returns: + scipy.sparse matrix, with rows corresponding to faces/components with + Neumann conditions eliminated. + + """ + if self.bc_type == "scalar": + exclude_neu = self._exclude_neu_dir * other + + elif self.bc_type == "vectorial": + raise NotImplementedError() + return exclude_neu + + def exclude_rob_dir(self, other): + """ Mapping to exclude faces/components with Robin and Dirichlet boundary + conditions from local systems. + + Parameters: + other (scipy.sparse matrix): Matrix of local equations for + continuity of flux and pressure. + + Returns: + scipy.sparse matrix, with rows corresponding to faces/components with + Neumann conditions eliminated. + + """ + if self.bc_type == "scalar": + exclude_rob = self._exclude_rob_dir* other + + elif self.bc_type == "vectorial": + raise NotImplementedError() + return exclude_rob + + def keep_robin(self, other): + """ Mapping to exclude faces/components that is not on the Robin boundary + conditions from local systems. + + Parameters: + other (scipy.sparse matrix): Matrix of local equations for + continuity of flux and pressure. + + Returns: + scipy.sparse matrix, with rows corresponding to faces/components with + all but robin conditions eliminated. + + """ + if self.bc_type == "scalar": + exclude_rob = self._keep_robin * other + + elif self.bc_type == "vectorial": + raise NotImplementedError() + return exclude_rob + + # -----------------End of class ExcludeBoundaries----------------------------- diff --git a/src/porepy/numerics/fv/mpfa.py b/src/porepy/numerics/fv/mpfa.py index 9e56f3fbb0..ee37f885dd 100644 --- a/src/porepy/numerics/fv/mpfa.py +++ b/src/porepy/numerics/fv/mpfa.py @@ -549,19 +549,19 @@ def _mpfa_local(g, k, bnd, eta=None, inverter="numba", apertures=None): # Obtain normal_vector * k, pairings of cells and nodes (which together # uniquely define sub-cells, and thus index for gradients. See comment # below for the ordering of elements in the subcell gradient. - nk_grad, cell_node_blocks, sub_cell_index = _tensor_vector_prod( + nk_grad_all, cell_node_blocks, sub_cell_index = _tensor_vector_prod( g, k, subcell_topology, apertures ) # Distance from cell centers to face centers, this will be the # contribution from gradient unknown to equations for pressure continuity - pr_cont_grad = fvutils.compute_dist_face_cell(g, subcell_topology, eta) + pr_cont_grad_all = fvutils.compute_dist_face_cell(g, subcell_topology, eta) # Darcy's law - darcy = -nk_grad[subcell_topology.unique_subfno] + darcy = -nk_grad_all[subcell_topology.unique_subfno] # Pair fluxes over subfaces, that is, enforce conservation - nk_grad = subcell_topology.pair_over_subfaces(nk_grad) + nk_grad_all = subcell_topology.pair_over_subfaces(nk_grad_all) # Contribution from cell center potentials to local systems # For pressure continuity, +-1 (Depending on whether the cell is on the @@ -569,7 +569,7 @@ def _mpfa_local(g, k, bnd, eta=None, inverter="numba", apertures=None): # The .A suffix is necessary to get a numpy array, instead of a scipy # matrix. sgn = g.cell_faces[subcell_topology.fno, subcell_topology.cno].A - pr_cont_cell = sps.coo_matrix( + pr_cont_cell_all = sps.coo_matrix( (sgn[0], (subcell_topology.subfno, subcell_topology.cno)) ).tocsr() # The cell centers give zero contribution to flux continuity @@ -577,6 +577,7 @@ def _mpfa_local(g, k, bnd, eta=None, inverter="numba", apertures=None): (np.zeros(1), (np.zeros(1), np.zeros(1))), shape=(subcell_topology.num_subfno, subcell_topology.num_cno), ).tocsr() + del sgn # Mapping from sub-faces to faces @@ -598,18 +599,25 @@ def _mpfa_local(g, k, bnd, eta=None, inverter="numba", apertures=None): bound_exclusion = fvutils.ExcludeBoundaries(subcell_topology, bnd, g.dim) # No flux conditions for Dirichlet boundary faces - nk_grad = bound_exclusion.exclude_dirichlet(nk_grad) - nk_cell = bound_exclusion.exclude_dirichlet(nk_cell) + nk_grad_n = bound_exclusion.exclude_rob_dir(nk_grad_all) + nk_cell = bound_exclusion.exclude_rob_dir(nk_cell) + + # Robin condition is only applied to Robin boundary faces + nk_grad_r = bound_exclusion.keep_robin(nk_grad_all) + pr_cont_grad_r = bound_exclusion.keep_robin(pr_cont_grad_all) + pr_cont_cell_r = bound_exclusion.keep_robin(pr_cont_cell_all) + + del nk_grad_all # No pressure condition for Neumann boundary faces - pr_cont_grad_all = pr_cont_grad - pr_cont_grad = bound_exclusion.exclude_neumann(pr_cont_grad) - pr_cont_cell = bound_exclusion.exclude_neumann(pr_cont_cell) + pr_cont_grad = bound_exclusion.exclude_neu_rob(pr_cont_grad_all) + pr_cont_cell = bound_exclusion.exclude_neu_rob(pr_cont_cell_all) # So far, the local numbering has been based on the numbering scheme # implemented in SubcellTopology (which treats one cell at a time). For # efficient inversion (below), it is desirable to get the system over to a # block-diagonal structure, with one block centered around each vertex. # Obtain the necessary mappings. + rows2blk_diag, cols2blk_diag, size_of_blocks = _block_diagonal_structure( sub_cell_index, cell_node_blocks, subcell_topology.nno_unique, bound_exclusion ) @@ -618,11 +626,14 @@ def _mpfa_local(g, k, bnd, eta=None, inverter="numba", apertures=None): # System of equations for the subcell gradient variables. On block diagonal # form. - grad_eqs = sps.vstack([nk_grad, pr_cont_grad]) + grad_eqs = sps.vstack([nk_grad_n, nk_grad_r + pr_cont_grad_r, pr_cont_grad]) num_nk_cell = nk_cell.shape[0] + num_nk_rob = nk_grad_r.shape[0] num_pr_cont_grad = pr_cont_grad.shape[0] - del nk_grad + import pdb; pdb.set_trace() + + del nk_grad_n, nk_grad_r, pr_cont_grad_r grad = rows2blk_diag * grad_eqs * cols2blk_diag @@ -666,7 +677,7 @@ def _mpfa_local(g, k, bnd, eta=None, inverter="numba", apertures=None): # it would require quite deep changes to the code. # Flux discretization: - flux = hf2f * darcy * igrad * (-sps.vstack([nk_cell, pr_cont_cell])) + flux = hf2f * darcy * igrad * (-sps.vstack([nk_cell, pr_cont_cell_r, pr_cont_cell])) #### # Boundary conditions @@ -676,7 +687,7 @@ def _mpfa_local(g, k, bnd, eta=None, inverter="numba", apertures=None): subcell_topology, sgn_unique, g, - num_nk_cell, + num_nk_cell + num_nk_rob, num_pr_cont_grad, ) # Discretization of boundary values @@ -695,7 +706,7 @@ def _mpfa_local(g, k, bnd, eta=None, inverter="numba", apertures=None): * hf2f * pr_cont_grad_all * igrad - * (-sps.vstack([nk_cell, pr_cont_cell])) + * (-sps.vstack([nk_cell, pr_cont_cell_r, pr_cont_cell])) ) # Internal faces, and boundary faces with a Dirichle condition do not need @@ -871,10 +882,13 @@ def _block_diagonal_structure(sub_cell_index, cell_node_blocks, nno, bound_exclu # Stack node numbers of equations on top of each other, and sort them to # get block-structure. First eliminate node numbers at the boundary, where - # the equations are either of flux or pressure continuity (not both) - nno_flux = bound_exclusion.exclude_dirichlet(nno) - nno_pressure = bound_exclusion.exclude_neumann(nno) - node_occ = np.hstack((nno_flux, nno_pressure)) + # the equations are either of flux, pressure continuity or robin + nno_flux = bound_exclusion.exclude_rob_dir(nno) + nno_pressure = bound_exclusion.exclude_neu_rob(nno) + # we have now eliminated all nodes related to robin, we therefore add them + nno_rob = bound_exclusion.keep_robin(nno) + + node_occ = np.hstack((nno_flux, nno_rob, nno_pressure)) sorted_ind = np.argsort(node_occ) sorted_nodes_rows = node_occ[sorted_ind] # Size of block systems @@ -919,14 +933,17 @@ def _create_bound_rhs(bnd, bound_exclusion, subcell_topology, sgn, g, num_flux, # are handled by assigning Neumann conditions. is_dir = np.logical_and(bnd.is_dir, np.logical_not(bnd.is_internal)) is_neu = np.logical_or(bnd.is_neu, bnd.is_internal) + is_rob = np.logical_and(bnd.is_rob, np.logical_not(bnd.is_internal)) + is_neu = np.logical_or(is_neu, is_rob) fno = subcell_topology.fno_unique num_neu = np.sum(is_neu[fno]) num_dir = np.sum(is_dir[fno]) num_bound = num_neu + num_dir - # Neumann boundary conditions - # Find Neumann faces, exclude Dirichlet faces (since these are excluded + # Neumann and Robin boundary conditions. We can handle Neumann and robin + # since they both are given in units of force. + # Find Neumann and Robin faces, exclude Dirichlet faces (since these are excluded # from the right hand side linear system), and do necessary formating. neu_ind = np.argwhere( bound_exclusion.exclude_dirichlet(is_neu[fno].astype("int64")) diff --git a/src/porepy/params/bc.py b/src/porepy/params/bc.py index 4f4a8bfab9..ec93111c8b 100644 --- a/src/porepy/params/bc.py +++ b/src/porepy/params/bc.py @@ -64,6 +64,7 @@ def __init__(self, g, faces=None, cond=None): self.is_neu = np.zeros(self.num_faces, dtype=bool) self.is_dir = np.zeros(self.num_faces, dtype=bool) + self.is_rob = np.zeros(self.num_faces, dtype=bool) # By default, all faces are Neumann. self.is_neu[bf] = True @@ -103,8 +104,13 @@ def __init__(self, g, faces=None, cond=None): elif s.lower() == "dir": self.is_dir[faces[l]] = True self.is_neu[faces[l]] = False + self.is_rob[faces[l]] = False + elif s.lower() == "rob": + self.is_dir[faces[l]] = False + self.is_neu[faces[l]] = False + self.is_rob[faces[l]] = True else: - raise ValueError("Boundary should be Dirichlet or Neumann") + raise ValueError("Boundary should be Dirichlet, Neumann or Robin") class BoundaryConditionNode(object): From f9e0a9de4da85a0e4c158cdb450046e3073aaa38 Mon Sep 17 00:00:00 2001 From: Runar Date: Fri, 7 Sep 2018 18:24:07 +0200 Subject: [PATCH 03/26] Added discretization of Robin condition in the local system for MPFA --- src/porepy/numerics/fv/mpfa.py | 125 +++++++++++++++++++++++---------- 1 file changed, 86 insertions(+), 39 deletions(-) diff --git a/src/porepy/numerics/fv/mpfa.py b/src/porepy/numerics/fv/mpfa.py index ee37f885dd..4a992b6e17 100644 --- a/src/porepy/numerics/fv/mpfa.py +++ b/src/porepy/numerics/fv/mpfa.py @@ -178,7 +178,7 @@ def discretize(self, g, data): # ------------------------------------------------------------------------------# -def mpfa(g, k, bnd, eta=None, inverter=None, apertures=None, max_memory=None, **kwargs): +def mpfa(g, k, bnd, eta=None, inverter=None, apertures=None, max_memory=None, alpha=None, **kwargs): """ Discretize the scalar elliptic equation by the multi-point flux approximation method. @@ -269,7 +269,7 @@ def mpfa(g, k, bnd, eta=None, inverter=None, apertures=None, max_memory=None, ** # TODO: We may want to estimate the memory need, and give a warning if # this seems excessive flux, bound_flux, bound_pressure_cell, bound_pressure_face = _mpfa_local( - g, k, bnd, eta=eta, inverter=inverter, apertures=apertures + g, k, bnd, eta=eta, inverter=inverter, apertures=apertures, alpha=alpha ) else: # Estimate number of partitions necessary based on prescribed memory @@ -456,7 +456,7 @@ def mpfa_partial( ) -def _mpfa_local(g, k, bnd, eta=None, inverter="numba", apertures=None): +def _mpfa_local(g, k, bnd, eta=None, inverter="numba", apertures=None, alpha=None): """ Actual implementation of the MPFA O-method. To calculate MPFA on a grid directly, either call this method, or, to respect the privacy of this @@ -515,6 +515,9 @@ def _mpfa_local(g, k, bnd, eta=None, inverter="numba", apertures=None): elif g.dim == 0: return sps.csr_matrix([0]), 0, 0, 0 + if alpha is None and np.sum(bnd.is_rob)>0: + raise ValueError('If applying Robin conditions you must supply an alpha') + # The grid coordinates are always three-dimensional, even if the grid is # really 2D. This means that there is not a 1-1 relation between the number # of coordinates of a point / vector and the real dimension. This again @@ -577,8 +580,23 @@ def _mpfa_local(g, k, bnd, eta=None, inverter="numba", apertures=None): (np.zeros(1), (np.zeros(1), np.zeros(1))), shape=(subcell_topology.num_subfno, subcell_topology.num_cno), ).tocsr() - - del sgn + + # For the Robin codition the distance from the cell centers to face centers + # will be the contribution from the gradients. We integrate over the subface + # and multiply by the area + num_nodes = np.diff(g.face_nodes.indptr) + scaled_sgn = ( + alpha * sgn[0] * g.face_areas[subcell_topology.fno] \ + / num_nodes[subcell_topology.fno] + ) + # pair_over_subfaces flips the sign so we flip it back + pr_trace_grad_all = pr_cont_grad_all * sps.diags(scaled_sgn) + pr_trace_cell_all = sps.coo_matrix( + (alpha * g.face_areas[subcell_topology.fno] / num_nodes[subcell_topology.fno], + (subcell_topology.subfno, subcell_topology.cno)) + ).tocsr() + + del sgn, scaled_sgn # Mapping from sub-faces to faces hf2f = sps.coo_matrix( @@ -604,11 +622,11 @@ def _mpfa_local(g, k, bnd, eta=None, inverter="numba", apertures=None): # Robin condition is only applied to Robin boundary faces nk_grad_r = bound_exclusion.keep_robin(nk_grad_all) - pr_cont_grad_r = bound_exclusion.keep_robin(pr_cont_grad_all) - pr_cont_cell_r = bound_exclusion.keep_robin(pr_cont_cell_all) + pr_trace_grad = bound_exclusion.keep_robin(pr_trace_grad_all) + pr_trace_cell = bound_exclusion.keep_robin(pr_trace_cell_all) del nk_grad_all - # No pressure condition for Neumann boundary faces + # No pressure condition for Neumann or Robin boundary faces pr_cont_grad = bound_exclusion.exclude_neu_rob(pr_cont_grad_all) pr_cont_cell = bound_exclusion.exclude_neu_rob(pr_cont_cell_all) @@ -626,14 +644,20 @@ def _mpfa_local(g, k, bnd, eta=None, inverter="numba", apertures=None): # System of equations for the subcell gradient variables. On block diagonal # form. - grad_eqs = sps.vstack([nk_grad_n, nk_grad_r + pr_cont_grad_r, pr_cont_grad]) + # NOTE: I think in the discretization for sub_cells a flow out of the cell is + # negative. This is a contradiction to what is done for the boundary conditions + # where we want to set dot(n, flux) where n is the normal pointing outwards. + # thats why we need +nk_grad_r - pr_trace_grad -pr_trace_cell instead of = rhs + # instead of how we would expect: -nk_grad_r + pr_trace_grad +pr_trace_cell= rhs. + # This is also why we multiply with -1 in scaled_sgn in _create_bound_rhs + grad_eqs = sps.vstack([nk_grad_n, nk_grad_r - pr_trace_grad, pr_cont_grad]) num_nk_cell = nk_cell.shape[0] num_nk_rob = nk_grad_r.shape[0] num_pr_cont_grad = pr_cont_grad.shape[0] - import pdb; pdb.set_trace() - del nk_grad_n, nk_grad_r, pr_cont_grad_r + + del nk_grad_n, nk_grad_r, pr_trace_grad grad = rows2blk_diag * grad_eqs * cols2blk_diag @@ -677,7 +701,8 @@ def _mpfa_local(g, k, bnd, eta=None, inverter="numba", apertures=None): # it would require quite deep changes to the code. # Flux discretization: - flux = hf2f * darcy * igrad * (-sps.vstack([nk_cell, pr_cont_cell_r, pr_cont_cell])) + # The negative in front of pr_trace_cell comes from the grad_egs + flux = hf2f * darcy * igrad * (-sps.vstack([nk_cell, -pr_trace_cell, pr_cont_cell])) #### # Boundary conditions @@ -687,7 +712,8 @@ def _mpfa_local(g, k, bnd, eta=None, inverter="numba", apertures=None): subcell_topology, sgn_unique, g, - num_nk_cell + num_nk_rob, + num_nk_cell, + num_nk_rob, num_pr_cont_grad, ) # Discretization of boundary values @@ -706,7 +732,7 @@ def _mpfa_local(g, k, bnd, eta=None, inverter="numba", apertures=None): * hf2f * pr_cont_grad_all * igrad - * (-sps.vstack([nk_cell, pr_cont_cell_r, pr_cont_cell])) + * (-sps.vstack([nk_cell, pr_trace_cell, pr_cont_cell])) ) # Internal faces, and boundary faces with a Dirichle condition do not need @@ -855,7 +881,7 @@ def _tensor_vector_prod(g, k, subcell_topology, apertures=None): ) nk = normals_mat * k_mat - + # Unique sub-cell indexes are pulled from column indices, we only need # every nd column (since nd faces of the cell meet at each vertex) sub_cell_ind = j[::, 0::nd] @@ -908,7 +934,7 @@ def _block_diagonal_structure(sub_cell_index, cell_node_blocks, nno, bound_exclu return rows2blk_diag, cols2blk_diag, size_of_blocks -def _create_bound_rhs(bnd, bound_exclusion, subcell_topology, sgn, g, num_flux, num_pr): +def _create_bound_rhs(bnd, bound_exclusion, subcell_topology, sgn, g, num_flux, num_rob, num_pr): """ Define rhs matrix to get basis functions for incorporates boundary conditions @@ -934,51 +960,71 @@ def _create_bound_rhs(bnd, bound_exclusion, subcell_topology, sgn, g, num_flux, is_dir = np.logical_and(bnd.is_dir, np.logical_not(bnd.is_internal)) is_neu = np.logical_or(bnd.is_neu, bnd.is_internal) is_rob = np.logical_and(bnd.is_rob, np.logical_not(bnd.is_internal)) - is_neu = np.logical_or(is_neu, is_rob) fno = subcell_topology.fno_unique num_neu = np.sum(is_neu[fno]) num_dir = np.sum(is_dir[fno]) - num_bound = num_neu + num_dir - - # Neumann and Robin boundary conditions. We can handle Neumann and robin - # since they both are given in units of force. + num_rob = np.sum(is_rob[fno]) + num_bound = num_neu + num_dir + num_rob + + # Neumann and Robin boundary conditions. Neumann and Robin conditions + # are essentially the same for the rhs (the rhs for both is a flux). + # However, we need to be carefull and get the indexing correct as seen + # from the local system, that is, first Neumann, then Robin and last + # Dirichlet. # Find Neumann and Robin faces, exclude Dirichlet faces (since these are excluded # from the right hand side linear system), and do necessary formating. neu_ind = np.argwhere( - bound_exclusion.exclude_dirichlet(is_neu[fno].astype("int64")) + bound_exclusion.exclude_rob_dir(is_neu[fno].astype("int64")) + ).ravel("F") + rob_ind = np.argwhere( + bound_exclusion.keep_robin(is_rob[fno].astype("int64")) + ).ravel("F") + neu_rob_ind = np.argwhere( + bound_exclusion.exclude_dirichlet((is_rob[fno] + is_neu[fno]).astype("int64")) ).ravel("F") - # We also need to map the respective Neumann and Dirichlet half-faces to + + # We also need to map the respective Neumann, Robin, and Dirichlet half-faces to # the global half-face numbering (also interior faces). The latter should # not have Dirichlet and Neumann excluded (respectively), and thus we need # new fields neu_ind_all = np.argwhere(is_neu[fno].astype("int")).ravel("F") + rob_ind_all = np.argwhere(is_rob[fno].astype("int")).ravel("F") dir_ind_all = np.argwhere(is_dir[fno].astype("int")).ravel("F") num_face_nodes = g.face_nodes.sum(axis=0).A.ravel(order="F") - # For the Neumann boundary conditions, we define the value as seen from + # We now merge the neuman and robin indices since they are treated equivalent + if rob_ind.size==0: + neu_rob_ind = neu_ind + elif neu_ind.size==0: + neu_rob_ind = rob_ind + else: + neu_rob_ind = np.hstack((neu_ind, rob_ind + num_flux)) + neu_rob_ind_all = np.hstack((neu_ind_all, rob_ind_all)) + + # For the Neumann/Robin boundary conditions, we define the value as seen from # the innside of the domain. E.g. outflow is defined to be positive. We # therefore set the matrix indices to -1. We also have to scale it with # the number of nodes per face because the flux of face is the sum of its # half-faces. - scaled_sgn = -1 / num_face_nodes[fno[neu_ind_all]] - if neu_ind.size > 0: - neu_cell = sps.coo_matrix( - (scaled_sgn, (neu_ind, np.arange(neu_ind.size))), - shape=(num_flux, num_bound), + scaled_sgn = -1 / num_face_nodes[fno[neu_rob_ind_all]] + if neu_rob_ind.size > 0: + neu_rob_cell = sps.coo_matrix( + (scaled_sgn, (neu_rob_ind, np.arange(neu_rob_ind.size))), + shape=(num_flux+num_rob, num_bound), ) else: # Special handling when no elements are found. Not sure if this is # necessary, or if it is me being stupid - neu_cell = sps.coo_matrix((num_flux, num_bound)) + neu_rob_cell = sps.coo_matrix((num_flux + num_rob, num_bound)) # Dirichlet boundary conditions dir_ind = np.argwhere( - bound_exclusion.exclude_neumann(is_dir[fno].astype("int64")) + bound_exclusion.exclude_neu_rob(is_dir[fno].astype("int64")) ).ravel("F") if dir_ind.size > 0: dir_cell = sps.coo_matrix( - (sgn[dir_ind_all], (dir_ind, num_neu + np.arange(dir_ind.size))), + (sgn[dir_ind_all], (dir_ind, num_neu + num_rob + np.arange(dir_ind.size))), shape=(num_pr, num_bound), ) else: @@ -988,12 +1034,12 @@ def _create_bound_rhs(bnd, bound_exclusion, subcell_topology, sgn, g, num_flux, # Number of elements in neu_ind and neu_ind_all are equal, we can test with # any of them. Same with dir. - if neu_ind.size > 0 and dir_ind.size > 0: - neu_dir_ind = np.hstack([neu_ind_all, dir_ind_all]).ravel("F") - elif neu_ind.size > 0: - neu_dir_ind = neu_ind_all + if neu_rob_ind.size > 0 and dir_ind.size > 0: + neu_rob_dir_ind = np.hstack([neu_rob_ind_all, dir_ind_all]).ravel("F") + elif neu_rob_ind.size > 0: + neu_rob_dir_ind = neu_rob_ind_all elif dir_ind.size > 0: - neu_dir_ind = dir_ind_all + neu_rob_dir_ind = dir_ind_all else: raise ValueError("Boundary values should be either Dirichlet or " "Neumann") @@ -1002,7 +1048,7 @@ def _create_bound_rhs(bnd, bound_exclusion, subcell_topology, sgn, g, num_flux, # The columns in neu_cell, dir_cell are ordered from 0 to num_bound-1. # Map these to all half-face indices bnd_2_all_hf = sps.coo_matrix( - (np.ones(num_bound), (np.arange(num_bound), neu_dir_ind)), + (np.ones(num_bound), (np.arange(num_bound), neu_rob_dir_ind)), shape=(num_bound, num_subfno), ) # The user of the discretization should now nothing about half faces, @@ -1014,5 +1060,6 @@ def _create_bound_rhs(bnd, bound_exclusion, subcell_topology, sgn, g, num_flux, ), shape=(num_subfno, g.num_faces), ) - rhs_bound = sps.vstack([neu_cell, dir_cell]) * bnd_2_all_hf * hf_2_f + rhs_bound = sps.vstack([neu_rob_cell, dir_cell]) * bnd_2_all_hf * hf_2_f + return rhs_bound From 622253cb763caf0ff782b3c242913ae18698eb14 Mon Sep 17 00:00:00 2001 From: Runar Date: Fri, 7 Sep 2018 18:24:56 +0200 Subject: [PATCH 04/26] added unittest for mpfa robin condition --- src/porepy/numerics/fv/mpfa.py | 8 ++- test/unit/test_mpfa_robin.py | 105 +++++++++++++++++++++++++++++++++ 2 files changed, 110 insertions(+), 3 deletions(-) create mode 100644 test/unit/test_mpfa_robin.py diff --git a/src/porepy/numerics/fv/mpfa.py b/src/porepy/numerics/fv/mpfa.py index 4a992b6e17..3f52bb89e8 100644 --- a/src/porepy/numerics/fv/mpfa.py +++ b/src/porepy/numerics/fv/mpfa.py @@ -515,9 +515,11 @@ def _mpfa_local(g, k, bnd, eta=None, inverter="numba", apertures=None, alpha=Non elif g.dim == 0: return sps.csr_matrix([0]), 0, 0, 0 - if alpha is None and np.sum(bnd.is_rob)>0: - raise ValueError('If applying Robin conditions you must supply an alpha') - + if alpha is None: + if np.sum(bnd.is_rob) !=0: + raise ValueError('If applying Robin conditions you must supply an alpha') + else: + alpha = 1 # The grid coordinates are always three-dimensional, even if the grid is # really 2D. This means that there is not a 1-1 relation between the number # of coordinates of a point / vector and the real dimension. This again diff --git a/test/unit/test_mpfa_robin.py b/test/unit/test_mpfa_robin.py new file mode 100644 index 0000000000..d87806b652 --- /dev/null +++ b/test/unit/test_mpfa_robin.py @@ -0,0 +1,105 @@ +import numpy as np +import scipy.sparse as sps +import unittest + +import porepy as pp + +class RobinBoundTest(unittest.TestCase): + def test_flow_left_right(self): + nxs = [1, 3] + nys = [1, 3] + for nx in nxs: + for ny in nys: + g = pp.CartGrid([nx, ny], physdims=[1, 1]) + g.compute_geometry() + k = pp.SecondOrderTensor(2, np.ones(g.num_cells)) + alpha = 1.5 + + left = g.face_centers[0] < 1e-10 + right = g.face_centers[0] > 1 - 1e-10 + + dir_ind = np.ravel(np.argwhere(left)) + rob_ind = np.ravel(np.argwhere(right)) + + names = ['dir']*len(dir_ind) + ['rob'] * len(rob_ind) + bnd_ind = np.hstack((dir_ind, rob_ind)) + bnd = pp.BoundaryCondition(g, bnd_ind, names) + + p_bound = 1 + rob_bound = 0 + C = (rob_bound - alpha*p_bound) / (alpha - p_bound) + u_ex = - C * g.face_normals[0] + p_ex = C * g.cell_centers[0] + p_bound + self.solve_robin(g, k, bnd, alpha, p_bound, rob_bound, dir_ind, rob_ind, p_ex, u_ex) + + def test_flow_nz_rhs(self): + nxs = [1, 3] + nys = [1, 3] + for nx in nxs: + for ny in nys: + g = pp.CartGrid([nx, ny], physdims=[1, 1]) + g.compute_geometry() + k = pp.SecondOrderTensor(2, np.ones(g.num_cells)) + alpha = 1.5 + + left = g.face_centers[0] < 1e-10 + right = g.face_centers[0] > 1 - 1e-10 + + dir_ind = np.ravel(np.argwhere(left)) + rob_ind = np.ravel(np.argwhere(right)) + + names = ['dir']*len(dir_ind) + ['rob'] * len(rob_ind) + bnd_ind = np.hstack((dir_ind, rob_ind)) + bnd = pp.BoundaryCondition(g, bnd_ind, names) + + p_bound = 1 + rob_bound = 1 + + C = (rob_bound - alpha*p_bound) / (alpha - p_bound) + u_ex = - C * g.face_normals[0] + p_ex = C * g.cell_centers[0] + p_bound + self.solve_robin(g, k, bnd, alpha, p_bound, rob_bound, dir_ind, rob_ind, p_ex, u_ex) + def test_flow_down(self): + nxs = [1, 3] + nys = [1, 3] + for nx in nxs: + for ny in nys: + g = pp.CartGrid([nx, ny], physdims=[1, 1]) + g.compute_geometry() + k = pp.SecondOrderTensor(2, np.ones(g.num_cells)) + alpha = 1.5 + + bot = g.face_centers[1] < 1e-10 + top = g.face_centers[1] > 1 - 1e-10 + + dir_ind = np.ravel(np.argwhere(top)) + rob_ind = np.ravel(np.argwhere(bot)) + + names = ['dir']*len(dir_ind) + ['rob'] * len(rob_ind) + bnd_ind = np.hstack((dir_ind, rob_ind)) + bnd = pp.BoundaryCondition(g, bnd_ind, names) + + p_bound = 1 + rob_bound = 1 + C = (rob_bound - alpha*p_bound) / (alpha - p_bound) + u_ex = C * g.face_normals[1] + p_ex = C *(1- g.cell_centers[1]) + p_bound + self.solve_robin(g, k, bnd, alpha, p_bound, rob_bound, dir_ind, rob_ind, p_ex, u_ex) + + + def solve_robin(self, g, k, bnd, alpha, p_bound, rob_bound, dir_ind, rob_ind, p_ex, u_ex): + flux, bound_flux, _, _ = pp.numerics.fv.mpfa._mpfa_local(g, k, bnd, alpha=alpha) + + div = pp.fvutils.scalar_divergence(g) + + u_bound = np.zeros(g.num_faces) + u_bound[dir_ind] = p_bound + u_bound[rob_ind] = rob_bound * g.face_areas[rob_ind] + + + a = div * flux + b = -div * bound_flux * u_bound + + p = np.linalg.solve(a.A, b) + assert np.allclose(p, p_ex) + assert np.allclose(flux * p + bound_flux * u_bound, u_ex) From 21b2fc707138e1a77589711330ae1f53d699f590 Mon Sep 17 00:00:00 2001 From: Runar Date: Fri, 7 Sep 2018 18:33:01 +0200 Subject: [PATCH 05/26] Added Robin condition to the convergence tests --- examples/example2/ConvergenceTestMPFA.ipynb | 141 +++++++++++++++++++- 1 file changed, 139 insertions(+), 2 deletions(-) diff --git a/examples/example2/ConvergenceTestMPFA.ipynb b/examples/example2/ConvergenceTestMPFA.ipynb index 061e056e18..646f1477b8 100644 --- a/examples/example2/ConvergenceTestMPFA.ipynb +++ b/examples/example2/ConvergenceTestMPFA.ipynb @@ -294,6 +294,143 @@ "assert np.abs(f_triang_nopert[2] - 0.0033976662367782447) < 1e-10\n" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Robin boundary condition" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Cartesian errors - no perturbations\n", + "[0.562515011909662, 0.14015471392164033, 0.034984342689359696] pressure error\n", + "[0.22376060444435023, 0.061834410597352465, 0.01641158802403413] flux error\n", + "Triangular errors - no perturbations\n", + "[1.1696464364631856, 0.8830933068427524, 0.49114177612455606] pressure error\n", + "[0.5558236600090738, 0.40733862448942615, 0.22657669056462687] flux error\n" + ] + } + ], + "source": [ + "# Analytical solution\n", + "x, y = sympy.symbols('x y')\n", + "u = x * (1-x) * sympy.sin(x) * sympy.cos(y)\n", + "u_f = sympy.lambdify((x, y), u, 'numpy')\n", + "dux = sympy.diff(u, x)\n", + "duy = sympy.diff(u, y)\n", + "dux_f = sympy.lambdify((x, y), dux, 'numpy')\n", + "duy_f = sympy.lambdify((x, y), duy, 'numpy')\n", + "rhs = -sympy.diff(dux, x) - sympy.diff(duy, y)\n", + "rhs_f = sympy.lambdify((x, y), rhs, 'numpy')\n", + "\n", + "###\n", + "# This is where parameters can be modified to alter the convergence test.\n", + "# The remaining lines should \n", + "np.random.seed(42)\n", + "base = 4\n", + "domain = np.array([1, 1])\n", + "basedim = np.array([base, base])\n", + "num_refs = 3\n", + "grid_type = 'cart'\n", + "\n", + "def run_convergence(grid_type, pert):\n", + " u_err = []\n", + " flux_err = []\n", + "\n", + " for g in setup_grids.grid_sequence(basedim, num_refs, grid_type, pert, domain=domain):\n", + " # Reset the random seed for every grid realization.\n", + " # This should make no difference for the convergence test, \n", + " # but it makes sure that we can run unit tests based on the values obtained\n", + " # here.\n", + " np.random.seed(42)\n", + " \n", + " # Permeability tensor\n", + " k = perm(2, np.ones(g.num_cells))\n", + " alpha = 1.2\n", + " # Exact solution\n", + " xf = g.face_centers\n", + " xc = g.cell_centers\n", + " u_ex = u_f(xc[0], xc[1])\n", + " du_ex_faces = np.vstack((dux_f(xf[0], xf[1]), duy_f(xf[0], xf[1])))\n", + " flux_ex = -np.sum(g.face_normals[:2] * du_ex_faces, axis=0)\n", + "\n", + " # Set type of boundary conditions - Dirichlet\n", + " bound_faces = g.tags['domain_boundary_faces'].nonzero()[0]\n", + " n = g.nodes\n", + " top = bc.face_on_side(g, 'ymax')\n", + " bot = bc.face_on_side(g, 'ymin')\n", + " left = bc.face_on_side(g, 'xmin')\n", + " right = bc.face_on_side(g, 'xmax')\n", + "\n", + " dir_faces = np.asarray(left[0])\n", + " rob_faces = np.hstack((right[0]))\n", + " neu_faces = np.hstack((bot[0], top[0]))\n", + " bc_faces = np.hstack([dir_faces, rob_faces])\n", + " bc_names = ['dir']*len(dir_faces) + ['rob'] * len(rob_faces)\n", + " bound_cond = bc.BoundaryCondition(g, bc_faces, bc_names)\n", + " \n", + " # MPFA discretization, and system matrix\n", + " flux, bound_flux, _, _ = mpfa.mpfa(g, k, bound_cond, alpha=alpha)\n", + " div = fvutils.scalar_divergence(g)\n", + " a = div * flux\n", + " \n", + " # Boundary conditions\n", + " nfi, _, sgn_n = sps.find(g.cell_faces[neu_faces,:])\n", + " rfi, _, sgn = sps.find(g.cell_faces[rob_faces,:])\n", + " u_bound = np.zeros(g.num_faces)\n", + " u_bound[dir_faces] = u_f(xf[0, dir_faces], xf[1, dir_faces])\n", + " # innflow is always negative, so need to flip flux according to normal direction.\n", + " u_bound[neu_faces[nfi]] = flux_ex[neu_faces[nfi]] * (sgn_n)\n", + " u_bound[rob_faces[rfi]] = flux_ex[rob_faces[rfi]] * (sgn) +\\\n", + " alpha * u_f(xf[0, rob_faces], xf[1, rob_faces]) * g.face_areas[rob_faces]\n", + "\n", + " # Right hand side - contribution from the solution and the boundary conditions\n", + " xc = g.cell_centers\n", + " b = rhs_f(xc[0], xc[1]) * g.cell_volumes - div * bound_flux * u_bound\n", + "\n", + " # Solve system, derive fluxes\n", + " u_num = spsolve(a, b)\n", + " flux_num = flux * u_num + bound_flux * u_bound\n", + "\n", + " # Exact solution\n", + " flux_diff = flux_num - flux_ex\n", + " \n", + " u_err.append(np.sqrt(np.sum(g.cell_volumes * (u_num - u_ex)**2)) /\n", + " np.sqrt(np.sum(g.cell_volumes * u_ex**2)))\n", + " flux_err.append(np.sqrt(np.sum((g.face_areas ** g.dim) * flux_diff**2))/\n", + " np.sqrt(np.sum((g.face_areas ** g.dim) * flux_ex**2)))\n", + " return u_err, flux_err\n", + "\n", + "grids = ['cart', 'triangular']\n", + "\n", + "## No tests on perturbed grids - I'm too lazy to find faces of the respective boundaries\n", + "\n", + "\n", + "u_cart_nopert, f_cart_nopert = run_convergence('cart', 0)\n", + "print('Cartesian errors - no perturbations')\n", + "print(u_cart_nopert,'pressure error')\n", + "print(f_cart_nopert,'flux error')\n", + "\n", + "u_triang_nopert, f_triang_nopert = run_convergence('triangular', 0)\n", + "print('Triangular errors - no perturbations')\n", + "print(u_triang_nopert,'pressure error')\n", + "print(f_triang_nopert,'flux error')\n", + "\n", + "# RB: These values were hard-coded 07.09.2018.\n", + "assert np.abs(u_cart_nopert[2] - 0.034984342689359696) < 1e-10\n", + "assert np.abs(f_cart_nopert[2] - 0.01641158802403413) < 1e-10\n", + "assert np.abs(u_triang_nopert[2] - 0.49114177612455606) < 1e-10\n", + "assert np.abs(f_triang_nopert[2] - 0.22657669056462687) < 1e-10\n" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -304,7 +441,7 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 5, "metadata": {}, "outputs": [ { @@ -471,7 +608,7 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 6, "metadata": {}, "outputs": [ { From c8eb5fe83f4518d12e8d52460c346d2eaac8aa28 Mon Sep 17 00:00:00 2001 From: Runar Date: Mon, 10 Sep 2018 16:18:09 +0200 Subject: [PATCH 06/26] Fixed bug in mpfa when applying only robin + dirichlet That is, applying Robin without any Neuman boundaries --- src/porepy/numerics/fv/mpfa.py | 7 +++--- test/unit/test_mpfa_robin.py | 41 ++++++++++++++++++++++++++++++++++ 2 files changed, 45 insertions(+), 3 deletions(-) diff --git a/src/porepy/numerics/fv/mpfa.py b/src/porepy/numerics/fv/mpfa.py index 3f52bb89e8..edf1fc1f81 100644 --- a/src/porepy/numerics/fv/mpfa.py +++ b/src/porepy/numerics/fv/mpfa.py @@ -591,7 +591,8 @@ def _mpfa_local(g, k, bnd, eta=None, inverter="numba", apertures=None, alpha=Non alpha * sgn[0] * g.face_areas[subcell_topology.fno] \ / num_nodes[subcell_topology.fno] ) - # pair_over_subfaces flips the sign so we flip it back + + # pair_over_subfaces flips the sign so we flip it backo pr_trace_grad_all = pr_cont_grad_all * sps.diags(scaled_sgn) pr_trace_cell_all = sps.coo_matrix( (alpha * g.face_areas[subcell_topology.fno] / num_nodes[subcell_topology.fno], @@ -966,7 +967,7 @@ def _create_bound_rhs(bnd, bound_exclusion, subcell_topology, sgn, g, num_flux, fno = subcell_topology.fno_unique num_neu = np.sum(is_neu[fno]) num_dir = np.sum(is_dir[fno]) - num_rob = np.sum(is_rob[fno]) + assert num_rob == np.sum(is_rob[fno]) num_bound = num_neu + num_dir + num_rob # Neumann and Robin boundary conditions. Neumann and Robin conditions @@ -999,7 +1000,7 @@ def _create_bound_rhs(bnd, bound_exclusion, subcell_topology, sgn, g, num_flux, if rob_ind.size==0: neu_rob_ind = neu_ind elif neu_ind.size==0: - neu_rob_ind = rob_ind + neu_rob_ind = rob_ind + num_flux else: neu_rob_ind = np.hstack((neu_ind, rob_ind + num_flux)) neu_rob_ind_all = np.hstack((neu_ind_all, rob_ind_all)) diff --git a/test/unit/test_mpfa_robin.py b/test/unit/test_mpfa_robin.py index d87806b652..a3adffdf7a 100644 --- a/test/unit/test_mpfa_robin.py +++ b/test/unit/test_mpfa_robin.py @@ -86,6 +86,47 @@ def test_flow_down(self): p_ex = C *(1- g.cell_centers[1]) + p_bound self.solve_robin(g, k, bnd, alpha, p_bound, rob_bound, dir_ind, rob_ind, p_ex, u_ex) + def test_no_neumann(self): + g = pp.CartGrid([2, 2], physdims=[1, 1]) + g.compute_geometry() + k = pp.SecondOrderTensor(2, np.ones(g.num_cells)) + alpha = 1.5 + + bot = g.face_centers[1] < 1e-10 + top = g.face_centers[1] > 1 - 1e-10 + left = g.face_centers[0] < 1e-10 + right = g.face_centers[0] > 1 - 1e-10 + + dir_ind = np.ravel(np.argwhere(top + left)) + rob_ind = np.ravel(np.argwhere(bot + right)) + + names = ['dir']*len(dir_ind) + ['rob'] * len(rob_ind) + bnd_ind = np.hstack((dir_ind, rob_ind)) + bnd = pp.BoundaryCondition(g, bnd_ind, names) + + flux, bound_flux, _, _ = pp.numerics.fv.mpfa._mpfa_local(g, k, bnd, alpha=alpha) + + div = pp.fvutils.scalar_divergence(g) + + rob_ex = [alpha * .25, alpha * .75, 1, 1] + u_bound = np.zeros(g.num_faces) + u_bound[dir_ind] = g.face_centers[1, dir_ind] + u_bound[rob_ind] = rob_ex * g.face_areas[rob_ind] + + + a = div * flux + b = -div * bound_flux * u_bound + + p = np.linalg.solve(a.A, b) + + u_ex = [np.dot(g.face_normals[:, f], np.array([0, -1, 0])) + for f in range(g.num_faces)] + u_ex = np.array(u_ex).ravel('F') + p_ex = g.cell_centers[1] + + assert np.allclose(p, p_ex) + assert np.allclose(flux * p + bound_flux * u_bound, u_ex) + def solve_robin(self, g, k, bnd, alpha, p_bound, rob_bound, dir_ind, rob_ind, p_ex, u_ex): flux, bound_flux, _, _ = pp.numerics.fv.mpfa._mpfa_local(g, k, bnd, alpha=alpha) From 4e1c573beee570af3a0d9ec1a6c91e259ae9304e Mon Sep 17 00:00:00 2001 From: Runar Date: Mon, 10 Sep 2018 16:53:36 +0200 Subject: [PATCH 07/26] Added Robin condition for the elasticity problem --- src/porepy/numerics/fv/fvutils.py | 51 +++++++++++ src/porepy/numerics/fv/mpsa.py | 143 +++++++++++++++++++++++------- test/unit/test_mpfa_robin.py | 1 - test/unit/test_mpsa_robin.py | 110 +++++++++++++++++++++++ 4 files changed, 271 insertions(+), 34 deletions(-) create mode 100644 test/unit/test_mpsa_robin.py diff --git a/src/porepy/numerics/fv/fvutils.py b/src/porepy/numerics/fv/fvutils.py index 1e9d12be5b..c825845cf4 100644 --- a/src/porepy/numerics/fv/fvutils.py +++ b/src/porepy/numerics/fv/fvutils.py @@ -995,6 +995,57 @@ def exclude_neumann_nd(self, other): return exclude_neumann_nd * other + def exclude_neu_rob_nd(self, other): + """ Exclusion of Neumann and robin conditions for vector equations (elasticity). + See above method without _nd suffix for description. + + """ + if self.bc_type == "scalar": + exclude_neumann_nd = sps.kron(sps.eye(self.nd), self._exclude_neu_rob) + + elif self.bc_type == "vectorial": + raise NotImplementedError() + + return exclude_neumann_nd * other + + def exclude_neu_dir_nd(self, other): + """ Exclusion of Neumann and Dirichlet conditions for vector equations + (elasticity). See above method without _nd suffix for description. + + """ + if self.bc_type == "scalar": + exclude_neu_dir_nd = sps.kron(sps.eye(self.nd), self._exclude_neu_dir) + + elif self.bc_type == "vectorial": + raise NotImplementedError() + + return exclude_neu_dir_nd * other + + def exclude_rob_dir_nd(self, other): + """ Exclusion of Roben and Dirichlet conditions for vector equations + (elasticity). See above method without _nd suffix for description. + + """ + if self.bc_type == "scalar": + exclude_neu_dir_nd = sps.kron(sps.eye(self.nd), self._exclude_rob_dir) + + elif self.bc_type == "vectorial": + raise NotImplementedError() + + return exclude_neu_dir_nd * other + + def keep_robin_nd(self, other): + """ Keep Roben conditions for vector equations (elasticity). + See above method without _nd suffix for description. + """ + if self.bc_type == "scalar": + exclude_neu_dir_nd = sps.kron(sps.eye(self.nd), self._keep_robin) + + elif self.bc_type == "vectorial": + raise NotImplementedError() + + return exclude_neu_dir_nd * other + def exclude_neu_rob(self, other): """ Mapping to exclude faces/components with Neumann and Robin boundary conditions from local systems. diff --git a/src/porepy/numerics/fv/mpsa.py b/src/porepy/numerics/fv/mpsa.py index f9d11eaed3..ac55796839 100644 --- a/src/porepy/numerics/fv/mpsa.py +++ b/src/porepy/numerics/fv/mpsa.py @@ -973,7 +973,7 @@ def mpsa_partial( return stress_glob, bound_stress_glob, active_faces -def _mpsa_local(g, constit, bound, eta=0, inverter="numba"): +def _mpsa_local(g, constit, bound, eta=0, alpha=None, inverter="numba"): """ Actual implementation of the MPSA W-method. To calculate the MPSA discretization on a grid, either call this method, or, to respect the @@ -1037,6 +1037,11 @@ def _mpsa_local(g, constit, bound, eta=0, inverter="numba"): nd = g.dim + if alpha is None: + if np.sum(bnd.is_rob) !=0: + raise ValueError('If applying Robin conditions you must supply an alpha') + else: + alpha = 1 # Define subcell topology subcell_topology = fvutils.SubcellTopology(g) # Obtain mappings to exclude boundary faces @@ -1045,7 +1050,7 @@ def _mpsa_local(g, constit, bound, eta=0, inverter="numba"): # elasticity and poro-elasticity). hook, igrad, rhs_cells, _, _ = mpsa_elasticity( - g, constit, subcell_topology, bound_exclusion, eta, inverter + g, constit, subcell_topology, bound_exclusion, eta, inverter, alpha=alpha ) hook_igrad = hook * igrad @@ -1073,7 +1078,7 @@ def _mpsa_local(g, constit, bound, eta=0, inverter="numba"): return stress, bound_stress -def mpsa_elasticity(g, constit, subcell_topology, bound_exclusion, eta, inverter): +def mpsa_elasticity(g, constit, subcell_topology, bound_exclusion, eta, inverter, alpha=None): """ This is the function where the real discretization takes place. It contains the parts that are common for elasticity and poro-elasticity, and was thus @@ -1108,7 +1113,7 @@ def mpsa_elasticity(g, constit, subcell_topology, bound_exclusion, eta, inverter nd = g.dim # Compute product between normal vectors and stiffness matrices - ncsym, ncasym, cell_node_blocks, sub_cell_index = _tensor_vector_prod( + ncsym_all, ncasym, cell_node_blocks, sub_cell_index = _tensor_vector_prod( g, constit, subcell_topology ) @@ -1120,20 +1125,20 @@ def mpsa_elasticity(g, constit, subcell_topology, bound_exclusion, eta, inverter hook_normal = sps.coo_matrix( (np.ones(ind_f.size), (np.arange(ind_f.size), ind_f)), shape=(ind_f.size, ind_f.size), - ) * (ncsym + ncasym) + ) * (ncsym_all + ncasym) del ind_f # The final expression of Hook's law will involve deformation gradients # on one side of the faces only; eliminate the other one. # Note that this must be done before we can pair forces from the two # sides of the faces. - hook = __unique_hooks_law(ncsym, ncasym, subcell_topology, nd) + hook = __unique_hooks_law(ncsym_all, ncasym, subcell_topology, nd) del ncasym # Pair the forces from each side - ncsym = subcell_topology.pair_over_subfaces_nd(ncsym) - ncsym = bound_exclusion.exclude_dirichlet_nd(ncsym) + ncsym_all = subcell_topology.pair_over_subfaces_nd(ncsym_all) + ncsym = bound_exclusion.exclude_rob_dir_nd(ncsym_all) num_subfno = subcell_topology.subfno.max() + 1 hook_cell = sps.coo_matrix( @@ -1141,16 +1146,24 @@ def mpsa_elasticity(g, constit, subcell_topology, bound_exclusion, eta, inverter shape=(num_subfno * nd, (np.max(subcell_topology.cno) + 1) * nd), ).tocsr() - hook_cell = bound_exclusion.exclude_dirichlet_nd(hook_cell) + hook_cell = bound_exclusion.exclude_rob_dir_nd(hook_cell) + # For the Robin boundary conditions we need to pair the forces with the + # displacement + ncsym_rob = bound_exclusion.keep_robin_nd(ncsym_all) # Book keeping num_sub_cells = cell_node_blocks[0].size + rob_grad, rob_cell = __get_displacement_submatrices_rob( + g, subcell_topology, eta, num_sub_cells, bound_exclusion, alpha + ) + # Matrices to enforce displacement continuity d_cont_grad, d_cont_cell = __get_displacement_submatrices( g, subcell_topology, eta, num_sub_cells, bound_exclusion ) - grad_eqs = sps.vstack([ncsym, d_cont_grad]) + grad_eqs = sps.vstack([ncsym, ncsym_rob + rob_grad, d_cont_grad]) + del ncsym, d_cont_grad igrad = _inverse_gradient( @@ -1164,7 +1177,7 @@ def mpsa_elasticity(g, constit, subcell_topology, bound_exclusion, eta, inverter ) # Right hand side for cell center variables - rhs_cells = -sps.vstack([hook_cell, d_cont_cell]) + rhs_cells = -sps.vstack([hook_cell, rob_cell, d_cont_cell]) return hook, igrad, rhs_cells, cell_node_blocks, hook_normal @@ -1229,8 +1242,8 @@ def __get_displacement_submatrices( # Expand equations for displacement balance, and eliminate rows # associated with neumann boundary conditions - d_cont_grad = bound_exclusion.exclude_neumann_nd(d_cont_grad) - d_cont_cell = bound_exclusion.exclude_neumann_nd(d_cont_cell) + d_cont_grad = bound_exclusion.exclude_neu_rob_nd(d_cont_grad) + d_cont_cell = bound_exclusion.exclude_neu_rob_nd(d_cont_cell) # The column ordering of the displacement equilibrium equations are # formed as a Kronecker product of scalar equations. Bring them to the @@ -1241,6 +1254,48 @@ def __get_displacement_submatrices( return d_cont_grad, d_cont_cell +def __get_displacement_submatrices_rob( + g, subcell_topology, eta, num_sub_cells, bound_exclusion, alpha +): + nd = g.dim + # Distance from cell centers to face centers, this will be the + # contribution from gradient unknown to equations for displacement + # at the boundary + rob_grad = fvutils.compute_dist_face_cell(g, subcell_topology, eta) + + # For the Robin codition the distance from the cell centers to face centers + # will be the contribution from the gradients. We integrate over the subface + # and multiply by the area + num_nodes = np.diff(g.face_nodes.indptr) + sgn = g.cell_faces[subcell_topology.fno, subcell_topology.cno].A + scaled_sgn = ( + alpha * sgn[0] * g.face_areas[subcell_topology.fno] \ + / num_nodes[subcell_topology.fno] + ) + # pair_over_subfaces flips the sign so we flip it back + rob_grad = sps.kron(sps.eye(nd), rob_grad * sps.diags(scaled_sgn)) + + # Contribution from cell center potentials to local systems + rob_cell = sps.coo_matrix( + (alpha * g.face_areas[subcell_topology.fno] / num_nodes[subcell_topology.fno], + (subcell_topology.subfno, subcell_topology.cno)) + ).tocsr() + rob_cell = sps.kron(sps.eye(nd), rob_cell) + + # Expand equations for displacement balance, and eliminate rows + # associated with neumann boundary conditions + rob_grad = bound_exclusion.keep_robin_nd(rob_grad) + rob_cell = bound_exclusion.keep_robin_nd(rob_cell) + + # The column ordering of the displacement equilibrium equations are + # formed as a Kronecker product of scalar equations. Bring them to the + # same form as that applied in the force balance equations + rob_grad, rob_cell = __rearange_columns_displacement_eqs( + rob_grad, rob_cell, num_sub_cells, nd + ) + + return rob_grad, rob_cell + def _split_stiffness_matrix(constit): """ @@ -1496,15 +1551,20 @@ def _block_diagonal_structure( # get block-structure. First eliminate node numbers at the boundary, where # the equations are either of flux or pressure continuity (not both) - nno_stress = bound_exclusion.exclude_dirichlet(nno) - nno_displacement = bound_exclusion.exclude_neumann(nno) + nno_stress = bound_exclusion.exclude_rob_dir(nno) + nno_displacement = bound_exclusion.exclude_neu_rob(nno) + # we have now eliminated all nodes related to robin, we therefore add them + nno_rob = bound_exclusion.keep_robin(nno) if bound_exclusion.bc_type == "scalar": - node_occ = np.hstack((np.tile(nno_stress, nd), np.tile(nno_displacement, nd))) + node_occ = np.hstack((np.tile(nno_stress, nd), + np.tile(nno_rob, nd), + np.tile(nno_displacement, nd))) elif bound_exclusion.bc_type == "vectorial": node_occ = np.hstack((nno_stress, nno_displacement)) - + if nno_rob.size >=0: + raise NotImplementedError() sorted_ind = np.argsort(node_occ, kind="mergesort") rows2blk_diag = sps.coo_matrix( (np.ones(sorted_ind.size), (np.arange(sorted_ind.size), sorted_ind)) @@ -1545,16 +1605,20 @@ def create_bound_rhs(bound, bound_exclusion, subcell_topology, g): basis functions for boundary values """ nd = g.dim - num_stress = bound_exclusion.exclude_dir.shape[0] * nd - num_displ = bound_exclusion.exclude_neu.shape[0] * nd + num_stress = bound_exclusion._exclude_rob_dir.shape[0] * nd + num_displ = bound_exclusion._exclude_neu_rob.shape[0] * nd + num_rob = bound_exclusion._keep_robin.shape[0] * nd + fno = subcell_topology.fno_unique subfno = subcell_topology.subfno_unique sgn = g.cell_faces[ subcell_topology.fno_unique, subcell_topology.cno_unique ].A.ravel("F") + num_neu = sum(bound.is_neu[fno]) * nd num_dir = sum(bound.is_dir[fno]) * nd - num_bound = num_neu + num_dir + assert num_rob == sum(bound.is_rob[fno]) * nd + num_bound = num_neu + num_dir + num_rob # Convenience method for duplicating a list, with a certain increment def expand_ind(ind, dim, increment): @@ -1568,18 +1632,31 @@ def expand_ind(ind, dim, increment): # Define right hand side for Neumann boundary conditions # First row indices in rhs matrix - is_neu = bound_exclusion.exclude_dirichlet(bound.is_neu[fno].astype("int64")) + is_neu = bound_exclusion.exclude_rob_dir(bound.is_neu[fno].astype("int64")) neu_ind_single = np.argwhere(is_neu).ravel("F") + is_rob = bound_exclusion.keep_robin(bound.is_rob[fno].astype("int64")) + rob_ind_single = np.argwhere(is_rob).ravel("F") - # There are is_neu.size Neumann conditions per dimension + # There are is_neu.size Neumann conditions per dimension and + # there are is_rob.size Robin conditions per dimension neu_ind = expand_ind(neu_ind_single, nd, is_neu.size) - + rob_ind = expand_ind(rob_ind_single, nd, is_rob.size) + # We now merge the neuman and robin indices since they are treated equivalent + if rob_ind.size==0: + neu_rob_ind = neu_ind + elif neu_ind.size==0: + neu_rob_ind = rob_ind + num_stress + else: + neu_rob_ind = np.hstack((neu_ind, rob_ind + num_stress)) + # We also need to account for all half faces, that is, do not exclude # Dirichlet and Neumann boundaries. neu_ind_single_all = np.argwhere(bound.is_neu[fno].astype("int")).ravel("F") + rob_ind_single_all = np.argwhere(bound.is_rob[fno].astype("int")).ravel("F") dir_ind_single_all = np.argwhere(bound.is_dir[fno].astype("int")).ravel("F") - neu_ind_all = np.tile(neu_ind_single_all, nd) + neu_rob_ind_single_all = np.hstack((neu_ind_single_all, rob_ind_single_all)) + neu_rob_ind_all = np.tile(neu_rob_ind_single_all, nd) # Some care is needed to compute coefficients in Neumann matrix: sgn is # already defined according to the subcell topology [fno], while areas @@ -1593,21 +1670,21 @@ def expand_ind(ind, dim, increment): # have to do # so, and we will flip the sign later. This means that a stress [1,1] on a # boundary face pushes(or pulls) the face to the top right corner. - neu_val = 1 / num_face_nodes[fno_ext[neu_ind_all]] - # The columns will be 0:neu_ind.size - if neu_ind.size > 0: + neu_val = 1 / num_face_nodes[fno_ext[neu_rob_ind_all]] + # The columns will be 0:neu_rob_ind.size + if neu_rob_ind.size > 0: neu_cell = sps.coo_matrix( - (neu_val.ravel("F"), (neu_ind, np.arange(neu_ind.size))), - shape=(num_stress, num_bound), + (neu_val.ravel("F"), (neu_rob_ind, np.arange(neu_rob_ind.size))), + shape=(num_stress+num_rob, num_bound), ).tocsr() else: # Special handling when no elements are found. Not sure if this is # necessary, or if it is me being stupid - neu_cell = sps.coo_matrix((num_stress, num_bound)).tocsr() + neu_cell = sps.coo_matrix((num_stress + num_rob, num_bound)).tocsr() # Dirichlet boundary conditions, procedure is similar to that for Neumann - is_dir = bound_exclusion.exclude_neumann(bound.is_dir[fno].astype("int64")) + is_dir = bound_exclusion.exclude_neu_rob(bound.is_dir[fno].astype("int64")) dir_ind_single = np.argwhere(is_dir).ravel("F") dir_ind = expand_ind(dir_ind_single, nd, is_dir.size) @@ -1622,7 +1699,7 @@ def expand_ind(ind, dim, increment): # for each face, then the y-component, etc. if dir_ind.size > 0: dir_cell = sps.coo_matrix( - (dir_val, (dir_ind, num_neu + np.arange(dir_ind.size))), + (dir_val, (dir_ind, num_neu + num_rob + np.arange(dir_ind.size))), shape=(num_displ, num_bound), ).tocsr() else: @@ -1635,7 +1712,7 @@ def expand_ind(ind, dim, increment): # The columns in neu_cell, dir_cell are ordered from 0 to num_bound-1. # Map these to all half-face indices - is_bnd = np.hstack((neu_ind_single_all, dir_ind_single_all)) + is_bnd = np.hstack((neu_rob_ind_single_all, dir_ind_single_all)) bnd_ind = fvutils.expand_indices_nd(is_bnd, nd) bnd_2_all_hf = sps.coo_matrix( diff --git a/test/unit/test_mpfa_robin.py b/test/unit/test_mpfa_robin.py index a3adffdf7a..909ff06c65 100644 --- a/test/unit/test_mpfa_robin.py +++ b/test/unit/test_mpfa_robin.py @@ -126,7 +126,6 @@ def test_no_neumann(self): assert np.allclose(p, p_ex) assert np.allclose(flux * p + bound_flux * u_bound, u_ex) - def solve_robin(self, g, k, bnd, alpha, p_bound, rob_bound, dir_ind, rob_ind, p_ex, u_ex): flux, bound_flux, _, _ = pp.numerics.fv.mpfa._mpfa_local(g, k, bnd, alpha=alpha) diff --git a/test/unit/test_mpsa_robin.py b/test/unit/test_mpsa_robin.py new file mode 100644 index 0000000000..24fb07efe1 --- /dev/null +++ b/test/unit/test_mpsa_robin.py @@ -0,0 +1,110 @@ +import numpy as np +import scipy.sparse as sps +import unittest + +import porepy as pp + +class RobinBoundTest(unittest.TestCase): + def test_dir_rob(self): + nx = 2 + ny = 2 + g = pp.CartGrid([nx, ny], physdims=[1, 1]) + g.compute_geometry() + c = pp.FourthOrderTensor(2, np.ones(g.num_cells), np.ones(g.num_cells)) + alpha = np.pi + + bot = g.face_centers[1] < 1e-10 + top = g.face_centers[1] > 1 - 1e-10 + left = g.face_centers[0] < 1e-10 + right = g.face_centers[0] > 1 - 1e-10 + + dir_ind = np.ravel(np.argwhere(left + bot + top)) + neu_ind = np.ravel(np.argwhere([])) + rob_ind = np.ravel(np.argwhere(right)) + + names = ['dir']*len(dir_ind) + ['rob'] * len(rob_ind) + bnd_ind = np.hstack((dir_ind, rob_ind)) + bnd = pp.BoundaryCondition(g, bnd_ind, names) + + def u_ex(x): + return np.vstack((x[0], 0 * x[1])) + + def T_ex(faces): + if np.size(faces) == 0: + return np.atleast_2d(np.array([])) + sigma = np.array([[3, 0], [0, 1]]) + T_r = [np.dot(sigma, g.face_normals[:2, f]) for f in faces] + return np.vstack(T_r).T + + u_bound = np.zeros((2, g.num_faces)) + + sgn_n = pp.numerics.fracture_deformation.sign_of_faces(g, neu_ind) + sgn_r = pp.numerics.fracture_deformation.sign_of_faces(g, rob_ind) + + u_bound[:, dir_ind] = u_ex(g.face_centers[:, dir_ind]) + u_bound[:, neu_ind] = T_ex(neu_ind) * sgn_n + u_bound[:, rob_ind] = ( + T_ex(rob_ind) * sgn_r + + alpha * u_ex(g.face_centers[:, rob_ind]) * g.face_areas[rob_ind] + ) + u, T = self.solve_mpsa(g, c, alpha, bnd, u_bound) + + assert np.allclose(u, u_ex(g.cell_centers).ravel('F')) + assert np.allclose(T, T_ex(np.arange(g.num_faces)).ravel('F')) + + def test_dir_neu_rob(self): + nx = 2 + ny = 3 + g = pp.CartGrid([nx, ny], physdims=[1, 1]) + g.compute_geometry() + c = pp.FourthOrderTensor(2, np.ones(g.num_cells), np.ones(g.num_cells)) + alpha = np.pi + + bot = g.face_centers[1] < 1e-10 + top = g.face_centers[1] > 1 - 1e-10 + left = g.face_centers[0] < 1e-10 + right = g.face_centers[0] > 1 - 1e-10 + + dir_ind = np.ravel(np.argwhere(left)) + neu_ind = np.ravel(np.argwhere(top)) + rob_ind = np.ravel(np.argwhere(right + bot)) + + names = ['dir']*len(dir_ind) + ['rob'] * len(rob_ind) + bnd_ind = np.hstack((dir_ind, rob_ind)) + bnd = pp.BoundaryCondition(g, bnd_ind, names) + + def u_ex(x): + return np.vstack((x[0], 0 * x[1])) + + def T_ex(faces): + if np.size(faces) == 0: + return np.atleast_2d(np.array([])) + sigma = np.array([[3, 0], [0, 1]]) + T_r = [np.dot(sigma, g.face_normals[:2, f]) for f in faces] + return np.vstack(T_r).T + + u_bound = np.zeros((2, g.num_faces)) + + sgn_n = pp.numerics.fracture_deformation.sign_of_faces(g, neu_ind) + sgn_r = pp.numerics.fracture_deformation.sign_of_faces(g, rob_ind) + + u_bound[:, dir_ind] = u_ex(g.face_centers[:, dir_ind]) + u_bound[:, neu_ind] = T_ex(neu_ind) * sgn_n + u_bound[:, rob_ind] = ( + T_ex(rob_ind) * sgn_r + + alpha * u_ex(g.face_centers[:, rob_ind]) * g.face_areas[rob_ind] + ) + u, T = self.solve_mpsa(g, c, alpha, bnd, u_bound) + + assert np.allclose(u, u_ex(g.cell_centers).ravel('F')) + assert np.allclose(T, T_ex(np.arange(g.num_faces)).ravel('F')) + + def solve_mpsa(self, g, c, alpha, bnd, u_bound): + stress, bound_stress = pp.numerics.fv.mpsa._mpsa_local(g, c, bnd, alpha=alpha) + div = pp.fvutils.vector_divergence(g) + a = div * stress + b = -div * bound_stress * u_bound.ravel('F') + + u = np.linalg.solve(a.A, b) + T = stress * u + bound_stress * u_bound.ravel('F') + return u, T From 3fab5c8f1d5f71e612226f6ade993d1442523d45 Mon Sep 17 00:00:00 2001 From: Runar Date: Tue, 11 Sep 2018 10:38:33 +0200 Subject: [PATCH 08/26] Fixed bug in mpsa where asymetric part was dropped for Neumann boundaries --- src/porepy/numerics/fv/mpsa.py | 32 ++++++++++++++++---------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/src/porepy/numerics/fv/mpsa.py b/src/porepy/numerics/fv/mpsa.py index ac55796839..9e04b6d95c 100644 --- a/src/porepy/numerics/fv/mpsa.py +++ b/src/porepy/numerics/fv/mpsa.py @@ -703,7 +703,7 @@ def given_slip_distance(self, g, stress, bound_stress, faces=None): # ------------------------------------------------------------------------------# -def mpsa(g, constit, bound, eta=None, inverter=None, max_memory=None, **kwargs): +def mpsa(g, constit, bound, eta=None, alpha=None, inverter=None, max_memory=None, **kwargs): """ Discretize the vector elliptic equation by the multi-point stress approximation method, specifically the weakly symmetric MPSA-W method. @@ -799,7 +799,7 @@ def mpsa(g, constit, bound, eta=None, inverter=None, max_memory=None, **kwargs): # TODO: We may want to estimate the memory need, and give a warning if # this seems excessive stress, bound_stress = _mpsa_local( - g, constit, bound, eta=eta, inverter=inverter + g, constit, bound, eta=eta, inverter=inverter, alpha=alpha ) else: # Estimate number of partitions necessary based on prescribed memory @@ -835,7 +835,8 @@ def mpsa(g, constit, bound, eta=None, inverter=None, max_memory=None, **kwargs): # Perform local discretization. loc_stress, loc_bound_stress, loc_faces = mpsa_partial( - g, constit, bound, eta=eta, inverter=inverter, nodes=active_nodes + g, constit, bound, eta=eta, inverter=inverter, nodes=active_nodes, + alpha=alpha ) # Eliminate contribution from faces already covered @@ -861,6 +862,7 @@ def mpsa_partial( faces=None, nodes=None, apertures=None, + alpha=None ): """ Run an MPFA discretization on subgrid, and return discretization in terms @@ -952,7 +954,7 @@ def mpsa_partial( # Discretization of sub-problem stress_loc, bound_stress_loc = _mpsa_local( - sub_g, loc_c, loc_bnd, eta=eta, inverter=inverter + sub_g, loc_c, loc_bnd, eta=eta, inverter=inverter, alpha=alpha ) face_map, cell_map = fvutils.map_subgrid_to_grid( @@ -1038,7 +1040,7 @@ def _mpsa_local(g, constit, bound, eta=0, alpha=None, inverter="numba"): nd = g.dim if alpha is None: - if np.sum(bnd.is_rob) !=0: + if np.sum(bound.is_rob) !=0: raise ValueError('If applying Robin conditions you must supply an alpha') else: alpha = 1 @@ -1134,10 +1136,10 @@ def mpsa_elasticity(g, constit, subcell_topology, bound_exclusion, eta, inverter # sides of the faces. hook = __unique_hooks_law(ncsym_all, ncasym, subcell_topology, nd) + # Pair the forces from each side + ncsym_all = subcell_topology.pair_over_subfaces_nd(ncsym_all + ncasym) del ncasym - # Pair the forces from each side - ncsym_all = subcell_topology.pair_over_subfaces_nd(ncsym_all) ncsym = bound_exclusion.exclude_rob_dir_nd(ncsym_all) num_subfno = subcell_topology.subfno.max() + 1 @@ -1155,7 +1157,7 @@ def mpsa_elasticity(g, constit, subcell_topology, bound_exclusion, eta, inverter num_sub_cells = cell_node_blocks[0].size rob_grad, rob_cell = __get_displacement_submatrices_rob( g, subcell_topology, eta, num_sub_cells, bound_exclusion, alpha - ) + ) # Matrices to enforce displacement continuity d_cont_grad, d_cont_cell = __get_displacement_submatrices( @@ -1163,8 +1165,7 @@ def mpsa_elasticity(g, constit, subcell_topology, bound_exclusion, eta, inverter ) grad_eqs = sps.vstack([ncsym, ncsym_rob + rob_grad, d_cont_grad]) - - del ncsym, d_cont_grad + del ncsym, d_cont_grad, ncsym_rob, rob_grad igrad = _inverse_gradient( grad_eqs, @@ -1262,19 +1263,18 @@ def __get_displacement_submatrices_rob( # contribution from gradient unknown to equations for displacement # at the boundary rob_grad = fvutils.compute_dist_face_cell(g, subcell_topology, eta) - + # For the Robin codition the distance from the cell centers to face centers # will be the contribution from the gradients. We integrate over the subface # and multiply by the area num_nodes = np.diff(g.face_nodes.indptr) - sgn = g.cell_faces[subcell_topology.fno, subcell_topology.cno].A + sgn = g.cell_faces[subcell_topology.fno_unique, subcell_topology.cno_unique].A scaled_sgn = ( - alpha * sgn[0] * g.face_areas[subcell_topology.fno] \ - / num_nodes[subcell_topology.fno] + alpha * sgn[0] * g.face_areas[subcell_topology.fno_unique] \ + / num_nodes[subcell_topology.fno_unique] ) # pair_over_subfaces flips the sign so we flip it back - rob_grad = sps.kron(sps.eye(nd), rob_grad * sps.diags(scaled_sgn)) - + rob_grad = sps.kron(sps.eye(nd), sps.diags(scaled_sgn)*rob_grad) # Contribution from cell center potentials to local systems rob_cell = sps.coo_matrix( (alpha * g.face_areas[subcell_topology.fno] / num_nodes[subcell_topology.fno], From ca106002d3114cfec0380bf632d0b957be8dfd14 Mon Sep 17 00:00:00 2001 From: Runar Date: Tue, 11 Sep 2018 10:39:09 +0200 Subject: [PATCH 09/26] Added unit-tests for Robin conditions in mpsa --- test/unit/test_mpsa_robin.py | 95 ++++++++++++++++++++++++++++++++++++ 1 file changed, 95 insertions(+) diff --git a/test/unit/test_mpsa_robin.py b/test/unit/test_mpsa_robin.py index 24fb07efe1..9cacdd9fb1 100644 --- a/test/unit/test_mpsa_robin.py +++ b/test/unit/test_mpsa_robin.py @@ -99,6 +99,101 @@ def T_ex(faces): assert np.allclose(u, u_ex(g.cell_centers).ravel('F')) assert np.allclose(T, T_ex(np.arange(g.num_faces)).ravel('F')) + def test_structured_triang(self): + nx = 1 + ny = 1 + g = pp.StructuredTriangleGrid([nx, ny], physdims=[1, 1]) + g.compute_geometry() + c = pp.FourthOrderTensor(2, np.ones(g.num_cells), np.ones(g.num_cells)) + alpha = np.pi + + bot = g.face_centers[1] < 1e-10 + top = g.face_centers[1] > 1 - 1e-10 + left = g.face_centers[0] < 1e-10 + right = g.face_centers[0] > 1 - 1e-10 + + dir_ind = np.ravel(np.argwhere(left)) + neu_ind = np.ravel(np.argwhere(top)) + rob_ind = np.ravel(np.argwhere(right + bot)) + + names = ['dir']*len(dir_ind) + ['rob'] * len(rob_ind) + bnd_ind = np.hstack((dir_ind, rob_ind)) + bnd = pp.BoundaryCondition(g, bnd_ind, names) + + def u_ex(x): + return np.vstack((x[1], x[0])) + + def T_ex(faces): + if np.size(faces) == 0: + return np.atleast_2d(np.array([])) + sigma = np.array([[0, 2], [2, 0]]) + T_r = [np.dot(sigma, g.face_normals[:2, f]) for f in faces] + return np.vstack(T_r).T + + u_bound = np.zeros((2, g.num_faces)) + + sgn_n = pp.numerics.fracture_deformation.sign_of_faces(g, neu_ind) + sgn_r = pp.numerics.fracture_deformation.sign_of_faces(g, rob_ind) + + u_bound[:, dir_ind] = u_ex(g.face_centers[:, dir_ind]) + u_bound[:, neu_ind] = T_ex(neu_ind) * sgn_n + u_bound[:, rob_ind] = ( + T_ex(rob_ind) * sgn_r + + alpha * u_ex(g.face_centers[:, rob_ind]) * g.face_areas[rob_ind] + ) + u, T = self.solve_mpsa(g, c, alpha, bnd, u_bound) + + assert np.allclose(u, u_ex(g.cell_centers).ravel('F')) + assert np.allclose(T, T_ex(np.arange(g.num_faces)).ravel('F')) + + def test_unstruct_triang(self): + corners = np.array([[0, 0, 1, 1], [0, 1, 1, 0]]) + points = np.random.rand(2, 2) + points = np.hstack((corners, points)) + g = pp.TriangleGrid(points) + g.compute_geometry() + c = pp.FourthOrderTensor(2, np.ones(g.num_cells), np.ones(g.num_cells)) + alpha = np.pi + + bot = g.face_centers[1] < 1e-10 + top = g.face_centers[1] > 1 - 1e-10 + left = g.face_centers[0] < 1e-10 + right = g.face_centers[0] > 1 - 1e-10 + + dir_ind = np.ravel(np.argwhere(left)) + neu_ind = np.ravel(np.argwhere(top)) + rob_ind = np.ravel(np.argwhere(right + bot)) + + names = ['dir']*len(dir_ind) + ['rob'] * len(rob_ind) + bnd_ind = np.hstack((dir_ind, rob_ind)) + bnd = pp.BoundaryCondition(g, bnd_ind, names) + + def u_ex(x): + return np.vstack((x[1], x[0])) + + def T_ex(faces): + if np.size(faces) == 0: + return np.atleast_2d(np.array([])) + sigma = np.array([[0, 2], [2, 0]]) + T_r = [np.dot(sigma, g.face_normals[:2, f]) for f in faces] + return np.vstack(T_r).T + + u_bound = np.zeros((2, g.num_faces)) + + sgn_n = pp.numerics.fracture_deformation.sign_of_faces(g, neu_ind) + sgn_r = pp.numerics.fracture_deformation.sign_of_faces(g, rob_ind) + + u_bound[:, dir_ind] = u_ex(g.face_centers[:, dir_ind]) + u_bound[:, neu_ind] = T_ex(neu_ind) * sgn_n + u_bound[:, rob_ind] = ( + T_ex(rob_ind) * sgn_r + + alpha * u_ex(g.face_centers[:, rob_ind]) * g.face_areas[rob_ind] + ) + u, T = self.solve_mpsa(g, c, alpha, bnd, u_bound) + + assert np.allclose(u, u_ex(g.cell_centers).ravel('F')) + assert np.allclose(T, T_ex(np.arange(g.num_faces)).ravel('F')) + def solve_mpsa(self, g, c, alpha, bnd, u_bound): stress, bound_stress = pp.numerics.fv.mpsa._mpsa_local(g, c, bnd, alpha=alpha) div = pp.fvutils.vector_divergence(g) From 8ddfde7fb6c271501271490ffa641555374dbfbf Mon Sep 17 00:00:00 2001 From: Runar Date: Tue, 11 Sep 2018 11:00:43 +0200 Subject: [PATCH 10/26] Added convergence test for MPSA with Robin --- examples/example2/ConvergenceTestMPSA.ipynb | 293 +++++++++++++++++--- 1 file changed, 252 insertions(+), 41 deletions(-) diff --git a/examples/example2/ConvergenceTestMPSA.ipynb b/examples/example2/ConvergenceTestMPSA.ipynb index 000d283961..d60ca47022 100644 --- a/examples/example2/ConvergenceTestMPSA.ipynb +++ b/examples/example2/ConvergenceTestMPSA.ipynb @@ -216,11 +216,18 @@ "output_type": "stream", "text": [ "Cartesian errors - no perturbations\n", - "[0.19341322457532645, 0.05398606173752483, 0.019785099166683015] displacement error\n", - "[0.03080448569971605, 0.022727907422006446, 0.017633923832997576] stress error\n", + "[0.18471065947993767, 0.04467625698150768, 0.01073709992570399] displacement error\n", + "[4.13442558 4.16092402] convergence rate u^n-1 / u^n\n", + "\n", + "[0.06476977708690358, 0.021191370180076885, 0.006402282045426171] stress error\n", + "[3.05642233 3.30997136] convergence rate f^n-1 / f^n\n", + "\n", "Triangular errors - no perturbations\n", - "[0.15803686415889875, 0.04310823056510322, 0.015475821210292553] displacement error\n", - "[0.04467366719714846, 0.02540025474312989, 0.016373893423034032] stress error\n" + "[0.16176733322360393, 0.03341692740958871, 0.007161359044861829] displacement error\n", + "[4.84087993 4.66628292] convergence rate u^n-1 / u^n\n", + "\n", + "[0.06486596857883836, 0.02141863302701911, 0.006979068019066017] stress error\n", + "[3.02848312 3.06898184] convergence rate f^n-1 / f^n\n" ] } ], @@ -297,8 +304,8 @@ " left = np.ravel(np.argwhere(g.face_centers[0, :] < 1e-10))\n", " right = np.ravel(np.argwhere(g.face_centers[0, :] > n - 1e-10))\n", "\n", - " dir_faces = np.hstack((bot, top))\n", - " neu_faces = np.hstack((left, right))\n", + " dir_faces = np.hstack((left, bot, top))\n", + " neu_faces = np.hstack((right))\n", " bound_cond = bc.BoundaryCondition(g, dir_faces, ['dir'] * dir_faces.size)\n", " \n", " # MPFA discretization, and system matrix\n", @@ -350,22 +357,226 @@ " return u_err, flux_err\n", "\n", "u_cart_nopert, f_cart_nopert = run_convergence('cart', 0)\n", + "u_cart_rate = np.array(u_cart_nopert)\n", + "u_cart_rate = u_cart_rate[:-1] / u_cart_rate[1:]\n", + "f_cart_rate = np.array(f_cart_nopert)\n", + "f_cart_rate = f_cart_rate[:-1] / f_cart_rate[1:]\n", "print('Cartesian errors - no perturbations')\n", "print(u_cart_nopert,'displacement error')\n", + "print(u_cart_rate, 'convergence rate u^n-1 / u^n\\n')\n", "print(f_cart_nopert,'stress error')\n", + "print(f_cart_rate, 'convergence rate f^n-1 / f^n\\n')\n", "\n", "u_triang_nopert, f_triang_nopert = run_convergence('triangular', 0)\n", + "u_triang_rate = np.array(u_triang_nopert)\n", + "u_triang_rate = u_triang_rate[:-1] / u_triang_rate[1:]\n", + "f_triang_rate = np.array(f_triang_nopert)\n", + "f_triang_rate = f_triang_rate[:-1] / f_triang_rate[1:]\n", + "\n", "print('Triangular errors - no perturbations')\n", "print(u_triang_nopert,'displacement error')\n", + "print(u_triang_rate, 'convergence rate u^n-1 / u^n\\n')\n", "print(f_triang_nopert,'stress error')\n", + "print(f_triang_rate, 'convergence rate f^n-1 / f^n')\n", + "\n", + "# RB: These values were hard-coded 11.09.2018.\n", + "assert np.abs(u_cart_nopert[2] - 0.01073709992570399) < 1e-10\n", + "assert np.abs(f_cart_nopert[2] - 0.006402282045426171) < 1e-10\n", + "assert np.abs(u_triang_nopert[2] - 0.007161359044861829) < 1e-10\n", + "assert np.abs(f_triang_nopert[2] - 0.006979068019066017) < 1e-10" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Robin Condition" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Cartesian errors - no perturbations\n", + "[0.13785825178140515, 0.033672908859893956, 0.0082926160732301] displacement error\n", + "[4.09404047 4.06058939] convergence rate u^n-1 / u^n\n", + "\n", + "[0.04932864932955223, 0.01582468790884128, 0.004760957079948613] stress error\n", + "[3.11719571 3.32384595] convergence rate f^n-1 / f^n\n", + "\n", + "Triangular errors - no perturbations\n", + "[0.11226477176063937, 0.025455746327683503, 0.006007054619672307] displacement error\n", + "[4.41019369 4.2376419 ] convergence rate u^n-1 / u^n\n", + "\n", + "[0.044251086212033885, 0.016626084791879308, 0.0058474670379793795] stress error\n", + "[2.6615458 2.84329688] convergence rate f^n-1 / f^n\n" + ] + } + ], + "source": [ + "# Permeability tensor, scalar for simplicity\n", + "mu = 1\n", + "lmbda = 1\n", "\n", - "# # EK: These values were hard-coded 12.06.2017.\n", - "assert np.abs(u_cart_nopert[2] - 0.01978509916668196) < 1e-10\n", - "assert np.abs(f_cart_nopert[2] - 0.017633923832997497) < 1e-10\n", + "###\n", + "# This is where parameters can be modified to alter the convergence test.\n", + "# The remaining lines should \n", + "np.random.seed(42)\n", + "base = 3\n", + "domain = np.array([1, 1])\n", + "basedim = np.array([base, base])\n", + "num_refs = 3\n", + "grid_type = 'cart'\n", "\n", - "# RB: These values were hard-coded 06.09.2018.\n", - "assert np.abs(u_triang_nopert[2] - 0.015475821210292553) < 1e-10\n", - "assert np.abs(f_triang_nopert[2] - 0.016373893423034032) < 1e-10" + "# Analytical solution\n", + "x, y = sympy.symbols('x y')\n", + "ux = (0.5-y)**2 * sympy.sin(y) * sympy.cos(y)\n", + "uy = y * (y-1) * sympy.cosh(x)\n", + "ux_f = sympy.lambdify((x, y), ux, 'numpy')\n", + "uy_f = sympy.lambdify((x, y), uy, 'numpy')\n", + "dux_x = sympy.diff(ux, x)\n", + "dux_y = sympy.diff(ux, y)\n", + "duy_x = sympy.diff(uy, x)\n", + "duy_y = sympy.diff(uy, y)\n", + "divu = dux_x + duy_y\n", + "\n", + "sxx = 2 * mu * dux_x + lmbda * divu\n", + "sxy = mu * (dux_y + duy_x)\n", + "syx = mu * (duy_x + dux_y)\n", + "syy = 2 * mu * duy_y + lmbda * divu\n", + "\n", + "sxx_f = sympy.lambdify((x, y), sxx, 'numpy')\n", + "sxy_f = sympy.lambdify((x, y), sxy, 'numpy')\n", + "syx_f = sympy.lambdify((x, y), syx, 'numpy')\n", + "syy_f = sympy.lambdify((x, y), syy, 'numpy')\n", + "\n", + "rhs_x = sympy.diff(sxx, x) + sympy.diff(syx, y)\n", + "rhs_y = sympy.diff(sxy, x) + sympy.diff(syy, y)\n", + "rhs_x_f = sympy.lambdify((x, y), rhs_x, 'numpy')\n", + "rhs_y_f = sympy.lambdify((x, y), rhs_y, 'numpy')\n", + "\n", + "def run_convergence(grid_type, pert):\n", + " u_err = []\n", + " flux_err = []\n", + "\n", + " for g in setup_grids.grid_sequence(basedim, num_refs, grid_type, pert, domain=domain):\n", + " # Reset the random seed for every grid realization.\n", + " # This should make no difference for the convergence test, \n", + " # but it makes sure that we can run unit tests based on the values obtained\n", + " # here.\n", + " np.random.seed(42)\n", + " \n", + " xf = g.face_centers\n", + " sx_ex_faces = np.vstack((sxx_f(xf[0], xf[1]), sxy_f(xf[0], xf[1])))\n", + " sy_ex_faces = np.vstack((syx_f(xf[0], xf[1]), syy_f(xf[0], xf[1])))\n", + " \n", + " stress_x_ex = np.sum(g.face_normals[:2] * sx_ex_faces, axis=0)\n", + " stress_y_ex = np.sum(g.face_normals[:2] * sy_ex_faces, axis=0)\n", + "\n", + " # Permeability tensor\n", + " mu_c = mu * np.ones(g.num_cells)\n", + " lmbda_c = lmbda * np.ones(g.num_cells)\n", + " k = stiffness(2, mu_c, lmbda_c)\n", + " alpha = 1.0\n", + " # Set type of boundary conditions - Dirichlet\n", + " n = g.nodes.max()\n", + " # Set type of boundary conditions - Dirichlet\n", + " top = np.ravel(np.argwhere(g.face_centers[1, :] > n - 1e-10))\n", + " bot = np.ravel(np.argwhere(g.face_centers[1, :] < 1e-10))\n", + " left = np.ravel(np.argwhere(g.face_centers[0, :] < 1e-10))\n", + " right = np.ravel(np.argwhere(g.face_centers[0, :] > n - 1e-10))\n", + "\n", + " dir_faces = np.hstack((top, left))\n", + " rob_faces = np.hstack((right, bot))\n", + " bc_faces = np.hstack([dir_faces, rob_faces])\n", + " bc_names = ['dir']*len(dir_faces) + ['rob'] * len(rob_faces)\n", + " bound_cond = bc.BoundaryCondition(g, bc_faces, bc_names)\n", + " \n", + " # MPFA discretization, and system matrix\n", + " stress, bound_stress = mpsa.mpsa(g, k, bound_cond,alpha=alpha)\n", + " div = fvutils.vector_divergence(g)\n", + " a = div * stress\n", + " xf = g.face_centers\n", + " \n", + " rfi, _, sgn = sps.find(g.cell_faces[rob_faces,:])\n", + " \n", + " u_bound = np.zeros((g.dim, g.num_faces))\n", + " u_bound[0, dir_faces] = ux_f(xf[0, dir_faces], xf[1, dir_faces])\n", + " u_bound[1, dir_faces] = uy_f(xf[0, dir_faces], xf[1, dir_faces])\n", + " \n", + " # When setting the Neumann boundary condition we need to flip the \n", + " # sign of the faces with a normal pointing inwards as it is assumed\n", + " # that the Neumann condition is the force from the boundary on to face\n", + " u_bound[0, rob_faces[rfi]] = stress_x_ex[rob_faces[rfi]] * (sgn) + \\\n", + " alpha * ux_f(xf[0, rob_faces[rfi]], xf[1, rob_faces[rfi]]) * g.face_areas[rob_faces[rfi]]\n", + "\n", + " u_bound[1, rob_faces[rfi]] = stress_y_ex[rob_faces[rfi]] * (sgn) +\\\n", + " alpha * uy_f(xf[0, rob_faces[rfi]], xf[1, rob_faces[rfi]]) * g.face_areas[rob_faces[rfi]]\n", + " # Right hand side - contribution from the solution and the boundary conditions\n", + " xc = g.cell_centers\n", + " rhs = np.vstack((rhs_x_f(xc[0], xc[1]), rhs_y_f(xc[0], xc[1]))) * g.cell_volumes\n", + " b = rhs.ravel('F') - div * bound_stress * u_bound.ravel('F')\n", + "\n", + " # Solve system, derive fluxes\n", + " u_num = spsolve(a, b)\n", + " stress_num = stress * u_num + bound_stress * u_bound.ravel('F')\n", + "\n", + " ux_num = u_num[::2]\n", + " uy_num = u_num[1::2]\n", + " \n", + " stress_x_num = stress_num[::2]\n", + " stress_y_num = stress_num[1::2]\n", + " # Exact solution\n", + " ux_ex = ux_f(xc[0], xc[1])\n", + " uy_ex = uy_f(xc[0], xc[1])\n", + " u_ex = np.vstack((ux_ex, uy_ex))\n", + " u_diff = np.vstack((ux_num - ux_ex, uy_num - uy_ex))\n", + " \n", + " stress_diff = np.vstack((stress_x_num - stress_x_ex, \n", + " stress_y_num - stress_y_ex))\n", + " stress_ex = np.vstack((stress_x_ex, stress_y_ex))\n", + " \n", + " u_err.append(np.sqrt(np.sum(g.cell_volumes * u_diff**2)) /\n", + " np.sqrt(np.sum(g.cell_volumes * u_ex**2)))\n", + " flux_err.append(np.sqrt(np.sum((g.face_areas ** g.dim) * stress_diff**2))/\n", + " np.sqrt(np.sum((g.face_areas ** g.dim) * stress_ex**2)))\n", + " \n", + " \n", + " return u_err, flux_err\n", + "\n", + "u_cart_nopert, f_cart_nopert = run_convergence('cart', 0)\n", + "u_cart_rate = np.array(u_cart_nopert)\n", + "u_cart_rate = u_cart_rate[:-1] / u_cart_rate[1:]\n", + "f_cart_rate = np.array(f_cart_nopert)\n", + "f_cart_rate = f_cart_rate[:-1] / f_cart_rate[1:]\n", + "print('Cartesian errors - no perturbations')\n", + "print(u_cart_nopert,'displacement error')\n", + "print(u_cart_rate, 'convergence rate u^n-1 / u^n\\n')\n", + "print(f_cart_nopert,'stress error')\n", + "print(f_cart_rate, 'convergence rate f^n-1 / f^n\\n')\n", + "\n", + "u_triang_nopert, f_triang_nopert = run_convergence('triangular', 0)\n", + "u_triang_rate = np.array(u_triang_nopert)\n", + "u_triang_rate = u_triang_rate[:-1] / u_triang_rate[1:]\n", + "f_triang_rate = np.array(f_triang_nopert)\n", + "f_triang_rate = f_triang_rate[:-1] / f_triang_rate[1:]\n", + "\n", + "print('Triangular errors - no perturbations')\n", + "print(u_triang_nopert,'displacement error')\n", + "print(u_triang_rate, 'convergence rate u^n-1 / u^n\\n')\n", + "print(f_triang_nopert,'stress error')\n", + "print(f_triang_rate, 'convergence rate f^n-1 / f^n')\n", + "\n", + "# RB: These values were hard-coded 11.09.2018.\n", + "assert np.abs(u_cart_nopert[2] - 0.0082926160732301) < 1e-10\n", + "assert np.abs(f_cart_nopert[2] - 0.004760957079948613) < 1e-10\n", + "assert np.abs(u_triang_nopert[2] - 0.006007054619672307) < 1e-10\n", + "assert np.abs(f_triang_nopert[2] - 0.0058474670379793795) < 1e-10\n" ] }, { @@ -377,7 +588,7 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 5, "metadata": {}, "outputs": [ { @@ -581,7 +792,7 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 6, "metadata": {}, "outputs": [ { @@ -793,7 +1004,7 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 7, "metadata": { "scrolled": false }, @@ -802,42 +1013,42 @@ "name": "stdout", "output_type": "stream", "text": [ - "Time for discretization 1.2574834823608398\n", - "Time for AMG 0.011420965194702148\n", - "Time for discretization 1.0321142673492432\n", - "Time for AMG 0.07628965377807617\n", - "Time for discretization 2.3645739555358887\n", - "Time for AMG 0.3568305969238281\n", + "Time for discretization 1.8826446533203125\n", + "Time for AMG 0.04243111610412598\n", + "Time for discretization 2.377439260482788\n", + "Time for AMG 0.1383042335510254\n", + "Time for discretization 2.542255401611328\n", + "Time for AMG 0.3778951168060303\n", "Cartesian errors - no perturbations\n", "[0.14582494480280667, 0.036685508143867125, 0.009623477148373078] displacement error\n", "[0.14173578960074407, 0.05453461058511372, 0.018777567804334296] stress error\n", - "Time for discretization 1.199927806854248\n", - "Time for AMG 0.010801315307617188\n", - "Time for discretization 1.252537727355957\n", - "Time for AMG 0.07427144050598145\n", - "Time for discretization 2.594367265701294\n", - "Time for AMG 0.3139197826385498\n", + "Time for discretization 1.409085750579834\n", + "Time for AMG 0.015498161315917969\n", + "Time for discretization 1.3807973861694336\n", + "Time for AMG 0.07768678665161133\n", + "Time for discretization 7.617006778717041\n", + "Time for AMG 0.3285984992980957\n", "Cartesian errors - perturbation 0.5\n", "[0.1823425575505196, 0.03669658584037131, 0.009785391062793553] displacement error\n", "[0.15394769153699012, 0.06655533666156926, 0.024831033744125167] stress error\n", - "Time for discretization 1.2206642627716064\n", - "Time for AMG 0.061220407485961914\n", - "Time for discretization 1.5830979347229004\n", - "Time for AMG 0.3124885559082031\n", - "Time for discretization 5.6347815990448\n", - "Time for AMG 4.207815170288086\n", + "Time for discretization 1.9705743789672852\n", + "Time for AMG 0.07087302207946777\n", + "Time for discretization 7.326463222503662\n", + "Time for AMG 0.4231894016265869\n", + "Time for discretization 10.036577224731445\n", + "Time for AMG 5.2540881633758545\n", "Triangular errors - no perturbations\n", "[0.06256710001786384, 0.015518298189089326, 0.004112646550453208] displacement error\n", "[0.09819533736107755, 0.040432049644100486, 0.01608810670167512] stress error\n", - "Time for discretization 1.2895686626434326\n", - "Time for AMG 0.05126047134399414\n", - "Time for discretization 1.8259963989257812\n", - "Time for AMG 0.30768799781799316\n", - "Time for discretization 6.248790264129639\n", - "Time for AMG 4.1906561851501465\n", + "Time for discretization 1.2782678604125977\n", + "Time for AMG 0.07084035873413086\n", + "Time for discretization 3.013735055923462\n", + "Time for AMG 0.4166250228881836\n", + "Time for discretization 10.38901948928833\n", + "Time for AMG 5.382190704345703\n", "Triangular errors - perturbation 0.5\n", "[0.07257865039340129, 0.01699696434861487, 0.004336430438721095] displacement error\n", - "[0.10603237914286788, 0.04292864828702439, 0.017041050867287796] stress error\n" + "[0.10603237914286788, 0.04292864828702438, 0.017041050867287796] stress error\n" ] } ], @@ -1049,7 +1260,7 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 8, "metadata": {}, "outputs": [], "source": [ From c6fa6e773bfcbe37daa167704daf43be3186495d Mon Sep 17 00:00:00 2001 From: Runar Date: Tue, 11 Sep 2018 11:10:32 +0200 Subject: [PATCH 11/26] Fixed an index bug in MPFA where the robin faces was given wrong area --- examples/example2/ConvergenceTestMPFA.ipynb | 53 ++++++++++++++++++--- src/porepy/numerics/fv/mpfa.py | 10 ++-- src/porepy/numerics/fv/mpsa.py | 11 +++-- src/porepy/params/bc.py | 1 + 4 files changed, 58 insertions(+), 17 deletions(-) diff --git a/examples/example2/ConvergenceTestMPFA.ipynb b/examples/example2/ConvergenceTestMPFA.ipynb index 646f1477b8..77efc77afe 100644 --- a/examples/example2/ConvergenceTestMPFA.ipynb +++ b/examples/example2/ConvergenceTestMPFA.ipynb @@ -179,10 +179,17 @@ "text": [ "Cartesian errors - no perturbations\n", "[0.2177630338617197, 0.053854028644456195, 0.013418192291060563] pressure error\n", + "[4.04357927 4.01350849] convergence rate u^n-1 / u^n\n", + "\n", "[0.020449627547754923, 0.00558959841388443, 0.0014594922987497814] flux error\n", + "[3.65851463 3.82982385] convergence rate f^n-1 / f^n\n", + "\n", "Triangular errors - no perturbations\n", "[0.1154963072658415, 0.028265747189434955, 0.007020068836548777] pressure error\n", - "[0.04223113257437278, 0.011882184577715638, 0.0033976662367782447] flux error\n" + "[4.08608718 4.02642023] convergence rate u^n-1 / u^n\n", + "\n", + "[0.04223113257437278, 0.011882184577715638, 0.0033976662367782447] flux error\n", + "[3.55415558 3.49716062] convergence rate f^n-1 / f^n\n" ] } ], @@ -276,15 +283,27 @@ "\n", "\n", "u_cart_nopert, f_cart_nopert = run_convergence('cart', 0)\n", + "u_cart_rate = np.array(u_cart_nopert)\n", + "u_cart_rate = u_cart_rate[:-1] / u_cart_rate[1:]\n", + "f_cart_rate = np.array(f_cart_nopert)\n", + "f_cart_rate = f_cart_rate[:-1] / f_cart_rate[1:]\n", "print('Cartesian errors - no perturbations')\n", "print(u_cart_nopert,'pressure error')\n", + "print(u_cart_rate, 'convergence rate u^n-1 / u^n\\n')\n", "print(f_cart_nopert,'flux error')\n", + "print(f_cart_rate, 'convergence rate f^n-1 / f^n\\n')\n", "\n", "u_triang_nopert, f_triang_nopert = run_convergence('triangular', 0)\n", + "u_triang_rate = np.array(u_triang_nopert)\n", + "u_triang_rate = u_triang_rate[:-1] / u_triang_rate[1:]\n", + "f_triang_rate = np.array(f_triang_nopert)\n", + "f_triang_rate = f_triang_rate[:-1] / f_triang_rate[1:]\n", + "\n", "print('Triangular errors - no perturbations')\n", "print(u_triang_nopert,'pressure error')\n", + "print(u_triang_rate, 'convergence rate u^n-1 / u^n\\n')\n", "print(f_triang_nopert,'flux error')\n", - "\n", + "print(f_triang_rate, 'convergence rate f^n-1 / f^n')\n", "# EK: These values were hard-coded 11.06.2017.\n", "assert np.abs(u_cart_nopert[2] - 0.01341819229106214) < 1e-10\n", "assert np.abs(f_cart_nopert[2] - 0.0014594922987503009) < 1e-10\n", @@ -312,10 +331,17 @@ "text": [ "Cartesian errors - no perturbations\n", "[0.562515011909662, 0.14015471392164033, 0.034984342689359696] pressure error\n", + "[4.01352902 4.00621258] convergence rate u^n-1 / u^n\n", + "\n", "[0.22376060444435023, 0.061834410597352465, 0.01641158802403413] flux error\n", + "[3.61870684 3.76772866] convergence rate f^n-1 / f^n\n", + "\n", "Triangular errors - no perturbations\n", - "[1.1696464364631856, 0.8830933068427524, 0.49114177612455606] pressure error\n", - "[0.5558236600090738, 0.40733862448942615, 0.22657669056462687] flux error\n" + "[0.11192090940414492, 0.02516849817177693, 0.005995729167675974] pressure error\n", + "[4.44686483 4.19773767] convergence rate u^n-1 / u^n\n", + "\n", + "[0.03311849158106908, 0.00866805360137708, 0.002734187490448117] flux error\n", + "[3.82075297 3.17024843] convergence rate f^n-1 / f^n\n" ] } ], @@ -415,20 +441,33 @@ "\n", "\n", "u_cart_nopert, f_cart_nopert = run_convergence('cart', 0)\n", + "u_cart_rate = np.array(u_cart_nopert)\n", + "u_cart_rate = u_cart_rate[:-1] / u_cart_rate[1:]\n", + "f_cart_rate = np.array(f_cart_nopert)\n", + "f_cart_rate = f_cart_rate[:-1] / f_cart_rate[1:]\n", "print('Cartesian errors - no perturbations')\n", "print(u_cart_nopert,'pressure error')\n", + "print(u_cart_rate, 'convergence rate u^n-1 / u^n\\n')\n", "print(f_cart_nopert,'flux error')\n", + "print(f_cart_rate, 'convergence rate f^n-1 / f^n\\n')\n", "\n", "u_triang_nopert, f_triang_nopert = run_convergence('triangular', 0)\n", + "u_triang_rate = np.array(u_triang_nopert)\n", + "u_triang_rate = u_triang_rate[:-1] / u_triang_rate[1:]\n", + "f_triang_rate = np.array(f_triang_nopert)\n", + "f_triang_rate = f_triang_rate[:-1] / f_triang_rate[1:]\n", + "\n", "print('Triangular errors - no perturbations')\n", "print(u_triang_nopert,'pressure error')\n", + "print(u_triang_rate, 'convergence rate u^n-1 / u^n\\n')\n", "print(f_triang_nopert,'flux error')\n", + "print(f_triang_rate, 'convergence rate f^n-1 / f^n')\n", "\n", - "# RB: These values were hard-coded 07.09.2018.\n", + "# RB: These values were hard-coded 11.09.2018.\n", "assert np.abs(u_cart_nopert[2] - 0.034984342689359696) < 1e-10\n", "assert np.abs(f_cart_nopert[2] - 0.01641158802403413) < 1e-10\n", - "assert np.abs(u_triang_nopert[2] - 0.49114177612455606) < 1e-10\n", - "assert np.abs(f_triang_nopert[2] - 0.22657669056462687) < 1e-10\n" + "assert np.abs(u_triang_nopert[2] - 0.005995729167675974) < 1e-10\n", + "assert np.abs(f_triang_nopert[2] - 0.002734187490448117) < 1e-10\n" ] }, { diff --git a/src/porepy/numerics/fv/mpfa.py b/src/porepy/numerics/fv/mpfa.py index edf1fc1f81..1d28e7733f 100644 --- a/src/porepy/numerics/fv/mpfa.py +++ b/src/porepy/numerics/fv/mpfa.py @@ -587,13 +587,13 @@ def _mpfa_local(g, k, bnd, eta=None, inverter="numba", apertures=None, alpha=Non # will be the contribution from the gradients. We integrate over the subface # and multiply by the area num_nodes = np.diff(g.face_nodes.indptr) + sgn = g.cell_faces[subcell_topology.fno_unique, subcell_topology.cno_unique].A scaled_sgn = ( - alpha * sgn[0] * g.face_areas[subcell_topology.fno] \ - / num_nodes[subcell_topology.fno] - ) - + alpha * sgn[0] * g.face_areas[subcell_topology.fno_unique] \ + / num_nodes[subcell_topology.fno_unique] + ) # pair_over_subfaces flips the sign so we flip it backo - pr_trace_grad_all = pr_cont_grad_all * sps.diags(scaled_sgn) + pr_trace_grad_all = sps.diags(scaled_sgn) * pr_cont_grad_all pr_trace_cell_all = sps.coo_matrix( (alpha * g.face_areas[subcell_topology.fno] / num_nodes[subcell_topology.fno], (subcell_topology.subfno, subcell_topology.cno)) diff --git a/src/porepy/numerics/fv/mpsa.py b/src/porepy/numerics/fv/mpsa.py index 9e04b6d95c..38b1a12c57 100644 --- a/src/porepy/numerics/fv/mpsa.py +++ b/src/porepy/numerics/fv/mpsa.py @@ -1039,11 +1039,6 @@ def _mpsa_local(g, constit, bound, eta=0, alpha=None, inverter="numba"): nd = g.dim - if alpha is None: - if np.sum(bound.is_rob) !=0: - raise ValueError('If applying Robin conditions you must supply an alpha') - else: - alpha = 1 # Define subcell topology subcell_topology = fvutils.SubcellTopology(g) # Obtain mappings to exclude boundary faces @@ -1114,6 +1109,12 @@ def mpsa_elasticity(g, constit, subcell_topology, bound_exclusion, eta, inverter nd = g.dim + if alpha is None: + if np.sum(bound.is_rob) !=0: + raise ValueError('If applying Robin conditions you must supply an alpha') + else: + alpha = 1 + # Compute product between normal vectors and stiffness matrices ncsym_all, ncasym, cell_node_blocks, sub_cell_index = _tensor_vector_prod( g, constit, subcell_topology diff --git a/src/porepy/params/bc.py b/src/porepy/params/bc.py index ec93111c8b..6347c02a38 100644 --- a/src/porepy/params/bc.py +++ b/src/porepy/params/bc.py @@ -233,6 +233,7 @@ def __init__(self, g, faces=None, cond=None): self.is_neu = np.zeros((g.dim, self.num_faces), dtype=bool) self.is_dir = np.zeros((g.dim, self.num_faces), dtype=bool) + self.is_rob = np.zeros((g.dim, self.num_faces), dtype=bool) self.is_neu[:, self.bf] = True self.set_bc(faces, cond) From 6b074042dc54cd1213e97af7cee384c20a0220b6 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Tue, 11 Sep 2018 14:32:27 +0200 Subject: [PATCH 12/26] Version bump --- setup.py | 2 +- src/porepy/__init__.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/setup.py b/setup.py index 9688d051b5..c2d98f1db0 100644 --- a/setup.py +++ b/setup.py @@ -60,7 +60,7 @@ def read(fname): setup( name="porepy", - version="0.4.1", + version="0.4.2", license="GPL", keywords=["porous media simulation fractures deformable"], author="Runar Berge, Alessio Fumagalli, Eirik Keilegavlen and Ivar Stefansson", diff --git a/src/porepy/__init__.py b/src/porepy/__init__.py index ecf99a92e7..56398dbc8b 100644 --- a/src/porepy/__init__.py +++ b/src/porepy/__init__.py @@ -16,7 +16,7 @@ """ -__version__ = "0.4.1" +__version__ = "0.4.2" # ------------------------------------ # Simplified namespaces. The rue of thumb is that classes and modules that a From fe1010eb71d70fd8dde467e9cd33a112b533a36c Mon Sep 17 00:00:00 2001 From: Runar Date: Tue, 11 Sep 2018 15:32:24 +0200 Subject: [PATCH 13/26] Added 3D test case for mpsa_robin --- test/unit/test_mpsa_robin.py | 56 ++++++++++++++++++++++++++++++++++++ 1 file changed, 56 insertions(+) diff --git a/test/unit/test_mpsa_robin.py b/test/unit/test_mpsa_robin.py index 9cacdd9fb1..d784708d84 100644 --- a/test/unit/test_mpsa_robin.py +++ b/test/unit/test_mpsa_robin.py @@ -194,6 +194,62 @@ def T_ex(faces): assert np.allclose(u, u_ex(g.cell_centers).ravel('F')) assert np.allclose(T, T_ex(np.arange(g.num_faces)).ravel('F')) + def test_unstruct_tetrahedron(self): + box = {'xmin': 0, 'xmax': 1, + 'ymin': 0, 'ymax': 1, + 'zmin': 0, 'zmax': 1} + + g = pp.meshing.simplex_grid([], box, mesh_size_min=3, mesh_size_frac=3) + g = g.grids_of_dimension(3)[0] + g = pp.CartGrid([4,4,4], [1,1,1]) + g.compute_geometry() + c = pp.FourthOrderTensor(3, np.ones(g.num_cells), np.ones(g.num_cells)) + alpha = 1.0 + + bot = g.face_centers[2] < 1e-10 + top = g.face_centers[2] > 1 - 1e-10 + west = g.face_centers[0] < 1e-10 + east = g.face_centers[0] > 1 - 1e-10 + north = g.face_centers[1] > 1 - 1e-10 + south = g.face_centers[1] < 1e-10 + + dir_ind = np.ravel(np.argwhere(west + top)) + neu_ind = np.ravel(np.argwhere(bot)) + rob_ind = np.ravel(np.argwhere(east + north + south)) + + names = ['dir']*len(dir_ind) + ['rob'] * len(rob_ind) + bnd_ind = np.hstack((dir_ind, rob_ind)) + bnd = pp.BoundaryCondition(g, bnd_ind, names) + + def u_ex(x): + return np.vstack((x[1], x[0], 0*x[2])) + + def T_ex(faces): + if np.size(faces) == 0: + return np.atleast_2d(np.array([])) + sigma = np.array([[0, 2, 0], + [2, 0, 0], + [0, 0, 0]]) + T_r = [np.dot(sigma, g.face_normals[:, f]) for f in faces] + return np.vstack(T_r).T + + u_bound = np.zeros((3, g.num_faces)) + + sgn_n = pp.numerics.fracture_deformation.sign_of_faces(g, neu_ind) + sgn_r = pp.numerics.fracture_deformation.sign_of_faces(g, rob_ind) + + u_bound[:, dir_ind] = u_ex(g.face_centers[:, dir_ind]) + u_bound[:, neu_ind] = T_ex(neu_ind) * sgn_n + u_bound[:, rob_ind] = ( + T_ex(rob_ind) * sgn_r + + alpha * u_ex(g.face_centers[:, rob_ind]) * g.face_areas[rob_ind] + ) + u, T = self.solve_mpsa(g, c, alpha, bnd, u_bound) + + assert np.allclose(u, u_ex(g.cell_centers).ravel('F')) + assert np.allclose(T, T_ex(np.arange(g.num_faces)).ravel('F')) + + def solve_mpsa(self, g, c, alpha, bnd, u_bound): stress, bound_stress = pp.numerics.fv.mpsa._mpsa_local(g, c, bnd, alpha=alpha) div = pp.fvutils.vector_divergence(g) From b5873e9681a9395d3adcc3f66adb18422aa2e85d Mon Sep 17 00:00:00 2001 From: Runar Date: Tue, 11 Sep 2018 15:32:40 +0200 Subject: [PATCH 14/26] Fixed some notation for robin condition --- src/porepy/numerics/fv/fvutils.py | 162 +++++++++++----------- src/porepy/numerics/fv/mpfa.py | 36 ++--- src/porepy/numerics/fv/mpsa.py | 218 +++--------------------------- 3 files changed, 121 insertions(+), 295 deletions(-) diff --git a/src/porepy/numerics/fv/fvutils.py b/src/porepy/numerics/fv/fvutils.py index c825845cf4..83b2ba6e3c 100644 --- a/src/porepy/numerics/fv/fvutils.py +++ b/src/porepy/numerics/fv/fvutils.py @@ -917,6 +917,86 @@ def exclude_robin(self, other): return exclude_rob + def exclude_neu_rob(self, other): + """ Mapping to exclude faces/components with Neumann and Robin boundary + conditions from local systems. + + Parameters: + other (scipy.sparse matrix): Matrix of local equations for + continuity of flux and pressure. + + Returns: + scipy.sparse matrix, with rows corresponding to faces/components with + Neumann conditions eliminated. + + """ + if self.bc_type == "scalar": + exclude_neu = self._exclude_neu_rob* other + + elif self.bc_type == "vectorial": + raise NotImplementedError() + return exclude_neu + + def exclude_neu_dir(self, other): + """ Mapping to exclude faces/components with Neumann and Dirichlet boundary + conditions from local systems. + + Parameters: + other (scipy.sparse matrix): Matrix of local equations for + continuity of flux and pressure. + + Returns: + scipy.sparse matrix, with rows corresponding to faces/components with + Neumann conditions eliminated. + + """ + if self.bc_type == "scalar": + exclude_neu = self._exclude_neu_dir * other + + elif self.bc_type == "vectorial": + raise NotImplementedError() + return exclude_neu + + def exclude_rob_dir(self, other): + """ Mapping to exclude faces/components with Robin and Dirichlet boundary + conditions from local systems. + + Parameters: + other (scipy.sparse matrix): Matrix of local equations for + continuity of flux and pressure. + + Returns: + scipy.sparse matrix, with rows corresponding to faces/components with + Neumann conditions eliminated. + + """ + if self.bc_type == "scalar": + exclude_rob = self._exclude_rob_dir* other + + elif self.bc_type == "vectorial": + raise NotImplementedError() + return exclude_rob + + def keep_robin(self, other): + """ Mapping to exclude faces/components that is not on the Robin boundary + conditions from local systems. + + Parameters: + other (scipy.sparse matrix): Matrix of local equations for + continuity of flux and pressure. + + Returns: + scipy.sparse matrix, with rows corresponding to faces/components with + all but robin conditions eliminated. + + """ + if self.bc_type == "scalar": + exclude_rob = self._keep_robin * other + + elif self.bc_type == "vectorial": + raise NotImplementedError() + return exclude_rob + def exclude_neumann_x(self, other): """ Mapping to exclude components in the x-direction with Neumann boundary conditions from local systems. @@ -1046,88 +1126,6 @@ def keep_robin_nd(self, other): return exclude_neu_dir_nd * other - def exclude_neu_rob(self, other): - """ Mapping to exclude faces/components with Neumann and Robin boundary - conditions from local systems. - - Parameters: - other (scipy.sparse matrix): Matrix of local equations for - continuity of flux and pressure. - - Returns: - scipy.sparse matrix, with rows corresponding to faces/components with - Neumann conditions eliminated. - - """ - if self.bc_type == "scalar": - exclude_neu = self._exclude_neu_rob* other - - elif self.bc_type == "vectorial": - raise NotImplementedError() - return exclude_neu - - def exclude_neu_dir(self, other): - """ Mapping to exclude faces/components with Neumann and Dirichlet boundary - conditions from local systems. - - Parameters: - other (scipy.sparse matrix): Matrix of local equations for - continuity of flux and pressure. - - Returns: - scipy.sparse matrix, with rows corresponding to faces/components with - Neumann conditions eliminated. - - """ - if self.bc_type == "scalar": - exclude_neu = self._exclude_neu_dir * other - - elif self.bc_type == "vectorial": - raise NotImplementedError() - return exclude_neu - - def exclude_rob_dir(self, other): - """ Mapping to exclude faces/components with Robin and Dirichlet boundary - conditions from local systems. - - Parameters: - other (scipy.sparse matrix): Matrix of local equations for - continuity of flux and pressure. - - Returns: - scipy.sparse matrix, with rows corresponding to faces/components with - Neumann conditions eliminated. - - """ - if self.bc_type == "scalar": - exclude_rob = self._exclude_rob_dir* other - - elif self.bc_type == "vectorial": - raise NotImplementedError() - return exclude_rob - - def keep_robin(self, other): - """ Mapping to exclude faces/components that is not on the Robin boundary - conditions from local systems. - - Parameters: - other (scipy.sparse matrix): Matrix of local equations for - continuity of flux and pressure. - - Returns: - scipy.sparse matrix, with rows corresponding to faces/components with - all but robin conditions eliminated. - - """ - if self.bc_type == "scalar": - exclude_rob = self._keep_robin * other - - elif self.bc_type == "vectorial": - raise NotImplementedError() - return exclude_rob - - - # -----------------End of class ExcludeBoundaries----------------------------- diff --git a/src/porepy/numerics/fv/mpfa.py b/src/porepy/numerics/fv/mpfa.py index 1d28e7733f..b7963ae9e9 100644 --- a/src/porepy/numerics/fv/mpfa.py +++ b/src/porepy/numerics/fv/mpfa.py @@ -178,7 +178,7 @@ def discretize(self, g, data): # ------------------------------------------------------------------------------# -def mpfa(g, k, bnd, eta=None, inverter=None, apertures=None, max_memory=None, alpha=None, **kwargs): +def mpfa(g, k, bnd, robin_weight=None, eta=None, inverter=None, apertures=None, max_memory=None, **kwargs): """ Discretize the scalar elliptic equation by the multi-point flux approximation method. @@ -203,6 +203,8 @@ def mpfa(g, k, bnd, eta=None, inverter=None, apertures=None, max_memory=None, al g (core.grids.grid): grid to be discretized k (core.constit.second_order_tensor) permeability tensor bnd (core.bc.bc) class for boundary values + robin_weight (float): Weight robin_weight for displacement term in Robin conditions + sigma*n + robin_weight * u = G eta Location of pressure continuity point. Defaults to 1/3 for simplex grids, 0 otherwise. On boundary faces with Dirichlet conditions, eta=0 will be enforced. @@ -269,7 +271,7 @@ def mpfa(g, k, bnd, eta=None, inverter=None, apertures=None, max_memory=None, al # TODO: We may want to estimate the memory need, and give a warning if # this seems excessive flux, bound_flux, bound_pressure_cell, bound_pressure_face = _mpfa_local( - g, k, bnd, eta=eta, inverter=inverter, apertures=apertures, alpha=alpha + g, k, bnd, eta=eta, inverter=inverter, apertures=apertures, robin_weight=robin_weight ) else: # Estimate number of partitions necessary based on prescribed memory @@ -456,7 +458,7 @@ def mpfa_partial( ) -def _mpfa_local(g, k, bnd, eta=None, inverter="numba", apertures=None, alpha=None): +def _mpfa_local(g, k, bnd, eta=None, inverter="numba", apertures=None, robin_weight=None): """ Actual implementation of the MPFA O-method. To calculate MPFA on a grid directly, either call this method, or, to respect the privacy of this @@ -471,15 +473,17 @@ def _mpfa_local(g, k, bnd, eta=None, inverter="numba", apertures=None, alpha=Non The method will give continuous fluxes over the faces, and pressure continuity for certain points (controlled by the parameter eta). This can be expressed as a linear system on the form - (i) A * grad_p = 0 - (ii) B * grad_p + C * p_cc = 0 - (iii) 0 D * p_cc = I + (i) A * grad_p = 0 + (ii) Ar * grad_p + Cr * P_cc = 0 + (iii) B * grad_p + C * p_cc = 0 + (iv) 0 D * p_cc = I Here, the first equation represents flux continuity, and involves only the - pressure gradients (grad_p). The second equation gives pressure continuity + pressure gradients (grad_p). The second equation gives the Robin conditions, + relating flux to the pressure. The third equation gives pressure continuity over cell faces, thus B will contain distances between cell centers and the face continuity points, while C consists of +- 1 (depending on which side - the cell is relative to the face normal vector). The third equation - enforces the pressure to be unity in one cell at a time. Thus (i)-(iii) can + the cell is relative to the face normal vector). The fourth equation + enforces the pressure to be unity in one cell at a time. Thus (i)-(iv) can be inverted to express the pressure gradients as in terms of the cell center variables, that is, we can compute the basis functions on the sub-cells. Because of the method construction (again see reference paper), @@ -489,7 +493,7 @@ def _mpfa_local(g, k, bnd, eta=None, inverter="numba", apertures=None, alpha=Non (i), that is, only consider contribution from one side of the face. Boundary values can be incorporated with appropriate modifications - Neumann conditions will have a non-zero right hand side for (i), while - Dirichlet gives a right hand side for (ii). + Dirichlet gives a right hand side for (iii). """ if eta is None: @@ -515,11 +519,11 @@ def _mpfa_local(g, k, bnd, eta=None, inverter="numba", apertures=None, alpha=Non elif g.dim == 0: return sps.csr_matrix([0]), 0, 0, 0 - if alpha is None: + if robin_weight is None: if np.sum(bnd.is_rob) !=0: - raise ValueError('If applying Robin conditions you must supply an alpha') + raise ValueError('If applying Robin conditions you must supply an robin_weight') else: - alpha = 1 + robin_weight = 1 # The grid coordinates are always three-dimensional, even if the grid is # really 2D. This means that there is not a 1-1 relation between the number # of coordinates of a point / vector and the real dimension. This again @@ -589,13 +593,13 @@ def _mpfa_local(g, k, bnd, eta=None, inverter="numba", apertures=None, alpha=Non num_nodes = np.diff(g.face_nodes.indptr) sgn = g.cell_faces[subcell_topology.fno_unique, subcell_topology.cno_unique].A scaled_sgn = ( - alpha * sgn[0] * g.face_areas[subcell_topology.fno_unique] \ + robin_weight * sgn[0] * g.face_areas[subcell_topology.fno_unique] \ / num_nodes[subcell_topology.fno_unique] ) - # pair_over_subfaces flips the sign so we flip it backo + # pair_over_subfaces flips the sign so we flip it back pr_trace_grad_all = sps.diags(scaled_sgn) * pr_cont_grad_all pr_trace_cell_all = sps.coo_matrix( - (alpha * g.face_areas[subcell_topology.fno] / num_nodes[subcell_topology.fno], + (robin_weight * g.face_areas[subcell_topology.fno] / num_nodes[subcell_topology.fno], (subcell_topology.subfno, subcell_topology.cno)) ).tocsr() diff --git a/src/porepy/numerics/fv/mpsa.py b/src/porepy/numerics/fv/mpsa.py index 38b1a12c57..f2ae477bdb 100644 --- a/src/porepy/numerics/fv/mpsa.py +++ b/src/porepy/numerics/fv/mpsa.py @@ -134,185 +134,6 @@ def rhs(self, g, bound_stress, bc_val, f): return -div * bound_stress * bc_val - f - -# ------------------------------------------------------------------------------# - - -class RobinMpsa(Mpsa): - """ - Subclass of MPSA for discretizing a domain with robin boundary conditinos. - Adds DOFs on each robin face. - """ - - def __init__(self, robin_faces, **kwargs): - """ - robin_faces is a function that takes a grid and return an array - of length grid.num_faces that are true on the robin faces and false otherwise - """ - Mpsa.__init__(self, **kwargs) - assert hasattr(self, "physics"), "Mpsa must assign physics" - self.robin_faces = robin_faces - - def ndof(self, g): - """ - Return the number of degrees of freedom associated to the method. - In this case number of cells pluss the number of robin face times - dimension (stress dof). - - Parameter - --------- - g: grid, or a subclass. - - Return - ------ - dof: the number of degrees of freedom. - - """ - return g.dim * (g.num_cells + np.sum(self.robin_faces(g))) - - def matrix_rhs(self, g, data, discretize=True): - """ - Return the matrix and right-hand side for a discretization of a second - order elliptic equation using a FV method with a multi-point stress - approximation with dofs added on the robin faces - - Parameters - ---------- - g : grid, or a subclass, with geometry fields computed. - data: dictionary to store the data. For details on necessary keywords, - see method discretize() - discretize (boolean, optional): default True. Whether to discetize - prior to matrix assembly. If False, data should already contain - discretization. - - Return - ------ - matrix: sparse csr (g.dim * (g_num_cells + {#of robin faces}, - g.dim * (g_num_cells + {#of robin faces})) - Discretization matrix. - rhs: array (g.dim * g_num_cells + {#of robin faces}) - Right-hand side which contains the boundary conditions and the scalar - source term. - """ - if discretize: - self.discretize_robin(g, data) - - b_e = data["b_e"] - A_e = data["A_e"] - - bc_val = data["param"].get_bc_val(self) - - rhs = b_e * bc_val - - return A_e, rhs - - def discretize_robin(self, g, data, faces=None, **kwargs): - """ - Discretize the vector elliptic equation by the multi-point stress and added - degrees of freedom on the robin faces - - The method computes fluxes over faces in terms of displacements in - adjacent cells (defined as the two cells sharing the face). - - The name of data in the input dictionary (data) are: - param : Parameter(Class). Contains the following parameters: - tensor : fourth_order_tensor - Permeability defined cell-wise. If not given a identity permeability - is assumed and a warning arised. - bc : boundary conditions (optional) - bc_val : dictionary (optional) - Values of the boundary conditions. The dictionary has at most the - following keys: 'dir' and 'neu', for Dirichlet and Neumann boundary - conditions, respectively. - apertures : (np.ndarray) (optional) apertures of the cells for scaling of - the face normals. - - Parameters - ---------- - g : grid, or a subclass, with geometry fields computed. - data: dictionary to store the data. - """ - - # dir_bound = g.get_all_boundary_faces() - # bound = bc.BoundaryCondition(g, dir_bound, ['dir'] * dir_bound.size) - - rob_faces = self.robin_faces(g) - - bound = data["param"].get_bc(self) - - if bound.bc_type == "scalar": - bound.is_dir[rob_faces] = True - bound.is_neu[rob_faces] = False - elif bound.bc_type == "vectorial": - bound.is_dir[:, rob_faces] = True - bound.is_neu[:, rob_faces] = False - else: - raise ValueError("Unknow boundary condition type: " + bound.bc_type) - if np.sum(bound.is_dir * bound.is_neu) !=0: - raise AssertionError('Found faces that are both dirichlet and neuman') - - # Discretize with normal mpsa - self.discretize(g, data, **kwargs) - stress, bound_stress = data["stress"], data["bound_stress"] - - # Create A and rhs - div = fvutils.vector_divergence(g) - a = div * stress - b = div * bound_stress - - # we find the matrix indices of the robin faces - rob_ind = mcolon.mcolon( - g.dim * rob_faces, g.dim * rob_faces + g.dim - ) - # We find the sign of the robin faces. - sgn_rob = _sign_matrix(g, rob_faces) - - # The displacement on the robin boundary face are considered unknowns, - # so we move them over to the lhs. The rhs now only consists of the - # external boundary faces - b_rob = b[:, rob_ind] - b_external = b.copy() - sparse_mat.zero_columns(b_external, rob_ind) - - bound_stress_external = bound_stress.copy().tocsc() - sparse_mat.zero_columns(bound_stress_external, rob_ind) - # We calculate the stress on the robin faces due to the internal cells - # and robin face displacements (the other boundaries are on the rhs) - roben_stress = sps.hstack( - ( - sgn_rob * stress[rob_ind, :], - (sgn_rob * bound_stress[rob_ind, :])[:, rob_ind] - ) - ) - # we add the displacement on the robin boundaries - alpha = data["param"].get_robin_factor()[rob_faces] - alpha = np.ravel([alpha]*g.dim, 'F') - roben_disp = sps.hstack( - ( - sps.csr_matrix((rob_ind.size, g.dim * g.num_cells)), - sps.diags(alpha, 0), - ) - ) - roben_condition = roben_stress + roben_disp - A = sps.vstack((sps.hstack((a, b_rob)), roben_condition), format="csr") - - # negative sign since we have moved b_external from lhs to rhs - d_b = -b_external - # sps.csr_matrix((int_b_left.size, g.num_faces * g.dim)) - d_t = ( - sps.csr_matrix( - (np.ones(rob_ind.size), (np.arange(rob_ind.size), rob_ind)), - (rob_ind.size, g.num_faces * g.dim), - ) - -sgn_rob * bound_stress_external[rob_ind] - ) - - b_matrix = sps.vstack((d_b, d_t), format="csr") - - data["b_e"] = b_matrix - data["A_e"] = A - - class FracturedMpsa(Mpsa): """ Subclass of MPSA for discretizing a fractured domain. Adds DOFs on each @@ -703,7 +524,7 @@ def given_slip_distance(self, g, stress, bound_stress, faces=None): # ------------------------------------------------------------------------------# -def mpsa(g, constit, bound, eta=None, alpha=None, inverter=None, max_memory=None, **kwargs): +def mpsa(g, constit, bound, robin_weight=None, eta=None, inverter=None, max_memory=None, **kwargs): """ Discretize the vector elliptic equation by the multi-point stress approximation method, specifically the weakly symmetric MPSA-W method. @@ -731,7 +552,10 @@ def mpsa(g, constit, bound, eta=None, alpha=None, inverter=None, max_memory=None Parameters: g (core.grids.grid): grid to be discretized - constit (core.bc.bc) class for boundary values + constit (pp.FourthOrderTensor) Constitutive law + bound (pp.BoundarCondition) Class for boundary condition + robin_weight (float): Weight alpha for displacement term in Robin conditions + sigma*n + alpha * u = G eta Location of pressure continuity point. Should be 1/3 for simplex grids, 0 otherwise. On boundary faces with Dirichlet conditions, eta=0 will be enforced. @@ -799,7 +623,7 @@ def mpsa(g, constit, bound, eta=None, alpha=None, inverter=None, max_memory=None # TODO: We may want to estimate the memory need, and give a warning if # this seems excessive stress, bound_stress = _mpsa_local( - g, constit, bound, eta=eta, inverter=inverter, alpha=alpha + g, constit, bound, eta=eta, inverter=inverter, robin_weight=robin_weight ) else: # Estimate number of partitions necessary based on prescribed memory @@ -836,7 +660,7 @@ def mpsa(g, constit, bound, eta=None, alpha=None, inverter=None, max_memory=None # Perform local discretization. loc_stress, loc_bound_stress, loc_faces = mpsa_partial( g, constit, bound, eta=eta, inverter=inverter, nodes=active_nodes, - alpha=alpha + robin_weight=robin_weight ) # Eliminate contribution from faces already covered @@ -862,7 +686,7 @@ def mpsa_partial( faces=None, nodes=None, apertures=None, - alpha=None + robin_weight=None ): """ Run an MPFA discretization on subgrid, and return discretization in terms @@ -954,7 +778,7 @@ def mpsa_partial( # Discretization of sub-problem stress_loc, bound_stress_loc = _mpsa_local( - sub_g, loc_c, loc_bnd, eta=eta, inverter=inverter, alpha=alpha + sub_g, loc_c, loc_bnd, eta=eta, inverter=inverter, robin_weight=robin_weight ) face_map, cell_map = fvutils.map_subgrid_to_grid( @@ -975,7 +799,7 @@ def mpsa_partial( return stress_glob, bound_stress_glob, active_faces -def _mpsa_local(g, constit, bound, eta=0, alpha=None, inverter="numba"): +def _mpsa_local(g, constit, bound, eta=0, robin_weight=None, inverter="numba"): """ Actual implementation of the MPSA W-method. To calculate the MPSA discretization on a grid, either call this method, or, to respect the @@ -1047,7 +871,7 @@ def _mpsa_local(g, constit, bound, eta=0, alpha=None, inverter="numba"): # elasticity and poro-elasticity). hook, igrad, rhs_cells, _, _ = mpsa_elasticity( - g, constit, subcell_topology, bound_exclusion, eta, inverter, alpha=alpha + g, constit, subcell_topology, bound_exclusion, eta, inverter, robin_weight=robin_weight ) hook_igrad = hook * igrad @@ -1075,7 +899,7 @@ def _mpsa_local(g, constit, bound, eta=0, alpha=None, inverter="numba"): return stress, bound_stress -def mpsa_elasticity(g, constit, subcell_topology, bound_exclusion, eta, inverter, alpha=None): +def mpsa_elasticity(g, constit, subcell_topology, bound_exclusion, eta, inverter, robin_weight=None): """ This is the function where the real discretization takes place. It contains the parts that are common for elasticity and poro-elasticity, and was thus @@ -1109,11 +933,11 @@ def mpsa_elasticity(g, constit, subcell_topology, bound_exclusion, eta, inverter nd = g.dim - if alpha is None: - if np.sum(bound.is_rob) !=0: - raise ValueError('If applying Robin conditions you must supply an alpha') + if robin_weight is None: + if np.sum(bound_exclusion._keep_robin.shape[0]) !=0: + raise ValueError('If applying Robin conditions you must supply an robin_weight') else: - alpha = 1 + robin_weight = 1 # Compute product between normal vectors and stiffness matrices ncsym_all, ncasym, cell_node_blocks, sub_cell_index = _tensor_vector_prod( @@ -1157,7 +981,7 @@ def mpsa_elasticity(g, constit, subcell_topology, bound_exclusion, eta, inverter # Book keeping num_sub_cells = cell_node_blocks[0].size rob_grad, rob_cell = __get_displacement_submatrices_rob( - g, subcell_topology, eta, num_sub_cells, bound_exclusion, alpha + g, subcell_topology, eta, num_sub_cells, bound_exclusion, robin_weight ) # Matrices to enforce displacement continuity @@ -1257,7 +1081,7 @@ def __get_displacement_submatrices( return d_cont_grad, d_cont_cell def __get_displacement_submatrices_rob( - g, subcell_topology, eta, num_sub_cells, bound_exclusion, alpha + g, subcell_topology, eta, num_sub_cells, bound_exclusion, robin_weight ): nd = g.dim # Distance from cell centers to face centers, this will be the @@ -1265,20 +1089,20 @@ def __get_displacement_submatrices_rob( # at the boundary rob_grad = fvutils.compute_dist_face_cell(g, subcell_topology, eta) - # For the Robin codition the distance from the cell centers to face centers + # For the Robin condition the distance from the cell centers to face centers # will be the contribution from the gradients. We integrate over the subface # and multiply by the area num_nodes = np.diff(g.face_nodes.indptr) sgn = g.cell_faces[subcell_topology.fno_unique, subcell_topology.cno_unique].A scaled_sgn = ( - alpha * sgn[0] * g.face_areas[subcell_topology.fno_unique] \ + robin_weight * sgn[0] * g.face_areas[subcell_topology.fno_unique] \ / num_nodes[subcell_topology.fno_unique] ) # pair_over_subfaces flips the sign so we flip it back rob_grad = sps.kron(sps.eye(nd), sps.diags(scaled_sgn)*rob_grad) # Contribution from cell center potentials to local systems rob_cell = sps.coo_matrix( - (alpha * g.face_areas[subcell_topology.fno] / num_nodes[subcell_topology.fno], + (robin_weight * g.face_areas[subcell_topology.fno] / num_nodes[subcell_topology.fno], (subcell_topology.subfno, subcell_topology.cno)) ).tocsr() rob_cell = sps.kron(sps.eye(nd), rob_cell) From 220d12233454cb2f916f18f031df49874c3253e4 Mon Sep 17 00:00:00 2001 From: Runar Date: Tue, 11 Sep 2018 15:34:10 +0200 Subject: [PATCH 15/26] Typo --- src/porepy/numerics/fv/mpfa.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/porepy/numerics/fv/mpfa.py b/src/porepy/numerics/fv/mpfa.py index b7963ae9e9..33b37a190b 100644 --- a/src/porepy/numerics/fv/mpfa.py +++ b/src/porepy/numerics/fv/mpfa.py @@ -587,7 +587,7 @@ def _mpfa_local(g, k, bnd, eta=None, inverter="numba", apertures=None, robin_wei shape=(subcell_topology.num_subfno, subcell_topology.num_cno), ).tocsr() - # For the Robin codition the distance from the cell centers to face centers + # For the Robin condition the distance from the cell centers to face centers # will be the contribution from the gradients. We integrate over the subface # and multiply by the area num_nodes = np.diff(g.face_nodes.indptr) From 9ff1451b55d9c1d83fb65238cdad8bc7da83f14e Mon Sep 17 00:00:00 2001 From: Runar Date: Tue, 11 Sep 2018 16:09:31 +0200 Subject: [PATCH 16/26] Fixed import mistake --- src/porepy/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/porepy/__init__.py b/src/porepy/__init__.py index 56398dbc8b..a78abd5384 100644 --- a/src/porepy/__init__.py +++ b/src/porepy/__init__.py @@ -27,7 +27,7 @@ # Numerics # Control volume, elliptic -from porepy.numerics.fv.mpsa import Mpsa, FracturedMpsa, RobinMpsa +from porepy.numerics.fv.mpsa import Mpsa, FracturedMpsa from porepy.numerics.fv.tpfa import Tpfa, TpfaMixedDim from porepy.numerics.fv.mpfa import Mpfa, MpfaMixedDim from porepy.numerics.fv.biot import Biot From 36dbc661913149cdb1f0ec4fa09c7359c69033d9 Mon Sep 17 00:00:00 2001 From: Runar Date: Wed, 12 Sep 2018 10:18:29 +0200 Subject: [PATCH 17/26] Fixed keyword error alpha->robin_weight --- src/porepy/numerics/fv/mpsa.py | 13 +++++++------ test/unit/test_mpfa_robin.py | 28 ++++++++++++++-------------- test/unit/test_mpsa_robin.py | 34 +++++++++++++++++----------------- 3 files changed, 38 insertions(+), 37 deletions(-) diff --git a/src/porepy/numerics/fv/mpsa.py b/src/porepy/numerics/fv/mpsa.py index f2ae477bdb..f45147c700 100644 --- a/src/porepy/numerics/fv/mpsa.py +++ b/src/porepy/numerics/fv/mpsa.py @@ -933,12 +933,6 @@ def mpsa_elasticity(g, constit, subcell_topology, bound_exclusion, eta, inverter nd = g.dim - if robin_weight is None: - if np.sum(bound_exclusion._keep_robin.shape[0]) !=0: - raise ValueError('If applying Robin conditions you must supply an robin_weight') - else: - robin_weight = 1 - # Compute product between normal vectors and stiffness matrices ncsym_all, ncasym, cell_node_blocks, sub_cell_index = _tensor_vector_prod( g, constit, subcell_topology @@ -978,6 +972,13 @@ def mpsa_elasticity(g, constit, subcell_topology, bound_exclusion, eta, inverter # For the Robin boundary conditions we need to pair the forces with the # displacement ncsym_rob = bound_exclusion.keep_robin_nd(ncsym_all) + + if robin_weight is None: + if ncsym_rob.shape[0] !=0: + raise ValueError('If applying Robin conditions you must supply an robin_weight') + else: + robin_weight = 1 + # Book keeping num_sub_cells = cell_node_blocks[0].size rob_grad, rob_cell = __get_displacement_submatrices_rob( diff --git a/test/unit/test_mpfa_robin.py b/test/unit/test_mpfa_robin.py index 909ff06c65..ba88272bbc 100644 --- a/test/unit/test_mpfa_robin.py +++ b/test/unit/test_mpfa_robin.py @@ -13,7 +13,7 @@ def test_flow_left_right(self): g = pp.CartGrid([nx, ny], physdims=[1, 1]) g.compute_geometry() k = pp.SecondOrderTensor(2, np.ones(g.num_cells)) - alpha = 1.5 + robin_weight = 1.5 left = g.face_centers[0] < 1e-10 right = g.face_centers[0] > 1 - 1e-10 @@ -27,10 +27,10 @@ def test_flow_left_right(self): p_bound = 1 rob_bound = 0 - C = (rob_bound - alpha*p_bound) / (alpha - p_bound) + C = (rob_bound - robin_weight*p_bound) / (robin_weight - p_bound) u_ex = - C * g.face_normals[0] p_ex = C * g.cell_centers[0] + p_bound - self.solve_robin(g, k, bnd, alpha, p_bound, rob_bound, dir_ind, rob_ind, p_ex, u_ex) + self.solve_robin(g, k, bnd, robin_weight, p_bound, rob_bound, dir_ind, rob_ind, p_ex, u_ex) def test_flow_nz_rhs(self): nxs = [1, 3] @@ -40,7 +40,7 @@ def test_flow_nz_rhs(self): g = pp.CartGrid([nx, ny], physdims=[1, 1]) g.compute_geometry() k = pp.SecondOrderTensor(2, np.ones(g.num_cells)) - alpha = 1.5 + robin_weight = 1.5 left = g.face_centers[0] < 1e-10 right = g.face_centers[0] > 1 - 1e-10 @@ -55,10 +55,10 @@ def test_flow_nz_rhs(self): p_bound = 1 rob_bound = 1 - C = (rob_bound - alpha*p_bound) / (alpha - p_bound) + C = (rob_bound - robin_weight*p_bound) / (robin_weight - p_bound) u_ex = - C * g.face_normals[0] p_ex = C * g.cell_centers[0] + p_bound - self.solve_robin(g, k, bnd, alpha, p_bound, rob_bound, dir_ind, rob_ind, p_ex, u_ex) + self.solve_robin(g, k, bnd, robin_weight, p_bound, rob_bound, dir_ind, rob_ind, p_ex, u_ex) def test_flow_down(self): nxs = [1, 3] nys = [1, 3] @@ -67,7 +67,7 @@ def test_flow_down(self): g = pp.CartGrid([nx, ny], physdims=[1, 1]) g.compute_geometry() k = pp.SecondOrderTensor(2, np.ones(g.num_cells)) - alpha = 1.5 + robin_weight = 1.5 bot = g.face_centers[1] < 1e-10 top = g.face_centers[1] > 1 - 1e-10 @@ -81,16 +81,16 @@ def test_flow_down(self): p_bound = 1 rob_bound = 1 - C = (rob_bound - alpha*p_bound) / (alpha - p_bound) + C = (rob_bound - robin_weight*p_bound) / (robin_weight - p_bound) u_ex = C * g.face_normals[1] p_ex = C *(1- g.cell_centers[1]) + p_bound - self.solve_robin(g, k, bnd, alpha, p_bound, rob_bound, dir_ind, rob_ind, p_ex, u_ex) + self.solve_robin(g, k, bnd, robin_weight, p_bound, rob_bound, dir_ind, rob_ind, p_ex, u_ex) def test_no_neumann(self): g = pp.CartGrid([2, 2], physdims=[1, 1]) g.compute_geometry() k = pp.SecondOrderTensor(2, np.ones(g.num_cells)) - alpha = 1.5 + robin_weight = 1.5 bot = g.face_centers[1] < 1e-10 top = g.face_centers[1] > 1 - 1e-10 @@ -104,11 +104,11 @@ def test_no_neumann(self): bnd_ind = np.hstack((dir_ind, rob_ind)) bnd = pp.BoundaryCondition(g, bnd_ind, names) - flux, bound_flux, _, _ = pp.numerics.fv.mpfa._mpfa_local(g, k, bnd, alpha=alpha) + flux, bound_flux, _, _ = pp.numerics.fv.mpfa._mpfa_local(g, k, bnd, robin_weight=robin_weight) div = pp.fvutils.scalar_divergence(g) - rob_ex = [alpha * .25, alpha * .75, 1, 1] + rob_ex = [robin_weight * .25, robin_weight * .75, 1, 1] u_bound = np.zeros(g.num_faces) u_bound[dir_ind] = g.face_centers[1, dir_ind] u_bound[rob_ind] = rob_ex * g.face_areas[rob_ind] @@ -127,8 +127,8 @@ def test_no_neumann(self): assert np.allclose(p, p_ex) assert np.allclose(flux * p + bound_flux * u_bound, u_ex) - def solve_robin(self, g, k, bnd, alpha, p_bound, rob_bound, dir_ind, rob_ind, p_ex, u_ex): - flux, bound_flux, _, _ = pp.numerics.fv.mpfa._mpfa_local(g, k, bnd, alpha=alpha) + def solve_robin(self, g, k, bnd, robin_weight, p_bound, rob_bound, dir_ind, rob_ind, p_ex, u_ex): + flux, bound_flux, _, _ = pp.numerics.fv.mpfa._mpfa_local(g, k, bnd, robin_weight=robin_weight) div = pp.fvutils.scalar_divergence(g) diff --git a/test/unit/test_mpsa_robin.py b/test/unit/test_mpsa_robin.py index d784708d84..aaf3fd3f35 100644 --- a/test/unit/test_mpsa_robin.py +++ b/test/unit/test_mpsa_robin.py @@ -11,7 +11,7 @@ def test_dir_rob(self): g = pp.CartGrid([nx, ny], physdims=[1, 1]) g.compute_geometry() c = pp.FourthOrderTensor(2, np.ones(g.num_cells), np.ones(g.num_cells)) - alpha = np.pi + robin_weight = np.pi bot = g.face_centers[1] < 1e-10 top = g.face_centers[1] > 1 - 1e-10 @@ -45,9 +45,9 @@ def T_ex(faces): u_bound[:, neu_ind] = T_ex(neu_ind) * sgn_n u_bound[:, rob_ind] = ( T_ex(rob_ind) * sgn_r - + alpha * u_ex(g.face_centers[:, rob_ind]) * g.face_areas[rob_ind] + + robin_weight * u_ex(g.face_centers[:, rob_ind]) * g.face_areas[rob_ind] ) - u, T = self.solve_mpsa(g, c, alpha, bnd, u_bound) + u, T = self.solve_mpsa(g, c, robin_weight, bnd, u_bound) assert np.allclose(u, u_ex(g.cell_centers).ravel('F')) assert np.allclose(T, T_ex(np.arange(g.num_faces)).ravel('F')) @@ -58,7 +58,7 @@ def test_dir_neu_rob(self): g = pp.CartGrid([nx, ny], physdims=[1, 1]) g.compute_geometry() c = pp.FourthOrderTensor(2, np.ones(g.num_cells), np.ones(g.num_cells)) - alpha = np.pi + robin_weight = np.pi bot = g.face_centers[1] < 1e-10 top = g.face_centers[1] > 1 - 1e-10 @@ -92,9 +92,9 @@ def T_ex(faces): u_bound[:, neu_ind] = T_ex(neu_ind) * sgn_n u_bound[:, rob_ind] = ( T_ex(rob_ind) * sgn_r - + alpha * u_ex(g.face_centers[:, rob_ind]) * g.face_areas[rob_ind] + + robin_weight * u_ex(g.face_centers[:, rob_ind]) * g.face_areas[rob_ind] ) - u, T = self.solve_mpsa(g, c, alpha, bnd, u_bound) + u, T = self.solve_mpsa(g, c, robin_weight, bnd, u_bound) assert np.allclose(u, u_ex(g.cell_centers).ravel('F')) assert np.allclose(T, T_ex(np.arange(g.num_faces)).ravel('F')) @@ -105,7 +105,7 @@ def test_structured_triang(self): g = pp.StructuredTriangleGrid([nx, ny], physdims=[1, 1]) g.compute_geometry() c = pp.FourthOrderTensor(2, np.ones(g.num_cells), np.ones(g.num_cells)) - alpha = np.pi + robin_weight = np.pi bot = g.face_centers[1] < 1e-10 top = g.face_centers[1] > 1 - 1e-10 @@ -139,9 +139,9 @@ def T_ex(faces): u_bound[:, neu_ind] = T_ex(neu_ind) * sgn_n u_bound[:, rob_ind] = ( T_ex(rob_ind) * sgn_r - + alpha * u_ex(g.face_centers[:, rob_ind]) * g.face_areas[rob_ind] + + robin_weight * u_ex(g.face_centers[:, rob_ind]) * g.face_areas[rob_ind] ) - u, T = self.solve_mpsa(g, c, alpha, bnd, u_bound) + u, T = self.solve_mpsa(g, c, robin_weight, bnd, u_bound) assert np.allclose(u, u_ex(g.cell_centers).ravel('F')) assert np.allclose(T, T_ex(np.arange(g.num_faces)).ravel('F')) @@ -153,7 +153,7 @@ def test_unstruct_triang(self): g = pp.TriangleGrid(points) g.compute_geometry() c = pp.FourthOrderTensor(2, np.ones(g.num_cells), np.ones(g.num_cells)) - alpha = np.pi + robin_weight = np.pi bot = g.face_centers[1] < 1e-10 top = g.face_centers[1] > 1 - 1e-10 @@ -187,9 +187,9 @@ def T_ex(faces): u_bound[:, neu_ind] = T_ex(neu_ind) * sgn_n u_bound[:, rob_ind] = ( T_ex(rob_ind) * sgn_r - + alpha * u_ex(g.face_centers[:, rob_ind]) * g.face_areas[rob_ind] + + robin_weight * u_ex(g.face_centers[:, rob_ind]) * g.face_areas[rob_ind] ) - u, T = self.solve_mpsa(g, c, alpha, bnd, u_bound) + u, T = self.solve_mpsa(g, c, robin_weight, bnd, u_bound) assert np.allclose(u, u_ex(g.cell_centers).ravel('F')) assert np.allclose(T, T_ex(np.arange(g.num_faces)).ravel('F')) @@ -204,7 +204,7 @@ def test_unstruct_tetrahedron(self): g = pp.CartGrid([4,4,4], [1,1,1]) g.compute_geometry() c = pp.FourthOrderTensor(3, np.ones(g.num_cells), np.ones(g.num_cells)) - alpha = 1.0 + robin_weight = 1.0 bot = g.face_centers[2] < 1e-10 top = g.face_centers[2] > 1 - 1e-10 @@ -242,16 +242,16 @@ def T_ex(faces): u_bound[:, neu_ind] = T_ex(neu_ind) * sgn_n u_bound[:, rob_ind] = ( T_ex(rob_ind) * sgn_r - + alpha * u_ex(g.face_centers[:, rob_ind]) * g.face_areas[rob_ind] + + robin_weight * u_ex(g.face_centers[:, rob_ind]) * g.face_areas[rob_ind] ) - u, T = self.solve_mpsa(g, c, alpha, bnd, u_bound) + u, T = self.solve_mpsa(g, c, robin_weight, bnd, u_bound) assert np.allclose(u, u_ex(g.cell_centers).ravel('F')) assert np.allclose(T, T_ex(np.arange(g.num_faces)).ravel('F')) - def solve_mpsa(self, g, c, alpha, bnd, u_bound): - stress, bound_stress = pp.numerics.fv.mpsa._mpsa_local(g, c, bnd, alpha=alpha) + def solve_mpsa(self, g, c, robin_weight, bnd, u_bound): + stress, bound_stress = pp.numerics.fv.mpsa._mpsa_local(g, c, bnd, robin_weight=robin_weight) div = pp.fvutils.vector_divergence(g) a = div * stress b = -div * bound_stress * u_bound.ravel('F') From 308419eb547f0c633cc2417c616b624d66ca3cca Mon Sep 17 00:00:00 2001 From: Runar Date: Wed, 12 Sep 2018 10:20:53 +0200 Subject: [PATCH 18/26] Cleaned up the ExcludeBoundaries --- src/porepy/numerics/fv/fvutils.py | 113 +++++++++++++----------------- src/porepy/numerics/fv/mpfa.py | 16 ++--- src/porepy/numerics/fv/mpsa.py | 23 +++--- src/porepy/params/bc.py | 4 ++ 4 files changed, 70 insertions(+), 86 deletions(-) diff --git a/src/porepy/numerics/fv/fvutils.py b/src/porepy/numerics/fv/fvutils.py index 83b2ba6e3c..1ab3a3688b 100644 --- a/src/porepy/numerics/fv/fvutils.py +++ b/src/porepy/numerics/fv/fvutils.py @@ -671,58 +671,20 @@ def __init__(self, subcell_topology, bound, nd): fno = subcell_topology.fno_unique num_subfno = subcell_topology.num_subfno_unique + self.fno = fno + self.num_subfno = num_subfno + # Define mappings to exclude boundary values if self.bc_type == "scalar": - col_neu = np.argwhere([not it for it in bound.is_neu[fno]]) - row_neu = np.arange(col_neu.size) - self.exclude_neu = sps.coo_matrix( - (np.ones(row_neu.size), (row_neu, col_neu.ravel("C"))), - shape=(row_neu.size, num_subfno), - ).tocsr() - - col_dir = np.argwhere([not it for it in bound.is_dir[fno]]) - row_dir = np.arange(col_dir.size) - self.exclude_dir = sps.coo_matrix( - (np.ones(row_dir.size), (row_dir, col_dir.ravel("C"))), - shape=(row_dir.size, num_subfno), - ).tocsr() - - col_dir = np.argwhere([not it for it in bound.is_rob[fno]]) - row_dir = np.arange(col_dir.size) - self.exclude_rob = sps.coo_matrix( - (np.ones(row_dir.size), (row_dir, col_dir.ravel("C"))), - shape=(row_dir.size, num_subfno), - ).tocsr() - - col_neu_dir = np.argwhere([not it for it in bound.is_neu[fno] | bound.is_dir[fno]]) - row_neu_dir = np.arange(col_neu_dir.size) - self._exclude_neu_dir = sps.coo_matrix( - (np.ones(row_neu_dir.size), (row_neu_dir, col_neu_dir.ravel("C"))), - shape=(row_neu_dir.size, num_subfno), - ).tocsr() - - col_neu_rob = np.argwhere([not it for it in bound.is_neu[fno] | bound.is_rob[fno]]) - row_neu_rob = np.arange(col_neu_rob.size) - self._exclude_neu_rob = sps.coo_matrix( - (np.ones(row_neu_rob.size), (row_neu_rob, col_neu_rob.ravel("C"))), - shape=(row_neu_rob.size, num_subfno), - ).tocsr() - - col_rob_dir = np.argwhere([not it for it in bound.is_dir[fno] | bound.is_rob[fno]]) - row_rob_dir = np.arange(col_rob_dir.size) - self._exclude_rob_dir = sps.coo_matrix( - (np.ones(row_rob_dir.size), (row_rob_dir, col_rob_dir.ravel("C"))), - shape=(row_rob_dir.size, num_subfno), - ).tocsr() - - col_rob = np.argwhere([it for it in bound.is_rob[fno]]) - row_rob = np.arange(col_rob.size) - self._keep_robin = sps.coo_matrix( - (np.ones(row_rob.size), (row_rob, col_rob.ravel("C"))), - shape=(row_rob.size, num_subfno), - ).tocsr() + self.exclude_neu = self.__exclude_matrix(bound.is_neu) + self.exclude_dir = self.__exclude_matrix(bound.is_dir) + self.exclude_rob = self.__exclude_matrix(bound.is_rob) + self.exclude_neu_dir = self.__exclude_matrix(bound.is_neu | bound.is_dir) + self.exclude_neu_rob = self.__exclude_matrix(bound.is_neu | bound.is_rob) + self.exclude_rob_dir = self.__exclude_matrix(bound.is_rob | bound.is_dir) + self.keep_rob = self.__exclude_matrix(~bound.is_rob) elif self.bc_type == "vectorial": # Neumann @@ -794,6 +756,25 @@ def __init__(self, subcell_topology, bound, nd): shape=(row_dir.size, nd * num_subfno), ).tocsr() + def __exclude_matrix(self, ids): + """ + creates an exclusion matrix. This is a mapping from sub-faces to + all sub-faces except those given by ids. + Example: + ids = [0, 2] + self.num_subfno = 4 + print(sef.__exclude_matrix(ids)) + [[0, 1, 0, 0], + [0, 0, 0, 1]] + """ + col = np.argwhere([not it for it in ids[self.fno]]) + row = np.arange(col.size) + return sps.coo_matrix( + (np.ones(row.size), (row, col.ravel("C"))), + shape=(row.size, self.num_subfno), + ).tocsr() + + def exclude_dirichlet(self, other): """ Mapping to exclude faces/components with Dirichlet boundary conditions from local systems. @@ -917,7 +898,7 @@ def exclude_robin(self, other): return exclude_rob - def exclude_neu_rob(self, other): + def exclude_neumann_robin(self, other): """ Mapping to exclude faces/components with Neumann and Robin boundary conditions from local systems. @@ -931,13 +912,13 @@ def exclude_neu_rob(self, other): """ if self.bc_type == "scalar": - exclude_neu = self._exclude_neu_rob* other + exclude_neu = self.exclude_neu_rob* other elif self.bc_type == "vectorial": raise NotImplementedError() return exclude_neu - def exclude_neu_dir(self, other): + def exclude_neumann_dirichlet(self, other): """ Mapping to exclude faces/components with Neumann and Dirichlet boundary conditions from local systems. @@ -951,13 +932,13 @@ def exclude_neu_dir(self, other): """ if self.bc_type == "scalar": - exclude_neu = self._exclude_neu_dir * other + exclude_neu = self.exclude_neu_dir * other elif self.bc_type == "vectorial": raise NotImplementedError() return exclude_neu - def exclude_rob_dir(self, other): + def exclude_robin_dirichlet(self, other): """ Mapping to exclude faces/components with Robin and Dirichlet boundary conditions from local systems. @@ -971,7 +952,7 @@ def exclude_rob_dir(self, other): """ if self.bc_type == "scalar": - exclude_rob = self._exclude_rob_dir* other + exclude_rob = self.exclude_rob_dir* other elif self.bc_type == "vectorial": raise NotImplementedError() @@ -991,7 +972,7 @@ def keep_robin(self, other): """ if self.bc_type == "scalar": - exclude_rob = self._keep_robin * other + exclude_rob = self.keep_rob * other elif self.bc_type == "vectorial": raise NotImplementedError() @@ -1075,51 +1056,51 @@ def exclude_neumann_nd(self, other): return exclude_neumann_nd * other - def exclude_neu_rob_nd(self, other): + def exclude_neumann_robin_nd(self, other): """ Exclusion of Neumann and robin conditions for vector equations (elasticity). See above method without _nd suffix for description. """ if self.bc_type == "scalar": - exclude_neumann_nd = sps.kron(sps.eye(self.nd), self._exclude_neu_rob) + exclude_neu_rob_nd = sps.kron(sps.eye(self.nd), self.exclude_neu_rob) elif self.bc_type == "vectorial": raise NotImplementedError() - return exclude_neumann_nd * other + return exclude_neu_rob_nd * other - def exclude_neu_dir_nd(self, other): + def exclude_neumann_dirichlet_nd(self, other): """ Exclusion of Neumann and Dirichlet conditions for vector equations (elasticity). See above method without _nd suffix for description. """ if self.bc_type == "scalar": - exclude_neu_dir_nd = sps.kron(sps.eye(self.nd), self._exclude_neu_dir) + exclude_neu_dir_nd = sps.kron(sps.eye(self.nd), self.exclude_neu_dir) elif self.bc_type == "vectorial": raise NotImplementedError() return exclude_neu_dir_nd * other - def exclude_rob_dir_nd(self, other): - """ Exclusion of Roben and Dirichlet conditions for vector equations + def exclude_robin_dirichlet_nd(self, other): + """ Exclusion of Robin and Dirichlet conditions for vector equations (elasticity). See above method without _nd suffix for description. """ if self.bc_type == "scalar": - exclude_neu_dir_nd = sps.kron(sps.eye(self.nd), self._exclude_rob_dir) + exclude_rob_dir_nd = sps.kron(sps.eye(self.nd), self.exclude_rob_dir) elif self.bc_type == "vectorial": raise NotImplementedError() - return exclude_neu_dir_nd * other + return exclude_rob_dir_nd * other def keep_robin_nd(self, other): - """ Keep Roben conditions for vector equations (elasticity). + """ Keep Robin conditions for vector equations (elasticity). See above method without _nd suffix for description. """ if self.bc_type == "scalar": - exclude_neu_dir_nd = sps.kron(sps.eye(self.nd), self._keep_robin) + exclude_neu_dir_nd = sps.kron(sps.eye(self.nd), self.keep_rob) elif self.bc_type == "vectorial": raise NotImplementedError() diff --git a/src/porepy/numerics/fv/mpfa.py b/src/porepy/numerics/fv/mpfa.py index 33b37a190b..33cb42c6a0 100644 --- a/src/porepy/numerics/fv/mpfa.py +++ b/src/porepy/numerics/fv/mpfa.py @@ -624,8 +624,8 @@ def _mpfa_local(g, k, bnd, eta=None, inverter="numba", apertures=None, robin_wei bound_exclusion = fvutils.ExcludeBoundaries(subcell_topology, bnd, g.dim) # No flux conditions for Dirichlet boundary faces - nk_grad_n = bound_exclusion.exclude_rob_dir(nk_grad_all) - nk_cell = bound_exclusion.exclude_rob_dir(nk_cell) + nk_grad_n = bound_exclusion.exclude_robin_dirichlet(nk_grad_all) + nk_cell = bound_exclusion.exclude_robin_dirichlet(nk_cell) # Robin condition is only applied to Robin boundary faces nk_grad_r = bound_exclusion.keep_robin(nk_grad_all) @@ -634,8 +634,8 @@ def _mpfa_local(g, k, bnd, eta=None, inverter="numba", apertures=None, robin_wei del nk_grad_all # No pressure condition for Neumann or Robin boundary faces - pr_cont_grad = bound_exclusion.exclude_neu_rob(pr_cont_grad_all) - pr_cont_cell = bound_exclusion.exclude_neu_rob(pr_cont_cell_all) + pr_cont_grad = bound_exclusion.exclude_neumann_robin(pr_cont_grad_all) + pr_cont_cell = bound_exclusion.exclude_neumann_robin(pr_cont_cell_all) # So far, the local numbering has been based on the numbering scheme # implemented in SubcellTopology (which treats one cell at a time). For @@ -916,8 +916,8 @@ def _block_diagonal_structure(sub_cell_index, cell_node_blocks, nno, bound_exclu # Stack node numbers of equations on top of each other, and sort them to # get block-structure. First eliminate node numbers at the boundary, where # the equations are either of flux, pressure continuity or robin - nno_flux = bound_exclusion.exclude_rob_dir(nno) - nno_pressure = bound_exclusion.exclude_neu_rob(nno) + nno_flux = bound_exclusion.exclude_robin_dirichlet(nno) + nno_pressure = bound_exclusion.exclude_neumann_robin(nno) # we have now eliminated all nodes related to robin, we therefore add them nno_rob = bound_exclusion.keep_robin(nno) @@ -982,7 +982,7 @@ def _create_bound_rhs(bnd, bound_exclusion, subcell_topology, sgn, g, num_flux, # Find Neumann and Robin faces, exclude Dirichlet faces (since these are excluded # from the right hand side linear system), and do necessary formating. neu_ind = np.argwhere( - bound_exclusion.exclude_rob_dir(is_neu[fno].astype("int64")) + bound_exclusion.exclude_robin_dirichlet(is_neu[fno].astype("int64")) ).ravel("F") rob_ind = np.argwhere( bound_exclusion.keep_robin(is_rob[fno].astype("int64")) @@ -1027,7 +1027,7 @@ def _create_bound_rhs(bnd, bound_exclusion, subcell_topology, sgn, g, num_flux, # Dirichlet boundary conditions dir_ind = np.argwhere( - bound_exclusion.exclude_neu_rob(is_dir[fno].astype("int64")) + bound_exclusion.exclude_neumann_robin(is_dir[fno].astype("int64")) ).ravel("F") if dir_ind.size > 0: dir_cell = sps.coo_matrix( diff --git a/src/porepy/numerics/fv/mpsa.py b/src/porepy/numerics/fv/mpsa.py index f45147c700..bd6476b6dd 100644 --- a/src/porepy/numerics/fv/mpsa.py +++ b/src/porepy/numerics/fv/mpsa.py @@ -1,5 +1,4 @@ """ - Implementation of the multi-point stress appoximation method, and also terms related to poro-elastic coupling. @@ -959,7 +958,7 @@ def mpsa_elasticity(g, constit, subcell_topology, bound_exclusion, eta, inverter ncsym_all = subcell_topology.pair_over_subfaces_nd(ncsym_all + ncasym) del ncasym - ncsym = bound_exclusion.exclude_rob_dir_nd(ncsym_all) + ncsym = bound_exclusion.exclude_robin_dirichlet_nd(ncsym_all) num_subfno = subcell_topology.subfno.max() + 1 hook_cell = sps.coo_matrix( @@ -967,7 +966,7 @@ def mpsa_elasticity(g, constit, subcell_topology, bound_exclusion, eta, inverter shape=(num_subfno * nd, (np.max(subcell_topology.cno) + 1) * nd), ).tocsr() - hook_cell = bound_exclusion.exclude_rob_dir_nd(hook_cell) + hook_cell = bound_exclusion.exclude_robin_dirichlet_nd(hook_cell) # For the Robin boundary conditions we need to pair the forces with the # displacement @@ -1069,8 +1068,8 @@ def __get_displacement_submatrices( # Expand equations for displacement balance, and eliminate rows # associated with neumann boundary conditions - d_cont_grad = bound_exclusion.exclude_neu_rob_nd(d_cont_grad) - d_cont_cell = bound_exclusion.exclude_neu_rob_nd(d_cont_cell) + d_cont_grad = bound_exclusion.exclude_neumann_robin_nd(d_cont_grad) + d_cont_cell = bound_exclusion.exclude_neumann_robin_nd(d_cont_cell) # The column ordering of the displacement equilibrium equations are # formed as a Kronecker product of scalar equations. Bring them to the @@ -1377,8 +1376,8 @@ def _block_diagonal_structure( # get block-structure. First eliminate node numbers at the boundary, where # the equations are either of flux or pressure continuity (not both) - nno_stress = bound_exclusion.exclude_rob_dir(nno) - nno_displacement = bound_exclusion.exclude_neu_rob(nno) + nno_stress = bound_exclusion.exclude_robin_dirichlet(nno) + nno_displacement = bound_exclusion.exclude_neumann_robin(nno) # we have now eliminated all nodes related to robin, we therefore add them nno_rob = bound_exclusion.keep_robin(nno) @@ -1431,9 +1430,9 @@ def create_bound_rhs(bound, bound_exclusion, subcell_topology, g): basis functions for boundary values """ nd = g.dim - num_stress = bound_exclusion._exclude_rob_dir.shape[0] * nd - num_displ = bound_exclusion._exclude_neu_rob.shape[0] * nd - num_rob = bound_exclusion._keep_robin.shape[0] * nd + num_stress = bound_exclusion.exclude_rob_dir.shape[0] * nd + num_displ = bound_exclusion.exclude_neu_rob.shape[0] * nd + num_rob = bound_exclusion.keep_rob.shape[0] * nd fno = subcell_topology.fno_unique subfno = subcell_topology.subfno_unique @@ -1458,7 +1457,7 @@ def expand_ind(ind, dim, increment): # Define right hand side for Neumann boundary conditions # First row indices in rhs matrix - is_neu = bound_exclusion.exclude_rob_dir(bound.is_neu[fno].astype("int64")) + is_neu = bound_exclusion.exclude_robin_dirichlet(bound.is_neu[fno].astype("int64")) neu_ind_single = np.argwhere(is_neu).ravel("F") is_rob = bound_exclusion.keep_robin(bound.is_rob[fno].astype("int64")) rob_ind_single = np.argwhere(is_rob).ravel("F") @@ -1510,7 +1509,7 @@ def expand_ind(ind, dim, increment): neu_cell = sps.coo_matrix((num_stress + num_rob, num_bound)).tocsr() # Dirichlet boundary conditions, procedure is similar to that for Neumann - is_dir = bound_exclusion.exclude_neu_rob(bound.is_dir[fno].astype("int64")) + is_dir = bound_exclusion.exclude_neumann_robin(bound.is_dir[fno].astype("int64")) dir_ind_single = np.argwhere(is_dir).ravel("F") dir_ind = expand_ind(dir_ind_single, nd, is_dir.size) diff --git a/src/porepy/params/bc.py b/src/porepy/params/bc.py index 6347c02a38..79feda5e65 100644 --- a/src/porepy/params/bc.py +++ b/src/porepy/params/bc.py @@ -265,6 +265,10 @@ def set_bc(self, faces, cond): elif s.lower() == "dir": self.is_dir[:, faces[j]] = True self.is_neu[:, faces[j]] = False + elif s.lower() == "rob": + self.is_rob[:, faces[j]] = True + self.is_neu[:, faces[j]] = False + self.is_dir[:, faces[j]] = False elif s.lower() == "dir_x": self.is_dir[0, faces[j]] = True self.is_neu[0, faces[j]] = False From 34da31a994e1d7578ba99302ad78f7fdfce6e82c Mon Sep 17 00:00:00 2001 From: Runar Date: Wed, 12 Sep 2018 11:10:36 +0200 Subject: [PATCH 19/26] Formated code with black --- .../papers/arXiv_1802_05961/example_1_data.py | 4 +- .../papers/arXiv_1802_05961/example_2_data.py | 9 +- examples/papers/arXiv_1802_05961/solvers.py | 1 + src/porepy/grids/gmsh/gmsh_interface.py | 2 +- src/porepy/grids/grid.py | 1 - src/porepy/numerics/fem/p1.py | 1 + src/porepy/numerics/fem/source.py | 1 + src/porepy/numerics/fv/fvutils.py | 12 +-- src/porepy/numerics/fv/mpfa.py | 61 +++++++++---- src/porepy/numerics/fv/mpsa.py | 88 +++++++++++++------ src/porepy/numerics/mixed_dim/coupler.py | 6 +- src/porepy/params/bc.py | 5 +- src/porepy/params/data.py | 1 - src/porepy/utils/tags.py | 1 + src/porepy/viz/exporter.py | 80 ++++++++++------- test/unit/test_mpfa_robin.py | 83 ++++++++++++----- test/unit/test_mpsa_robin.py | 82 +++++++++-------- test/unit/test_vtk.py | 19 ++-- 18 files changed, 287 insertions(+), 170 deletions(-) diff --git a/examples/papers/arXiv_1802_05961/example_1_data.py b/examples/papers/arXiv_1802_05961/example_1_data.py index d0b6c4c418..0a9d9740ad 100644 --- a/examples/papers/arXiv_1802_05961/example_1_data.py +++ b/examples/papers/arXiv_1802_05961/example_1_data.py @@ -29,7 +29,7 @@ def create_gb(cells_2d, alpha_1d=None, alpha_mortar=None): mesh_kwargs = { "tol": tol(), "mesh_size_frac": mesh_size, - "mesh_size_min": mesh_size / 20 + "mesh_size_min": mesh_size / 20, } mesh_kwargs = {"mesh_size_frac": mesh_size, "mesh_size_min": mesh_size / 20} @@ -99,7 +99,7 @@ def add_data(gb, solver): param.set_source("flow", np.zeros(g.num_cells)) - if if_p1: # for P1 a different handling of the boundary conditions + if if_p1: # for P1 a different handling of the boundary conditions bound_nodes = g.get_boundary_nodes() if bound_nodes.size == 0: bc = pp.BoundaryConditionNode(g, np.empty(0), np.empty(0)) diff --git a/examples/papers/arXiv_1802_05961/example_2_data.py b/examples/papers/arXiv_1802_05961/example_2_data.py index ee066a598b..76e3914b19 100644 --- a/examples/papers/arXiv_1802_05961/example_2_data.py +++ b/examples/papers/arXiv_1802_05961/example_2_data.py @@ -33,7 +33,7 @@ def str_of_nodes(n): file_csv = "geiger_3d.csv" file_geo = geo_file_name(h, False) - print("Create grid "+file_geo) + print("Create grid " + file_geo) _, _, domain = importer.network_3d_from_csv(file_csv, tol=tol()) network = pickle.load(open("geiger_3d_network", "rb")) gb = importer.dfm_from_gmsh(file_geo, 3, network=network, tol=tol()) @@ -146,7 +146,7 @@ def add_data(gb, domain, solver, case): param.set_source("flow", np.zeros(g.num_cells)) - if if_p1: # for P1 a different handling of the boundary conditions + if if_p1: # for P1 a different handling of the boundary conditions bound_nodes = g.get_boundary_nodes() if bound_nodes.size == 0: bc = pp.BoundaryConditionNode(g, np.empty(0), np.empty(0)) @@ -226,6 +226,7 @@ def b_pressure(g): # ------------------------------------------------------------------------------# + def b_pressure_node(g): b_nodes = g.get_boundary_nodes() @@ -236,9 +237,7 @@ def b_pressure_node(g): b_node_coords = g.nodes[:, b_nodes] val = 0.4 + tol() - b_in = np.logical_and.reduce( - tuple(b_node_coords[i, :] < val for i in range(3)) - ) + b_in = np.logical_and.reduce(tuple(b_node_coords[i, :] < val for i in range(3))) val = 0.8 - tol() b_out = np.logical_and.reduce( diff --git a/examples/papers/arXiv_1802_05961/solvers.py b/examples/papers/arXiv_1802_05961/solvers.py index 216a3fd81b..ee301792a8 100644 --- a/examples/papers/arXiv_1802_05961/solvers.py +++ b/examples/papers/arXiv_1802_05961/solvers.py @@ -98,6 +98,7 @@ def solve_p1(gb, folder, return_only_matrix=False): save = Exporter(gb, "sol", folder=folder, simplicial=True) save.write_vtk("pressure", point_data=True) + # ------------------------------------------------------------------------------# diff --git a/src/porepy/grids/gmsh/gmsh_interface.py b/src/porepy/grids/gmsh/gmsh_interface.py index fde5cb3636..7600dc1d18 100644 --- a/src/porepy/grids/gmsh/gmsh_interface.py +++ b/src/porepy/grids/gmsh/gmsh_interface.py @@ -557,7 +557,7 @@ def run_gmsh(in_file, out_file, dims, **kwargs): """ if not os.path.isfile(in_file): - raise FileNotFoundError("file "+in_file+" not found") + raise FileNotFoundError("file " + in_file + " not found") # Import config file to get location of gmsh executable. config = read_config.read() diff --git a/src/porepy/grids/grid.py b/src/porepy/grids/grid.py index 63982a1d7b..36025b0093 100644 --- a/src/porepy/grids/grid.py +++ b/src/porepy/grids/grid.py @@ -649,7 +649,6 @@ def update_boundary_node_tag(self): nodes = self.face_nodes.indices[mcolon.mcolon(first, second)] self.tags[node_tag][nodes] = True - def cell_diameters(self, cn=None): """ Compute the cell diameters. If self.dim == 0, return 0 diff --git a/src/porepy/numerics/fem/p1.py b/src/porepy/numerics/fem/p1.py index 363e049616..03beb6c6d4 100644 --- a/src/porepy/numerics/fem/p1.py +++ b/src/porepy/numerics/fem/p1.py @@ -267,6 +267,7 @@ def stiffH1(self, K, c_volume, coord, dim): # ------------------------------------------------------------------------------# + class P1Coupling(AbstractCoupling): # ------------------------------------------------------------------------------# diff --git a/src/porepy/numerics/fem/source.py b/src/porepy/numerics/fem/source.py index ccf7b1d171..c26a1bd048 100644 --- a/src/porepy/numerics/fem/source.py +++ b/src/porepy/numerics/fem/source.py @@ -55,4 +55,5 @@ def matrix_rhs(self, g, data): return lhs, M.dot(param.get_source(self)) + # ------------------------------------------------------------------------------# diff --git a/src/porepy/numerics/fv/fvutils.py b/src/porepy/numerics/fv/fvutils.py index 1ab3a3688b..254b3f50fa 100644 --- a/src/porepy/numerics/fv/fvutils.py +++ b/src/porepy/numerics/fv/fvutils.py @@ -685,7 +685,7 @@ def __init__(self, subcell_topology, bound, nd): self.exclude_neu_rob = self.__exclude_matrix(bound.is_neu | bound.is_rob) self.exclude_rob_dir = self.__exclude_matrix(bound.is_rob | bound.is_dir) self.keep_rob = self.__exclude_matrix(~bound.is_rob) - + elif self.bc_type == "vectorial": # Neumann col_neu_x = np.argwhere([not it for it in bound.is_neu[0, fno]]) @@ -774,7 +774,6 @@ def __exclude_matrix(self, ids): shape=(row.size, self.num_subfno), ).tocsr() - def exclude_dirichlet(self, other): """ Mapping to exclude faces/components with Dirichlet boundary conditions from local systems. @@ -894,10 +893,10 @@ def exclude_robin(self, other): exclude_rob = self.exclude_rob * other elif self.bc_type == "vectorial": - raise NotImplementedError('can not exclude robin for vectorial bc') + raise NotImplementedError("can not exclude robin for vectorial bc") return exclude_rob - + def exclude_neumann_robin(self, other): """ Mapping to exclude faces/components with Neumann and Robin boundary conditions from local systems. @@ -912,7 +911,7 @@ def exclude_neumann_robin(self, other): """ if self.bc_type == "scalar": - exclude_neu = self.exclude_neu_rob* other + exclude_neu = self.exclude_neu_rob * other elif self.bc_type == "vectorial": raise NotImplementedError() @@ -952,7 +951,7 @@ def exclude_robin_dirichlet(self, other): """ if self.bc_type == "scalar": - exclude_rob = self.exclude_rob_dir* other + exclude_rob = self.exclude_rob_dir * other elif self.bc_type == "vectorial": raise NotImplementedError() @@ -1107,6 +1106,7 @@ def keep_robin_nd(self, other): return exclude_neu_dir_nd * other + # -----------------End of class ExcludeBoundaries----------------------------- diff --git a/src/porepy/numerics/fv/mpfa.py b/src/porepy/numerics/fv/mpfa.py index 33cb42c6a0..be526d51b1 100644 --- a/src/porepy/numerics/fv/mpfa.py +++ b/src/porepy/numerics/fv/mpfa.py @@ -178,7 +178,17 @@ def discretize(self, g, data): # ------------------------------------------------------------------------------# -def mpfa(g, k, bnd, robin_weight=None, eta=None, inverter=None, apertures=None, max_memory=None, **kwargs): +def mpfa( + g, + k, + bnd, + robin_weight=None, + eta=None, + inverter=None, + apertures=None, + max_memory=None, + **kwargs +): """ Discretize the scalar elliptic equation by the multi-point flux approximation method. @@ -271,7 +281,13 @@ def mpfa(g, k, bnd, robin_weight=None, eta=None, inverter=None, apertures=None, # TODO: We may want to estimate the memory need, and give a warning if # this seems excessive flux, bound_flux, bound_pressure_cell, bound_pressure_face = _mpfa_local( - g, k, bnd, eta=eta, inverter=inverter, apertures=apertures, robin_weight=robin_weight + g, + k, + bnd, + eta=eta, + inverter=inverter, + apertures=apertures, + robin_weight=robin_weight, ) else: # Estimate number of partitions necessary based on prescribed memory @@ -458,7 +474,9 @@ def mpfa_partial( ) -def _mpfa_local(g, k, bnd, eta=None, inverter="numba", apertures=None, robin_weight=None): +def _mpfa_local( + g, k, bnd, eta=None, inverter="numba", apertures=None, robin_weight=None +): """ Actual implementation of the MPFA O-method. To calculate MPFA on a grid directly, either call this method, or, to respect the privacy of this @@ -520,8 +538,10 @@ def _mpfa_local(g, k, bnd, eta=None, inverter="numba", apertures=None, robin_wei return sps.csr_matrix([0]), 0, 0, 0 if robin_weight is None: - if np.sum(bnd.is_rob) !=0: - raise ValueError('If applying Robin conditions you must supply an robin_weight') + if np.sum(bnd.is_rob) != 0: + raise ValueError( + "If applying Robin conditions you must supply an robin_weight" + ) else: robin_weight = 1 # The grid coordinates are always three-dimensional, even if the grid is @@ -593,14 +613,20 @@ def _mpfa_local(g, k, bnd, eta=None, inverter="numba", apertures=None, robin_wei num_nodes = np.diff(g.face_nodes.indptr) sgn = g.cell_faces[subcell_topology.fno_unique, subcell_topology.cno_unique].A scaled_sgn = ( - robin_weight * sgn[0] * g.face_areas[subcell_topology.fno_unique] \ + robin_weight + * sgn[0] + * g.face_areas[subcell_topology.fno_unique] / num_nodes[subcell_topology.fno_unique] ) # pair_over_subfaces flips the sign so we flip it back pr_trace_grad_all = sps.diags(scaled_sgn) * pr_cont_grad_all pr_trace_cell_all = sps.coo_matrix( - (robin_weight * g.face_areas[subcell_topology.fno] / num_nodes[subcell_topology.fno], - (subcell_topology.subfno, subcell_topology.cno)) + ( + robin_weight + * g.face_areas[subcell_topology.fno] + / num_nodes[subcell_topology.fno], + (subcell_topology.subfno, subcell_topology.cno), + ) ).tocsr() del sgn, scaled_sgn @@ -654,7 +680,7 @@ def _mpfa_local(g, k, bnd, eta=None, inverter="numba", apertures=None, robin_wei # NOTE: I think in the discretization for sub_cells a flow out of the cell is # negative. This is a contradiction to what is done for the boundary conditions # where we want to set dot(n, flux) where n is the normal pointing outwards. - # thats why we need +nk_grad_r - pr_trace_grad -pr_trace_cell instead of = rhs + # thats why we need +nk_grad_r - pr_trace_grad -pr_trace_cell instead of = rhs # instead of how we would expect: -nk_grad_r + pr_trace_grad +pr_trace_cell= rhs. # This is also why we multiply with -1 in scaled_sgn in _create_bound_rhs grad_eqs = sps.vstack([nk_grad_n, nk_grad_r - pr_trace_grad, pr_cont_grad]) @@ -662,8 +688,7 @@ def _mpfa_local(g, k, bnd, eta=None, inverter="numba", apertures=None, robin_wei num_nk_cell = nk_cell.shape[0] num_nk_rob = nk_grad_r.shape[0] num_pr_cont_grad = pr_cont_grad.shape[0] - - + del nk_grad_n, nk_grad_r, pr_trace_grad grad = rows2blk_diag * grad_eqs * cols2blk_diag @@ -888,7 +913,7 @@ def _tensor_vector_prod(g, k, subcell_topology, apertures=None): ) nk = normals_mat * k_mat - + # Unique sub-cell indexes are pulled from column indices, we only need # every nd column (since nd faces of the cell meet at each vertex) sub_cell_ind = j[::, 0::nd] @@ -941,7 +966,9 @@ def _block_diagonal_structure(sub_cell_index, cell_node_blocks, nno, bound_exclu return rows2blk_diag, cols2blk_diag, size_of_blocks -def _create_bound_rhs(bnd, bound_exclusion, subcell_topology, sgn, g, num_flux, num_rob, num_pr): +def _create_bound_rhs( + bnd, bound_exclusion, subcell_topology, sgn, g, num_flux, num_rob, num_pr +): """ Define rhs matrix to get basis functions for incorporates boundary conditions @@ -996,14 +1023,14 @@ def _create_bound_rhs(bnd, bound_exclusion, subcell_topology, sgn, g, num_flux, # not have Dirichlet and Neumann excluded (respectively), and thus we need # new fields neu_ind_all = np.argwhere(is_neu[fno].astype("int")).ravel("F") - rob_ind_all = np.argwhere(is_rob[fno].astype("int")).ravel("F") + rob_ind_all = np.argwhere(is_rob[fno].astype("int")).ravel("F") dir_ind_all = np.argwhere(is_dir[fno].astype("int")).ravel("F") num_face_nodes = g.face_nodes.sum(axis=0).A.ravel(order="F") # We now merge the neuman and robin indices since they are treated equivalent - if rob_ind.size==0: + if rob_ind.size == 0: neu_rob_ind = neu_ind - elif neu_ind.size==0: + elif neu_ind.size == 0: neu_rob_ind = rob_ind + num_flux else: neu_rob_ind = np.hstack((neu_ind, rob_ind + num_flux)) @@ -1018,7 +1045,7 @@ def _create_bound_rhs(bnd, bound_exclusion, subcell_topology, sgn, g, num_flux, if neu_rob_ind.size > 0: neu_rob_cell = sps.coo_matrix( (scaled_sgn, (neu_rob_ind, np.arange(neu_rob_ind.size))), - shape=(num_flux+num_rob, num_bound), + shape=(num_flux + num_rob, num_bound), ) else: # Special handling when no elements are found. Not sure if this is diff --git a/src/porepy/numerics/fv/mpsa.py b/src/porepy/numerics/fv/mpsa.py index bd6476b6dd..399da77f82 100644 --- a/src/porepy/numerics/fv/mpsa.py +++ b/src/porepy/numerics/fv/mpsa.py @@ -133,6 +133,7 @@ def rhs(self, g, bound_stress, bc_val, f): return -div * bound_stress * bc_val - f + class FracturedMpsa(Mpsa): """ Subclass of MPSA for discretizing a fractured domain. Adds DOFs on each @@ -375,8 +376,8 @@ def discretize_fractures(self, g, data, faces=None, **kwargs): bound.is_neu[:, frac_faces] = False else: raise ValueError("Unknow boundary condition type: " + bound.bc_type) - if np.sum(bound.is_dir * bound.is_neu) !=0: - raise AssertionError('Found faces that are both dirichlet and neuman') + if np.sum(bound.is_dir * bound.is_neu) != 0: + raise AssertionError("Found faces that are both dirichlet and neuman") # Discretize with normal mpsa self.discretize(g, data, **kwargs) stress, bound_stress = data["stress"], data["bound_stress"] @@ -523,7 +524,16 @@ def given_slip_distance(self, g, stress, bound_stress, faces=None): # ------------------------------------------------------------------------------# -def mpsa(g, constit, bound, robin_weight=None, eta=None, inverter=None, max_memory=None, **kwargs): +def mpsa( + g, + constit, + bound, + robin_weight=None, + eta=None, + inverter=None, + max_memory=None, + **kwargs +): """ Discretize the vector elliptic equation by the multi-point stress approximation method, specifically the weakly symmetric MPSA-W method. @@ -658,8 +668,13 @@ def mpsa(g, constit, bound, robin_weight=None, eta=None, inverter=None, max_memo # Perform local discretization. loc_stress, loc_bound_stress, loc_faces = mpsa_partial( - g, constit, bound, eta=eta, inverter=inverter, nodes=active_nodes, - robin_weight=robin_weight + g, + constit, + bound, + eta=eta, + inverter=inverter, + nodes=active_nodes, + robin_weight=robin_weight, ) # Eliminate contribution from faces already covered @@ -685,7 +700,7 @@ def mpsa_partial( faces=None, nodes=None, apertures=None, - robin_weight=None + robin_weight=None, ): """ Run an MPFA discretization on subgrid, and return discretization in terms @@ -870,7 +885,13 @@ def _mpsa_local(g, constit, bound, eta=0, robin_weight=None, inverter="numba"): # elasticity and poro-elasticity). hook, igrad, rhs_cells, _, _ = mpsa_elasticity( - g, constit, subcell_topology, bound_exclusion, eta, inverter, robin_weight=robin_weight + g, + constit, + subcell_topology, + bound_exclusion, + eta, + inverter, + robin_weight=robin_weight, ) hook_igrad = hook * igrad @@ -898,7 +919,9 @@ def _mpsa_local(g, constit, bound, eta=0, robin_weight=None, inverter="numba"): return stress, bound_stress -def mpsa_elasticity(g, constit, subcell_topology, bound_exclusion, eta, inverter, robin_weight=None): +def mpsa_elasticity( + g, constit, subcell_topology, bound_exclusion, eta, inverter, robin_weight=None +): """ This is the function where the real discretization takes place. It contains the parts that are common for elasticity and poro-elasticity, and was thus @@ -973,8 +996,10 @@ def mpsa_elasticity(g, constit, subcell_topology, bound_exclusion, eta, inverter ncsym_rob = bound_exclusion.keep_robin_nd(ncsym_all) if robin_weight is None: - if ncsym_rob.shape[0] !=0: - raise ValueError('If applying Robin conditions you must supply an robin_weight') + if ncsym_rob.shape[0] != 0: + raise ValueError( + "If applying Robin conditions you must supply an robin_weight" + ) else: robin_weight = 1 @@ -1080,6 +1105,7 @@ def __get_displacement_submatrices( return d_cont_grad, d_cont_cell + def __get_displacement_submatrices_rob( g, subcell_topology, eta, num_sub_cells, bound_exclusion, robin_weight ): @@ -1088,22 +1114,28 @@ def __get_displacement_submatrices_rob( # contribution from gradient unknown to equations for displacement # at the boundary rob_grad = fvutils.compute_dist_face_cell(g, subcell_topology, eta) - + # For the Robin condition the distance from the cell centers to face centers # will be the contribution from the gradients. We integrate over the subface # and multiply by the area num_nodes = np.diff(g.face_nodes.indptr) sgn = g.cell_faces[subcell_topology.fno_unique, subcell_topology.cno_unique].A scaled_sgn = ( - robin_weight * sgn[0] * g.face_areas[subcell_topology.fno_unique] \ + robin_weight + * sgn[0] + * g.face_areas[subcell_topology.fno_unique] / num_nodes[subcell_topology.fno_unique] - ) + ) # pair_over_subfaces flips the sign so we flip it back - rob_grad = sps.kron(sps.eye(nd), sps.diags(scaled_sgn)*rob_grad) + rob_grad = sps.kron(sps.eye(nd), sps.diags(scaled_sgn) * rob_grad) # Contribution from cell center potentials to local systems rob_cell = sps.coo_matrix( - (robin_weight * g.face_areas[subcell_topology.fno] / num_nodes[subcell_topology.fno], - (subcell_topology.subfno, subcell_topology.cno)) + ( + robin_weight + * g.face_areas[subcell_topology.fno] + / num_nodes[subcell_topology.fno], + (subcell_topology.subfno, subcell_topology.cno), + ) ).tocsr() rob_cell = sps.kron(sps.eye(nd), rob_cell) @@ -1382,13 +1414,17 @@ def _block_diagonal_structure( nno_rob = bound_exclusion.keep_robin(nno) if bound_exclusion.bc_type == "scalar": - node_occ = np.hstack((np.tile(nno_stress, nd), - np.tile(nno_rob, nd), - np.tile(nno_displacement, nd))) + node_occ = np.hstack( + ( + np.tile(nno_stress, nd), + np.tile(nno_rob, nd), + np.tile(nno_displacement, nd), + ) + ) elif bound_exclusion.bc_type == "vectorial": node_occ = np.hstack((nno_stress, nno_displacement)) - if nno_rob.size >=0: + if nno_rob.size >= 0: raise NotImplementedError() sorted_ind = np.argsort(node_occ, kind="mergesort") rows2blk_diag = sps.coo_matrix( @@ -1460,24 +1496,24 @@ def expand_ind(ind, dim, increment): is_neu = bound_exclusion.exclude_robin_dirichlet(bound.is_neu[fno].astype("int64")) neu_ind_single = np.argwhere(is_neu).ravel("F") is_rob = bound_exclusion.keep_robin(bound.is_rob[fno].astype("int64")) - rob_ind_single = np.argwhere(is_rob).ravel("F") + rob_ind_single = np.argwhere(is_rob).ravel("F") # There are is_neu.size Neumann conditions per dimension and # there are is_rob.size Robin conditions per dimension neu_ind = expand_ind(neu_ind_single, nd, is_neu.size) rob_ind = expand_ind(rob_ind_single, nd, is_rob.size) # We now merge the neuman and robin indices since they are treated equivalent - if rob_ind.size==0: + if rob_ind.size == 0: neu_rob_ind = neu_ind - elif neu_ind.size==0: + elif neu_ind.size == 0: neu_rob_ind = rob_ind + num_stress else: neu_rob_ind = np.hstack((neu_ind, rob_ind + num_stress)) - + # We also need to account for all half faces, that is, do not exclude # Dirichlet and Neumann boundaries. neu_ind_single_all = np.argwhere(bound.is_neu[fno].astype("int")).ravel("F") - rob_ind_single_all = np.argwhere(bound.is_rob[fno].astype("int")).ravel("F") + rob_ind_single_all = np.argwhere(bound.is_rob[fno].astype("int")).ravel("F") dir_ind_single_all = np.argwhere(bound.is_dir[fno].astype("int")).ravel("F") neu_rob_ind_single_all = np.hstack((neu_ind_single_all, rob_ind_single_all)) @@ -1500,7 +1536,7 @@ def expand_ind(ind, dim, increment): if neu_rob_ind.size > 0: neu_cell = sps.coo_matrix( (neu_val.ravel("F"), (neu_rob_ind, np.arange(neu_rob_ind.size))), - shape=(num_stress+num_rob, num_bound), + shape=(num_stress + num_rob, num_bound), ).tocsr() else: diff --git a/src/porepy/numerics/mixed_dim/coupler.py b/src/porepy/numerics/mixed_dim/coupler.py index ac3aed3c42..357985d835 100644 --- a/src/porepy/numerics/mixed_dim/coupler.py +++ b/src/porepy/numerics/mixed_dim/coupler.py @@ -45,15 +45,15 @@ def __init__(self, discr=None, coupling=None, **kwargs): if not isinstance(coupling_fct, list): coupling_fct = [coupling_fct] self.coupling_fct = [c for c in coupling_fct] - else: # we assign coupling as coupling function + else: # we assign coupling as coupling function # Make sure coupling is list if not isinstance(coupling, list): coupling = [coupling] coupling_matrix_rhs = [] for c in coupling: - if c is None: # we assign None coupling + if c is None: # we assign None coupling coupling_matrix_rhs.append(c) - else: # Assign the matrix rhs coupling function + else: # Assign the matrix rhs coupling function coupling_matrix_rhs.append(c.matrix_rhs) self.coupling_fct = [c for c in coupling_matrix_rhs] diff --git a/src/porepy/params/bc.py b/src/porepy/params/bc.py index 79feda5e65..493b8a54fd 100644 --- a/src/porepy/params/bc.py +++ b/src/porepy/params/bc.py @@ -64,7 +64,7 @@ def __init__(self, g, faces=None, cond=None): self.is_neu = np.zeros(self.num_faces, dtype=bool) self.is_dir = np.zeros(self.num_faces, dtype=bool) - self.is_rob = np.zeros(self.num_faces, dtype=bool) + self.is_rob = np.zeros(self.num_faces, dtype=bool) # By default, all faces are Neumann. self.is_neu[bf] = True @@ -108,10 +108,11 @@ def __init__(self, g, faces=None, cond=None): elif s.lower() == "rob": self.is_dir[faces[l]] = False self.is_neu[faces[l]] = False - self.is_rob[faces[l]] = True + self.is_rob[faces[l]] = True else: raise ValueError("Boundary should be Dirichlet, Neumann or Robin") + class BoundaryConditionNode(object): """ Class to store information on boundary conditions for nodal numerical diff --git a/src/porepy/params/data.py b/src/porepy/params/data.py index 5505268867..5fef51112c 100644 --- a/src/porepy/params/data.py +++ b/src/porepy/params/data.py @@ -766,7 +766,6 @@ def get_bc_val_mechanics(self): bc_val_mechanics = property(get_bc_val_mechanics) - # ------------------------------------ def get_robin_factor(self, default=1): """ double or array-like diff --git a/src/porepy/utils/tags.py b/src/porepy/utils/tags.py index 677fd8b54c..ec471ade08 100644 --- a/src/porepy/utils/tags.py +++ b/src/porepy/utils/tags.py @@ -39,6 +39,7 @@ def standard_face_tags(): """ return ["fracture_faces", "tip_faces", "domain_boundary_faces"] + def standard_node_tags(): """ Returns the standard node tag key. diff --git a/src/porepy/viz/exporter.py b/src/porepy/viz/exporter.py index 5b1efdad5d..4ea47e91a6 100644 --- a/src/porepy/viz/exporter.py +++ b/src/porepy/viz/exporter.py @@ -26,11 +26,13 @@ # ------------------------------------------------------------------------------# + class Field(object): """ Internal class to store information for the data to export. """ - def __init__(self, name, cell_data = False, point_data = False, values = None): + + def __init__(self, name, cell_data=False, point_data=False, values=None): assert np.logical_xor(cell_data, point_data) # name of the field self.name = name @@ -82,12 +84,7 @@ def check(self, values, g): ) num_elem = g.num_cells if self.cell_data else g.num_nodes if np.atleast_2d(values).shape[1] != num_elem: - raise ValueError( - "Field " - + str(self.name) - + " has wrong dimension." - ) - + raise ValueError("Field " + str(self.name) + " has wrong dimension.") def set_values(self, values): """ @@ -100,12 +97,15 @@ def _check_values(self): if self.values is None: raise ValueError("Field " + str(self.name) + " values not valid") + # ------------------------------------------------------------------------------# + class Fields(object): """ Internal class to store a list of field. """ + def __init__(self): self.fields = None @@ -147,8 +147,10 @@ def names(self): """ return [f.name for f in self] + # ------------------------------------------------------------------------------# + class Exporter: def __init__(self, grid, name, folder=None, **kwargs): """ @@ -325,18 +327,24 @@ def _export_vtk_single(self, data, time_step, point_data): fields = Fields() if len(data) > 0: if point_data: - fields.extend([Field(n, point_data = True, values = v) - for n, v in data.items()]) + fields.extend( + [Field(n, point_data=True, values=v) for n, v in data.items()] + ) else: - fields.extend([Field(n, cell_data = True, values = v) - for n, v in data.items()]) + fields.extend( + [Field(n, cell_data=True, values=v) for n, v in data.items()] + ) - fields.extend([ - Field("grid_dim", cell_data = True, - values = self.gb.dim * np.ones(self.gb.num_cells)), - Field("cell_id", cell_data = True, - values = np.arange(self.gb.num_cells)) - ]) + fields.extend( + [ + Field( + "grid_dim", + cell_data=True, + values=self.gb.dim * np.ones(self.gb.num_cells), + ), + Field("cell_id", cell_data=True, values=np.arange(self.gb.num_cells)), + ] + ) self._write_vtk(fields, name, self.gb_VTK) @@ -358,13 +366,15 @@ def _export_vtk_gb(self, data, time_step, point_data): # consider the grid_bucket node data extra_fields = Fields() - extra_fields.extend([ - Field("grid_dim", cell_data=True), - Field("cell_id", cell_data=True), - Field("grid_node_number", cell_data=True), - Field("is_mortar", cell_data=True), - Field("mortar_side", cell_data=True), - ]) + extra_fields.extend( + [ + Field("grid_dim", cell_data=True), + Field("cell_id", cell_data=True), + Field("grid_node_number", cell_data=True), + Field("is_mortar", cell_data=True), + Field("mortar_side", cell_data=True), + ] + ) fields.extend(extra_fields) self.gb.assign_node_ordering(overwrite_existing=False) @@ -397,13 +407,15 @@ def _export_vtk_gb(self, data, time_step, point_data): # consider the grid_bucket edge data extra_fields = Fields() - extra_fields.extend([ - Field("grid_dim", cell_data=True), - Field("cell_id", cell_data=True), - Field("grid_edge_number", cell_data=True), - Field("is_mortar", cell_data=True), - Field("mortar_side", cell_data=True), - ]) + extra_fields.extend( + [ + Field("grid_dim", cell_data=True), + Field("cell_id", cell_data=True), + Field("grid_edge_number", cell_data=True), + Field("is_mortar", cell_data=True), + Field("mortar_side", cell_data=True), + ] + ) self.gb.add_edge_props(extra_fields.names()) for _, d in self.gb.edges(): d["grid_dim"] = {} @@ -566,9 +578,9 @@ def _write_vtk(self, fields, name, g_VTK): for field in fields: if field.values is None: continue - dataVTK = ns.numpy_to_vtk(field.values, - deep = True, - array_type=field.dtype()) + dataVTK = ns.numpy_to_vtk( + field.values, deep=True, array_type=field.dtype() + ) dataVTK.SetName(field.name) dataVTK.SetNumberOfComponents(field.num_components) diff --git a/test/unit/test_mpfa_robin.py b/test/unit/test_mpfa_robin.py index ba88272bbc..73ca7500a6 100644 --- a/test/unit/test_mpfa_robin.py +++ b/test/unit/test_mpfa_robin.py @@ -4,6 +4,7 @@ import porepy as pp + class RobinBoundTest(unittest.TestCase): def test_flow_left_right(self): nxs = [1, 3] @@ -21,16 +22,27 @@ def test_flow_left_right(self): dir_ind = np.ravel(np.argwhere(left)) rob_ind = np.ravel(np.argwhere(right)) - names = ['dir']*len(dir_ind) + ['rob'] * len(rob_ind) + names = ["dir"] * len(dir_ind) + ["rob"] * len(rob_ind) bnd_ind = np.hstack((dir_ind, rob_ind)) bnd = pp.BoundaryCondition(g, bnd_ind, names) p_bound = 1 rob_bound = 0 - C = (rob_bound - robin_weight*p_bound) / (robin_weight - p_bound) - u_ex = - C * g.face_normals[0] + C = (rob_bound - robin_weight * p_bound) / (robin_weight - p_bound) + u_ex = -C * g.face_normals[0] p_ex = C * g.cell_centers[0] + p_bound - self.solve_robin(g, k, bnd, robin_weight, p_bound, rob_bound, dir_ind, rob_ind, p_ex, u_ex) + self.solve_robin( + g, + k, + bnd, + robin_weight, + p_bound, + rob_bound, + dir_ind, + rob_ind, + p_ex, + u_ex, + ) def test_flow_nz_rhs(self): nxs = [1, 3] @@ -48,17 +60,29 @@ def test_flow_nz_rhs(self): dir_ind = np.ravel(np.argwhere(left)) rob_ind = np.ravel(np.argwhere(right)) - names = ['dir']*len(dir_ind) + ['rob'] * len(rob_ind) + names = ["dir"] * len(dir_ind) + ["rob"] * len(rob_ind) bnd_ind = np.hstack((dir_ind, rob_ind)) bnd = pp.BoundaryCondition(g, bnd_ind, names) p_bound = 1 rob_bound = 1 - C = (rob_bound - robin_weight*p_bound) / (robin_weight - p_bound) - u_ex = - C * g.face_normals[0] + C = (rob_bound - robin_weight * p_bound) / (robin_weight - p_bound) + u_ex = -C * g.face_normals[0] p_ex = C * g.cell_centers[0] + p_bound - self.solve_robin(g, k, bnd, robin_weight, p_bound, rob_bound, dir_ind, rob_ind, p_ex, u_ex) + self.solve_robin( + g, + k, + bnd, + robin_weight, + p_bound, + rob_bound, + dir_ind, + rob_ind, + p_ex, + u_ex, + ) + def test_flow_down(self): nxs = [1, 3] nys = [1, 3] @@ -75,16 +99,27 @@ def test_flow_down(self): dir_ind = np.ravel(np.argwhere(top)) rob_ind = np.ravel(np.argwhere(bot)) - names = ['dir']*len(dir_ind) + ['rob'] * len(rob_ind) + names = ["dir"] * len(dir_ind) + ["rob"] * len(rob_ind) bnd_ind = np.hstack((dir_ind, rob_ind)) bnd = pp.BoundaryCondition(g, bnd_ind, names) p_bound = 1 rob_bound = 1 - C = (rob_bound - robin_weight*p_bound) / (robin_weight - p_bound) + C = (rob_bound - robin_weight * p_bound) / (robin_weight - p_bound) u_ex = C * g.face_normals[1] - p_ex = C *(1- g.cell_centers[1]) + p_bound - self.solve_robin(g, k, bnd, robin_weight, p_bound, rob_bound, dir_ind, rob_ind, p_ex, u_ex) + p_ex = C * (1 - g.cell_centers[1]) + p_bound + self.solve_robin( + g, + k, + bnd, + robin_weight, + p_bound, + rob_bound, + dir_ind, + rob_ind, + p_ex, + u_ex, + ) def test_no_neumann(self): g = pp.CartGrid([2, 2], physdims=[1, 1]) @@ -100,11 +135,13 @@ def test_no_neumann(self): dir_ind = np.ravel(np.argwhere(top + left)) rob_ind = np.ravel(np.argwhere(bot + right)) - names = ['dir']*len(dir_ind) + ['rob'] * len(rob_ind) + names = ["dir"] * len(dir_ind) + ["rob"] * len(rob_ind) bnd_ind = np.hstack((dir_ind, rob_ind)) bnd = pp.BoundaryCondition(g, bnd_ind, names) - flux, bound_flux, _, _ = pp.numerics.fv.mpfa._mpfa_local(g, k, bnd, robin_weight=robin_weight) + flux, bound_flux, _, _ = pp.numerics.fv.mpfa._mpfa_local( + g, k, bnd, robin_weight=robin_weight + ) div = pp.fvutils.scalar_divergence(g) @@ -113,22 +150,27 @@ def test_no_neumann(self): u_bound[dir_ind] = g.face_centers[1, dir_ind] u_bound[rob_ind] = rob_ex * g.face_areas[rob_ind] - a = div * flux b = -div * bound_flux * u_bound p = np.linalg.solve(a.A, b) - u_ex = [np.dot(g.face_normals[:, f], np.array([0, -1, 0])) - for f in range(g.num_faces)] - u_ex = np.array(u_ex).ravel('F') + u_ex = [ + np.dot(g.face_normals[:, f], np.array([0, -1, 0])) + for f in range(g.num_faces) + ] + u_ex = np.array(u_ex).ravel("F") p_ex = g.cell_centers[1] assert np.allclose(p, p_ex) assert np.allclose(flux * p + bound_flux * u_bound, u_ex) - def solve_robin(self, g, k, bnd, robin_weight, p_bound, rob_bound, dir_ind, rob_ind, p_ex, u_ex): - flux, bound_flux, _, _ = pp.numerics.fv.mpfa._mpfa_local(g, k, bnd, robin_weight=robin_weight) + def solve_robin( + self, g, k, bnd, robin_weight, p_bound, rob_bound, dir_ind, rob_ind, p_ex, u_ex + ): + flux, bound_flux, _, _ = pp.numerics.fv.mpfa._mpfa_local( + g, k, bnd, robin_weight=robin_weight + ) div = pp.fvutils.scalar_divergence(g) @@ -136,7 +178,6 @@ def solve_robin(self, g, k, bnd, robin_weight, p_bound, rob_bound, dir_ind, rob_ u_bound[dir_ind] = p_bound u_bound[rob_ind] = rob_bound * g.face_areas[rob_ind] - a = div * flux b = -div * bound_flux * u_bound diff --git a/test/unit/test_mpsa_robin.py b/test/unit/test_mpsa_robin.py index aaf3fd3f35..559503785e 100644 --- a/test/unit/test_mpsa_robin.py +++ b/test/unit/test_mpsa_robin.py @@ -4,6 +4,7 @@ import porepy as pp + class RobinBoundTest(unittest.TestCase): def test_dir_rob(self): nx = 2 @@ -22,7 +23,7 @@ def test_dir_rob(self): neu_ind = np.ravel(np.argwhere([])) rob_ind = np.ravel(np.argwhere(right)) - names = ['dir']*len(dir_ind) + ['rob'] * len(rob_ind) + names = ["dir"] * len(dir_ind) + ["rob"] * len(rob_ind) bnd_ind = np.hstack((dir_ind, rob_ind)) bnd = pp.BoundaryCondition(g, bnd_ind, names) @@ -40,17 +41,17 @@ def T_ex(faces): sgn_n = pp.numerics.fracture_deformation.sign_of_faces(g, neu_ind) sgn_r = pp.numerics.fracture_deformation.sign_of_faces(g, rob_ind) - + u_bound[:, dir_ind] = u_ex(g.face_centers[:, dir_ind]) u_bound[:, neu_ind] = T_ex(neu_ind) * sgn_n u_bound[:, rob_ind] = ( - T_ex(rob_ind) * sgn_r - + robin_weight * u_ex(g.face_centers[:, rob_ind]) * g.face_areas[rob_ind] + T_ex(rob_ind) * sgn_r + + robin_weight * u_ex(g.face_centers[:, rob_ind]) * g.face_areas[rob_ind] ) u, T = self.solve_mpsa(g, c, robin_weight, bnd, u_bound) - assert np.allclose(u, u_ex(g.cell_centers).ravel('F')) - assert np.allclose(T, T_ex(np.arange(g.num_faces)).ravel('F')) + assert np.allclose(u, u_ex(g.cell_centers).ravel("F")) + assert np.allclose(T, T_ex(np.arange(g.num_faces)).ravel("F")) def test_dir_neu_rob(self): nx = 2 @@ -69,7 +70,7 @@ def test_dir_neu_rob(self): neu_ind = np.ravel(np.argwhere(top)) rob_ind = np.ravel(np.argwhere(right + bot)) - names = ['dir']*len(dir_ind) + ['rob'] * len(rob_ind) + names = ["dir"] * len(dir_ind) + ["rob"] * len(rob_ind) bnd_ind = np.hstack((dir_ind, rob_ind)) bnd = pp.BoundaryCondition(g, bnd_ind, names) @@ -87,17 +88,17 @@ def T_ex(faces): sgn_n = pp.numerics.fracture_deformation.sign_of_faces(g, neu_ind) sgn_r = pp.numerics.fracture_deformation.sign_of_faces(g, rob_ind) - + u_bound[:, dir_ind] = u_ex(g.face_centers[:, dir_ind]) u_bound[:, neu_ind] = T_ex(neu_ind) * sgn_n u_bound[:, rob_ind] = ( - T_ex(rob_ind) * sgn_r - + robin_weight * u_ex(g.face_centers[:, rob_ind]) * g.face_areas[rob_ind] + T_ex(rob_ind) * sgn_r + + robin_weight * u_ex(g.face_centers[:, rob_ind]) * g.face_areas[rob_ind] ) u, T = self.solve_mpsa(g, c, robin_weight, bnd, u_bound) - assert np.allclose(u, u_ex(g.cell_centers).ravel('F')) - assert np.allclose(T, T_ex(np.arange(g.num_faces)).ravel('F')) + assert np.allclose(u, u_ex(g.cell_centers).ravel("F")) + assert np.allclose(T, T_ex(np.arange(g.num_faces)).ravel("F")) def test_structured_triang(self): nx = 1 @@ -116,7 +117,7 @@ def test_structured_triang(self): neu_ind = np.ravel(np.argwhere(top)) rob_ind = np.ravel(np.argwhere(right + bot)) - names = ['dir']*len(dir_ind) + ['rob'] * len(rob_ind) + names = ["dir"] * len(dir_ind) + ["rob"] * len(rob_ind) bnd_ind = np.hstack((dir_ind, rob_ind)) bnd = pp.BoundaryCondition(g, bnd_ind, names) @@ -134,17 +135,17 @@ def T_ex(faces): sgn_n = pp.numerics.fracture_deformation.sign_of_faces(g, neu_ind) sgn_r = pp.numerics.fracture_deformation.sign_of_faces(g, rob_ind) - + u_bound[:, dir_ind] = u_ex(g.face_centers[:, dir_ind]) u_bound[:, neu_ind] = T_ex(neu_ind) * sgn_n u_bound[:, rob_ind] = ( - T_ex(rob_ind) * sgn_r - + robin_weight * u_ex(g.face_centers[:, rob_ind]) * g.face_areas[rob_ind] + T_ex(rob_ind) * sgn_r + + robin_weight * u_ex(g.face_centers[:, rob_ind]) * g.face_areas[rob_ind] ) u, T = self.solve_mpsa(g, c, robin_weight, bnd, u_bound) - assert np.allclose(u, u_ex(g.cell_centers).ravel('F')) - assert np.allclose(T, T_ex(np.arange(g.num_faces)).ravel('F')) + assert np.allclose(u, u_ex(g.cell_centers).ravel("F")) + assert np.allclose(T, T_ex(np.arange(g.num_faces)).ravel("F")) def test_unstruct_triang(self): corners = np.array([[0, 0, 1, 1], [0, 1, 1, 0]]) @@ -164,7 +165,7 @@ def test_unstruct_triang(self): neu_ind = np.ravel(np.argwhere(top)) rob_ind = np.ravel(np.argwhere(right + bot)) - names = ['dir']*len(dir_ind) + ['rob'] * len(rob_ind) + names = ["dir"] * len(dir_ind) + ["rob"] * len(rob_ind) bnd_ind = np.hstack((dir_ind, rob_ind)) bnd = pp.BoundaryCondition(g, bnd_ind, names) @@ -182,26 +183,24 @@ def T_ex(faces): sgn_n = pp.numerics.fracture_deformation.sign_of_faces(g, neu_ind) sgn_r = pp.numerics.fracture_deformation.sign_of_faces(g, rob_ind) - + u_bound[:, dir_ind] = u_ex(g.face_centers[:, dir_ind]) u_bound[:, neu_ind] = T_ex(neu_ind) * sgn_n u_bound[:, rob_ind] = ( - T_ex(rob_ind) * sgn_r - + robin_weight * u_ex(g.face_centers[:, rob_ind]) * g.face_areas[rob_ind] + T_ex(rob_ind) * sgn_r + + robin_weight * u_ex(g.face_centers[:, rob_ind]) * g.face_areas[rob_ind] ) u, T = self.solve_mpsa(g, c, robin_weight, bnd, u_bound) - assert np.allclose(u, u_ex(g.cell_centers).ravel('F')) - assert np.allclose(T, T_ex(np.arange(g.num_faces)).ravel('F')) + assert np.allclose(u, u_ex(g.cell_centers).ravel("F")) + assert np.allclose(T, T_ex(np.arange(g.num_faces)).ravel("F")) def test_unstruct_tetrahedron(self): - box = {'xmin': 0, 'xmax': 1, - 'ymin': 0, 'ymax': 1, - 'zmin': 0, 'zmax': 1} + box = {"xmin": 0, "xmax": 1, "ymin": 0, "ymax": 1, "zmin": 0, "zmax": 1} g = pp.meshing.simplex_grid([], box, mesh_size_min=3, mesh_size_frac=3) g = g.grids_of_dimension(3)[0] - g = pp.CartGrid([4,4,4], [1,1,1]) + g = pp.CartGrid([4, 4, 4], [1, 1, 1]) g.compute_geometry() c = pp.FourthOrderTensor(3, np.ones(g.num_cells), np.ones(g.num_cells)) robin_weight = 1.0 @@ -217,19 +216,17 @@ def test_unstruct_tetrahedron(self): neu_ind = np.ravel(np.argwhere(bot)) rob_ind = np.ravel(np.argwhere(east + north + south)) - names = ['dir']*len(dir_ind) + ['rob'] * len(rob_ind) + names = ["dir"] * len(dir_ind) + ["rob"] * len(rob_ind) bnd_ind = np.hstack((dir_ind, rob_ind)) bnd = pp.BoundaryCondition(g, bnd_ind, names) def u_ex(x): - return np.vstack((x[1], x[0], 0*x[2])) + return np.vstack((x[1], x[0], 0 * x[2])) def T_ex(faces): if np.size(faces) == 0: return np.atleast_2d(np.array([])) - sigma = np.array([[0, 2, 0], - [2, 0, 0], - [0, 0, 0]]) + sigma = np.array([[0, 2, 0], [2, 0, 0], [0, 0, 0]]) T_r = [np.dot(sigma, g.face_normals[:, f]) for f in faces] return np.vstack(T_r).T @@ -237,25 +234,26 @@ def T_ex(faces): sgn_n = pp.numerics.fracture_deformation.sign_of_faces(g, neu_ind) sgn_r = pp.numerics.fracture_deformation.sign_of_faces(g, rob_ind) - + u_bound[:, dir_ind] = u_ex(g.face_centers[:, dir_ind]) u_bound[:, neu_ind] = T_ex(neu_ind) * sgn_n u_bound[:, rob_ind] = ( - T_ex(rob_ind) * sgn_r - + robin_weight * u_ex(g.face_centers[:, rob_ind]) * g.face_areas[rob_ind] + T_ex(rob_ind) * sgn_r + + robin_weight * u_ex(g.face_centers[:, rob_ind]) * g.face_areas[rob_ind] ) u, T = self.solve_mpsa(g, c, robin_weight, bnd, u_bound) - assert np.allclose(u, u_ex(g.cell_centers).ravel('F')) - assert np.allclose(T, T_ex(np.arange(g.num_faces)).ravel('F')) - + assert np.allclose(u, u_ex(g.cell_centers).ravel("F")) + assert np.allclose(T, T_ex(np.arange(g.num_faces)).ravel("F")) def solve_mpsa(self, g, c, robin_weight, bnd, u_bound): - stress, bound_stress = pp.numerics.fv.mpsa._mpsa_local(g, c, bnd, robin_weight=robin_weight) + stress, bound_stress = pp.numerics.fv.mpsa._mpsa_local( + g, c, bnd, robin_weight=robin_weight + ) div = pp.fvutils.vector_divergence(g) a = div * stress - b = -div * bound_stress * u_bound.ravel('F') + b = -div * bound_stress * u_bound.ravel("F") u = np.linalg.solve(a.A, b) - T = stress * u + bound_stress * u_bound.ravel('F') + T = stress * u + bound_stress * u_bound.ravel("F") return u, T diff --git a/test/unit/test_vtk.py b/test/unit/test_vtk.py index f152b2ce03..66fc0f1570 100644 --- a/test/unit/test_vtk.py +++ b/test/unit/test_vtk.py @@ -206,11 +206,11 @@ def test_gb_1(self): content = content_file.read() assert content == self._gb_1_grid_2_vtu() - with open(folder+"grid_mortar_1.vtu", "r") as content_file: + with open(folder + "grid_mortar_1.vtu", "r") as content_file: content = content_file.read() assert content == self._gb_1_mortar_grid_vtu() -#------------------------------------------------------------------------------# + # ------------------------------------------------------------------------------# def test_gb_2(self): if not if_vtk: @@ -245,11 +245,11 @@ def test_gb_2(self): print(content) assert content == self._gb_2_grid_2_vtu() - with open(folder+"grid_mortar_1.vtu", "r") as content_file: + with open(folder + "grid_mortar_1.vtu", "r") as content_file: content = content_file.read() assert content == self._gb_2_mortar_grid_1_vtu() -#------------------------------------------------------------------------------# + # ------------------------------------------------------------------------------# def _single_grid_1d_grid_vtu(self): return """ @@ -2388,7 +2388,7 @@ def _gb_1_mortar_grid_vtu(self): """ -#------------------------------------------------------------------------------# + # ------------------------------------------------------------------------------# def _gb_1_mortar_grid_vtu(self): return """ @@ -2448,7 +2448,7 @@ def _gb_1_mortar_grid_vtu(self): """ -#------------------------------------------------------------------------------# + # ------------------------------------------------------------------------------# def _gb_2_grid_pvd(self): return """ @@ -2622,7 +2622,7 @@ def _gb_2_grid_2_vtu(self): """ -#------------------------------------------------------------------------------# + # ------------------------------------------------------------------------------# def _gb_2_mortar_grid_1_vtu(self): return """ @@ -2688,8 +2688,9 @@ def _gb_2_mortar_grid_1_vtu(self): """ -#------------------------------------------------------------------------------# -if __name__ == '__main__': +# ------------------------------------------------------------------------------# + +if __name__ == "__main__": BasicsTest().test_gb_2() unittest.main() From 98dc22598b7b3e5402ba2ea77755865fb13c0f67 Mon Sep 17 00:00:00 2001 From: Runar Date: Wed, 12 Sep 2018 11:33:23 +0200 Subject: [PATCH 20/26] changed assert -> Raise AssertionError() --- src/porepy/numerics/fv/mpfa.py | 7 +++++-- src/porepy/numerics/fv/mpsa.py | 16 ++++++++++------ 2 files changed, 15 insertions(+), 8 deletions(-) diff --git a/src/porepy/numerics/fv/mpfa.py b/src/porepy/numerics/fv/mpfa.py index be526d51b1..95d7bb1133 100644 --- a/src/porepy/numerics/fv/mpfa.py +++ b/src/porepy/numerics/fv/mpfa.py @@ -889,7 +889,8 @@ def _tensor_vector_prod(g, k, subcell_topology, apertures=None): # Duplicates in [cno, nno] corresponds to different faces meeting at the # same node. There should be exactly nd of these. This test will fail # for pyramids in 3D - assert np.all(blocksz == nd) + if not np.all(blocksz == nd): + raise AssertionError() # Define row and column indices to be used for normal_vectors * perm. # Rows are based on sub-face numbers. @@ -998,7 +999,9 @@ def _create_bound_rhs( fno = subcell_topology.fno_unique num_neu = np.sum(is_neu[fno]) num_dir = np.sum(is_dir[fno]) - assert num_rob == np.sum(is_rob[fno]) + if not num_rob == np.sum(is_rob[fno]): + raise AssertionError() + num_bound = num_neu + num_dir + num_rob # Neumann and Robin boundary conditions. Neumann and Robin conditions diff --git a/src/porepy/numerics/fv/mpsa.py b/src/porepy/numerics/fv/mpsa.py index 399da77f82..25d4b0a04f 100644 --- a/src/porepy/numerics/fv/mpsa.py +++ b/src/porepy/numerics/fv/mpsa.py @@ -142,7 +142,8 @@ class FracturedMpsa(Mpsa): def __init__(self, given_traction=False, **kwargs): Mpsa.__init__(self, **kwargs) - assert hasattr(self, "physics"), "Mpsa must assign physics" + if not hasattr(self, "physics"): + raise AttributeError("Mpsa must assign physics") self.given_traction_flag = given_traction def ndof(self, g): @@ -1255,7 +1256,8 @@ def _tensor_vector_prod(g, constit, subcell_topology): # Duplicates in [cno, nno] corresponds to different faces meeting at the # same node. There should be exactly nd of these. This test will fail # for pyramids in 3D - assert np.all(blocksz == nd) + if not np.all(blocksz == nd): + raise AssertionError() # Define row and column indices to be used for normal vector matrix # Rows are based on sub-face numbers. @@ -1478,7 +1480,9 @@ def create_bound_rhs(bound, bound_exclusion, subcell_topology, g): num_neu = sum(bound.is_neu[fno]) * nd num_dir = sum(bound.is_dir[fno]) * nd - assert num_rob == sum(bound.is_rob[fno]) * nd + if not num_rob == sum(bound.is_rob[fno]) * nd: + raise AssertionError() + num_bound = num_neu + num_dir + num_rob # Convenience method for duplicating a list, with a certain increment @@ -1936,9 +1940,9 @@ def __rearange_columns_displacement_eqs(d_cont_grad, d_cont_cell, num_sub_cells, def _neu_face_sgn(g, neu_ind): neu_sgn = (g.cell_faces[neu_ind, :]).data - assert ( - neu_sgn.size == neu_ind.size - ), "A normal sign is only well defined for a boundary face" + if not neu_sgn.size == neu_ind.size: + raise AssertionError("A normal sign is only well defined for a boundary face") + sort_id = np.argsort(g.cell_faces[neu_ind, :].indices) return neu_sgn[sort_id] From e69b85eb901f75a4d656fd10f8f1a37f004ae790 Mon Sep 17 00:00:00 2001 From: Runar Date: Wed, 12 Sep 2018 17:14:22 +0200 Subject: [PATCH 21/26] Added robin condition to vectorial bc NOTE: not tested --- src/porepy/numerics/fv/fvutils.py | 476 +++++++++++++++++++++++------- src/porepy/numerics/fv/mpsa.py | 6 +- 2 files changed, 375 insertions(+), 107 deletions(-) diff --git a/src/porepy/numerics/fv/fvutils.py b/src/porepy/numerics/fv/fvutils.py index 254b3f50fa..2f7f24224e 100644 --- a/src/porepy/numerics/fv/fvutils.py +++ b/src/porepy/numerics/fv/fvutils.py @@ -687,74 +687,27 @@ def __init__(self, subcell_topology, bound, nd): self.keep_rob = self.__exclude_matrix(~bound.is_rob) elif self.bc_type == "vectorial": - # Neumann - col_neu_x = np.argwhere([not it for it in bound.is_neu[0, fno]]) - row_neu_x = np.arange(col_neu_x.size) - self.exclude_neu_x = sps.coo_matrix( - (np.ones(row_neu_x.size), (row_neu_x, col_neu_x.ravel("C"))), - shape=(row_neu_x.size, num_subfno), - ).tocsr() - - col_neu_y = np.argwhere([not it for it in bound.is_neu[1, fno]]) - row_neu_y = np.arange(col_neu_y.size) - self.exclude_neu_y = sps.coo_matrix( - (np.ones(row_neu_y.size), (row_neu_y, col_neu_y.ravel("C"))), - shape=(row_neu_y.size, num_subfno), - ).tocsr() - col_neu_y += num_subfno - col_neu = np.append(col_neu_x, [col_neu_y]) - - if nd == 3: - col_neu_z = np.argwhere([not it for it in bound.is_neu[2, fno]]) - row_neu_z = np.arange(col_neu_z.size) - self.exclude_neu_z = sps.coo_matrix( - (np.ones(row_neu_z.size), (row_neu_z, col_neu_z.ravel("C"))), - shape=(row_neu_z.size, num_subfno), - ).tocsr() - - col_neu_z += 2 * num_subfno - col_neu = np.append(col_neu, [col_neu_z]) - - row_neu = np.arange(col_neu.size) - self.exclude_neu_nd = sps.coo_matrix( - (np.ones(row_neu.size), (row_neu, col_neu.ravel("C"))), - shape=(row_neu.size, nd * num_subfno), - ).tocsr() - - # Dirichlet, same procedure - col_dir_x = np.argwhere([not it for it in bound.is_dir[0, fno]]) - row_dir_x = np.arange(col_dir_x.size) - self.exclude_dir_x = sps.coo_matrix( - (np.ones(row_dir_x.size), (row_dir_x, col_dir_x.ravel("C"))), - shape=(row_dir_x.size, num_subfno), - ).tocsr() - - col_dir_y = np.argwhere([not it for it in bound.is_dir[1, fno]]) - row_dir_y = np.arange(col_dir_y.size) - self.exclude_dir_y = sps.coo_matrix( - (np.ones(row_dir_y.size), (row_dir_y, col_dir_y.ravel("C"))), - shape=(row_dir_y.size, num_subfno), - ).tocsr() - - col_dir_y += num_subfno - col_dir = np.append(col_dir_x, [col_dir_y]) - - if nd == 3: - col_dir_z = np.argwhere([not it for it in bound.is_dir[2, fno]]) - row_dir_z = np.arange(col_dir_z.size) - self.exclude_dir_z = sps.coo_matrix( - (np.ones(row_dir_z.size), (row_dir_z, col_dir_z.ravel("C"))), - shape=(row_dir_z.size, num_subfno), - ).tocsr() - - col_dir_z += 2 * num_subfno - col_dir = np.append(col_dir, [col_dir_z]) - - row_dir = np.arange(col_dir.size) - self.exclude_dir_nd = sps.coo_matrix( - (np.ones(row_dir.size), (row_dir, col_dir.ravel("C"))), - shape=(row_dir.size, nd * num_subfno), - ).tocsr() + self.exclude_neu_nd, self.exclude_neu_x, self.exclude_neu_y, self.exclude_neu_z = self.__exclude_matrix_xyz( + bound.is_neu + ) + self.exclude_dir_nd, self.exclude_dir_x, self.exclude_dir_y, self.exclude_dir_z = self.__exclude_matrix_xyz( + bound.is_dir + ) + self.exclude_rob_nd, self.exclude_rob_x, self.exclude_rob_y, self.exclude_rob_z = self.__exclude_matrix_xyz( + bound.is_rob + ) + self.exclude_neu_dir_nd, self.exclude_neu_dir_x, self.exclude_neu_dir_y, self.exclude_neu_dir_z = self.__exclude_matrix_xyz( + bound.is_neu | bound.is_dir + ) + self.exclude_neu_rob_nd, self.exclude_neu_rob_x, self.exclude_neu_rob_y, self.exclude_neu_rob_z = self.__exclude_matrix_xyz( + bound.is_neu | bound.is_rob + ) + self.exclude_rob_dir_nd, self.exclude_rob_dir_x, self.exclude_rob_dir_y, self.exclude_rob_dir_z = self.__exclude_matrix_xyz( + bound.is_rob | bound.is_dir + ) + self.keep_rob_nd, self.keep_rob_x, self.keep_rob_y, self.keep_rob_z = self.__exclude_matrix_xyz( + ~bound.is_rob + ) def __exclude_matrix(self, ids): """ @@ -774,6 +727,43 @@ def __exclude_matrix(self, ids): shape=(row.size, self.num_subfno), ).tocsr() + def __exclude_matrix_xyz(self, ids): + col_x = np.argwhere([not it for it in ids[0, self.fno]]) + row_x = np.arange(col_x.size) + exclude_x = sps.coo_matrix( + (np.ones(row_x.size), (row_x, col_x.ravel("C"))), + shape=(row_x.size, self.num_subfno), + ).tocsr() + + col_y = np.argwhere([not it for it in ids[1, self.fno]]) + row_y = np.arange(col_y.size) + exclude_y = sps.coo_matrix( + (np.ones(row_y.size), (row_y, col_y.ravel("C"))), + shape=(row_y.size, self.num_subfno), + ).tocsr() + col_y += self.num_subfno + col_neu = np.append(col_x, [col_y]) + + if self.nd == 3: + col_z = np.argwhere([not it for it in ids[2, self.fno]]) + row_z = np.arange(col_z.size) + exclude_z = sps.coo_matrix( + (np.ones(row_z.size), (row_z, col_z.ravel("C"))), + shape=(row_z.size, self.num_subfno), + ).tocsr() + + col_z += 2 * self.num_subfno + col_neu = np.append(col_neu, [col_z]) + + else: + exclude_z = sps.coo_matrix(shape=(0, self.num_subfno)) + row_neu = np.arange(col_neu.size) + exclude_nd = sps.coo_matrix( + (np.ones(row_neu.size), (row_neu, col_neu.ravel("C"))), + shape=(row_neu.size, self.nd * self.num_subfno), + ).tocsr() + return exclude_nd, exclude_x, exclude_y, exclude_z + def exclude_dirichlet(self, other): """ Mapping to exclude faces/components with Dirichlet boundary conditions from local systems. @@ -876,6 +866,57 @@ def exclude_neumann(self, other): return exclude_neu + def exclude_neumann_x(self, other): + """ Mapping to exclude components in the x-direction with Neumann boundary conditions from + local systems. + NOTE: Currently works for boundary faces aligned with the coordinate system. + + Parameters: + other (scipy.sparse matrix): Matrix of local equations for + continuity of flux and pressure. + + Returns: + scipy.sparse matrix, with components in the x-direction corresponding to faces with + Neumann conditions eliminated. + + """ + + return self.exclude_neu_x * other + + def exclude_neumann_y(self, other): + """ Mapping to exclude components in the y-direction with Neumann boundary conditions from + local systems. + NOTE: Currently works for boundary faces aligned with the coordinate system. + + Parameters: + other (scipy.sparse matrix): Matrix of local equations for + continuity of flux and pressure. + + Returns: + scipy.sparse matrix, with components in the y-direction corresponding to faces with + Neumann conditions eliminated. + + """ + + return self.exclude_neu_y * other + + def exclude_neumann_z(self, other): + """ Mapping to exclude components in the z-direction with Neumann boundary conditions from + local systems. + NOTE: Currently works for boundary faces aligned with the coordinate system. + + Parameters: + other (scipy.sparse matrix): Matrix of local equations for + continuity of flux and pressure. + + Returns: + scipy.sparse matrix, with components in the z-direction corresponding to faces with + Neumann conditions eliminated. + + """ + + return self.exclude_neu_z * other + def exclude_robin(self, other): """ Mapping to exclude faces/components with robin boundary conditions from local systems. @@ -893,10 +934,67 @@ def exclude_robin(self, other): exclude_rob = self.exclude_rob * other elif self.bc_type == "vectorial": + exclude_rob = np.append( + self.exclude_rob_x * other, [self.exclude_rob_y * other] + ) + if self.nd == 3: + exclude_rob = np.append(exclude_rob, [self.exclude_rob_z * other]) + raise NotImplementedError("can not exclude robin for vectorial bc") return exclude_rob + def exclude_robin_x(self, other): + """ Mapping to exclude components in the x-direction with Robin boundary conditions from + local systems. + NOTE: Currently works for boundary faces aligned with the coordinate system. + + Parameters: + other (scipy.sparse matrix): Matrix of local equations for + continuity of flux and pressure. + + Returns: + scipy.sparse matrix, with components in the x-direction corresponding to faces with + Robin conditions eliminated. + + """ + + return self.exclude_rob_x * other + + def exclude_robin_y(self, other): + """ Mapping to exclude components in the y-direction with Robin boundary conditions from + local systems. + NOTE: Currently works for boundary faces aligned with the coordinate system. + + Parameters: + other (scipy.sparse matrix): Matrix of local equations for + continuity of flux and pressure. + + Returns: + scipy.sparse matrix, with components in the y-direction corresponding to faces with + Robin conditions eliminated. + + """ + + return self.exclude_rob_y * other + + def exclude_robin_z(self, other): + """ Mapping to exclude components in the z-direction with Robin boundary conditions from + local systems. + NOTE: Currently works for boundary faces aligned with the coordinate system. + + Parameters: + other (scipy.sparse matrix): Matrix of local equations for + continuity of flux and pressure. + + Returns: + scipy.sparse matrix, with components in the z-direction corresponding to faces with + Robin conditions eliminated. + + """ + + return self.exclude_rob_z * other + def exclude_neumann_robin(self, other): """ Mapping to exclude faces/components with Neumann and Robin boundary conditions from local systems. @@ -911,11 +1009,69 @@ def exclude_neumann_robin(self, other): """ if self.bc_type == "scalar": - exclude_neu = self.exclude_neu_rob * other + exclude_neu_rob = self.exclude_neu_rob * other elif self.bc_type == "vectorial": - raise NotImplementedError() - return exclude_neu + exclude_neu_rob = np.append( + self.exclude_neu_rob_x * other, [self.exclude_neu_rob_y * other] + ) + if self.nd == 3: + exclude_neu_rob = np.append( + exclude_neu_rob, [self.exclude_neu_rob_z * other] + ) + + return exclude_neu_rob + + def exclude_neumann_robin_x(self, other): + """ Mapping to exclude components in the x-direction with Neumann or Robin boundary conditions from + local systems. + NOTE: Currently works for boundary faces aligned with the coordinate system. + + Parameters: + other (scipy.sparse matrix): Matrix of local equations for + continuity of flux and pressure. + + Returns: + scipy.sparse matrix, with components in the x-direction corresponding to faces with + Neumann and Robin conditions eliminated. + + """ + + return self.exclude_neu_rob_x * other + + def exclude_neumann_robin_y(self, other): + """ Mapping to exclude components in the y-direction with Neumann or Robin boundary conditions from + local systems. + NOTE: Currently works for boundary faces aligned with the coordinate system. + + Parameters: + other (scipy.sparse matrix): Matrix of local equations for + continuity of flux and pressure. + + Returns: + scipy.sparse matrix, with components in the y-direction corresponding to faces with + Neumann and Robin conditions eliminated. + + """ + + return self.exclude_neu_rob_y * other + + def exclude_neumann_robin_z(self, other): + """ Mapping to exclude components in the z-direction with Neumann or Robin boundary conditions from + local systems. + NOTE: Currently works for boundary faces aligned with the coordinate system. + + Parameters: + other (scipy.sparse matrix): Matrix of local equations for + continuity of flux and pressure. + + Returns: + scipy.sparse matrix, with components in the z-direction corresponding to faces with + Neuman and Robin conditions eliminated. + + """ + + return self.exclude_neu_rob_z * other def exclude_neumann_dirichlet(self, other): """ Mapping to exclude faces/components with Neumann and Dirichlet boundary @@ -931,11 +1087,69 @@ def exclude_neumann_dirichlet(self, other): """ if self.bc_type == "scalar": - exclude_neu = self.exclude_neu_dir * other + exclude_neu_dir = self.exclude_neu_dir * other elif self.bc_type == "vectorial": - raise NotImplementedError() - return exclude_neu + exclude_neu_dir = np.append( + self.exclude_neu_dir_x * other, [self.exclude_neu_dir_y * other] + ) + if self.nd == 3: + exclude_neu_dir = np.append( + exclude_neu_dir, [self.exclude_neu_dir_z * other] + ) + + return exclude_neu_dir + + def exclude_neumann_dirichlet_x(self, other): + """ Mapping to exclude components in the x-direction with Neumann or Dirichlet boundary conditions from + local systems. + NOTE: Currently works for boundary faces aligned with the coordinate system. + + Parameters: + other (scipy.sparse matrix): Matrix of local equations for + continuity of flux and pressure. + + Returns: + scipy.sparse matrix, with components in the x-direction corresponding to faces with + Neumann and Dirichlet conditions eliminated. + + """ + + return self.exclude_neu_dir_x * other + + def exclude_neumann_dirichlet_y(self, other): + """ Mapping to exclude components in the y-direction with Neumann or Dirichlet boundary conditions from + local systems. + NOTE: Currently works for boundary faces aligned with the coordinate system. + + Parameters: + other (scipy.sparse matrix): Matrix of local equations for + continuity of flux and pressure. + + Returns: + scipy.sparse matrix, with components in the y-direction corresponding to faces with + Neumann and Dirichlet conditions eliminated. + + """ + + return self.exclude_neu_dir_y * other + + def exclude_neumann_dirichlet_z(self, other): + """ Mapping to exclude components in the z-direction with Neumann or Dirichlet boundary conditions from + local systems. + NOTE: Currently works for boundary faces aligned with the coordinate system. + + Parameters: + other (scipy.sparse matrix): Matrix of local equations for + continuity of flux and pressure. + + Returns: + scipy.sparse matrix, with components in the z-direction corresponding to faces with + Neumann and Dirichlet conditions eliminated. + + """ + + return self.exclude_neu_dir_z * other def exclude_robin_dirichlet(self, other): """ Mapping to exclude faces/components with Robin and Dirichlet boundary @@ -951,11 +1165,69 @@ def exclude_robin_dirichlet(self, other): """ if self.bc_type == "scalar": - exclude_rob = self.exclude_rob_dir * other + exclude_rob_dir = self.exclude_rob_dir * other elif self.bc_type == "vectorial": - raise NotImplementedError() - return exclude_rob + exclude_rob_dir = np.append( + self.exclude_rob_dir_x * other, [self.exclude_rob_dir_y * other] + ) + if self.nd == 3: + exclude_rob_dir = np.append( + exclude_rob_dir, [self.exclude_rob_dir_z * other] + ) + + return exclude_rob_dir + + def exclude_robin_dirichlet_x(self, other): + """ Mapping to exclude components in the x-direction with Robin or Dirichlet boundary conditions from + local systems. + NOTE: Currently works for boundary faces aligned with the coordinate system. + + Parameters: + other (scipy.sparse matrix): Matrix of local equations for + continuity of flux and pressure. + + Returns: + scipy.sparse matrix, with components in the x-direction corresponding to faces with + Robin and Dirichlet conditions eliminated. + + """ + + return self.exclude_rob_dir_x * other + + def exclude_robin_dirichlet_y(self, other): + """ Mapping to exclude components in the y-direction with Robin or Dirichlet boundary conditions from + local systems. + NOTE: Currently works for boundary faces aligned with the coordinate system. + + Parameters: + other (scipy.sparse matrix): Matrix of local equations for + continuity of flux and pressure. + + Returns: + scipy.sparse matrix, with components in the y-direction corresponding to faces with + Robin and Dirichlet conditions eliminated. + + """ + + return self.exclude_rob_dir_y * other + + def exclude_robin_dirichlet_z(self, other): + """ Mapping to exclude components in the z-direction with Robin or Dirichlet boundary conditions from + local systems. + NOTE: Currently works for boundary faces aligned with the coordinate system. + + Parameters: + other (scipy.sparse matrix): Matrix of local equations for + continuity of flux and pressure. + + Returns: + scipy.sparse matrix, with components in the z-direction corresponding to faces with + Robin and Dirichlet conditions eliminated. + + """ + + return self.exclude_rob_dir_z * other def keep_robin(self, other): """ Mapping to exclude faces/components that is not on the Robin boundary @@ -971,14 +1243,17 @@ def keep_robin(self, other): """ if self.bc_type == "scalar": - exclude_rob = self.keep_rob * other + keep_rob = self.keep_rob * other elif self.bc_type == "vectorial": - raise NotImplementedError() - return exclude_rob + keep_rob = np.append(self.keep_rob_x * other, [self.keep_rob_y * other]) + if self.nd == 3: + keep_rob = np.append(keep_rob, [self.keep_rob_z * other]) - def exclude_neumann_x(self, other): - """ Mapping to exclude components in the x-direction with Neumann boundary conditions from + return keep_rob + + def keep_robin_x(self, other): + """ Mapping to exclude components in the x-direction that is not boundary conditions from local systems. NOTE: Currently works for boundary faces aligned with the coordinate system. @@ -988,14 +1263,13 @@ def exclude_neumann_x(self, other): Returns: scipy.sparse matrix, with components in the x-direction corresponding to faces with - Neumann conditions eliminated. - + all but Robin conditions eliminated. """ - return self.exclude_neu_x * other + return self.keep_rob_x * other - def exclude_neumann_y(self, other): - """ Mapping to exclude components in the y-direction with Neumann boundary conditions from + def keep_robin_y(self, other): + """ Mapping to exclude components in the y-direction that is not boundary conditions from local systems. NOTE: Currently works for boundary faces aligned with the coordinate system. @@ -1005,14 +1279,13 @@ def exclude_neumann_y(self, other): Returns: scipy.sparse matrix, with components in the y-direction corresponding to faces with - Neumann conditions eliminated. - + all but Robin conditions eliminated. """ - return self.exclude_neu_y * other + return self.keep_rob_y * other - def exclude_neumann_z(self, other): - """ Mapping to exclude components in the z-direction with Neumann boundary conditions from + def keep_robin_z(self, other): + """ Mapping to exclude components in the z-direction that is not boundary conditions from local systems. NOTE: Currently works for boundary faces aligned with the coordinate system. @@ -1022,18 +1295,15 @@ def exclude_neumann_z(self, other): Returns: scipy.sparse matrix, with components in the z-direction corresponding to faces with - Neumann conditions eliminated. - + all but Robin conditions eliminated. """ - return self.exclude_neu_z * other + return self.keep_rob_z * other def exclude_dirichlet_nd(self, other): """ Exclusion of Dirichlet conditions for vector equations (elasticity). See above method without _nd suffix for description. - """ - if self.bc_type == "scalar": exclude_dirichlet_nd = sps.kron(sps.eye(self.nd), self.exclude_dir) @@ -1064,7 +1334,7 @@ def exclude_neumann_robin_nd(self, other): exclude_neu_rob_nd = sps.kron(sps.eye(self.nd), self.exclude_neu_rob) elif self.bc_type == "vectorial": - raise NotImplementedError() + exclude_neu_rob_nd = self.exclude_neu_rob_nd return exclude_neu_rob_nd * other @@ -1077,7 +1347,7 @@ def exclude_neumann_dirichlet_nd(self, other): exclude_neu_dir_nd = sps.kron(sps.eye(self.nd), self.exclude_neu_dir) elif self.bc_type == "vectorial": - raise NotImplementedError() + exclude_neu_dir_nd = self.exclude_neu_dir_nd return exclude_neu_dir_nd * other @@ -1090,7 +1360,7 @@ def exclude_robin_dirichlet_nd(self, other): exclude_rob_dir_nd = sps.kron(sps.eye(self.nd), self.exclude_rob_dir) elif self.bc_type == "vectorial": - raise NotImplementedError() + exclude_rob_dir_nd = self.exclude_rob_dir_nd return exclude_rob_dir_nd * other @@ -1099,12 +1369,12 @@ def keep_robin_nd(self, other): See above method without _nd suffix for description. """ if self.bc_type == "scalar": - exclude_neu_dir_nd = sps.kron(sps.eye(self.nd), self.keep_rob) + keep_rob_nd = sps.kron(sps.eye(self.nd), self.keep_rob) elif self.bc_type == "vectorial": - raise NotImplementedError() + keep_rob_nd = self.keep_rob_nd - return exclude_neu_dir_nd * other + return keep_rob_nd * other # -----------------End of class ExcludeBoundaries----------------------------- diff --git a/src/porepy/numerics/fv/mpsa.py b/src/porepy/numerics/fv/mpsa.py index 25d4b0a04f..eb0222ac9e 100644 --- a/src/porepy/numerics/fv/mpsa.py +++ b/src/porepy/numerics/fv/mpsa.py @@ -1368,7 +1368,6 @@ def _inverse_gradient( nd, inverter, ): - # Mappings to convert linear system to block diagonal form rows2blk_diag, cols2blk_diag, size_of_blocks = _block_diagonal_structure( sub_cell_index, cell_node_blocks, nno_unique, bound_exclusion, nd @@ -1425,9 +1424,8 @@ def _block_diagonal_structure( ) elif bound_exclusion.bc_type == "vectorial": - node_occ = np.hstack((nno_stress, nno_displacement)) - if nno_rob.size >= 0: - raise NotImplementedError() + node_occ = np.hstack((nno_stress, nno_rob, nno_displacement)) + sorted_ind = np.argsort(node_occ, kind="mergesort") rows2blk_diag = sps.coo_matrix( (np.ones(sorted_ind.size), (np.arange(sorted_ind.size), sorted_ind)) From 38d20569ffbe238d42fab29e06bfd80d790c1501 Mon Sep 17 00:00:00 2001 From: Runar Date: Fri, 14 Sep 2018 13:15:48 +0200 Subject: [PATCH 22/26] Added the robin_weight to the parameter class --- src/porepy/numerics/fv/mpfa.py | 32 ++++----- src/porepy/params/data.py | 120 +++++++++++++++++++++++++-------- test/unit/test_parameters.py | 25 +++++++ 3 files changed, 131 insertions(+), 46 deletions(-) diff --git a/src/porepy/numerics/fv/mpfa.py b/src/porepy/numerics/fv/mpfa.py index 95d7bb1133..46a9fdd828 100644 --- a/src/porepy/numerics/fv/mpfa.py +++ b/src/porepy/numerics/fv/mpfa.py @@ -89,14 +89,12 @@ def matrix_rhs(self, g, data, discretize=True): order elliptic equation using a FV method with a multi-point flux approximation. - The name of data in the input dictionary (data) are: - k : second_order_tensor - Permeability defined cell-wise. - bc : boundary conditions (optional) - bc_val : dictionary (optional) - Values of the boundary conditions. The dictionary has at most the - following keys: 'dir' and 'neu', for Dirichlet and Neumann boundary - conditions, respectively. + The data should contain a parameter class under the field "param". + The following parameters will be accessed: + get_tensor : SecondOrderTensor. Permeability defined cell-wise. + get_bc : boundary conditions + get_bc_val : boundary values + get_robin_weight : float. Weight for pressure in Robin condition Parameters ---------- @@ -147,16 +145,11 @@ def rhs(self, g, bound_flux, bc_val): def discretize(self, g, data): """ - The name of data in the input dictionary (data) are: - k : second_order_tensor - Permeability defined cell-wise. If not given a identity permeability - is assumed and a warning arised. - bc : boundary conditions (optional) - bc_val : dictionary (optional) - Values of the boundary conditions. The dictionary has at most the - following keys: 'dir' and 'neu', for Dirichlet and Neumann boundary - conditions, respectively. - + The data should contain a parameter class under the field "param". + The following parameters will be accessed: + get_tensor : SecondOrderTensor. Permeability defined cell-wise. + get_bc : boundary conditions + get_robin_weight : float. Weight for pressure in Robin condition Parameters ---------- g : grid, or a subclass, with geometry fields computed. @@ -167,8 +160,9 @@ def discretize(self, g, data): k = param.get_tensor(self) bnd = param.get_bc(self) a = param.aperture + robin_weight = param.get_robin_weight(self) - trm, bound_flux, bp_cell, bp_face = mpfa(g, k, bnd, apertures=a) + 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 diff --git a/src/porepy/params/data.py b/src/porepy/params/data.py index 5fef51112c..a9eeedf75b 100644 --- a/src/porepy/params/data.py +++ b/src/porepy/params/data.py @@ -569,6 +569,99 @@ def get_stiffness(self): stiffness = property(get_stiffness) + def get_robin_weight(self, obj): + """ Pick out physics-specific weight for the Robin bc-condition. + + Discretization methods considering Robin conditions should access this method. + + Parameters: + + obj : Solver + Discretization object. Should have attribute 'physics'. + + Returns: + + tensor, representing + Permeability if obj.physics equals 'flow' + conductivity if obj.physics equals 'transport' + stiffness if physics equals 'mechanics' + + """ + physics = self._get_physics(obj) + + if physics == "flow": + return self.get_robin_weight_flow() + elif physics == "transport": + return self.get_robin_weight_transport() + elif physics == "mechanics": + return self.get_robin_weight_mechanics() + else: + raise ValueError( + 'Unknown physics "%s".\n Possible physics are: %s' + % (physics, self.known_physics) + ) + + def get_robin_weight_flow(self): + """ BoundaryCondition object + float giving the Robin factor + Solvers should rather access get_tensor(). + """ + if hasattr(self, "_robin_weight_flow"): + return self._robin_weight_flow + else: + return 1. + + robin_weight_flow = property(get_robin_weight_flow) + + def get_robin_weight_transport(self): + """ BoundaryCondition object + float giving the Robin factor + Solvers should rather access get_tensor(). + """ + if hasattr(self, "_robin_weight_transport"): + return self._robin_weight_transport + else: + return 1. + + robin_weight_transport = property(get_robin_weight_transport) + + def get_robin_weight_mechanics(self): + """ BoundaryCondition object + float giving the Robin factor + Solvers should rather access get_tensor(). + """ + if hasattr(self, "_robin_weight_mechanics"): + return self._robin_weight_mechanics + else: + return 1. + + robin_weight_mechanics = property(get_robin_weight_mechanics) + + def set_robin_weight(self, obj, val): + """ Set physics-specific Robin Weight + + Parameters: + + obj: Solver or str + Identification of physical regime. Either discretization object + with attribute 'physics' or a str. + + val : float, representing the weigth in the Robin boundary condition + """ + physics = self._get_physics(obj) + + if physics == "flow": + self._robin_weight_flow = val + elif physics == "transport": + self._robin_weight_transport = val + elif physics == "mechanics": + self._robin_weight_mechanics = val + else: + raise ValueError( + 'Unknown physics "%s".\n Possible physics are: %s' + % (physics, self.known_physics) + ) + # --------------------- Boundary conditions and values ------------------------ # Boundary condition @@ -765,30 +858,3 @@ def get_bc_val_mechanics(self): return np.zeros(self._num_faces * self.dim) bc_val_mechanics = property(get_bc_val_mechanics) - - # ------------------------------------ - def get_robin_factor(self, default=1): - """ double or array-like - face-wise representation of robin_factor. Set as either a np.ndarary, or a - scalar (uniform) value. Always returned as np.ndarray. - """ - if not hasattr(self, "_robin_factor"): - return default * np.ones(self._num_faces * self.dim) - - if isinstance(self._robin_factor, np.ndarray): - # Hope that the user did not initialize as array with wrong size - return self._robin_factor - else: - return self._robin_factor * np.ones(self._num_faces) - - def set_robin_factor(self, val): - """ Set physics-specific robin factor. The robin factor is given by - T + robin_factor * u = g - where T is the traction, and g is given as a boundary value - - Parameters: - val : robin_factor, representing the robin factor - """ - self._robin_factor = val - - robin_factor = property(get_robin_factor, set_robin_factor) diff --git a/test/unit/test_parameters.py b/test/unit/test_parameters.py index b0b7476d20..fb9c393fa2 100644 --- a/test/unit/test_parameters.py +++ b/test/unit/test_parameters.py @@ -189,5 +189,30 @@ def test_porosity_assertion_vector(self): ##### + def _validate_robin_weight(self, p, physics, val=1): + if physics=='flow': + self.assertTrue(np.allclose(p.robin_weight_flow, val)) + elif physics=='transport': + self.assertTrue(np.allclose(p.robin_weight_transport, val)) + elif physics=='mechanics': + self.assertTrue(np.allclose(p.robin_weight_mechanics, val)) + else: + self.assertTrue(False) + self.assertTrue(np.allclose(p.get_robin_weight(physics), val)) + + def test_robin_weight_default(self, val=1): + p = Parameters(self.g) + for name in p.known_physics: + 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: + p.set_robin_weight(name, .1) + self._validate_robin_weight(p, name, .1) + + ##### + if __name__ == "__main__": unittest.main() From 17c122995ae6c15f8ffca9b4411fe8101ec98771 Mon Sep 17 00:00:00 2001 From: Runar Date: Fri, 14 Sep 2018 13:16:24 +0200 Subject: [PATCH 23/26] Changed how the assymetric tensor is used. Now only the Robin-faces uses it --- src/porepy/numerics/fv/mpsa.py | 43 ++++++++++++++++++++-------------- test/unit/test_mpsa_robin.py | 10 ++++---- 2 files changed, 30 insertions(+), 23 deletions(-) diff --git a/src/porepy/numerics/fv/mpsa.py b/src/porepy/numerics/fv/mpsa.py index eb0222ac9e..ba61197a72 100644 --- a/src/porepy/numerics/fv/mpsa.py +++ b/src/porepy/numerics/fv/mpsa.py @@ -978,38 +978,45 @@ def mpsa_elasticity( # sides of the faces. hook = __unique_hooks_law(ncsym_all, ncasym, subcell_topology, nd) - # Pair the forces from each side - ncsym_all = subcell_topology.pair_over_subfaces_nd(ncsym_all + ncasym) - del ncasym - - ncsym = bound_exclusion.exclude_robin_dirichlet_nd(ncsym_all) - - num_subfno = subcell_topology.subfno.max() + 1 - hook_cell = sps.coo_matrix( - (np.zeros(1), (np.zeros(1), np.zeros(1))), - shape=(num_subfno * nd, (np.max(subcell_topology.cno) + 1) * nd), - ).tocsr() - - hook_cell = bound_exclusion.exclude_robin_dirichlet_nd(hook_cell) - # For the Robin boundary conditions we need to pair the forces with the - # displacement - ncsym_rob = bound_exclusion.keep_robin_nd(ncsym_all) - + # displacement. + # The contribution of the displacement at the Robin boundary is + # rob_grad * G + rob_cell * u (this is the displacement at the boundary scaled + # with the Robin weight times the area of the faces), + # where G are the gradients and u the cell center displacement. The Robin + # condtion is then ncsym_rob * G + rob_grad * G + rob_cell * u + 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( - "If applying Robin conditions you must supply an robin_weight" + "If applying Robin conditions you must supply a robin_weight" ) else: robin_weight = 1 # Book keeping num_sub_cells = cell_node_blocks[0].size + rob_grad, rob_cell = __get_displacement_submatrices_rob( g, subcell_topology, eta, num_sub_cells, bound_exclusion, robin_weight ) + # Pair the forces from each side + ncsym = subcell_topology.pair_over_subfaces_nd(ncsym_all) + del ncsym_all + ncsym = bound_exclusion.exclude_robin_dirichlet_nd(ncsym) + + num_subfno = subcell_topology.subfno.max() + 1 + hook_cell = sps.coo_matrix( + (np.zeros(1), (np.zeros(1), np.zeros(1))), + shape=(num_subfno * nd, (np.max(subcell_topology.cno) + 1) * nd), + ).tocsr() + + hook_cell = bound_exclusion.exclude_robin_dirichlet_nd(hook_cell) + # Matrices to enforce displacement continuity d_cont_grad, d_cont_cell = __get_displacement_submatrices( g, subcell_topology, eta, num_sub_cells, bound_exclusion diff --git a/test/unit/test_mpsa_robin.py b/test/unit/test_mpsa_robin.py index 559503785e..be2ca3b641 100644 --- a/test/unit/test_mpsa_robin.py +++ b/test/unit/test_mpsa_robin.py @@ -114,8 +114,8 @@ def test_structured_triang(self): right = g.face_centers[0] > 1 - 1e-10 dir_ind = np.ravel(np.argwhere(left)) - neu_ind = np.ravel(np.argwhere(top)) - rob_ind = np.ravel(np.argwhere(right + bot)) + neu_ind = np.ravel(np.argwhere(())) + rob_ind = np.ravel(np.argwhere(right + top + bot)) names = ["dir"] * len(dir_ind) + ["rob"] * len(rob_ind) bnd_ind = np.hstack((dir_ind, rob_ind)) @@ -161,9 +161,9 @@ def test_unstruct_triang(self): left = g.face_centers[0] < 1e-10 right = g.face_centers[0] > 1 - 1e-10 - dir_ind = np.ravel(np.argwhere(left)) - neu_ind = np.ravel(np.argwhere(top)) - rob_ind = np.ravel(np.argwhere(right + bot)) + dir_ind = np.ravel(np.argwhere(())) + neu_ind = np.ravel(np.argwhere(())) + rob_ind = np.ravel(np.argwhere(right + left + top + bot)) names = ["dir"] * len(dir_ind) + ["rob"] * len(rob_ind) bnd_ind = np.hstack((dir_ind, rob_ind)) From 1c4ad9394c711eda325efd61ddd6bc3ae618d17a Mon Sep 17 00:00:00 2001 From: Runar Date: Fri, 14 Sep 2018 14:11:22 +0200 Subject: [PATCH 24/26] Fixed the convergence test example for mpsa and mpfa --- examples/example2/ConvergenceTestMPFA.ipynb | 6 +- examples/example2/ConvergenceTestMPSA.ipynb | 91 +++++++++++---------- 2 files changed, 50 insertions(+), 47 deletions(-) diff --git a/examples/example2/ConvergenceTestMPFA.ipynb b/examples/example2/ConvergenceTestMPFA.ipynb index 77efc77afe..4c601f01cb 100644 --- a/examples/example2/ConvergenceTestMPFA.ipynb +++ b/examples/example2/ConvergenceTestMPFA.ipynb @@ -380,7 +380,7 @@ " \n", " # Permeability tensor\n", " k = perm(2, np.ones(g.num_cells))\n", - " alpha = 1.2\n", + " robin_weight = 1.2\n", " # Exact solution\n", " xf = g.face_centers\n", " xc = g.cell_centers\n", @@ -404,7 +404,7 @@ " bound_cond = bc.BoundaryCondition(g, bc_faces, bc_names)\n", " \n", " # MPFA discretization, and system matrix\n", - " flux, bound_flux, _, _ = mpfa.mpfa(g, k, bound_cond, alpha=alpha)\n", + " flux, bound_flux, _, _ = mpfa.mpfa(g, k, bound_cond, robin_weight=robin_weight)\n", " div = fvutils.scalar_divergence(g)\n", " a = div * flux\n", " \n", @@ -416,7 +416,7 @@ " # innflow is always negative, so need to flip flux according to normal direction.\n", " u_bound[neu_faces[nfi]] = flux_ex[neu_faces[nfi]] * (sgn_n)\n", " u_bound[rob_faces[rfi]] = flux_ex[rob_faces[rfi]] * (sgn) +\\\n", - " alpha * u_f(xf[0, rob_faces], xf[1, rob_faces]) * g.face_areas[rob_faces]\n", + " robin_weight * u_f(xf[0, rob_faces], xf[1, rob_faces]) * g.face_areas[rob_faces]\n", "\n", " # Right hand side - contribution from the solution and the boundary conditions\n", " xc = g.cell_centers\n", diff --git a/examples/example2/ConvergenceTestMPSA.ipynb b/examples/example2/ConvergenceTestMPSA.ipynb index d60ca47022..b39b629736 100644 --- a/examples/example2/ConvergenceTestMPSA.ipynb +++ b/examples/example2/ConvergenceTestMPSA.ipynb @@ -190,7 +190,8 @@ "assert np.abs(f_cart_pert[2] - 0.0053879362435559708) < 1e-10\n", "\n", "\n", - "# RB: These values were hard-coded 06.09.2018.\n", + "# RB: These values were hard-coded 06.09.2018. The values changed because\n", + "# the default value for eta for triangular grids was changed to 1/3\n", "assert np.abs(u_triang_nopert[2] - 0.0005339430679082081) < 1e-10\n", "assert np.abs(f_triang_nopert[2] - 0.005076956913423936) < 1e-10\n", "assert np.abs(u_triang_pert[2] - 0.0006930802706380639) < 1e-10\n", @@ -216,18 +217,18 @@ "output_type": "stream", "text": [ "Cartesian errors - no perturbations\n", - "[0.18471065947993767, 0.04467625698150768, 0.01073709992570399] displacement error\n", - "[4.13442558 4.16092402] convergence rate u^n-1 / u^n\n", + "[0.19341322457532645, 0.05398606173752483, 0.019785099166683015] displacement error\n", + "[3.58265112 2.72862225] convergence rate u^n-1 / u^n\n", "\n", - "[0.06476977708690358, 0.021191370180076885, 0.006402282045426171] stress error\n", - "[3.05642233 3.30997136] convergence rate f^n-1 / f^n\n", + "[0.03080448569971605, 0.022727907422006446, 0.017633923832997576] stress error\n", + "[1.35535952 1.28887408] convergence rate f^n-1 / f^n\n", "\n", "Triangular errors - no perturbations\n", - "[0.16176733322360393, 0.03341692740958871, 0.007161359044861829] displacement error\n", - "[4.84087993 4.66628292] convergence rate u^n-1 / u^n\n", + "[0.15803686415889875, 0.04310823056510322, 0.015475821210292553] displacement error\n", + "[3.6660485 2.78552136] convergence rate u^n-1 / u^n\n", "\n", - "[0.06486596857883836, 0.02141863302701911, 0.006979068019066017] stress error\n", - "[3.02848312 3.06898184] convergence rate f^n-1 / f^n\n" + "[0.04467366719714846, 0.02540025474312989, 0.016373893423034032] stress error\n", + "[1.75878816 1.55126543] convergence rate f^n-1 / f^n\n" ] } ], @@ -304,8 +305,8 @@ " left = np.ravel(np.argwhere(g.face_centers[0, :] < 1e-10))\n", " right = np.ravel(np.argwhere(g.face_centers[0, :] > n - 1e-10))\n", "\n", - " dir_faces = np.hstack((left, bot, top))\n", - " neu_faces = np.hstack((right))\n", + " dir_faces = np.hstack((bot, top))\n", + " neu_faces = np.hstack((left, right)) \n", " bound_cond = bc.BoundaryCondition(g, dir_faces, ['dir'] * dir_faces.size)\n", " \n", " # MPFA discretization, and system matrix\n", @@ -379,11 +380,13 @@ "print(f_triang_nopert,'stress error')\n", "print(f_triang_rate, 'convergence rate f^n-1 / f^n')\n", "\n", - "# RB: These values were hard-coded 11.09.2018.\n", - "assert np.abs(u_cart_nopert[2] - 0.01073709992570399) < 1e-10\n", - "assert np.abs(f_cart_nopert[2] - 0.006402282045426171) < 1e-10\n", - "assert np.abs(u_triang_nopert[2] - 0.007161359044861829) < 1e-10\n", - "assert np.abs(f_triang_nopert[2] - 0.006979068019066017) < 1e-10" + "# # EK: These values were hard-coded 12.06.2017.\n", + "assert np.abs(u_cart_nopert[2] - 0.01978509916668196) < 1e-10\n", + "assert np.abs(f_cart_nopert[2] - 0.017633923832997497) < 1e-10\n", + "\n", + "# RB: These values were hard-coded 06.09.2018.\n", + "assert np.abs(u_triang_nopert[2] - 0.015475821210292553) < 1e-10\n", + "assert np.abs(f_triang_nopert[2] - 0.016373893423034032) < 1e-10" ] }, { @@ -482,7 +485,7 @@ " mu_c = mu * np.ones(g.num_cells)\n", " lmbda_c = lmbda * np.ones(g.num_cells)\n", " k = stiffness(2, mu_c, lmbda_c)\n", - " alpha = 1.0\n", + " robin_weight = 1.0\n", " # Set type of boundary conditions - Dirichlet\n", " n = g.nodes.max()\n", " # Set type of boundary conditions - Dirichlet\n", @@ -498,7 +501,7 @@ " bound_cond = bc.BoundaryCondition(g, bc_faces, bc_names)\n", " \n", " # MPFA discretization, and system matrix\n", - " stress, bound_stress = mpsa.mpsa(g, k, bound_cond,alpha=alpha)\n", + " stress, bound_stress = mpsa.mpsa(g, k, bound_cond,robin_weight=robin_weight)\n", " div = fvutils.vector_divergence(g)\n", " a = div * stress\n", " xf = g.face_centers\n", @@ -513,10 +516,10 @@ " # sign of the faces with a normal pointing inwards as it is assumed\n", " # that the Neumann condition is the force from the boundary on to face\n", " u_bound[0, rob_faces[rfi]] = stress_x_ex[rob_faces[rfi]] * (sgn) + \\\n", - " alpha * ux_f(xf[0, rob_faces[rfi]], xf[1, rob_faces[rfi]]) * g.face_areas[rob_faces[rfi]]\n", + " robin_weight * ux_f(xf[0, rob_faces[rfi]], xf[1, rob_faces[rfi]]) * g.face_areas[rob_faces[rfi]]\n", "\n", " u_bound[1, rob_faces[rfi]] = stress_y_ex[rob_faces[rfi]] * (sgn) +\\\n", - " alpha * uy_f(xf[0, rob_faces[rfi]], xf[1, rob_faces[rfi]]) * g.face_areas[rob_faces[rfi]]\n", + " robin_weight * uy_f(xf[0, rob_faces[rfi]], xf[1, rob_faces[rfi]]) * g.face_areas[rob_faces[rfi]]\n", " # Right hand side - contribution from the solution and the boundary conditions\n", " xc = g.cell_centers\n", " rhs = np.vstack((rhs_x_f(xc[0], xc[1]), rhs_y_f(xc[0], xc[1]))) * g.cell_volumes\n", @@ -1013,39 +1016,39 @@ "name": "stdout", "output_type": "stream", "text": [ - "Time for discretization 1.8826446533203125\n", - "Time for AMG 0.04243111610412598\n", - "Time for discretization 2.377439260482788\n", - "Time for AMG 0.1383042335510254\n", - "Time for discretization 2.542255401611328\n", - "Time for AMG 0.3778951168060303\n", + "Time for discretization 1.08636474609375\n", + "Time for AMG 0.020662784576416016\n", + "Time for discretization 1.6143503189086914\n", + "Time for AMG 0.07893586158752441\n", + "Time for discretization 1.8830769062042236\n", + "Time for AMG 0.4144134521484375\n", "Cartesian errors - no perturbations\n", "[0.14582494480280667, 0.036685508143867125, 0.009623477148373078] displacement error\n", "[0.14173578960074407, 0.05453461058511372, 0.018777567804334296] stress error\n", - "Time for discretization 1.409085750579834\n", - "Time for AMG 0.015498161315917969\n", - "Time for discretization 1.3807973861694336\n", - "Time for AMG 0.07768678665161133\n", - "Time for discretization 7.617006778717041\n", - "Time for AMG 0.3285984992980957\n", + "Time for discretization 1.1395857334136963\n", + "Time for AMG 0.010153055191040039\n", + "Time for discretization 1.394442081451416\n", + "Time for AMG 0.0783224105834961\n", + "Time for discretization 1.727994680404663\n", + "Time for AMG 0.37769436836242676\n", "Cartesian errors - perturbation 0.5\n", "[0.1823425575505196, 0.03669658584037131, 0.009785391062793553] displacement error\n", "[0.15394769153699012, 0.06655533666156926, 0.024831033744125167] stress error\n", - "Time for discretization 1.9705743789672852\n", - "Time for AMG 0.07087302207946777\n", - "Time for discretization 7.326463222503662\n", - "Time for AMG 0.4231894016265869\n", - "Time for discretization 10.036577224731445\n", - "Time for AMG 5.2540881633758545\n", + "Time for discretization 1.5092995166778564\n", + "Time for AMG 0.04783058166503906\n", + "Time for discretization 2.5991227626800537\n", + "Time for AMG 0.35640907287597656\n", + "Time for discretization 5.721246719360352\n", + "Time for AMG 4.348304986953735\n", "Triangular errors - no perturbations\n", "[0.06256710001786384, 0.015518298189089326, 0.004112646550453208] displacement error\n", "[0.09819533736107755, 0.040432049644100486, 0.01608810670167512] stress error\n", - "Time for discretization 1.2782678604125977\n", - "Time for AMG 0.07084035873413086\n", - "Time for discretization 3.013735055923462\n", - "Time for AMG 0.4166250228881836\n", - "Time for discretization 10.38901948928833\n", - "Time for AMG 5.382190704345703\n", + "Time for discretization 1.1475048065185547\n", + "Time for AMG 0.04819488525390625\n", + "Time for discretization 1.789534568786621\n", + "Time for AMG 0.342240571975708\n", + "Time for discretization 6.363281965255737\n", + "Time for AMG 4.3170764446258545\n", "Triangular errors - perturbation 0.5\n", "[0.07257865039340129, 0.01699696434861487, 0.004336430438721095] displacement error\n", "[0.10603237914286788, 0.04292864828702438, 0.017041050867287796] stress error\n" From 589c225cead814e1b0988be0dbbaaa28d8ed211f Mon Sep 17 00:00:00 2001 From: Runar Date: Fri, 14 Sep 2018 14:41:37 +0200 Subject: [PATCH 25/26] Wrote some comments for mpsa --- src/porepy/numerics/fv/mpfa.py | 4 +++- src/porepy/numerics/fv/mpsa.py | 12 ++++++++++-- test/unit/test_parameters.py | 7 ++++--- 3 files changed, 17 insertions(+), 6 deletions(-) diff --git a/src/porepy/numerics/fv/mpfa.py b/src/porepy/numerics/fv/mpfa.py index 46a9fdd828..1ac8cde3b2 100644 --- a/src/porepy/numerics/fv/mpfa.py +++ b/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 diff --git a/src/porepy/numerics/fv/mpsa.py b/src/porepy/numerics/fv/mpsa.py index ba61197a72..e3d8560134 100644 --- a/src/porepy/numerics/fv/mpsa.py +++ b/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) diff --git a/test/unit/test_parameters.py b/test/unit/test_parameters.py index fb9c393fa2..b970496879 100644 --- a/test/unit/test_parameters.py +++ b/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: From 173beffe937fa7fc6727f77875bc769da934137a Mon Sep 17 00:00:00 2001 From: Runar Date: Fri, 14 Sep 2018 18:54:52 +0200 Subject: [PATCH 26/26] Checkouted test_bc_vectorial from develop and fixed a bug in fvutils --- src/porepy/numerics/fv/fvutils.py | 2 +- test/unit/test_bc_vectorial.py | 127 ++++++++++++++++++++++++++++++ 2 files changed, 128 insertions(+), 1 deletion(-) create mode 100644 test/unit/test_bc_vectorial.py diff --git a/src/porepy/numerics/fv/fvutils.py b/src/porepy/numerics/fv/fvutils.py index 2f7f24224e..e9c7b5a92b 100644 --- a/src/porepy/numerics/fv/fvutils.py +++ b/src/porepy/numerics/fv/fvutils.py @@ -756,7 +756,7 @@ def __exclude_matrix_xyz(self, ids): col_neu = np.append(col_neu, [col_z]) else: - exclude_z = sps.coo_matrix(shape=(0, self.num_subfno)) + exclude_z = sps.coo_matrix([], shape=(0, self.num_subfno)) row_neu = np.arange(col_neu.size) exclude_nd = sps.coo_matrix( (np.ones(row_neu.size), (row_neu, col_neu.ravel("C"))), diff --git a/test/unit/test_bc_vectorial.py b/test/unit/test_bc_vectorial.py new file mode 100644 index 0000000000..080032947e --- /dev/null +++ b/test/unit/test_bc_vectorial.py @@ -0,0 +1,127 @@ +import unittest +import numpy as np + +import math + +from porepy.grids import structured +from porepy.params.bc import BoundaryConditionVectorial +from porepy.numerics.fv import fvutils + +""" +Checks the actions done in porepy.numerics.fv.mpsa.create_bound_rhs_nd +for handling boundary conditions expressed in a vectorial form +""" + + +class testBoundaryConditionsVectorial(unittest.TestCase): + + def test_2d(self): + + """ + The domain consists of a 3x3 mesh. + Bottom faces are dirichlet + Left and right faces are rolling along y (dir_x and neu_y) + Top mid face is rolling along x (neu_x and dir_y) + Top left and right faces are neumann + """ + + g = structured.CartGrid([3, 3]) + g.compute_geometry() + nd = g.dim + + boundary_faces = np.array([0, 3, 4, 7, 8, 11, 12, 13, 14, 22]) + boundary_faces_type = ['dir_x'] * 6 + ['dir'] * 3 + ['dir_y'] * 1 + + bound = BoundaryConditionVectorial(g, boundary_faces, boundary_faces_type) + + subcell_topology = fvutils.SubcellTopology(g) + fno = subcell_topology.fno_unique + bound_exclusion = fvutils.ExcludeBoundaries(subcell_topology, bound, nd) + + # Neumann + is_neu_x = bound_exclusion.exclude_dirichlet_x(bound.is_neu[0, fno].astype( + 'int64')) + neu_ind_single_x = np.argwhere(is_neu_x).ravel('F') + + is_neu_y = bound_exclusion.exclude_dirichlet_y(bound.is_neu[1, fno].astype( + 'int64')) + neu_ind_single_y = np.argwhere(is_neu_y).ravel('F') + neu_ind_single_y += is_neu_x.size + + # We also need to account for all half faces, that is, do not exclude + # Dirichlet and Neumann boundaries. + neu_ind_single_all_x = np.argwhere(bound.is_neu[0, fno].astype('int'))\ + .ravel('F') + neu_ind_single_all_y = np.argwhere(bound.is_neu[1, fno].astype('int'))\ + .ravel('F') + + # expand the indices + # this procedure replaces the method 'expand_ind' in the above + # method 'create_bound_rhs' + + # 1 - stack and sort indices + + is_bnd_neu_x = nd * neu_ind_single_all_x + is_bnd_neu_y = nd * neu_ind_single_all_y + 1 + + is_bnd_neu = np.sort(np.append(is_bnd_neu_x, [is_bnd_neu_y])) + + # 2 - find the indices corresponding to the boundary components + # having Neumann condtion + + ind_is_bnd_neu_x = np.argwhere(np.isin(is_bnd_neu, is_bnd_neu_x)).ravel('F') + ind_is_bnd_neu_y = np.argwhere(np.isin(is_bnd_neu, is_bnd_neu_y)).ravel('F') + + neu_ind_sz = ind_is_bnd_neu_x.size + ind_is_bnd_neu_y.size + + # 3 - create the expanded neu_ind array + + neu_ind = np.zeros(neu_ind_sz, dtype='int') + + neu_ind[ind_is_bnd_neu_x] = neu_ind_single_x + neu_ind[ind_is_bnd_neu_y] = neu_ind_single_y + + assert np.alltrue(neu_ind == [30, 31, 36, 37, 38, 39, 44, 45, 46, 47, 52, 53, 24, 66, 25, 67, 26, 27, 28, 68, 29, 69]) + + # Dirichlet, same procedure + is_dir_x = bound_exclusion.exclude_neumann_x(bound.is_dir[0, fno].astype( + 'int64')) + dir_ind_single_x = np.argwhere(is_dir_x).ravel('F') + + is_dir_y = bound_exclusion.exclude_neumann_y(bound.is_dir[1, fno].astype( + 'int64')) + dir_ind_single_y = np.argwhere(is_dir_y).ravel('F') + dir_ind_single_y += is_dir_x.size + + dir_ind_single_all_x = np.argwhere(bound.is_dir[0, fno].astype('int'))\ + .ravel('F') + dir_ind_single_all_y = np.argwhere(bound.is_dir[1, fno].astype('int'))\ + .ravel('F') + + #expand indices + + is_bnd_dir_x = nd * dir_ind_single_all_x + is_bnd_dir_y = nd * dir_ind_single_all_y + 1 + + is_bnd_dir = np.sort(np.append(is_bnd_dir_x, [is_bnd_dir_y])) + + ind_is_bnd_dir_x = np.argwhere(np.isin(is_bnd_dir, is_bnd_dir_x)).ravel('F') + ind_is_bnd_dir_y = np.argwhere(np.isin(is_bnd_dir, is_bnd_dir_y)).ravel('F') + + dir_ind_sz = ind_is_bnd_dir_x.size + ind_is_bnd_dir_y.size + + dir_ind = np.zeros(dir_ind_sz, dtype='int') + + dir_ind[ind_is_bnd_dir_x] = dir_ind_single_x + dir_ind[ind_is_bnd_dir_y] = dir_ind_single_y + + assert np.alltrue(dir_ind == [0, 1, 6, 7, 8, 9, 14, 15, 16, 17, 22, 23, 24, 54, 25, 55, 26, 56, 27, 57, 28, 58, 29, 59, 72, 73]) + +if __name__ == '__main__': + unittest.main() + + + + + +