Skip to content

Commit

Permalink
Merge branch 'develop' into primal_contact
Browse files Browse the repository at this point in the history
  • Loading branch information
@keileg
keileg committed Jul 2, 2019
2 parents af5f74b + 0f04e8e commit 1c9cffa
Show file tree
Hide file tree
Showing 2 changed files with 138 additions and 35 deletions.
126 changes: 91 additions & 35 deletions src/porepy/numerics/fv/upwind.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,12 +40,12 @@ class Upwind:
"""

# ------------------------------------------------------------------------------#

def __init__(self, keyword="transport"):
self.keyword = keyword

# ------------------------------------------------------------------------------#
# Keywords used to store matrix and right hand side in the matrix_dictionary
self.matrix_keyword = "transport"
self.rhs_keyword = "rhs"

def ndof(self, g):
"""
Expand All @@ -63,9 +63,63 @@ def ndof(self, g):
"""
return g.num_cells

def assemble_matrix_rhs(self, g, data):
""" Return the matrix for an upwind discretization of a linear transport
problem.
Parameters:
g (Grid): Computational grid, with geometry fields computed.
data (dictionary): With data stored.
Returns:
scipy.sparse.csr_matrix: System matrix of this discretization.
Size: g.num_cells x g.num_cells.
np.ndarray: Right hand side vector with representation of boundary
conditions. The size of the vector will depend on the discretization.
"""
matrix_dictionary = data[pp.DISCRETIZATION_MATRICES][self.keyword]
if not self.matrix_keyword in matrix_dictionary.keys():
self.discretize(g, data)

return self.assemble_matrix(g, data), self.assemble_rhs(g, data)

def assemble_matrix(self, g, data):
""" Return the matrix for an upwind discretization of a linear transport
problem.
Parameters:
g (Grid): Computational grid, with geometry fields computed.
data (dictionary): With data stored.
Returns:
scipy.sparse.csr_matrix: System matrix of this discretization.
Size: g.num_cells x g.num_cells.
"""
matrix_dictionary = data[pp.DISCRETIZATION_MATRICES][self.keyword]
return matrix_dictionary[self.matrix_keyword]

# ------------------------------------------------------------------------------#

def assemble_matrix_rhs(self, g, data, d_name="darcy_flux"):
def assemble_rhs(self, g, data):
""" Return the right-hand side for an upwind discretization of a linear
transport problem.
Parameters:
g (Grid): Computational grid, with geometry fields computed.
data (dictionary): With data stored.
Returns:
np.ndarray: Right hand side vector with representation of boundary
conditions. The size of the vector will depend on the discretization.
"""
matrix_dictionary = data[pp.DISCRETIZATION_MATRICES][self.keyword]

return matrix_dictionary[self.rhs_keyword]

def discretize(self, g, data, d_name="darcy_flux"):
"""
Return the matrix and righ-hand side for a discretization of a scalar
linear transport problem using the upwind scheme.
Expand All @@ -83,6 +137,7 @@ def assemble_matrix_rhs(self, g, data, d_name="darcy_flux"):
following keys: 'dir' and 'neu', for Dirichlet and Neumann boundary
conditions, respectively.
source : array (g.num_cells) of source (positive) or sink (negative) terms.
Parameters
----------
g : grid, or a subclass, with geometry fields computed.
Expand Down Expand Up @@ -112,11 +167,16 @@ def assemble_matrix_rhs(self, g, data, d_name="darcy_flux"):
for i in np.arange( N ):
conc = invM.dot((M_minus_U).dot(conc) + rhs)
"""
if g.dim == 0:
data["flow_faces"] = sps.csr_matrix([0.0])
return sps.csr_matrix([0.0]), np.array([0.0])

parameter_dictionary = data[pp.PARAMETERS][self.keyword]
matrix_dictionary = data[pp.DISCRETIZATION_MATRICES][self.keyword]

# Shortcut for point grids
if g.dim == 0:
matrix_dictionary[self.matrix_keyword] = sps.csr_matrix([0.0])
matrix_dictionary[self.rhs_keyword] = np.array([0.0])
return

darcy_flux = parameter_dictionary[d_name]
bc = parameter_dictionary["bc"]
bc_val = parameter_dictionary["bc_values"]
Expand Down Expand Up @@ -165,31 +225,31 @@ def assemble_matrix_rhs(self, g, data, d_name="darcy_flux"):
flow_cells.tocsr()
flow_cells = flow_cells.astype(np.float)

data["flow_faces"] = flow_faces
if not has_bc:
return flow_cells, np.zeros(g.num_cells)

# Impose the boundary conditions
bc_val_dir = np.zeros(g.num_faces)
if np.any(bc.is_dir):
is_dir = np.where(bc.is_dir)[0]
bc_val_dir[is_dir] = bc_val[is_dir]

# We assume that for Neumann boundary condition a positive 'bc_val'
# represents an outflow for the domain. A negative 'bc_val' represents
# an inflow for the domain.
bc_val_neu = np.zeros(g.num_faces)
if np.any(bc.is_neu):
is_neu = np.where(bc.is_neu)[0]
bc_val_neu[is_neu] = bc_val[is_neu]

return (
flow_cells,
-inflow.transpose() * bc_val_dir
- np.abs(g.cell_faces.transpose()) * bc_val_neu,
)
# Store disrcetization matrix
matrix_dictionary[self.matrix_keyword] = flow_cells

# ------------------------------------------------------------------------------#
if not has_bc:
# Short cut if there are no trivial boundary conditions
matrix_dictionary[self.rhs_keyword] = np.zeros(g.num_cells)
else:
# Impose the boundary conditions
bc_val_dir = np.zeros(g.num_faces)
if np.any(bc.is_dir):
is_dir = np.where(bc.is_dir)[0]
bc_val_dir[is_dir] = bc_val[is_dir]

# We assume that for Neumann boundary condition a positive 'bc_val'
# represents an outflow for the domain. A negative 'bc_val' represents
# an inflow for the domain.
bc_val_neu = np.zeros(g.num_faces)
if np.any(bc.is_neu):
is_neu = np.where(bc.is_neu)[0]
bc_val_neu[is_neu] = bc_val[is_neu]

matrix_dictionary[self.rhs_keyword] = (
-inflow.transpose() * bc_val_dir
- np.abs(g.cell_faces.transpose()) * bc_val_neu
)

def cfl(self, g, data, d_name="darcy_flux"):
"""
Expand Down Expand Up @@ -241,8 +301,6 @@ def cfl(self, g, data, d_name="darcy_flux"):
# deltaT is deltaX/darcy_flux with coefficient
return np.amin(np.abs(np.divide(dist, darcy_flux[faces])) * coeff)

# ------------------------------------------------------------------------------#

def darcy_flux(self, g, beta, cell_apertures=None):
"""
Return the normal component of the velocity, for each face, weighted by
Expand Down Expand Up @@ -278,8 +336,6 @@ def darcy_flux(self, g, beta, cell_apertures=None):
[np.dot(n, a * beta) for n, a in zip(g.face_normals.T, face_apertures)]
)

# ------------------------------------------------------------------------------#

def outflow(self, g, data, d_name="darcy_flux"):
if g.dim == 0:
return sps.csr_matrix([0])
Expand Down
47 changes: 47 additions & 0 deletions test/unit/test_upwind.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,8 @@ def test_upwind_1d_darcy_flux_positive(self):
specified_parameters = {"bc": bc, "darcy_flux": dis}
data = pp.initialize_default_data(g, {}, "transport", specified_parameters)

solver.discretize(g, data)

M = solver.assemble_matrix_rhs(g, data)[0].todense()
deltaT = solver.cfl(g, data)

Expand All @@ -50,6 +52,9 @@ def test_upwind_1d_darcy_flux_negative(self):
bc = BoundaryCondition(g, bf, bf.size * ["neu"])
specified_parameters = {"bc": bc, "darcy_flux": dis}
data = pp.initialize_default_data(g, {}, "transport", specified_parameters)

solver.discretize(g, data)

M = solver.assemble_matrix_rhs(g, data)[0].todense()
deltaT = solver.cfl(g, data)

Expand All @@ -74,6 +79,9 @@ def test_upwind_2d_cart_darcy_flux_positive(self):
bc = BoundaryCondition(g, bf, bf.size * ["neu"])
specified_parameters = {"bc": bc, "darcy_flux": dis}
data = pp.initialize_default_data(g, {}, "transport", specified_parameters)

solver.discretize(g, data)

M = solver.assemble_matrix_rhs(g, data)[0].todense()
deltaT = solver.cfl(g, data)

Expand Down Expand Up @@ -108,6 +116,9 @@ def test_upwind_2d_cart_darcy_flux_negative(self):
bc = BoundaryCondition(g, bf, bf.size * ["neu"])
specified_parameters = {"bc": bc, "darcy_flux": dis}
data = pp.initialize_default_data(g, {}, "transport", specified_parameters)

solver.discretize(g, data)

M = solver.assemble_matrix_rhs(g, data)[0].todense()
deltaT = solver.cfl(g, data)

Expand Down Expand Up @@ -141,6 +152,9 @@ def test_upwind_2d_simplex_darcy_flux_positive(self):
bc = BoundaryCondition(g, bf, bf.size * ["neu"])
specified_parameters = {"bc": bc, "darcy_flux": dis}
data = pp.initialize_default_data(g, {}, "transport", specified_parameters)

solver.discretize(g, data)

M = solver.assemble_matrix_rhs(g, data)[0].todense()
deltaT = solver.cfl(g, data)

Expand All @@ -165,6 +179,9 @@ def test_upwind_2d_simplex_darcy_flux_negative(self):
bc = BoundaryCondition(g, bf, bf.size * ["neu"])
specified_parameters = {"bc": bc, "darcy_flux": dis}
data = pp.initialize_default_data(g, {}, "transport", specified_parameters)

solver.discretize(g, data)

M = solver.assemble_matrix_rhs(g, data)[0].todense()
deltaT = solver.cfl(g, data)

Expand All @@ -189,6 +206,9 @@ def test_upwind_3d_cart_darcy_flux_negative(self):
bc = BoundaryCondition(g, bf, bf.size * ["neu"])
specified_parameters = {"bc": bc, "darcy_flux": dis}
data = pp.initialize_default_data(g, {}, "transport", specified_parameters)

solver.discretize(g, data)

M = solver.assemble_matrix_rhs(g, data)[0].todense()
deltaT = solver.cfl(g, data)

Expand Down Expand Up @@ -225,6 +245,9 @@ def test_upwind_3d_cart_darcy_flux_positive(self):
bc = BoundaryCondition(g, bf, bf.size * ["neu"])
specified_parameters = {"bc": bc, "darcy_flux": dis}
data = pp.initialize_default_data(g, {}, "transport", specified_parameters)

solver.discretize(g, data)

M = solver.assemble_matrix_rhs(g, data)[0].todense()
deltaT = solver.cfl(g, data)

Expand Down Expand Up @@ -263,6 +286,9 @@ def test_upwind_1d_surf_darcy_flux_positive(self):
bc = BoundaryCondition(g, bf, bf.size * ["neu"])
specified_parameters = {"bc": bc, "darcy_flux": dis}
data = pp.initialize_default_data(g, {}, "transport", specified_parameters)

solver.discretize(g, data)

M = solver.assemble_matrix_rhs(g, data)[0].todense()
deltaT = solver.cfl(g, data)

Expand All @@ -289,6 +315,9 @@ def test_upwind_1d_surf_darcy_flux_negative(self):
bc = BoundaryCondition(g, bf, bf.size * ["neu"])
specified_parameters = {"bc": bc, "darcy_flux": dis}
data = pp.initialize_default_data(g, {}, "transport", specified_parameters)

solver.discretize(g, data)

M = solver.assemble_matrix_rhs(g, data)[0].todense()
deltaT = solver.cfl(g, data)

Expand All @@ -315,6 +344,9 @@ def test_upwind_2d_cart_surf_darcy_flux_positive(self):
bc = BoundaryCondition(g, bf, bf.size * ["neu"])
specified_parameters = {"bc": bc, "darcy_flux": dis}
data = pp.initialize_default_data(g, {}, "transport", specified_parameters)

solver.discretize(g, data)

M = solver.assemble_matrix_rhs(g, data)[0].todense()
deltaT = solver.cfl(g, data)

Expand Down Expand Up @@ -351,6 +383,9 @@ def test_upwind_2d_cart_surf_darcy_flux_negative(self):
bc = BoundaryCondition(g, bf, bf.size * ["neu"])
specified_parameters = {"bc": bc, "darcy_flux": dis}
data = pp.initialize_default_data(g, {}, "transport", specified_parameters)

solver.discretize(g, data)

M = solver.assemble_matrix_rhs(g, data)[0].todense()
deltaT = solver.cfl(g, data)

Expand Down Expand Up @@ -387,6 +422,9 @@ def test_upwind_2d_simplex_surf_darcy_flux_positive(self):
bc = BoundaryCondition(g, bf, bf.size * ["neu"])
specified_parameters = {"bc": bc, "darcy_flux": dis}
data = pp.initialize_default_data(g, {}, "transport", specified_parameters)

solver.discretize(g, data)

M = solver.assemble_matrix_rhs(g, data)[0].todense()
deltaT = solver.cfl(g, data)

Expand All @@ -413,6 +451,9 @@ def test_upwind_2d_simplex_surf_darcy_flux_negative(self):
bc = BoundaryCondition(g, bf, bf.size * ["neu"])
specified_parameters = {"bc": bc, "darcy_flux": dis}
data = pp.initialize_default_data(g, {}, "transport", specified_parameters)

solver.discretize(g, data)

M = solver.assemble_matrix_rhs(g, data)[0].todense()
deltaT = solver.cfl(g, data)

Expand All @@ -438,6 +479,9 @@ def test_upwind_1d_darcy_flux_negative_bc_dir(self):
bc_val = 3 * np.ones(g.num_faces).ravel("F")
specified_parameters = {"bc": bc, "bc_values": bc_val, "darcy_flux": dis}
data = pp.initialize_default_data(g, {}, "transport", specified_parameters)

solver.discretize(g, data)

M, rhs = solver.assemble_matrix_rhs(g, data)
deltaT = solver.cfl(g, data)

Expand Down Expand Up @@ -465,6 +509,9 @@ def test_upwind_1d_darcy_flux_negative_bc_neu(self):
bc_val = np.array([2, 0, 0, -2]).ravel("F")
specified_parameters = {"bc": bc, "bc_values": bc_val, "darcy_flux": dis}
data = pp.initialize_default_data(g, {}, "transport", specified_parameters)

solver.discretize(g, data)

M, rhs = solver.assemble_matrix_rhs(g, data)
deltaT = solver.cfl(g, data)

Expand Down

0 comments on commit 1c9cffa

Please sign in to comment.