Skip to content

Commit

Permalink
Merge c834d5b into 908d7a0
Browse files Browse the repository at this point in the history
  • Loading branch information
rbe051 committed Sep 17, 2018
2 parents 908d7a0 + c834d5b commit a892127
Show file tree
Hide file tree
Showing 5 changed files with 834 additions and 200 deletions.
203 changes: 59 additions & 144 deletions src/porepy/numerics/fv/mpsa.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1592,7 +1592,6 @@ def expand_ind(ind, dim, increment):

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(
(np.ones(num_bound), (np.arange(num_bound), bnd_ind)),
shape=(num_bound, num_subfno * nd),
Expand All @@ -1619,150 +1618,77 @@ def create_bound_rhs_nd(bound, bound_exclusion, subcell_topology, g):
"""
nd = g.dim

num_stress = (
bound_exclusion.exclude_dir_x.shape[0] + bound_exclusion.exclude_dir_y.shape[0]
)
num_displ = (
bound_exclusion.exclude_neu_x.shape[0] + bound_exclusion.exclude_neu_y.shape[0]
)
if nd == 3:
num_stress += bound_exclusion.exclude_dir_z.shape[0]
num_displ += bound_exclusion.exclude_neu_z.shape[0]
num_stress = bound_exclusion.exclude_rob_dir_nd.shape[0]
num_displ = bound_exclusion.exclude_neu_rob_nd.shape[0]

num_rob = bound_exclusion.keep_rob_nd.shape[0]

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[0, fno]) + sum(bound.is_neu[1, fno])
num_dir = sum(bound.is_dir[0, fno]) + sum(bound.is_dir[1, fno])
if nd == 3:
num_neu += sum(bound.is_neu[2, fno])
num_dir += sum(bound.is_dir[2, fno])
num_neu = np.sum(bound.is_neu[:, fno])
num_dir = np.sum(bound.is_dir[:, fno])
if not num_rob == np.sum(bound.is_rob[:, fno]):
raise AssertionError()

num_bound = num_neu + num_dir
num_bound = num_neu + num_dir + num_rob

# Define right hand side for Neumann boundary conditions
# First row indices in rhs matrix
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")

neu_ind_all = np.append(neu_ind_single_all_x, [neu_ind_single_all_y])

# 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]))

if nd == 3:
is_neu_z = bound_exclusion.exclude_dirichlet_z(
bound.is_neu[2, fno].astype("int64")
)
neu_ind_single_z = np.argwhere(is_neu_z).ravel("F")
neu_ind_single_z += is_neu_x.size + is_neu_y.size

neu_ind_single_all_z = np.argwhere(bound.is_neu[2, fno].astype("int")).ravel(
"F"
)

neu_ind_all = np.append(neu_ind_all, [neu_ind_single_all_z])

is_bnd_neu_z = nd * neu_ind_single_all_z + 2
is_bnd_neu = np.sort(np.append(is_bnd_neu, [is_bnd_neu_z]))

# 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

if nd == 3:
ind_is_bnd_neu_z = np.argwhere(np.isin(is_bnd_neu, is_bnd_neu_z)).ravel("F")
neu_ind_sz += ind_is_bnd_neu_z.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
if nd == 3:
neu_ind[ind_is_bnd_neu_z] = neu_ind_single_z
is_neu_nd = bound_exclusion.exclude_robin_dirichlet_nd(
bound.is_neu[:, fno].ravel("C")
).ravel("F")
neu_ind_nd = np.argwhere(is_neu_nd).ravel("F")
neu_ind = np.reshape(neu_ind_nd, (nd, -1), order="C").ravel("F")

# Robin, same procedure
is_rob_nd = bound_exclusion.keep_robin_nd(bound.is_rob[:, fno].ravel("C")).ravel(
"F"
)
rob_ind_nd = np.argwhere(is_rob_nd).ravel("F")
rob_ind = np.reshape(rob_ind_nd, (nd, -1), order="C").ravel("F")

# 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_dir_nd = bound_exclusion.exclude_neumann_robin_nd(
bound.is_dir[:, fno].ravel("C")
).ravel("F")
dir_ind_nd = np.argwhere(is_dir_nd).ravel("F")
dir_ind = np.reshape(dir_ind_nd, (nd, -1), order="C").ravel("F")

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]))

if nd == 3:
is_dir_z = bound_exclusion.exclude_neumann_z(
bound.is_dir[2, fno].astype("int64")
)
dir_ind_single_z = np.argwhere(is_dir_z).ravel("F")
dir_ind_single_z += is_dir_x.size + is_dir_y.size

dir_ind_single_all_z = np.argwhere(bound.is_dir[2, fno].astype("int")).ravel(
"F"
)

is_bnd_dir_z = nd * dir_ind_single_all_z + 2
is_bnd_dir = np.sort(np.append(is_bnd_dir, [is_bnd_dir_z]))

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

if nd == 3:
ind_is_bnd_dir_z = np.argwhere(np.isin(is_bnd_dir, is_bnd_dir_z)).ravel("F")
dir_ind_sz += ind_is_bnd_dir_z.size

dir_ind = np.zeros(dir_ind_sz, dtype="int")
# We also need to account for all half faces, that is, do not exclude
# Dirichlet and Neumann boundaries. This is the global indexing.
is_neu_all = bound.is_neu[:, fno].ravel("C")
neu_ind_all = np.argwhere(
np.reshape(is_neu_all, (nd, -1), order="C").ravel("F")
).ravel("F")
is_dir_all = bound.is_dir[:, fno].ravel("C")
dir_ind_all = np.argwhere(
np.reshape(is_dir_all, (nd, -1), order="C").ravel("F")
).ravel("F")
is_rob_all = bound.is_rob[:, fno].ravel("C")
rob_ind_all = np.argwhere(
np.reshape(is_rob_all, (nd, -1), order="C").ravel("F")
).ravel("F")

dir_ind[ind_is_bnd_dir_x] = dir_ind_single_x
dir_ind[ind_is_bnd_dir_y] = dir_ind_single_y
# 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))

if nd == 3:
dir_ind[ind_is_bnd_dir_z] = dir_ind_single_z
neu_rob_ind_all = np.hstack((neu_ind_all, rob_ind_all))

# stack together
bnd_ind = np.hstack((is_bnd_neu, is_bnd_dir))
bnd_ind = np.hstack((neu_rob_ind_all, dir_ind_all))

# Some care is needed to compute coefficients in Neumann matrix: sgn is
# already defined according to the subcell topology [fno], while areas
# must be drawn from the grid structure, and thus go through fno

fno_ext = np.tile(fno, nd)
num_face_nodes = g.face_nodes.sum(axis=0).A.ravel("F")

Expand All @@ -1771,41 +1697,31 @@ def create_bound_rhs_nd(bound, bound_exclusion, subcell_topology, g):
# 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]]
neu_val = 1 / num_face_nodes[fno_ext[neu_rob_ind_all]]

# The columns will be 0:neu_ind.size
if neu_ind.size > 0:
# 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()

# For Dirichlet, the coefficients in the matrix should be duplicated the same way as
# the row indices, but with no increment

dir_val_x = sgn[dir_ind_single_all_x]
dir_val_y = sgn[dir_ind_single_all_y]

dir_val = np.zeros(dir_ind_sz)

dir_val[ind_is_bnd_dir_x] = dir_val_x
dir_val[ind_is_bnd_dir_y] = dir_val_y

if nd == 3:
dir_val_z = sgn[dir_ind_single_all_z]
dir_val[ind_is_bnd_dir_z] = dir_val_z

sgn_nd = np.tile(sgn, (nd, 1)).ravel("F")
dir_val = sgn_nd[dir_ind_all]
del sgn_nd
# Column numbering starts right after the last Neumann column. dir_val
# is ordered [u_x_1, u_y_1, u_x_2, u_y_2, ...], and dir_ind shuffles this
# ordering. The final matrix will first have the x-coponent of the displacement
# 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:
Expand All @@ -1817,7 +1733,6 @@ def create_bound_rhs_nd(bound, bound_exclusion, subcell_topology, g):

# 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), bnd_ind)),
shape=(num_bound, num_subfno * nd),
Expand Down
27 changes: 15 additions & 12 deletions src/porepy/params/rock.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@
moduli etc.
"""
from porepy.params import units
import porepy as pp


def poisson_from_lame(mu, lmbda):
Expand Down Expand Up @@ -71,10 +71,10 @@ class SandStone(UnitRock):
def __init__(self):

# Fairly permeable rock.
self.PERMEABILITY = 1 * units.DARCY
self.PERMEABILITY = 1 * pp.DARCY
self.POROSITY = 0.2
# Reported range for Young's modulus is 0.5-8.6
self.YOUNG_MODULUS = 5 * units.KILOGRAM / units.CENTI ** 2 * 1e5
self.YOUNG_MODULUS = 5 * pp.KILOGRAM / pp.CENTI ** 2 * 1e5
# Reported range for Poisson's ratio is 0.066-0.125
self.POISSON_RATIO = 0.1

Expand All @@ -94,10 +94,10 @@ class Shale(UnitRock):

def __init__(self):
# No source for permeability and porosity.
self.PERMEABILITY = 1e-5 * units.DARCY
self.PERMEABILITY = 1e-5 * pp.DARCY
self.POROSITY = 0.01
# Reported range for Young's modulus is 0.8-3.0
self.YOUNG_MODULUS = 1.5 * units.KILOGRAM / units.CENTI ** 2 * 1e5
self.YOUNG_MODULUS = 1.5 * pp.KILOGRAM / pp.CENTI ** 2 * 1e5
# Reported range for Poisson's ratio is 0.11-0.54 (the latter is strange)
self.POISSON_RATIO = 0.3

Expand All @@ -108,22 +108,25 @@ def __init__(self):

class Granite(UnitRock):
""" Generic values for granite.
Data partially from:
http://civilblog.org/2015/02/13/what-are-the-values-of-modulus-of-elasticity-poissons-ratio-for-different-rocks/
And:
https://www.jsg.utexas.edu/tyzhu/files/Some-Useful-Numbers.pdf
"""

def __init__(self):
# No source for permeability and porosity
self.PERMEABILITY = 1e-8 * units.DARCY
UnitRock.__init__(self)

self.PERMEABILITY = 1e-8 * pp.DARCY
# jsg range: 0.005 - .015
self.POROSITY = 0.01
# Reported range for Young's modulus is 2.6-7.0
self.YOUNG_MODULUS = 5 * units.KILOGRAM / units.CENTI ** 2 * 1e5
# Reported range for Young's modulus by jsg is 10-70GPa
self.YOUNG_MODULUS = 4 * pp.GIGA * pp.PASCAL
# Reported range for Poisson's ratio is 0.125-0.25
self.POISSON_RATIO = 0.2

# Reported density
self.DENSITY = 2700 * pp.KILOGRAM / pp.METER ** 3
self.LAMBDA, self.MU = lame_from_young_poisson(
self.YOUNG_MODULUS, self.POISSON_RATIO
)

0 comments on commit a892127

Please sign in to comment.