diff --git a/autotest/test_gwf_csub_distypes.py b/autotest/test_gwf_csub_distypes.py
new file mode 100644
index 00000000000..bc5267dfdc0
--- /dev/null
+++ b/autotest/test_gwf_csub_distypes.py
@@ -0,0 +1,448 @@
+import pathlib as pl
+
+import flopy
+import numpy as np
+import pytest
+from flopy.utils.gridgen import Gridgen
+
+from framework import TestFramework
+from simulation import TestSimulation
+
+ex = ["csub_dis", "csub_disv", "csub_disu", "csub_disu01", "csub_disu02"]
+# ex = ["csub_dis"]
+ex_dict = {name: None for name in ex}
+ex_dict["csub_disu01"] = 0
+ex_dict["csub_disu02"] = 2
+
+paktest = "csub"
+
+# temporal discretization
+nper = 2
+perlen = [1.0, 100.0]
+nstp = [1, 10]
+tsmult = [1.0] * nper
+tdis_rc = []
+for idx in range(nper):
+ tdis_rc.append((perlen[idx], nstp[idx], tsmult[idx]))
+
+# base spatial discretization
+nlay, nrow, ncol = 3, 9, 9
+refinement_level = 2
+nrow_refined, ncol_refined = nrow * refinement_level, ncol * refinement_level
+shape3d = (nlay, nrow, ncol)
+shape3d_refined = (nlay, nrow_refined, ncol_refined)
+shape2d = (nrow, ncol)
+shape2d_refined = (nrow_refined, ncol_refined)
+size3d = nlay * nrow * ncol
+size3d_refined = nlay * nrow_refined * ncol_refined
+size2d = nrow * ncol
+size2d_refined = nrow_refined * ncol_refined
+
+delr = delc = 1000.0
+top = 0.0
+bot = -100.0
+dz = (top - bot) / nlay
+botm = [top - k * dz for k in range(1, nlay + 1)]
+z_node = [z + 0.5 * dz for z in botm]
+
+delr_refined = delr / refinement_level
+delc_refined = delc / refinement_level
+
+hk = [1.0, 0.001, 1.0]
+sy = [0.25, 0.45, 0.25]
+ss = [5e-5, 5e-4, 5e-5]
+
+well_coordinates = (
+ np.array([(4.25, 4.25), (4.25, 4.75), (4.75, 4.75), (4.75, 4.25)]) * delr
+)
+wellq = -1000.0
+
+nouter, ninner = 100, 300
+dvclose, rclose, relax = 1e-6, 0.01, 0.97
+
+# subwt data
+cc = 0.25
+cr = 0.25
+void = 0.82
+theta = void / (1.0 + void)
+kv = 999.0
+sgm = 1.7
+sgs = 2.0
+
+beta = 0.0
+# beta = 4.65120000e-10
+gammaw = 9806.65000000
+
+
+def get_interbed(modelgrid):
+ grid_type = modelgrid.grid_type
+ ia = []
+ x0, x1, y0, y1 = modelgrid.extent
+ for k in range(1, nlay, 1):
+ cellid = modelgrid.intersect(x0 + 0.1, y1 - 0.1, z=z_node[k])
+ ia.append(get_node_number(modelgrid, cellid))
+
+ package_data = []
+ ifno = 0
+ ini_stress = 0.0
+
+ nodes = [node for node in range(ia[0], ia[1])]
+ if grid_type == "structured":
+ cellids = modelgrid.get_lrc(nodes)
+ elif grid_type == "vertex":
+ cellids = modelgrid.get_lni(nodes)
+ else:
+ cellids = [(node,) for node in nodes]
+
+ for cellid in cellids:
+ rnb = 1.0
+ vk = 999.0
+ package_data.append(
+ (
+ ifno,
+ cellid, # will need to be detuplaized with *cellid - does not work for dis
+ "nodelay",
+ ini_stress,
+ modelgrid.cell_thickness[cellid],
+ rnb,
+ ss[1],
+ ss[1] * 1000.0,
+ theta,
+ vk,
+ top,
+ )
+ )
+ ifno += 1
+ return package_data
+
+
+def build_dis(gwf):
+ return flopy.mf6.ModflowGwfdis(
+ gwf,
+ nlay=nlay,
+ nrow=nrow,
+ ncol=ncol,
+ delr=delr,
+ delc=delc,
+ top=top,
+ botm=botm,
+ )
+
+
+def get_gridgen_ws(ws):
+ gridgen_ws = ws / "gridgen"
+ gridgen_ws.mkdir(parents=True, exist_ok=True)
+ return gridgen_ws
+
+
+def build_temp_gwf(ws):
+ gridgen_ws = get_gridgen_ws(ws)
+ gridgen_sim = flopy.mf6.MFSimulation(
+ sim_name="gridgen", sim_ws=gridgen_ws, exe_name="mf6"
+ )
+ gridgen_gwf = flopy.mf6.ModflowGwf(gridgen_sim, modelname="gridgen")
+ return gridgen_gwf
+
+
+def build_disv(ws, gwf, gridgen):
+ temp_gwf = build_temp_gwf(ws)
+ dis = build_dis(temp_gwf)
+ g = Gridgen(
+ temp_gwf.modelgrid,
+ model_ws=get_gridgen_ws(ws),
+ exe_name=gridgen,
+ )
+ g.build()
+ gridprops = g.get_gridprops_disv()
+ return flopy.mf6.ModflowGwfdisv(gwf, **gridprops)
+
+
+def build_disu(ws, gwf, refinement_layer, gridgen):
+ temp_gwf = build_temp_gwf(ws)
+ dis = build_dis(temp_gwf)
+ g = Gridgen(
+ temp_gwf.modelgrid,
+ model_ws=get_gridgen_ws(ws),
+ exe_name=gridgen,
+ )
+ if refinement_layer is not None:
+ x0, x1, y0, y1 = temp_gwf.modelgrid.extent
+ polys = [[[(x0, y0), (x1, y0), (x1, y1), (x0, y1), (x0, y0)]]]
+ g.add_refinement_features(
+ polys,
+ "polygon",
+ 1,
+ layers=[refinement_layer],
+ )
+ g.build()
+ gridprops = g.get_gridprops_disu6()
+ return flopy.mf6.ModflowGwfdisu(gwf, **gridprops)
+
+
+def get_node_number(modelgrid, cellid):
+ if modelgrid.grid_type == "unstructured":
+ node = cellid
+ elif modelgrid.grid_type == "vertex":
+ node = modelgrid.ncpl * cellid[0] + cellid[1]
+ else:
+ node = (
+ modelgrid.nrow * modelgrid.ncol * cellid[0]
+ + modelgrid.ncol * cellid[1]
+ + cellid[2]
+ )
+ return node
+
+
+def build_3d_array(modelgrid, values, dtype=float):
+ if isinstance(values, dtype):
+ arr = np.full(modelgrid.nnodes, values, dtype=dtype)
+ else:
+ arr = np.zeros(modelgrid.nnodes, dtype=dtype)
+ ia = []
+ x0, x1, y0, y1 = modelgrid.extent
+ for k in range(nlay):
+ cellid = modelgrid.intersect(x0 + 0.1, y1 - 0.1, z=z_node[k])
+ ia.append(get_node_number(modelgrid, cellid))
+ ia.append(modelgrid.nnodes + 1)
+ for k in range(nlay):
+ arr[ia[k] : ia[k + 1]] = values[k]
+ return arr.reshape(modelgrid.shape)
+
+
+def build_well_data(modelgrid):
+ well_spd = []
+ for x, y in well_coordinates:
+ cellid = modelgrid.intersect(x, y, z=z_node[-1])
+ if isinstance(cellid, tuple):
+ well_spd.append((*cellid, wellq))
+ else:
+ well_spd.append((cellid, wellq))
+ return {1: well_spd}
+
+
+def build_model(idx, ws, gridgen):
+ return build_mf6(idx, ws, gridgen), None
+
+
+# build MODFLOW 6 files
+def build_mf6(idx, ws, gridgen):
+ name = ex[idx]
+ sim = flopy.mf6.MFSimulation(
+ sim_name=name,
+ version="mf6",
+ exe_name="mf6",
+ sim_ws=ws,
+ )
+ # create tdis package
+ tdis = flopy.mf6.ModflowTdis(
+ sim,
+ time_units="DAYS",
+ nper=nper,
+ perioddata=tdis_rc,
+ )
+
+ # create iterative model solution and register the gwf model with it
+ ims = flopy.mf6.ModflowIms(
+ sim,
+ print_option="SUMMARY",
+ outer_dvclose=dvclose,
+ outer_maximum=nouter,
+ under_relaxation="NONE",
+ inner_maximum=ninner,
+ inner_dvclose=dvclose,
+ rcloserecord=rclose,
+ linear_acceleration="CG",
+ scaling_method="NONE",
+ reordering_method="NONE",
+ relaxation_factor=relax,
+ )
+
+ # create gwf model
+ gwf = flopy.mf6.ModflowGwf(
+ sim,
+ modelname=name,
+ print_input=True,
+ save_flows=True,
+ )
+
+ if "disv" in name:
+ dis = build_disv(ws, gwf, gridgen)
+ elif "disu" in name:
+ dis = build_disu(ws, gwf, ex_dict[name], gridgen)
+ else:
+ dis = build_dis(gwf)
+
+ # initial conditions
+ ic = flopy.mf6.ModflowGwfic(
+ gwf,
+ strt=top,
+ )
+
+ k11 = build_3d_array(gwf.modelgrid, hk)
+ icelltype = build_3d_array(gwf.modelgrid, [1, 0, 0], dtype=int)
+
+ # node property flow
+ npf = flopy.mf6.ModflowGwfnpf(
+ gwf,
+ icelltype=icelltype,
+ k=k11,
+ k33=k11,
+ )
+
+ # storage
+ sy_arr = build_3d_array(gwf.modelgrid, sy)
+ sto = flopy.mf6.ModflowGwfsto(
+ gwf,
+ iconvert=icelltype,
+ ss=0.0,
+ sy=sy_arr,
+ transient={0: True},
+ )
+
+ # well
+ wel = flopy.mf6.ModflowGwfwel(
+ gwf,
+ stress_period_data=build_well_data(gwf.modelgrid),
+ save_flows=False,
+ )
+
+ # csub
+ cg_ske_cr = build_3d_array(gwf.modelgrid, ss)
+ packagedata = get_interbed(gwf.modelgrid)
+
+ csub = flopy.mf6.ModflowGwfcsub(
+ gwf,
+ zdisplacement_filerecord=f"{name}.csub.zdis.bin",
+ compaction_filerecord=f"{name}.csub.comp.bin",
+ ninterbeds=len(packagedata),
+ sgs=sgs,
+ sgm=sgm,
+ beta=beta,
+ gammaw=gammaw,
+ cg_ske_cr=cg_ske_cr,
+ cg_theta=theta,
+ packagedata=packagedata,
+ )
+ # orecarray = {}
+ # orecarray["csub_obs.csv"] = [
+ # ("wc01", "compaction-cell", (1, 5, 8)),
+ # ("wc02", "compaction-cell", (3, 6, 11)),
+ # ]
+ # csub_obs_package = csub.obs.initialize(
+ # filename=opth, digits=10, print_input=True, continuous=orecarray
+ # )
+
+ # output control
+ oc = flopy.mf6.ModflowGwfoc(
+ gwf,
+ budget_filerecord=f"{name}.cbc",
+ head_filerecord=f"{name}.hds",
+ headprintrecord=[("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL")],
+ saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")],
+ printrecord=[("HEAD", "ALL"), ("BUDGET", "ALL")],
+ )
+ return sim
+
+
+def eval_zdis(sim):
+ print("evaluating z-displacement...")
+
+ name = ex[sim.idxsim]
+ ws = pl.Path(sim.simpath)
+ sim = flopy.mf6.MFSimulation.load(sim_name=name, sim_ws=ws)
+ gwf = sim.get_model()
+ x0, x1, y0, y1 = gwf.modelgrid.extent
+
+ comp_obj = flopy.utils.HeadFile(
+ ws / f"{name}.csub.comp.bin",
+ text="CSUB-COMPACTION",
+ precision="double",
+ )
+ zdis_obj = flopy.utils.HeadFile(
+ ws / f"{name}.csub.zdis.bin",
+ text="CSUB-ZDISPLACE",
+ precision="double",
+ )
+
+ layer_refinement = ex_dict[name]
+
+ # create reusable mapping dictionary so it can be used for all time step
+ # with refined disu grids - which do not have naturally ordered node
+ # numbers in refined layers
+ map_dict = {}
+ if layer_refinement is not None:
+ z = z_node[layer_refinement]
+ icnt = 0
+ for i in range(nrow_refined):
+ y = y1 - delc_refined * (i + 0.5)
+ for j in range(ncol_refined):
+ x = x0 + delr_refined * (j + 0.5)
+ node = gwf.modelgrid.intersect(x, y, z=z)
+ map_dict[icnt] = {"node": node, "cellid": (i, j)}
+ icnt += 1
+
+ for totim in comp_obj.get_times():
+ if layer_refinement is None:
+ comp = comp_obj.get_data(totim=totim).flatten().reshape(shape3d)
+ zdis = zdis_obj.get_data(totim=totim).flatten().reshape(shape3d)
+ else:
+ comp1d = comp_obj.get_data(totim=totim).squeeze()
+ zdis1d = zdis_obj.get_data(totim=totim).squeeze()
+ ia = [0]
+ for k in range(nlay):
+ if k == layer_refinement:
+ ia.append(ia[k] + size2d_refined)
+ else:
+ ia.append(ia[k] + size2d)
+
+ comp = np.zeros(shape3d, dtype=float)
+ zdis = np.zeros(shape3d, dtype=float)
+ for k in range(nlay):
+ ia0 = ia[k]
+ ia1 = ia[k + 1]
+ comp_slice = comp1d[ia0:ia1].copy()
+ zdis_slice = zdis1d[ia0:ia1].copy()
+ if k == layer_refinement:
+ comp_temp = np.zeros(shape2d_refined, dtype=float)
+ zdis_temp = np.zeros(shape2d_refined, dtype=float)
+ for key, value in map_dict.items():
+ comp_temp[value["cellid"]] = comp1d[value["node"]]
+ zdis_temp[value["cellid"]] = zdis1d[value["node"]]
+ comp[k] = comp_temp.reshape(
+ nrow_refined // 2, 2, ncol_refined // 2, 2
+ ).mean(axis=(1, -1))
+ zdis[k] = zdis_temp.reshape(
+ nrow_refined // 2, 2, ncol_refined // 2, 2
+ ).mean(axis=(1, -1))
+ else:
+ comp[k] = comp_slice.reshape(shape2d)
+ zdis[k] = zdis_slice.reshape(shape2d)
+
+ comp = comp.sum(axis=0)
+ zdis = zdis[0]
+ assert np.allclose(comp, zdis), (
+ "sum of compaction is not equal to the "
+ + f"z-displacement at time {totim}"
+ )
+
+ return
+
+
+@pytest.mark.parametrize(
+ "idx, name",
+ list(enumerate(ex)),
+)
+def test_mf6model(idx, name, function_tmpdir, targets):
+ ws = function_tmpdir
+ test = TestFramework()
+ test.build(lambda i, w: build_model(i, w, targets.gridgen), idx, ws)
+ test.run(
+ TestSimulation(
+ name=name,
+ exe_dict=targets,
+ exfunc=eval_zdis,
+ cmp_verbose=False,
+ idxsim=idx,
+ ),
+ ws,
+ )
diff --git a/environment.yml b/environment.yml
index c989e3bc197..61176675aae 100644
--- a/environment.yml
+++ b/environment.yml
@@ -13,6 +13,8 @@ dependencies:
- meson>=1.1.0
- ninja
- numpy
+ - pyshp
+ - shapely
- pip
- pip:
- git+https://github.com/modflowpy/flopy.git
diff --git a/msvs/mf6core.vfproj b/msvs/mf6core.vfproj
index ee13f7fa521..8270d990e13 100644
--- a/msvs/mf6core.vfproj
+++ b/msvs/mf6core.vfproj
@@ -2,50 +2,44 @@
-
-
+
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
@@ -78,13 +72,13 @@
-
+
-
+
-
+
-
+
@@ -181,7 +175,8 @@
-
+
+
@@ -232,13 +227,13 @@
-
+
-
+
-
+
-
+
@@ -247,13 +242,13 @@
-
+
-
+
-
+
-
+
@@ -351,23 +346,23 @@
-
+
-
+
-
+
-
+
-
+
-
+
-
+
-
+
@@ -397,15 +392,13 @@
-
+
-
+
-
+
-
+
-
-
-
-
+
+
diff --git a/src/Model/GroundWaterFlow/gwf3csub8.f90 b/src/Model/GroundWaterFlow/gwf3csub8.f90
index 8b1c8621e26..d7de64d3635 100644
--- a/src/Model/GroundWaterFlow/gwf3csub8.f90
+++ b/src/Model/GroundWaterFlow/gwf3csub8.f90
@@ -247,7 +247,7 @@ module GwfCsubModule
procedure, private :: csub_read_packagedata
!
! -- helper methods
- procedure, private :: csub_calc_void
+ procedure, private :: csub_calc_void_ratio
procedure, private :: csub_calc_theta
procedure, private :: csub_calc_znode
procedure, private :: csub_calc_adjes
@@ -3691,10 +3691,14 @@ subroutine csub_ot_dv(this, idvfl, idvprint)
integer(I4B) :: nodem
integer(I4B) :: nodeu
integer(I4B) :: i
+ integer(I4B) :: ii
+ integer(I4B) :: idx_conn
integer(I4B) :: k
integer(I4B) :: ncpl
integer(I4B) :: nlay
+ integer(I4B) :: ihc
real(DP) :: dinact
+ real(DP) :: va_scale
! -- formats
character(len=*), parameter :: fmtnconv = &
"(/4x, 'DELAY INTERBED CELL HEADS IN ', i0, ' INTERBEDS IN', &
@@ -3752,7 +3756,26 @@ subroutine csub_ot_dv(this, idvfl, idvprint)
!
! -- disu
if (this%dis%ndim == 1) then
- ! TO DO -
+ do node = this%dis%nodes, 1, -1
+ do ii = this%dis%con%ia(node) + 1, this%dis%con%ia(node + 1) - 1
+ !
+ ! -- Set the m cell number
+ nodem = this%dis%con%ja(ii)
+ idx_conn = this%dis%con%jas(ii)
+ !
+ ! -- vertical connection
+ ihc = this%dis%con%ihc(idx_conn)
+ if (ihc == 0) then
+ !
+ ! -- node has an underlying cell
+ if (node < nodem) then
+ va_scale = this%dis%get_area_factor(node, idx_conn)
+ this%buffusr(node) = this%buffusr(node) + &
+ va_scale * this%buffusr(nodem)
+ end if
+ end if
+ end do
+ end do
! -- disv or dis
else
nlay = this%dis%nodesuser / ncpl
@@ -3929,7 +3952,7 @@ subroutine csub_cg_calc_stress(this, nodes, hnew)
integer(I4B) :: ii
integer(I4B) :: nn
integer(I4B) :: m
- integer(I4B) :: iis
+ integer(I4B) :: idx_conn
real(DP) :: gs
real(DP) :: top
real(DP) :: bot
@@ -3938,11 +3961,8 @@ subroutine csub_cg_calc_stress(this, nodes, hnew)
real(DP) :: hcell
real(DP) :: hbar
real(DP) :: gs_conn
- real(DP) :: area_node
- real(DP) :: area_conn
real(DP) :: es
real(DP) :: phead
- real(DP) :: hwva
real(DP) :: sadd
!
! -- calculate geostatic stress if necessary
@@ -3986,9 +4006,6 @@ subroutine csub_cg_calc_stress(this, nodes, hnew)
!
! -- calculate geostatic stress above cell
do node = 1, this%dis%nodes
- !
- ! -- area of cell
- area_node = this%dis%get_area(node)
!
! -- geostatic stress of cell
gs = this%cg_gs(node)
@@ -3999,10 +4016,10 @@ subroutine csub_cg_calc_stress(this, nodes, hnew)
!
! -- Set the m cell number
m = this%dis%con%ja(ii)
- iis = this%dis%con%jas(ii)
+ idx_conn = this%dis%con%jas(ii)
!
! -- vertical connection
- if (this%dis%con%ihc(iis) == 0) then
+ if (this%dis%con%ihc(idx_conn) == 0) then
!
! -- node has an overlying cell
if (m < node) then
@@ -4012,15 +4029,11 @@ subroutine csub_cg_calc_stress(this, nodes, hnew)
gs = gs + this%cg_gs(m)
!
! -- disu discretization
- ! *** this needs to be checked ***
else
- area_conn = this%dis%get_area(m)
- hwva = this%dis%con%hwva(iis)
- va_scale = this%dis%con%hwva(iis) / this%dis%get_area(m)
+ va_scale = this%dis%get_area_factor(node, idx_conn)
gs_conn = this%cg_gs(m)
gs = gs + (gs_conn * va_scale)
end if
-
end if
end if
end do
@@ -4356,7 +4369,7 @@ subroutine csub_set_initial_state(this, nodes, hnew)
real(DP) :: fact
real(DP) :: top
real(DP) :: bot
- real(DP) :: void
+ real(DP) :: void_ratio
real(DP) :: es
real(DP) :: znode
real(DP) :: hcell
@@ -4450,7 +4463,7 @@ subroutine csub_set_initial_state(this, nodes, hnew)
! -- convert specific storage values since they are simulated to
! be a function of the average effective stress
else
- void = this%csub_calc_void(this%cg_theta(node))
+ void_ratio = this%csub_calc_void_ratio(this%cg_theta(node))
es = this%cg_es(node)
hcell = hnew(node)
!
@@ -4460,7 +4473,7 @@ subroutine csub_set_initial_state(this, nodes, hnew)
! -- calculate znode and factor
znode = this%csub_calc_znode(top, bot, hbar)
fact = this%csub_calc_adjes(node, es, bot, znode)
- fact = fact * (DONE + void)
+ fact = fact * (DONE + void_ratio)
end if
!
! -- user-specified compression indices - multiply by dlog10es
@@ -4496,7 +4509,7 @@ subroutine csub_set_initial_state(this, nodes, hnew)
! -- convert specific storage values since they are simulated to
! be a function of the average effective stress
else
- void = this%csub_calc_void(this%theta(ib))
+ void_ratio = this%csub_calc_void_ratio(this%theta(ib))
es = this%cg_es(node)
hcell = hnew(node)
!
@@ -4506,7 +4519,7 @@ subroutine csub_set_initial_state(this, nodes, hnew)
! -- calculate zone and factor
znode = this%csub_calc_znode(top, bot, hbar)
fact = this%csub_calc_adjes(node, es, bot, znode)
- fact = fact * (DONE + void)
+ fact = fact * (DONE + void_ratio)
end if
!
! -- user-specified compression indices - multiply by dlog10es
@@ -5426,18 +5439,18 @@ end subroutine csub_nodelay_wcomp_fn
!!
!! @return void void ratio
!<
- function csub_calc_void(this, theta) result(void)
+ function csub_calc_void_ratio(this, theta) result(void_ratio)
! -- dummy variables
class(GwfCsubType), intent(inout) :: this
real(DP), intent(in) :: theta !< porosity
! -- local variables
- real(DP) :: void
+ real(DP) :: void_ratio
! -- calculate void ratio
- void = theta / (DONE - theta)
+ void_ratio = theta / (DONE - theta)
!
! -- return
return
- end function csub_calc_void
+ end function csub_calc_void_ratio
!> @brief Calculate the porosity
!!
@@ -5445,15 +5458,15 @@ end function csub_calc_void
!!
!! @return theta porosity
!<
- function csub_calc_theta(this, void) result(theta)
+ function csub_calc_theta(this, void_ratio) result(theta)
! -- dummy variables
class(GwfCsubType), intent(inout) :: this
- real(DP), intent(in) :: void
+ real(DP), intent(in) :: void_ratio
! -- local variables
real(DP) :: theta
!
! -- calculate theta
- theta = void / (DONE + void)
+ theta = void_ratio / (DONE + void_ratio)
!
! -- return
return
@@ -5678,10 +5691,10 @@ subroutine csub_calc_sfacts(this, node, bot, znode, theta, es, es0, fact)
real(DP), intent(in) :: theta !< porosity
real(DP), intent(in) :: es !< current effective stress
real(DP), intent(in) :: es0 !< previous effective stress
- real(DP), intent(inout) :: fact !< skeletal storage coefficient factor (1/((1+void)*bar(es)))
+ real(DP), intent(inout) :: fact !< skeletal storage coefficient factor (1/((1+void_ratio)*bar(es)))
! -- local variables
real(DP) :: esv
- real(DP) :: void
+ real(DP) :: void_ratio
real(DP) :: denom
!
! -- initialize variables
@@ -5693,9 +5706,9 @@ subroutine csub_calc_sfacts(this, node, bot, znode, theta, es, es0, fact)
end if
!
! -- calculate storage factors for the effective stress case
- void = this%csub_calc_void(theta)
+ void_ratio = this%csub_calc_void_ratio(theta)
denom = this%csub_calc_adjes(node, esv, bot, znode)
- denom = denom * (DONE + void)
+ denom = denom * (DONE + void_ratio)
if (denom /= DZERO) then
fact = DONE / denom
end if
@@ -5720,18 +5733,18 @@ subroutine csub_adj_matprop(this, comp, thick, theta)
real(DP), intent(inout) :: theta !< porosity
! -- local variables
real(DP) :: strain
- real(DP) :: void
+ real(DP) :: void_ratio
!
! -- initialize variables
strain = DZERO
- void = this%csub_calc_void(theta)
+ void_ratio = this%csub_calc_void_ratio(theta)
!
! -- calculate strain
if (thick > DZERO) strain = -comp / thick
!
! -- update void ratio, theta, and thickness
- void = void + strain * (DONE + void)
- theta = this%csub_calc_theta(void)
+ void_ratio = void_ratio + strain * (DONE + void_ratio)
+ theta = this%csub_calc_theta(void_ratio)
thick = thick - comp
!
! -- return
diff --git a/src/Model/ModelUtilities/DiscretizationBase.f90 b/src/Model/ModelUtilities/DiscretizationBase.f90
index 6c43a150d01..03c807dc1fd 100644
--- a/src/Model/ModelUtilities/DiscretizationBase.f90
+++ b/src/Model/ModelUtilities/DiscretizationBase.f90
@@ -107,6 +107,7 @@ module BaseDisModule
procedure, public :: nlarray_to_nodelist
procedure, public :: highest_active
procedure, public :: get_area
+ procedure, public :: get_area_factor
end type DisBaseType
@@ -1505,4 +1506,37 @@ function get_area(this, node) result(area)
return
end function get_area
+ !> @ brief Calculate the area factor for the cell connection
+ !!
+ !! Function calculates the area factor for the cell connection. The sum of
+ !! all area factors for all cell connections to overlying or underlying
+ !! cells cells will be 1.
+ !!
+ !! TODO: confirm that this works for cells that are only partially covered
+ !! by overlying or underlying cells.
+ !!
+ !<
+ function get_area_factor(this, node, idx_conn) result(area_factor)
+ ! -- return
+ real(DP) :: area_factor !< connection cell area factor
+ ! -- dummy
+ class(DisBaseType) :: this
+ integer(I4B), intent(in) :: node !< cell node number
+ integer(I4B), intent(in) :: idx_conn !< connection index
+ ! -- local
+ real(DP) :: area_node
+ real(DP) :: area_conn
+ ! ------------------------------------------------------------------------------
+ !
+ ! -- calculate the cell area fraction
+ area_node = this%area(node)
+ area_conn = this%con%hwva(idx_conn)
+ !
+ ! -- return the cell area factor
+ area_factor = area_conn / area_node
+ !
+ ! -- return
+ return
+ end function get_area_factor
+
end module BaseDisModule