From af26abf950c993301c59ea2dfec0a81802321361 Mon Sep 17 00:00:00 2001 From: langevin-usgs Date: Fri, 23 Aug 2024 17:19:02 -0500 Subject: [PATCH] fix(swf-dis): fixes for swf dis functionality (#2000) * fix(swf-dis): fixes for swf dis functionality * rename botm to bottom, allow idomain for disv1d, add tests * format * ruff check . --fix * implement fractional cell distance calculation to offset disv1d nodes anywhere along the cell * add more rigorous test of manning flow calculation * for the disv1d channel grid, rename cell2d to cell1d * remove length and input parameter in disv1d griddata; it can be calculated from vertices --- autotest/test_chf_dfw.py | 9 +- autotest/test_chf_dfw_beg2022.py | 9 +- autotest/test_chf_dfw_bowl.py | 9 +- autotest/test_chf_dfw_loop.py | 20 +- autotest/test_chf_dfw_swrt2.py | 10 +- autotest/test_chf_dfw_swrt2b.py | 9 +- autotest/test_chf_dis.py | 231 ++++++++++++++++ autotest/test_chf_dis_fdc.py | 275 +++++++++++++++++++ autotest/test_chf_gwf.py | 9 +- autotest/test_olf_dfw_swrt2dis.py | 8 +- autotest/test_olf_dis.py | 327 +++++++++++++++++++++++ autotest/test_olf_gwf.py | 2 +- autotest/test_olf_vcatch.py | 8 +- doc/mf6io/mf6ivar/dfn/chf-dis2d.dfn | 2 +- doc/mf6io/mf6ivar/dfn/chf-disv1d.dfn | 35 +-- doc/mf6io/mf6ivar/dfn/olf-dis2d.dfn | 2 +- doc/mf6io/mf6ivar/dfn/swf-dis2d.dfn | 2 +- doc/mf6io/mf6ivar/dfn/swf-disv1d.dfn | 35 +-- src/Idm/chf-dis2didm.f90 | 10 +- src/Idm/chf-disv1didm.f90 | 56 ++-- src/Idm/olf-dis2didm.f90 | 10 +- src/Idm/swf-dis2didm.f90 | 10 +- src/Model/Discretization/Dis2d.f90 | 25 +- src/Model/Discretization/Disv1d.f90 | 181 +++++++------ src/Model/Discretization/Disv2d.f90 | 114 ++++---- src/Model/ModelUtilities/Connections.f90 | 240 +++++++++++------ 26 files changed, 1257 insertions(+), 391 deletions(-) create mode 100644 autotest/test_chf_dis.py create mode 100644 autotest/test_chf_dis_fdc.py create mode 100644 autotest/test_olf_dis.py diff --git a/autotest/test_chf_dfw.py b/autotest/test_chf_dfw.py index 6dd2bd6aeb6..a4cc1459d04 100644 --- a/autotest/test_chf_dfw.py +++ b/autotest/test_chf_dfw.py @@ -57,22 +57,21 @@ def build_models(idx, test): total_length = dx * nreach vertices = [] vertices = [[j, j * dx, 0.0] for j in range(nreach + 1)] - cell2d = [] + cell1d = [] for j in range(nreach): - cell2d.append([j, 0.5, 2, j, j + 1]) - nodes = len(cell2d) + cell1d.append([j, 0.5, 2, j, j + 1]) + nodes = len(cell1d) nvert = len(vertices) disv1d = flopy.mf6.ModflowChfdisv1D( chf, nodes=nodes, nvert=nvert, - length=dx, width=50.0, bottom=0.0, idomain=1, vertices=vertices, - cell2d=cell2d, + cell1d=cell1d, ) dfw = flopy.mf6.ModflowChfdfw( diff --git a/autotest/test_chf_dfw_beg2022.py b/autotest/test_chf_dfw_beg2022.py index c1caf5acacc..d122ee2b76c 100644 --- a/autotest/test_chf_dfw_beg2022.py +++ b/autotest/test_chf_dfw_beg2022.py @@ -66,10 +66,10 @@ def build_models(idx, test): nreach = int(total_length / dx) vertices = [] vertices = [[j, j * dx, 0.0] for j in range(nreach + 1)] - cell2d = [] + cell1d = [] for j in range(nreach): - cell2d.append([j, 0.5, 2, j, j + 1]) - nodes = len(cell2d) + cell1d.append([j, 0.5, 2, j, j + 1]) + nodes = len(cell1d) nvert = len(vertices) slope = 1.0 / 10000.0 @@ -80,12 +80,11 @@ def build_models(idx, test): chf, nodes=nodes, nvert=nvert, - length=dx, width=40.0, bottom=z, idomain=1, vertices=vertices, - cell2d=cell2d, + cell1d=cell1d, ) dfw = flopy.mf6.ModflowChfdfw( diff --git a/autotest/test_chf_dfw_bowl.py b/autotest/test_chf_dfw_bowl.py index 073affc38d3..b7f6422ce56 100644 --- a/autotest/test_chf_dfw_bowl.py +++ b/autotest/test_chf_dfw_bowl.py @@ -91,22 +91,21 @@ def build_models(idx, test): vertices = [] vertices = [[j, j * dx, 0.0] for j in range(nreach + 1)] - cell2d = [] + cell1d = [] for j in range(nreach): - cell2d.append([j, 0.5, 2, j, j + 1]) - nodes = len(cell2d) + cell1d.append([j, 0.5, 2, j, j + 1]) + nodes = len(cell1d) nvert = len(vertices) disv1d = flopy.mf6.ModflowChfdisv1D( chf, nodes=nodes, nvert=nvert, - length=dx, width=1.0, bottom=reach_bottom, idomain=1, vertices=vertices, - cell2d=cell2d, + cell1d=cell1d, ) dfw = flopy.mf6.ModflowChfdfw( diff --git a/autotest/test_chf_dfw_loop.py b/autotest/test_chf_dfw_loop.py index 60b1eeea978..426887823e7 100644 --- a/autotest/test_chf_dfw_loop.py +++ b/autotest/test_chf_dfw_loop.py @@ -82,7 +82,7 @@ def build_models(idx, test): [19, 3000.0, 2500.0], ] - cell2d = [ + cell1d = [ [0, 0.5, 2, 0, 1], [1, 0.5, 2, 1, 2], [2, 0.5, 2, 2, 3], @@ -218,19 +218,18 @@ def build_models(idx, test): (608400, 0, 0), ] - nodes = len(cell2d) + nodes = len(cell1d) nvert = len(vertices) disv1d = flopy.mf6.ModflowChfdisv1D( chf, nodes=nodes, nvert=nvert, - length=reach_length, width=1.0, bottom=reach_bottom, idomain=1, vertices=vertices, - cell2d=cell2d, + cell1d=cell1d, ) stage0 = np.array(14 * [3] + 4 * [2]) @@ -413,6 +412,19 @@ def make_plot(test): fname = ws / f"{name}.obs.2.png" plt.savefig(fname) + fig = plt.figure(figsize=(6, 6)) + ax = fig.add_subplot(1, 1, 1) + ax.set_xlim(-100, 6100) + ax.set_ylim(-100, 6100) + pmv = flopy.plot.PlotMapView(model=chf, ax=ax) + pmv.plot_grid() + vertices = chf.disv1d.vertices.get_data() + ax.plot(vertices["xv"], vertices["yv"], "bo") + for iv, x, y in vertices: + ax.text(x, y, f"{iv + 1}") + fname = ws / "grid.png" + plt.savefig(fname) + return diff --git a/autotest/test_chf_dfw_swrt2.py b/autotest/test_chf_dfw_swrt2.py index 8e3c99c3bc0..bd3883b4f97 100644 --- a/autotest/test_chf_dfw_swrt2.py +++ b/autotest/test_chf_dfw_swrt2.py @@ -76,10 +76,10 @@ def build_models(idx, test): vertices = [] vertices = [[j, j * dx, 0.0] for j in range(nreach + 1)] - cell2d = [] + cell1d = [] for j in range(nreach): - cell2d.append([j, 0.5, 2, j, j + 1]) - nodes = len(cell2d) + cell1d.append([j, 0.5, 2, j, j + 1]) + nodes = len(cell1d) nvert = len(vertices) reach_bottom = np.linspace(1.05, 0.05, nreach) @@ -89,12 +89,11 @@ def build_models(idx, test): export_array_ascii=True, nodes=nodes, nvert=nvert, - length=dx, width=dx, bottom=reach_bottom, idomain=1, vertices=vertices, - cell2d=cell2d, + cell1d=cell1d, ) dfw = flopy.mf6.ModflowChfdfw( @@ -226,7 +225,6 @@ def check_output(idx, test): # ensure export array is working properly flist = [ - "disv1d.length", "disv1d.width", "disv1d.bottom", "disv1d.idomain", diff --git a/autotest/test_chf_dfw_swrt2b.py b/autotest/test_chf_dfw_swrt2b.py index 3907e9bb544..95e556f0185 100644 --- a/autotest/test_chf_dfw_swrt2b.py +++ b/autotest/test_chf_dfw_swrt2b.py @@ -86,10 +86,10 @@ def build_models(idx, test): vertices = [] vertices = [[j, j * dx, 0.0] for j in range(nreach + 1)] - cell2d = [] + cell1d = [] for j in range(nreach): - cell2d.append([j, 0.5, 2, j, j + 1]) - nodes = len(cell2d) + cell1d.append([j, 0.5, 2, j, j + 1]) + nodes = len(cell1d) nvert = len(vertices) reach_bottom = np.linspace(1.05, 0.05, nreach) @@ -98,12 +98,11 @@ def build_models(idx, test): chf, nodes=nodes, nvert=nvert, - length=dx, width=dx, bottom=reach_bottom, idomain=1, vertices=vertices, - cell2d=cell2d, + cell1d=cell1d, ) dfw = flopy.mf6.ModflowChfdfw( diff --git a/autotest/test_chf_dis.py b/autotest/test_chf_dis.py new file mode 100644 index 00000000000..e70d784a4b3 --- /dev/null +++ b/autotest/test_chf_dis.py @@ -0,0 +1,231 @@ +""" + +Test the disv1d to ensure that it supports IDOMAIN and writes a valid binary +grid file. + +""" + +import flopy +import numpy as np +import pytest +from framework import TestFramework + +cases = [ + "chf-disv1d", +] + +# grid size +dx = 1000.0 +nreach = 9 +total_length = dx * nreach +vertices = [] +vertices = [[j, j * dx, 0.0] for j in range(nreach + 1)] +cell1d = [] +for j in range(nreach): + cell1d.append([j, 0.5, 2, j, j + 1]) +nodes = len(cell1d) +nvert = len(vertices) +xorigin = 100.0 +yorigin = 200.0 +angrot = 25.0 + +idomain = np.ones((nodes,), dtype=int) +idomain[3:6] = 0 + + +def build_models(idx, test): + perlen = [1] # 1 second + nstp = [1] + tsmult = [1.0] + nper = len(perlen) + + tdis_rc = [] + for i in range(nper): + tdis_rc.append((perlen[i], nstp[i], tsmult[i])) + + name = "chf" + + # build MODFLOW 6 files + ws = test.workspace + sim = flopy.mf6.MFSimulation( + sim_name=f"{name}_sim", version="mf6", exe_name="mf6", sim_ws=ws + ) + # create tdis package + tdis = flopy.mf6.ModflowTdis( + sim, time_units="SECONDS", nper=nper, perioddata=tdis_rc + ) + + nouter, ninner = 100, 50 + hclose, rclose, relax = 1e-6, 1e-8, 1.0 + ims = flopy.mf6.ModflowIms( + sim, + print_option="SUMMARY", + outer_dvclose=hclose, + outer_maximum=nouter, + under_relaxation="DBD", + under_relaxation_theta=0.9, + under_relaxation_kappa=0.0001, + under_relaxation_gamma=0.0, + inner_maximum=ninner, + inner_dvclose=hclose, + linear_acceleration="BICGSTAB", + scaling_method="NONE", + reordering_method="NONE", + relaxation_factor=relax, + # backtracking_number=5, + # backtracking_tolerance=1.0, + # backtracking_reduction_factor=0.3, + # backtracking_residual_limit=100.0, + ) + + add_chf_model_disv1d(sim) + + return sim, None + + +def add_chf_model_disv1d(sim): + name = "channel" + chf = flopy.mf6.ModflowChf(sim, modelname=name, save_flows=True) + + disv1d = flopy.mf6.ModflowChfdisv1D( + chf, + nodes=nodes, + nvert=nvert, + width=50.0, + bottom=0.0, + idomain=idomain, + vertices=vertices, + cell1d=cell1d, + xorigin=xorigin, + yorigin=yorigin, + angrot=angrot, + ) + + dfw = flopy.mf6.ModflowChfdfw( + chf, + print_flows=True, + save_flows=True, + length_conversion=1.0, + time_conversion=86400.0, + manningsn=0.035, + idcxs=None, + ) + + ic = flopy.mf6.ModflowChfic(chf, strt=1.0) + + # output control + oc = flopy.mf6.ModflowChfoc( + chf, + budget_filerecord=f"{name}.bud", + stage_filerecord=f"{name}.stage", + saverecord=[ + ("STAGE", "ALL"), + ("BUDGET", "ALL"), + ], + printrecord=[ + ("STAGE", "LAST"), + ("BUDGET", "ALL"), + ], + ) + + chd = flopy.mf6.ModflowChfchd( + chf, + maxbound=1, + print_input=True, + print_flows=True, + stress_period_data=[(0, 0.5), (nodes - 1, 1.0)], + ) + + return + + +def make_plot(test, mfsim, stage, idx): + print("making plots...") + import matplotlib.pyplot as plt + + chf = mfsim.chf[0] + pmv = flopy.plot.PlotMapView(model=chf) + pmv.plot_array(stage, masked_values=[3e30]) + pmv.plot_grid() + + fname = test.workspace / "results.png" + plt.savefig(fname) + + +def check_grb_disv1d(fpth): + grb = flopy.mf6.utils.MfGrdFile(fpth) + assert grb.grid_type == "DISV1D", "grb grid type not DISV1D" + assert grb.ncells == nodes, "grb ncells is incorrect" + assert grb._datadict["NCELLS"] == nodes, "grb nodes is incorrect" + assert grb.verts.shape == ( + nodes + 1, + 2, + ), "vertices shape is incorrect" + assert grb.nja == 14, "nja in grb file is not 14" + assert grb.xorigin == xorigin, "xorigin in grb file is not correct" + assert grb.yorigin == yorigin, "yorigin in grb file is not correct" + assert grb.angrot == angrot, "angrot in grb file is not correct" + assert np.allclose( + grb.bot.reshape((nodes,)), np.zeros((nodes,)) + ), "grb botm not correct" + cellx = np.linspace(dx / 2, nreach * dx - dx / 2, nreach) + celly = np.zeros(nreach) + assert np.allclose( + grb._datadict["CELLX"], cellx.flatten() + ), "cellx is not right" + assert np.allclose( + grb._datadict["CELLY"], celly.flatten() + ), "celly is not right" + assert ( + grb._datadict["IAVERT"].shape[0] == nodes + 1 + ), "iavert size not right" + assert ( + grb._datadict["IAVERT"][-1] - 1 == grb._datadict["JAVERT"].shape[0] + ), "javert size not right" + assert ( + grb.ia.shape[0] == grb.ncells + 1 + ), "ia in grb file is not correct size" + assert grb.ja.shape[0] == grb.nja, "ja in grb file is not corect size" + assert np.allclose( + grb.idomain.reshape((nodes,)), idomain.reshape((nodes,)) + ), "grb idomain not correct" + + +def check_output(idx, test): + print(f"evaluating model for case {idx}...") + + modelname = "channel" + ws = test.workspace + mfsim = flopy.mf6.MFSimulation.load(sim_ws=ws) + + # read the binary grid file + fpth = test.workspace / f"{modelname}.disv1d.grb" + check_grb_disv1d(fpth) + + # read binary stage file + fpth = test.workspace / f"{modelname}.stage" + sobj = flopy.utils.HeadFile(fpth, precision="double", text="STAGE") + stage = sobj.get_data().reshape((nodes,)) + assert np.allclose( + stage[idomain == 0], 3.0e30 + ), "stage should have nodata values where idomain is zero" + assert stage[idomain == 1].max() == 1.0, "maximum stage should be 1.0" + assert stage[idomain == 1].min() == 0.5, "minimum stage should be 0.5" + + makeplot = False + if makeplot: + # PlotMapView not working yet for disv1d + make_plot(test, mfsim, stage.flatten(), idx) + + +@pytest.mark.developmode +@pytest.mark.parametrize("idx, name", enumerate(cases)) +def test_mf6model(idx, name, function_tmpdir, targets): + test = TestFramework( + name=name, + workspace=function_tmpdir, + build=lambda t: build_models(idx, t), + check=lambda t: check_output(idx, t), + targets=targets, + ) + test.run() diff --git a/autotest/test_chf_dis_fdc.py b/autotest/test_chf_dis_fdc.py new file mode 100644 index 00000000000..6bc076def68 --- /dev/null +++ b/autotest/test_chf_dis_fdc.py @@ -0,0 +1,275 @@ +""" + +Test the fdc input parameter for disv1d to shift the node location along the +cell. Use a 2 cell model and assign different roughness values for cell 1 +and cell 2. Fix the head in both cells and make a hand calculation of the flow. +Compare the hand calculation of flow with the simulated flow. + +""" + +import flopy +import numpy as np +import pytest +from framework import TestFramework + +cases = [ + "chf-dis-fdc", +] + +# grid size +dx = 1000.0 +nreach = 2 +total_length = dx * nreach +vertices = [] +vertices = [[j, j * dx, 0.0] for j in range(nreach + 1)] +cell1d = [] +fdc = [0.0, 1.0] +for j in range(nreach): + cell1d.append([j, fdc[j], 2, j, j + 1]) +nodes = len(cell1d) +nvert = len(vertices) +xorigin = 100.0 +yorigin = 200.0 +angrot = 25.0 + +idomain = np.ones((nodes,), dtype=int) +idomain[3:6] = 0 +h0 = 0.5 +h1 = 1.0 +rough = [0.035, 0.35] + + +def build_models(idx, test): + perlen = [1] # 1 second + nstp = [1] + tsmult = [1.0] + nper = len(perlen) + + tdis_rc = [] + for i in range(nper): + tdis_rc.append((perlen[i], nstp[i], tsmult[i])) + + name = "chf" + + # build MODFLOW 6 files + ws = test.workspace + sim = flopy.mf6.MFSimulation( + sim_name=f"{name}_sim", version="mf6", exe_name="mf6", sim_ws=ws + ) + # create tdis package + tdis = flopy.mf6.ModflowTdis( + sim, time_units="SECONDS", nper=nper, perioddata=tdis_rc + ) + + nouter, ninner = 100, 50 + hclose, rclose, relax = 1e-6, 1e-8, 1.0 + ims = flopy.mf6.ModflowIms( + sim, + print_option="SUMMARY", + outer_dvclose=hclose, + outer_maximum=nouter, + under_relaxation="DBD", + under_relaxation_theta=0.9, + under_relaxation_kappa=0.0001, + under_relaxation_gamma=0.0, + inner_maximum=ninner, + inner_dvclose=hclose, + linear_acceleration="BICGSTAB", + scaling_method="NONE", + reordering_method="NONE", + relaxation_factor=relax, + # backtracking_number=5, + # backtracking_tolerance=1.0, + # backtracking_reduction_factor=0.3, + # backtracking_residual_limit=100.0, + ) + + add_chf_model_disv1d(sim) + + return sim, None + + +def add_chf_model_disv1d(sim): + name = "channel" + chf = flopy.mf6.ModflowChf(sim, modelname=name, save_flows=True) + + disv1d = flopy.mf6.ModflowChfdisv1D( + chf, + nodes=nodes, + nvert=nvert, + width=50.0, + bottom=0.0, + idomain=idomain, + vertices=vertices, + cell1d=cell1d, + xorigin=xorigin, + yorigin=yorigin, + angrot=angrot, + ) + + dfw = flopy.mf6.ModflowChfdfw( + chf, + print_flows=True, + save_flows=True, + length_conversion=1.0, + time_conversion=86400.0, + manningsn=rough, + idcxs=None, + ) + + ic = flopy.mf6.ModflowChfic(chf, strt=1.0) + + # output control + oc = flopy.mf6.ModflowChfoc( + chf, + budget_filerecord=f"{name}.bud", + stage_filerecord=f"{name}.stage", + saverecord=[ + ("STAGE", "ALL"), + ("BUDGET", "ALL"), + ], + printrecord=[ + ("STAGE", "LAST"), + ("BUDGET", "ALL"), + ], + ) + + chd = flopy.mf6.ModflowChfchd( + chf, + maxbound=1, + print_input=True, + print_flows=True, + stress_period_data=[(0, h0), (nodes - 1, h1)], + ) + + return + + +def make_plot(test, mfsim, stage, idx): + print("making plots...") + import matplotlib.pyplot as plt + + chf = mfsim.chf[0] + pmv = flopy.plot.PlotMapView(model=chf) + pmv.plot_array(stage, masked_values=[3e30]) # not working yet + pmv.plot_grid() + + fname = test.workspace / "results.png" + plt.savefig(fname) + + +def check_grb_disv1d(fpth): + grb = flopy.mf6.utils.MfGrdFile(fpth) + assert grb.grid_type == "DISV1D", "grb grid type not DISV1D" + assert grb.ncells == nodes, "grb ncells is incorrect" + assert grb._datadict["NCELLS"] == nodes, "grb nodes is incorrect" + assert grb.verts.shape == ( + nodes + 1, + 2, + ), "vertices shape is incorrect" + assert grb.nja == 4, "nja in grb file is not 4" + assert grb.xorigin == xorigin, "xorigin in grb file is not correct" + assert grb.yorigin == yorigin, "yorigin in grb file is not correct" + assert grb.angrot == angrot, "angrot in grb file is not correct" + assert np.allclose( + grb.bot.reshape((nodes,)), np.zeros((nodes,)) + ), "grb botm not correct" + cellx = np.array( + [0.0, 2 * dx] + ) # node centers pushed all the way to left and right + celly = np.zeros(nreach) + assert np.allclose( + grb._datadict["CELLX"], cellx.flatten() + ), "cellx is not right" + assert np.allclose( + grb._datadict["CELLY"], celly.flatten() + ), "celly is not right" + assert ( + grb._datadict["IAVERT"].shape[0] == nodes + 1 + ), "iavert size not right" + assert ( + grb._datadict["IAVERT"][-1] - 1 == grb._datadict["JAVERT"].shape[0] + ), "javert size not right" + assert ( + grb.ia.shape[0] == grb.ncells + 1 + ), "ia in grb file is not correct size" + assert grb.ja.shape[0] == grb.nja, "ja in grb file is not corect size" + assert np.allclose( + grb.idomain.reshape((nodes,)), idomain.reshape((nodes,)) + ), "grb idomain not correct" + + +def check_output(idx, test): + print(f"evaluating model for case {idx}...") + + modelname = "channel" + ws = test.workspace + mfsim = flopy.mf6.MFSimulation.load(sim_ws=ws) + chf = mfsim.get_model(modelname) + + # check binary grid file + fpth = test.workspace / f"{modelname}.disv1d.grb" + check_grb_disv1d(fpth) + + # read binary stage file + stage = chf.output.stage().get_data().reshape((nodes,)) + assert stage[idomain == 1].max() == 1.0, "maximum stage should be 1.0" + assert stage[idomain == 1].min() == 0.5, "minimum stage should be 0.5" + + # extract the simulated flow from the budget file + fpth = test.workspace / f"{modelname}.bud" + sobj = flopy.utils.CellBudgetFile(fpth, precision="double") + flowjaface = sobj.get_data(text="FLOW-JA-FACE")[-1].flatten() + flow_sim = flowjaface[1] + + # make an independent calculate of the flow between two constant head + # cells + def get_cond_n(depth, width, rough, dhds): + unitconv = 86400.0 + a = depth * width + rh = depth + conveyance = a * rh ** (2.0 / 3.0) / rough + dhds_sqr = dhds**0.5 + c = unitconv * conveyance / dx / dhds_sqr + return c + + cln = 1000.0 + clm = 1000.0 + dhds = abs(h1 - h0) / (cln + clm) + depth_n = h0 - 0.0 + depth_m = h1 - 0.0 + if h0 > h1: + depth_m = depth_n + else: + depth_n = depth_m + + width = 50.0 + rough_n, rough_m = rough + cn = get_cond_n(depth_n, width, rough_n, dhds) + cm = get_cond_n(depth_m, width, rough_m, dhds) + cond = cn * cm / (cn + cm) + flow = cond * (h1 - h0) + print(f"{cn=} {cm=} {cond=}") + print(f"Known flow is {flow} cubic meters per seconds") + print(f"Simulated flow is {flow_sim} cubic meters per seconds") + assert np.allclose( + flow, flow_sim + ), "known flow and simulated flow not the same" + + makeplot = False + if makeplot: + # PlotMapView not working yet for disv1d + make_plot(test, mfsim, stage.flatten(), idx) + + +@pytest.mark.developmode +@pytest.mark.parametrize("idx, name", enumerate(cases)) +def test_mf6model(idx, name, function_tmpdir, targets): + test = TestFramework( + name=name, + workspace=function_tmpdir, + build=lambda t: build_models(idx, t), + check=lambda t: check_output(idx, t), + targets=targets, + ) + test.run() diff --git a/autotest/test_chf_gwf.py b/autotest/test_chf_gwf.py index 08bf3d0bd47..611890efb15 100644 --- a/autotest/test_chf_gwf.py +++ b/autotest/test_chf_gwf.py @@ -82,22 +82,21 @@ def add_chf_model(sim): total_length = dx * nreach vertices = [] vertices = [[j, j * dx, 0.0] for j in range(nreach + 1)] - cell2d = [] + cell1d = [] for j in range(nreach): - cell2d.append([j, 0.5, 2, j, j + 1]) - nodes = len(cell2d) + cell1d.append([j, 0.5, 2, j, j + 1]) + nodes = len(cell1d) nvert = len(vertices) disv1d = flopy.mf6.ModflowChfdisv1D( chf, nodes=nodes, nvert=nvert, - length=dx, width=50.0, bottom=0.0, idomain=1, vertices=vertices, - cell2d=cell2d, + cell1d=cell1d, ) dfw = flopy.mf6.ModflowChfdfw( diff --git a/autotest/test_olf_dfw_swrt2dis.py b/autotest/test_olf_dfw_swrt2dis.py index 63634c38343..b9f381a4fe2 100644 --- a/autotest/test_olf_dfw_swrt2dis.py +++ b/autotest/test_olf_dfw_swrt2dis.py @@ -73,9 +73,9 @@ def build_models(idx, test): nrow = 11 ncol = 11 - botm = np.empty((nrow, ncol), dtype=float) + bottom = np.empty((nrow, ncol), dtype=float) for i in range(nrow): - botm[i, :] = np.linspace(1.05, 0.05, nrow) + bottom[i, :] = np.linspace(1.05, 0.05, nrow) dis = flopy.mf6.ModflowOlfdis2D( olf, @@ -83,7 +83,7 @@ def build_models(idx, test): ncol=ncol, delr=dx, delc=dx, - botm=botm, + bottom=bottom, ) dfw = flopy.mf6.ModflowOlfdfw( @@ -206,7 +206,7 @@ def check_output(idx, test): # at end of simulation, water depth should be 1.0 for all reaches olf = mfsim.get_model(olfname) - depth = stage_all[-1] - olf.dis.botm.array + depth = stage_all[-1] - olf.dis.bottom.array assert np.allclose( depth, 1.0 ), f"Simulated depth at end should be 1, but found {depth}" diff --git a/autotest/test_olf_dis.py b/autotest/test_olf_dis.py new file mode 100644 index 00000000000..f189e2967ea --- /dev/null +++ b/autotest/test_olf_dis.py @@ -0,0 +1,327 @@ +""" + +Test the dis2d and disv2d discretization packages so that they support +IDOMAIN and have valid binary grid files. + +""" + +import flopy +import numpy as np +import pytest +from flopy.utils.gridutil import get_disv_kwargs +from framework import TestFramework + +cases = [ + "olf-dis2d", + "olf-disv2d", +] + +# grid size +dx = 1000.0 +nrow = 10 +ncol = 10 +ncpl = nrow * ncol +xorigin = 75.0 +yorigin = 500.0 +angrot = 32.0 +botm = np.zeros((nrow, ncol), dtype=float) +idomain = np.ones((nrow, ncol), dtype=int) +idomain[4:6, 4:6] = 0 + + +def build_models(idx, test): + perlen = [1] # 1 second + nstp = [1] + tsmult = [1.0] + nper = len(perlen) + + tdis_rc = [] + for i in range(nper): + tdis_rc.append((perlen[i], nstp[i], tsmult[i])) + + name = "olf" + + # build MODFLOW 6 files + ws = test.workspace + sim = flopy.mf6.MFSimulation( + sim_name=f"{name}_sim", version="mf6", exe_name="mf6", sim_ws=ws + ) + # create tdis package + tdis = flopy.mf6.ModflowTdis( + sim, time_units="SECONDS", nper=nper, perioddata=tdis_rc + ) + + nouter, ninner = 100, 50 + hclose, rclose, relax = 1e-6, 1e-8, 1.0 + ims = flopy.mf6.ModflowIms( + sim, + print_option="SUMMARY", + outer_dvclose=hclose, + outer_maximum=nouter, + under_relaxation="DBD", + under_relaxation_theta=0.9, + under_relaxation_kappa=0.0001, + under_relaxation_gamma=0.0, + inner_maximum=ninner, + inner_dvclose=hclose, + linear_acceleration="BICGSTAB", + scaling_method="NONE", + reordering_method="NONE", + relaxation_factor=relax, + # backtracking_number=5, + # backtracking_tolerance=1.0, + # backtracking_reduction_factor=0.3, + # backtracking_residual_limit=100.0, + ) + + if idx == 0: + add_olf_model_dis2d(sim) + elif idx == 1: + add_olf_model_disv2d(sim) + + return sim, None + + +def add_olf_model_dis2d(sim): + name = "overland" + olf = flopy.mf6.ModflowOlf(sim, modelname=name, save_flows=True) + + dis = flopy.mf6.ModflowOlfdis2D( + olf, + nrow=nrow, + ncol=ncol, + delr=dx, + delc=dx, + bottom=botm, + idomain=idomain, + xorigin=xorigin, + yorigin=yorigin, + angrot=angrot, + ) + + dfw = flopy.mf6.ModflowOlfdfw( + olf, + print_flows=True, + save_flows=True, + length_conversion=1.0, + time_conversion=86400.0, + manningsn=3.5, + idcxs=None, + ) + + ic = flopy.mf6.ModflowOlfic(olf, strt=1.0) + + # output control + oc = flopy.mf6.ModflowOlfoc( + olf, + budget_filerecord=f"{name}.bud", + stage_filerecord=f"{name}.stage", + saverecord=[ + ("STAGE", "ALL"), + ("BUDGET", "ALL"), + ], + printrecord=[ + ("STAGE", "LAST"), + ("BUDGET", "ALL"), + ], + ) + + chd = flopy.mf6.ModflowOlfchd( + olf, + maxbound=1, + print_input=True, + print_flows=True, + stress_period_data=[(0, 0, 1.0), (nrow - 1, ncol - 1, 0.5)], + ) + + return + + +def add_olf_model_disv2d(sim): + name = "overland" + olf = flopy.mf6.ModflowOlf(sim, modelname=name, save_flows=True) + + disvkwargs = get_disv_kwargs( + 1, + nrow, + ncol, + dx, + dx, + 1.0, # top + 0.0, # botm + 0.0, # xoffset + 0.0, # yoffset + ) + + _ = disvkwargs.pop("top") + _ = disvkwargs.pop("nlay") + bottom = disvkwargs.pop("botm").reshape((ncpl,)) + disvkwargs["nodes"] = disvkwargs.pop("ncpl") + + dis = flopy.mf6.ModflowOlfdisv2D( + olf, + bottom=bottom, + idomain=idomain.reshape((ncpl,)), + xorigin=xorigin, + yorigin=yorigin, + angrot=angrot, + **disvkwargs, + ) + + dfw = flopy.mf6.ModflowOlfdfw( + olf, + print_flows=True, + save_flows=True, + length_conversion=1.0, + time_conversion=86400.0, + manningsn=3.5, + idcxs=None, + ) + + ic = flopy.mf6.ModflowOlfic(olf, strt=1.0) + + # output control + oc = flopy.mf6.ModflowOlfoc( + olf, + budget_filerecord=f"{name}.bud", + stage_filerecord=f"{name}.stage", + saverecord=[ + ("STAGE", "ALL"), + ("BUDGET", "ALL"), + ], + printrecord=[ + ("STAGE", "LAST"), + ("BUDGET", "ALL"), + ], + ) + + chd = flopy.mf6.ModflowOlfchd( + olf, + maxbound=1, + print_input=True, + print_flows=True, + stress_period_data=[(0, 1.0), (ncpl - 1, 0.5)], + ) + + return + + +def make_plot(test, mfsim, stage, idx): + print("making plots...") + import matplotlib.pyplot as plt + + olf = mfsim.olf[0] + pmv = flopy.plot.PlotMapView(model=olf) + pmv.plot_array(stage, masked_values=[3e30]) + pmv.plot_grid() + + fname = test.workspace / "results.png" + plt.savefig(fname) + + +def check_grb_dis2d(fpth): + grb = flopy.mf6.utils.MfGrdFile(fpth) + assert grb.grid_type == "DIS2D", "grb grid type not DIS2D" + assert grb.ncells == nrow * ncol, "grb ncells is incorrect" + assert grb.nrow == nrow, "nrow in grb file is not 10" + assert grb.ncol == ncol, "ncol in grb file is not 10" + assert grb.nja == 432, "nja in grb file is not 432" + assert grb.xorigin == xorigin, "xorigin in grb file is not correct" + assert grb.yorigin == yorigin, "yorigin in grb file is not correct" + assert grb.angrot == angrot, "angrot in grb file is not correct" + assert np.allclose(grb.delr, dx * np.ones((ncol))), "grb delr not correct" + assert np.allclose(grb.delc, dx * np.ones((nrow))), "grb delc not correct" + assert np.allclose( + grb.bot.reshape((nrow, ncol)), np.zeros((nrow, ncol)) + ), "grb botm not correct" + assert ( + grb.ia.shape[0] == grb.ncells + 1 + ), "ia in grb file is not correct size" + assert grb.ja.shape[0] == grb.nja, "ja in grb file is not corect size" + assert np.allclose( + grb.idomain.reshape((nrow, ncol)), idomain + ), "grb idomain not correct" + + +def check_grb_disv2d(fpth): + grb = flopy.mf6.utils.MfGrdFile(fpth) + assert grb.grid_type == "DISV2D", "grb grid type not DISV2D" + assert grb.ncells == ncpl, "grb ncells is incorrect" + assert grb._datadict["NODES"] == 96, "grb nodes is incorrect" + assert grb.verts.shape == ( + (nrow + 1) * (ncol + 1), + 2, + ), "vertices shape is incorrect" + assert grb.nja == 432, "nja in grb file is not 432" + assert grb.xorigin == xorigin, "xorigin in grb file is not correct" + assert grb.yorigin == yorigin, "yorigin in grb file is not correct" + assert grb.angrot == angrot, "angrot in grb file is not correct" + assert np.allclose( + grb.bot.reshape((nrow, ncol)), np.zeros((nrow, ncol)) + ), "grb botm not correct" + cellx, celly = np.meshgrid( + np.linspace(dx / 2, ncol * dx - dx / 2, ncol), + np.linspace(dx * nrow - dx / 2.0, dx / 2, nrow), + ) + assert np.allclose( + grb._datadict["CELLX"], cellx.flatten() + ), "cellx is not right" + assert np.allclose( + grb._datadict["CELLY"], celly.flatten() + ), "celly is not right" + assert ( + grb._datadict["IAVERT"].shape[0] == ncpl + 1 + ), "iavert size not right" + assert ( + grb._datadict["IAVERT"][-1] - 1 == grb._datadict["JAVERT"].shape[0] + ), "javert size not right" + assert ( + grb.ia.shape[0] == grb.ncells + 1 + ), "ia in grb file is not correct size" + assert grb.ja.shape[0] == grb.nja, "ja in grb file is not corect size" + assert np.allclose( + grb.idomain.reshape((ncpl,)), idomain.reshape((ncpl,)) + ), "grb idomain not correct" + + +def check_output(idx, test): + print(f"evaluating model for case {idx}...") + + modelname = "overland" + ws = test.workspace + mfsim = flopy.mf6.MFSimulation.load(sim_ws=ws) + + # read the binary grid file + if idx == 0: + fpth = test.workspace / f"{modelname}.dis2d.grb" + check_grb_dis2d(fpth) + elif idx == 1: + fpth = test.workspace / f"{modelname}.disv2d.grb" + check_grb_disv2d(fpth) + + # read binary stage file + fpth = test.workspace / f"{modelname}.stage" + sobj = flopy.utils.HeadFile(fpth, precision="double", text="STAGE") + stage = sobj.get_data().reshape((nrow, ncol)) + assert np.allclose( + stage[idomain == 0], 3.0e30 + ), "stage should have nodata values where idomain is zero" + assert stage[idomain == 1].max() == 1.0, "maximum stage should be 1.0" + assert stage[idomain == 1].min() == 0.5, "minimum stage should be 0.5" + + makeplot = False + if makeplot: + make_plot(test, mfsim, stage.flatten(), idx) + + +@pytest.mark.developmode +@pytest.mark.parametrize("idx, name", enumerate(cases)) +def test_mf6model(idx, name, function_tmpdir, targets): + test = TestFramework( + name=name, + workspace=function_tmpdir, + build=lambda t: build_models(idx, t), + check=lambda t: check_output(idx, t), + targets=targets, + ) + test.run() diff --git a/autotest/test_olf_gwf.py b/autotest/test_olf_gwf.py index b17a2fe4d7b..918bf1ccc23 100644 --- a/autotest/test_olf_gwf.py +++ b/autotest/test_olf_gwf.py @@ -88,7 +88,7 @@ def add_olf_model(sim): ncol=ncol, delr=dx, delc=dx, - botm=botm, + bottom=botm, ) dfw = flopy.mf6.ModflowOlfdfw( diff --git a/autotest/test_olf_vcatch.py b/autotest/test_olf_vcatch.py index e29a66e0540..a18794ad650 100644 --- a/autotest/test_olf_vcatch.py +++ b/autotest/test_olf_vcatch.py @@ -119,7 +119,7 @@ def build_models(idx, test): ) sim.register_ims_package(imsolf, [olf.name]) - botm = land_surface.reshape((nrow, ncol)) + bottom = land_surface.reshape((nrow, ncol)) dis = flopy.mf6.ModflowOlfdis2D( olf, export_array_ascii=True, @@ -127,7 +127,7 @@ def build_models(idx, test): ncol=ncol, delr=dx, delc=dx, - botm=botm, + bottom=bottom, xorigin=-810, ) @@ -153,7 +153,7 @@ def build_models(idx, test): ic = flopy.mf6.ModflowOlfic( olf, export_array_ascii=True, - strt=botm, + strt=bottom, ) # output control @@ -272,7 +272,7 @@ def check_output(idx, test): # ensure export array is working properly flist = [ - "dis2d.botm", + "dis2d.bottom", "dis2d.delc", "dis2d.delr", "dfw.manningsn", diff --git a/doc/mf6io/mf6ivar/dfn/chf-dis2d.dfn b/doc/mf6io/mf6ivar/dfn/chf-dis2d.dfn index 96b1837b572..439cdf91328 100644 --- a/doc/mf6io/mf6ivar/dfn/chf-dis2d.dfn +++ b/doc/mf6io/mf6ivar/dfn/chf-dis2d.dfn @@ -90,7 +90,7 @@ description is the row spacing in the column direction. default_value 1.0 block griddata -name botm +name bottom type double precision shape (ncol, nrow) reader readarray diff --git a/doc/mf6io/mf6ivar/dfn/chf-disv1d.dfn b/doc/mf6io/mf6ivar/dfn/chf-disv1d.dfn index 8e6b3c10f69..040add764ba 100644 --- a/doc/mf6io/mf6ivar/dfn/chf-disv1d.dfn +++ b/doc/mf6io/mf6ivar/dfn/chf-disv1d.dfn @@ -69,17 +69,6 @@ description is the total number of (x, y, z) vertex pairs used to characterize t # --------------------- swf disv1d griddata --------------------- -block griddata -name length -type double precision -shape (nodes) -valid -reader readarray -layered false -optional -longname length -description length for each one-dimensional cell - block griddata name width type double precision @@ -155,29 +144,29 @@ optional false longname y-coordinate for vertex description is the y-coordinate for the vertex. -# --------------------- swf disv1d cell2d --------------------- +# --------------------- swf disv1d cell1d --------------------- -block cell2d -name cell2d -type recarray icell2d fdc ncvert icvert +block cell1d +name cell1d +type recarray icell1d fdc ncvert icvert shape (nodes) reader urword optional false -longname cell2d data +longname cell1d data description -block cell2d -name icell2d +block cell1d +name icell1d type integer in_record true tagged false reader urword optional false -longname cell2d number -description is the cell2d number. Records in the cell2d block must be listed in consecutive order from the first to the last. +longname cell1d number +description is the cell1d number. Records in the cell1d block must be listed in consecutive order from the first to the last. numeric_index true -block cell2d +block cell1d name fdc type double precision in_record true @@ -187,7 +176,7 @@ optional false longname fractional distance to the cell center description is the fractional distance to the cell center. FDC is relative to the first vertex in the ICVERT array. In most cases FDC should be 0.5, which would place the center of the line segment that defines the cell. If the value of FDC is 1, the cell center would located at the last vertex. FDC values of 0 and 1 can be used to place the node at either end of the cell which can be useful for cells with boundary conditions. -block cell2d +block cell1d name ncvert type integer in_record true @@ -197,7 +186,7 @@ optional false longname number of cell vertices description is the number of vertices required to define the cell. There may be a different number of vertices for each cell. -block cell2d +block cell1d name icvert type integer shape (ncvert) diff --git a/doc/mf6io/mf6ivar/dfn/olf-dis2d.dfn b/doc/mf6io/mf6ivar/dfn/olf-dis2d.dfn index 96b1837b572..439cdf91328 100644 --- a/doc/mf6io/mf6ivar/dfn/olf-dis2d.dfn +++ b/doc/mf6io/mf6ivar/dfn/olf-dis2d.dfn @@ -90,7 +90,7 @@ description is the row spacing in the column direction. default_value 1.0 block griddata -name botm +name bottom type double precision shape (ncol, nrow) reader readarray diff --git a/doc/mf6io/mf6ivar/dfn/swf-dis2d.dfn b/doc/mf6io/mf6ivar/dfn/swf-dis2d.dfn index 96b1837b572..439cdf91328 100644 --- a/doc/mf6io/mf6ivar/dfn/swf-dis2d.dfn +++ b/doc/mf6io/mf6ivar/dfn/swf-dis2d.dfn @@ -90,7 +90,7 @@ description is the row spacing in the column direction. default_value 1.0 block griddata -name botm +name bottom type double precision shape (ncol, nrow) reader readarray diff --git a/doc/mf6io/mf6ivar/dfn/swf-disv1d.dfn b/doc/mf6io/mf6ivar/dfn/swf-disv1d.dfn index 8e6b3c10f69..040add764ba 100644 --- a/doc/mf6io/mf6ivar/dfn/swf-disv1d.dfn +++ b/doc/mf6io/mf6ivar/dfn/swf-disv1d.dfn @@ -69,17 +69,6 @@ description is the total number of (x, y, z) vertex pairs used to characterize t # --------------------- swf disv1d griddata --------------------- -block griddata -name length -type double precision -shape (nodes) -valid -reader readarray -layered false -optional -longname length -description length for each one-dimensional cell - block griddata name width type double precision @@ -155,29 +144,29 @@ optional false longname y-coordinate for vertex description is the y-coordinate for the vertex. -# --------------------- swf disv1d cell2d --------------------- +# --------------------- swf disv1d cell1d --------------------- -block cell2d -name cell2d -type recarray icell2d fdc ncvert icvert +block cell1d +name cell1d +type recarray icell1d fdc ncvert icvert shape (nodes) reader urword optional false -longname cell2d data +longname cell1d data description -block cell2d -name icell2d +block cell1d +name icell1d type integer in_record true tagged false reader urword optional false -longname cell2d number -description is the cell2d number. Records in the cell2d block must be listed in consecutive order from the first to the last. +longname cell1d number +description is the cell1d number. Records in the cell1d block must be listed in consecutive order from the first to the last. numeric_index true -block cell2d +block cell1d name fdc type double precision in_record true @@ -187,7 +176,7 @@ optional false longname fractional distance to the cell center description is the fractional distance to the cell center. FDC is relative to the first vertex in the ICVERT array. In most cases FDC should be 0.5, which would place the center of the line segment that defines the cell. If the value of FDC is 1, the cell center would located at the last vertex. FDC values of 0 and 1 can be used to place the node at either end of the cell which can be useful for cells with boundary conditions. -block cell2d +block cell1d name ncvert type integer in_record true @@ -197,7 +186,7 @@ optional false longname number of cell vertices description is the number of vertices required to define the cell. There may be a different number of vertices for each cell. -block cell2d +block cell1d name icvert type integer shape (ncvert) diff --git a/src/Idm/chf-dis2didm.f90 b/src/Idm/chf-dis2didm.f90 index 8f4eabec9c8..47c2051c179 100644 --- a/src/Idm/chf-dis2didm.f90 +++ b/src/Idm/chf-dis2didm.f90 @@ -22,7 +22,7 @@ module ChfDis2DInputModule logical :: ncol = .false. logical :: delr = .false. logical :: delc = .false. - logical :: botm = .false. + logical :: bottom = .false. logical :: idomain = .false. end type ChfDis2dParamFoundType @@ -215,13 +215,13 @@ module ChfDis2DInputModule ) type(InputParamDefinitionType), parameter :: & - chfdis2d_botm = InputParamDefinitionType & + chfdis2d_bottom = InputParamDefinitionType & ( & 'CHF', & ! component 'DIS2D', & ! subcomponent 'GRIDDATA', & ! block - 'BOTM', & ! tag name - 'BOTM', & ! fortran variable + 'BOTTOM', & ! tag name + 'BOTTOM', & ! fortran variable 'DOUBLE2D', & ! type 'NCOL NROW', & ! shape 'cell bottom elevation', & ! longname @@ -263,7 +263,7 @@ module ChfDis2DInputModule chfdis2d_ncol, & chfdis2d_delr, & chfdis2d_delc, & - chfdis2d_botm, & + chfdis2d_bottom, & chfdis2d_idomain & ] diff --git a/src/Idm/chf-disv1didm.f90 b/src/Idm/chf-disv1didm.f90 index e678128bf8d..331cd46438e 100644 --- a/src/Idm/chf-disv1didm.f90 +++ b/src/Idm/chf-disv1didm.f90 @@ -20,14 +20,13 @@ module ChfDisv1DInputModule logical :: export_ascii = .false. logical :: nodes = .false. logical :: nvert = .false. - logical :: length = .false. logical :: width = .false. logical :: bottom = .false. logical :: idomain = .false. logical :: iv = .false. logical :: xv = .false. logical :: yv = .false. - logical :: icell2d = .false. + logical :: icell1d = .false. logical :: fdc = .false. logical :: ncvert = .false. logical :: icvert = .false. @@ -185,24 +184,6 @@ module ChfDisv1DInputModule .false. & ! timeseries ) - type(InputParamDefinitionType), parameter :: & - chfdisv1d_length = InputParamDefinitionType & - ( & - 'CHF', & ! component - 'DISV1D', & ! subcomponent - 'GRIDDATA', & ! block - 'LENGTH', & ! tag name - 'LENGTH', & ! fortran variable - 'DOUBLE1D', & ! type - 'NODES', & ! shape - 'length', & ! longname - .true., & ! required - .false., & ! multi-record - .false., & ! preserve case - .false., & ! layered - .false. & ! timeseries - ) - type(InputParamDefinitionType), parameter :: & chfdisv1d_width = InputParamDefinitionType & ( & @@ -312,16 +293,16 @@ module ChfDisv1DInputModule ) type(InputParamDefinitionType), parameter :: & - chfdisv1d_icell2d = InputParamDefinitionType & + chfdisv1d_icell1d = InputParamDefinitionType & ( & 'CHF', & ! component 'DISV1D', & ! subcomponent - 'CELL2D', & ! block - 'ICELL2D', & ! tag name - 'ICELL2D', & ! fortran variable + 'CELL1D', & ! block + 'ICELL1D', & ! tag name + 'ICELL1D', & ! fortran variable 'INTEGER', & ! type '', & ! shape - 'cell2d number', & ! longname + 'cell1d number', & ! longname .true., & ! required .true., & ! multi-record .false., & ! preserve case @@ -334,7 +315,7 @@ module ChfDisv1DInputModule ( & 'CHF', & ! component 'DISV1D', & ! subcomponent - 'CELL2D', & ! block + 'CELL1D', & ! block 'FDC', & ! tag name 'FDC', & ! fortran variable 'DOUBLE', & ! type @@ -352,7 +333,7 @@ module ChfDisv1DInputModule ( & 'CHF', & ! component 'DISV1D', & ! subcomponent - 'CELL2D', & ! block + 'CELL1D', & ! block 'NCVERT', & ! tag name 'NCVERT', & ! fortran variable 'INTEGER', & ! type @@ -370,7 +351,7 @@ module ChfDisv1DInputModule ( & 'CHF', & ! component 'DISV1D', & ! subcomponent - 'CELL2D', & ! block + 'CELL1D', & ! block 'ICVERT', & ! tag name 'ICVERT', & ! fortran variable 'INTEGER1D', & ! type @@ -394,14 +375,13 @@ module ChfDisv1DInputModule chfdisv1d_export_ascii, & chfdisv1d_nodes, & chfdisv1d_nvert, & - chfdisv1d_length, & chfdisv1d_width, & chfdisv1d_bottom, & chfdisv1d_idomain, & chfdisv1d_iv, & chfdisv1d_xv, & chfdisv1d_yv, & - chfdisv1d_icell2d, & + chfdisv1d_icell1d, & chfdisv1d_fdc, & chfdisv1d_ncvert, & chfdisv1d_icvert & @@ -426,16 +406,16 @@ module ChfDisv1DInputModule ) type(InputParamDefinitionType), parameter :: & - chfdisv1d_cell2d = InputParamDefinitionType & + chfdisv1d_cell1d = InputParamDefinitionType & ( & 'CHF', & ! component 'DISV1D', & ! subcomponent - 'CELL2D', & ! block - 'CELL2D', & ! tag name - 'CELL2D', & ! fortran variable - 'RECARRAY ICELL2D FDC NCVERT ICVERT', & ! type + 'CELL1D', & ! block + 'CELL1D', & ! tag name + 'CELL1D', & ! fortran variable + 'RECARRAY ICELL1D FDC NCVERT ICVERT', & ! type 'NODES', & ! shape - 'cell2d data', & ! longname + 'cell1d data', & ! longname .true., & ! required .false., & ! multi-record .false., & ! preserve case @@ -447,7 +427,7 @@ module ChfDisv1DInputModule chf_disv1d_aggregate_definitions(*) = & [ & chfdisv1d_vertices, & - chfdisv1d_cell2d & + chfdisv1d_cell1d & ] type(InputBlockDefinitionType), parameter :: & @@ -478,7 +458,7 @@ module ChfDisv1DInputModule .false. & ! block_variable ), & InputBlockDefinitionType( & - 'CELL2D', & ! blockname + 'CELL1D', & ! blockname .true., & ! required .true., & ! aggregate .false. & ! block_variable diff --git a/src/Idm/olf-dis2didm.f90 b/src/Idm/olf-dis2didm.f90 index 4f3b669f83a..16aa1095ab3 100644 --- a/src/Idm/olf-dis2didm.f90 +++ b/src/Idm/olf-dis2didm.f90 @@ -22,7 +22,7 @@ module OlfDis2DInputModule logical :: ncol = .false. logical :: delr = .false. logical :: delc = .false. - logical :: botm = .false. + logical :: bottom = .false. logical :: idomain = .false. end type OlfDis2dParamFoundType @@ -215,13 +215,13 @@ module OlfDis2DInputModule ) type(InputParamDefinitionType), parameter :: & - olfdis2d_botm = InputParamDefinitionType & + olfdis2d_bottom = InputParamDefinitionType & ( & 'OLF', & ! component 'DIS2D', & ! subcomponent 'GRIDDATA', & ! block - 'BOTM', & ! tag name - 'BOTM', & ! fortran variable + 'BOTTOM', & ! tag name + 'BOTTOM', & ! fortran variable 'DOUBLE2D', & ! type 'NCOL NROW', & ! shape 'cell bottom elevation', & ! longname @@ -263,7 +263,7 @@ module OlfDis2DInputModule olfdis2d_ncol, & olfdis2d_delr, & olfdis2d_delc, & - olfdis2d_botm, & + olfdis2d_bottom, & olfdis2d_idomain & ] diff --git a/src/Idm/swf-dis2didm.f90 b/src/Idm/swf-dis2didm.f90 index f46547a178c..d0419084dd3 100644 --- a/src/Idm/swf-dis2didm.f90 +++ b/src/Idm/swf-dis2didm.f90 @@ -22,7 +22,7 @@ module SwfDis2DInputModule logical :: ncol = .false. logical :: delr = .false. logical :: delc = .false. - logical :: botm = .false. + logical :: bottom = .false. logical :: idomain = .false. end type SwfDis2dParamFoundType @@ -215,13 +215,13 @@ module SwfDis2DInputModule ) type(InputParamDefinitionType), parameter :: & - swfdis2d_botm = InputParamDefinitionType & + swfdis2d_bottom = InputParamDefinitionType & ( & 'SWF', & ! component 'DIS2D', & ! subcomponent 'GRIDDATA', & ! block - 'BOTM', & ! tag name - 'BOTM', & ! fortran variable + 'BOTTOM', & ! tag name + 'BOTTOM', & ! fortran variable 'DOUBLE2D', & ! type 'NCOL NROW', & ! shape 'cell bottom elevation', & ! longname @@ -263,7 +263,7 @@ module SwfDis2DInputModule swfdis2d_ncol, & swfdis2d_delr, & swfdis2d_delc, & - swfdis2d_botm, & + swfdis2d_bottom, & swfdis2d_idomain & ] diff --git a/src/Model/Discretization/Dis2d.f90 b/src/Model/Discretization/Dis2d.f90 index 7c37d73e224..6e1446e2c8a 100644 --- a/src/Model/Discretization/Dis2d.f90 +++ b/src/Model/Discretization/Dis2d.f90 @@ -25,7 +25,7 @@ module Dis2dModule integer(I4B), pointer :: ncol => null() !< number of columns real(DP), dimension(:), pointer, contiguous :: delr => null() !< spacing along a row real(DP), dimension(:), pointer, contiguous :: delc => null() !< spacing along a column - real(DP), dimension(:, :), pointer, contiguous :: botm => null() !< bottom elevations for each cell (ncol, nrow) + real(DP), dimension(:, :), pointer, contiguous :: bottom => null() !< bottom elevations for each cell (ncol, nrow) integer(I4B), dimension(:, :), pointer, contiguous :: idomain => null() !< idomain (ncol, nrow) real(DP), dimension(:), pointer, contiguous :: cellx => null() !< cell center x coordinate for column j real(DP), dimension(:), pointer, contiguous :: celly => null() !< cell center y coordinate for row i @@ -80,8 +80,7 @@ module Dis2dModule logical :: ncol = .false. logical :: delr = .false. logical :: delc = .false. - logical :: top = .false. - logical :: botm = .false. + logical :: bottom = .false. logical :: idomain = .false. end type DisFoundtype @@ -167,7 +166,7 @@ subroutine dis3d_da(this) ! -- Deallocate Arrays call mem_deallocate(this%nodereduced) call mem_deallocate(this%nodeuser) - call mem_deallocate(this%botm) + call mem_deallocate(this%bottom) call mem_deallocate(this%idomain) ! end subroutine dis3d_da @@ -270,7 +269,7 @@ subroutine source_dimensions(this) call mem_allocate(this%delc, this%nrow, 'DELC', this%memoryPath) call mem_allocate(this%idomain, this%ncol, this%nrow, 'IDOMAIN', & this%memoryPath) - call mem_allocate(this%botm, this%ncol, this%nrow, 'BOTM', & + call mem_allocate(this%bottom, this%ncol, this%nrow, 'BOTTOM', & this%memoryPath) call mem_allocate(this%cellx, this%ncol, 'CELLX', this%memoryPath) call mem_allocate(this%celly, this%nrow, 'CELLY', this%memoryPath) @@ -315,7 +314,7 @@ subroutine source_griddata(this) ! -- update defaults with idm sourced values call mem_set_value(this%delr, 'DELR', this%input_mempath, found%delr) call mem_set_value(this%delc, 'DELC', this%input_mempath, found%delc) - call mem_set_value(this%botm, 'BOTM', this%input_mempath, found%botm) + call mem_set_value(this%bottom, 'BOTTOM', this%input_mempath, found%bottom) call mem_set_value(this%idomain, 'IDOMAIN', this%input_mempath, found%idomain) ! ! -- log simulation values @@ -342,12 +341,8 @@ subroutine log_griddata(this, found) write (this%iout, '(4x,a)') 'DELC set from input file' end if ! - if (found%top) then - write (this%iout, '(4x,a)') 'TOP set from input file' - end if - ! - if (found%botm) then - write (this%iout, '(4x,a)') 'BOTM set from input file' + if (found%bottom) then + write (this%iout, '(4x,a)') 'BOTTOM set from input file' end if ! if (found%idomain) then @@ -453,7 +448,7 @@ subroutine grid_finalize(this) DHALF * this%delc(i) end do ! - ! -- Move botm into bot, and calculate area + ! -- Move bottom into bot, and calculate area node = 0 do i = 1, this%nrow do j = 1, this%ncol @@ -461,7 +456,7 @@ subroutine grid_finalize(this) noder = node if (this%nodes < this%nodesuser) noder = this%nodereduced(node) if (noder <= 0) cycle - this%bot(noder) = this%botm(j, i) + this%bot(noder) = this%bottom(j, i) this%area(noder) = this%delr(j) * this%delc(i) this%xc(noder) = this%cellx(j) this%yc(noder) = this%celly(i) @@ -578,7 +573,7 @@ subroutine write_grb(this, icelltype) write (iunit) this%angrot ! angrot write (iunit) this%delr ! delr write (iunit) this%delc ! delc - write (iunit) this%botm ! botm + write (iunit) this%bottom ! bottom write (iunit) this%con%iausr ! iausr write (iunit) this%con%jausr ! jausr write (iunit) this%idomain ! idomain diff --git a/src/Model/Discretization/Disv1d.f90 b/src/Model/Discretization/Disv1d.f90 index 4d9596827b0..497a9623bc4 100644 --- a/src/Model/Discretization/Disv1d.f90 +++ b/src/Model/Discretization/Disv1d.f90 @@ -20,10 +20,10 @@ module Disv1dModule type, extends(DisBaseType) :: Disv1dType integer(I4B), pointer :: nvert => null() !< number of x,y vertices - real(DP), dimension(:), pointer, contiguous :: length => null() !< length of each reach + real(DP), dimension(:), pointer, contiguous :: length => null() !< length of each reach (of size nodesuser) real(DP), dimension(:), pointer, contiguous :: width => null() !< reach width real(DP), dimension(:), pointer, contiguous :: bottom => null() !< reach bottom elevation - integer(I4B), dimension(:), pointer, contiguous :: idomain => null() !< idomain (nodes) + integer(I4B), dimension(:), pointer, contiguous :: idomain => null() !< idomain (of size nodesuser) real(DP), dimension(:, :), pointer, contiguous :: vertices => null() !< cell vertices stored as 2d array with columns of x, y real(DP), dimension(:, :), pointer, contiguous :: cellxy => null() !< reach midpoints stored as 2d array with columns of x, y real(DP), dimension(:), pointer, contiguous :: fdc => null() !< fdc stored as array @@ -45,7 +45,7 @@ module Disv1dModule procedure :: source_dimensions procedure :: source_griddata procedure :: source_vertices - procedure :: source_cell2d + procedure :: source_cell1d procedure :: log_options procedure :: log_dimensions procedure :: log_griddata @@ -71,14 +71,13 @@ module Disv1dModule logical :: angrot = .false. logical :: nodes = .false. logical :: nvert = .false. - logical :: length = .false. logical :: width = .false. logical :: bottom = .false. logical :: idomain = .false. logical :: iv = .false. logical :: xv = .false. logical :: yv = .false. - logical :: icell2d = .false. + logical :: icell1d = .false. logical :: fdc = .false. logical :: ncvert = .false. logical :: icvert = .false. @@ -138,9 +137,6 @@ subroutine disv1d_df(this) call this%disv1d_load() end if - ! create connectivity using vertices and cell2d - call this%create_connections() - ! finalize the grid call this%grid_finalize() @@ -271,7 +267,7 @@ subroutine disv1d_load(this) ! If vertices provided by user, read and store vertices if (this%nvert > 0) then call this%source_vertices() - call this%source_cell2d() + call this%source_cell1d() end if end subroutine disv1d_load @@ -386,7 +382,7 @@ subroutine source_dimensions(this) if (this%nvert < 1) then call store_warning( & 'NVERT was not specified or was specified as zero. The & - &VERTICES and CELL2D blocks will not be read for the DISV1D6 & + &VERTICES and CELL1D blocks will not be read for the DISV1D6 & &Package in model '//trim(this%memoryPath)//'.') end if ! @@ -402,12 +398,12 @@ subroutine source_dimensions(this) ! ! -- Allocate vertices array if (this%nvert > 0) then - call mem_allocate(this%vertices, 3, this%nvert, & + call mem_allocate(this%vertices, 2, this%nvert, & 'VERTICES', this%memoryPath) call mem_allocate(this%fdc, this%nodesuser, & 'FDC', this%memoryPath) - call mem_allocate(this%cellxy, 3, this%nodesuser, & - 'CELLXYZ', this%memoryPath) + call mem_allocate(this%cellxy, 2, this%nodesuser, & + 'CELLXY', this%memoryPath) end if ! ! -- initialize all cells to be active (idomain = 1) @@ -443,33 +439,25 @@ subroutine log_dimensions(this, found) end subroutine log_dimensions subroutine source_griddata(this) - ! -- modules + ! modules use MemoryManagerExtModule, only: mem_set_value use SimVariablesModule, only: idm_context - ! -- dummy + ! dummy class(Disv1dType) :: this - ! -- locals + ! locals character(len=LENMEMPATH) :: idmMemoryPath type(DisFoundType) :: found - ! -- formats - ! - ! -- set memory path + ! formats + + ! set memory path idmMemoryPath = create_mem_path(this%name_model, 'DISV1D', idm_context) - ! - ! -- update defaults with idm sourced values - call mem_set_value(this%length, 'LENGTH', idmMemoryPath, & - found%length) + call mem_set_value(this%width, 'WIDTH', idmMemoryPath, & found%width) call mem_set_value(this%bottom, 'BOTTOM', idmMemoryPath, & found%bottom) call mem_set_value(this%idomain, 'IDOMAIN', idmMemoryPath, found%idomain) - if (.not. found%length) then - write (errmsg, '(a)') 'Error in GRIDDATA block: LENGTH not found.' - call store_error(errmsg) - end if - if (.not. found%width) then write (errmsg, '(a)') 'Error in GRIDDATA block: WIDTH not found.' call store_error(errmsg) @@ -484,13 +472,11 @@ subroutine source_griddata(this) call store_error_filename(this%input_fname) end if - ! -- log simulation values + ! log simulation values if (this%iout > 0) then call this%log_griddata(found) end if - ! - ! -- Return - return + end subroutine source_griddata !> @brief Write griddata found to list file @@ -501,10 +487,6 @@ subroutine log_griddata(this, found) write (this%iout, '(1x,a)') 'Setting Discretization Griddata' - if (found%length) then - write (this%iout, '(4x,a)') 'LENGTH set from input file' - end if - if (found%width) then write (this%iout, '(4x,a)') 'WIDTH set from input file' end if @@ -565,10 +547,10 @@ subroutine source_vertices(this) return end subroutine source_vertices - !> @brief Copy cell2d information from input data context + !> @brief Copy cell1d information from input data context !! to model context !< - subroutine source_cell2d(this) + subroutine source_cell1d(this) ! -- modules use MemoryHelperModule, only: create_mem_path use MemoryManagerModule, only: mem_setptr @@ -578,7 +560,7 @@ subroutine source_cell2d(this) class(Disv1dType) :: this ! -- locals character(len=LENMEMPATH) :: idmMemoryPath - integer(I4B), dimension(:), contiguous, pointer :: icell2d => null() + integer(I4B), dimension(:), contiguous, pointer :: icell1d => null() integer(I4B), dimension(:), contiguous, pointer :: ncvert => null() integer(I4B), dimension(:), contiguous, pointer :: icvert => null() real(DP), dimension(:), contiguous, pointer :: fdc => null() @@ -589,14 +571,14 @@ subroutine source_cell2d(this) idmMemoryPath = create_mem_path(this%name_model, 'DISV1D', idm_context) ! ! -- set pointers to input path ncvert and icvert - call mem_setptr(icell2d, 'ICELL2D', idmMemoryPath) + call mem_setptr(icell1d, 'ICELL1D', idmMemoryPath) call mem_setptr(ncvert, 'NCVERT', idmMemoryPath) call mem_setptr(icvert, 'ICVERT', idmMemoryPath) ! ! -- - if (associated(icell2d) .and. associated(ncvert) & + if (associated(icell1d) .and. associated(ncvert) & .and. associated(icvert)) then - call this%define_cellverts(icell2d, ncvert, icvert) + call this%define_cellverts(icell1d, ncvert, icvert) else call store_error('Required cell vertex arrays not found.') end if @@ -609,32 +591,38 @@ subroutine source_cell2d(this) do i = 1, this%nodesuser this%fdc(i) = fdc(i) end do - call calculate_cellxy(this%vertices, this%fdc, this%iavert, & - this%javert, this%cellxy) else call store_error('Required fdc array not found.') end if - ! - ! -- log + + ! calculate length from vertices + call calculate_cell_length(this%vertices, this%iavert, this%javert, & + this%length) + + ! calculate cellxy from vertices and fdc + call calculate_cellxy(this%vertices, this%fdc, this%iavert, & + this%javert, this%length, this%cellxy) + + ! log if (this%iout > 0) then - write (this%iout, '(1x,a)') 'Setting Discretization CELL2D' - write (this%iout, '(1x,a,/)') 'End Setting Discretization CELL2D' + write (this%iout, '(1x,a)') 'Setting Discretization CELL1D' + write (this%iout, '(1x,a,/)') 'End Setting Discretization CELL1D' end if ! ! -- Return return - end subroutine source_cell2d + end subroutine source_cell1d !> @brief Construct the iavert and javert integer vectors which !! are compressed sparse row index arrays that relate the vertices !! to reaches !< - subroutine define_cellverts(this, icell2d, ncvert, icvert) + subroutine define_cellverts(this, icell1d, ncvert, icvert) ! -- modules use SparseModule, only: sparsematrix ! -- dummy class(Disv1dType) :: this - integer(I4B), dimension(:), contiguous, pointer, intent(in) :: icell2d + integer(I4B), dimension(:), contiguous, pointer, intent(in) :: icell1d integer(I4B), dimension(:), contiguous, pointer, intent(in) :: ncvert integer(I4B), dimension(:), contiguous, pointer, intent(in) :: icvert ! -- locals @@ -648,7 +636,7 @@ subroutine define_cellverts(this, icell2d, ncvert, icvert) ! -- add sparse matrix connections from input memory paths icv_idx = 1 do i = 1, this%nodesuser - if (icell2d(i) /= i) call store_error('ICELL2D input sequence violation.') + if (icell1d(i) /= i) call store_error('ICELL1D input sequence violation.') do j = 1, ncvert(i) call vert_spm%addconnection(i, icvert(icv_idx), 0) if (j == 1) then @@ -670,12 +658,13 @@ end subroutine define_cellverts !> @brief Calculate x, y, coordinates of reach midpoint !< - subroutine calculate_cellxy(vertices, fdc, iavert, javert, cellxy) + subroutine calculate_cellxy(vertices, fdc, iavert, javert, length, cellxy) ! -- dummy real(DP), dimension(:, :), intent(in) :: vertices !< 2d array of vertices with x, y as columns real(DP), dimension(:), intent(in) :: fdc !< fractional distance to reach midpoint (normally 0.5) integer(I4B), dimension(:), intent(in) :: iavert !< csr mapping of vertices to cell reaches integer(I4B), dimension(:), intent(in) :: javert !< csr mapping of vertices to cell reaches + real(DP), dimension(:), intent(in) :: length !< vector of cell lengths real(DP), dimension(:, :), intent(inout) :: cellxy !< 2d array of reach midpoint with x, y as columns ! -- local integer(I4B) :: nodes !< number of nodes @@ -684,7 +673,6 @@ subroutine calculate_cellxy(vertices, fdc, iavert, javert, cellxy) integer(I4B) :: iv0 !< index for line reach start integer(I4B) :: iv1 !< index for linen reach end integer(I4B) :: ixy !< x, y column index - real(DP) :: length !< reach length = sum of individual line reaches real(DP) :: fd0 !< fractional distance to start of this line reach real(DP) :: fd1 !< fractional distance to end of this line reach real(DP) :: fd !< fractional distance where midpoint (defined by fdc) is located @@ -693,20 +681,13 @@ subroutine calculate_cellxy(vertices, fdc, iavert, javert, cellxy) nodes = size(iavert) - 1 do n = 1, nodes - ! calculate length of this reach - length = DZERO - do j = iavert(n), iavert(n + 1) - 2 - length = length + & - calcdist(vertices, javert(j), javert(j + 1)) - end do - ! find vertices that span midpoint iv0 = 0 iv1 = 0 fd0 = DZERO do j = iavert(n), iavert(n + 1) - 2 d = calcdist(vertices, javert(j), javert(j + 1)) - fd1 = fd0 + d / length + fd1 = fd0 + d / length(n) ! if true, then we know the midpoint is some fractional distance (fd) ! from vertex j to vertex j + 1 @@ -720,7 +701,7 @@ subroutine calculate_cellxy(vertices, fdc, iavert, javert, cellxy) end do ! find x, y position of point on line - do ixy = 1, 3 + do ixy = 1, 2 cellxy(ixy, n) = (DONE - fd) * vertices(ixy, iv0) + & fd * vertices(ixy, iv1) end do @@ -728,6 +709,33 @@ subroutine calculate_cellxy(vertices, fdc, iavert, javert, cellxy) end do end subroutine calculate_cellxy + !> @brief Calculate x, y, coordinates of reach midpoint + !< + subroutine calculate_cell_length(vertices, iavert, javert, length) + ! -- dummy + real(DP), dimension(:, :), intent(in) :: vertices !< 2d array of vertices with x, y as columns + integer(I4B), dimension(:), intent(in) :: iavert !< csr mapping of vertices to cell reaches + integer(I4B), dimension(:), intent(in) :: javert !< csr mapping of vertices to cell reaches + real(DP), dimension(:), intent(inout) :: length !< 2d array of reach midpoint with x, y as columns + ! -- local + integer(I4B) :: nodes !< number of nodes + integer(I4B) :: n !< node index + integer(I4B) :: j !< vertex index + real(DP) :: dlen !< length + + nodes = size(iavert) - 1 + do n = 1, nodes + + ! calculate length of this reach + dlen = DZERO + do j = iavert(n), iavert(n + 1) - 2 + dlen = dlen + calcdist(vertices, javert(j), javert(j + 1)) + end do + length(n) = dlen + + end do + end subroutine calculate_cell_length + !> @brief Finalize grid construction !< subroutine grid_finalize(this) @@ -738,6 +746,11 @@ subroutine grid_finalize(this) class(Disv1dType) :: this ! local integer(I4B) :: node, noder, k + ! format + character(len=*), parameter :: fmtnr = & + "(/1x, 'The specified IDOMAIN results in a reduced number of cells.',& + &/1x, 'Number of user nodes: ',I0,& + &/1X, 'Number of nodes in solution: ', I0, //)" ! count active cells this%nodes = 0 @@ -750,12 +763,11 @@ subroutine grid_finalize(this) call store_error('Model does not have any active nodes. Make sure & &IDOMAIN has some values greater than zero.') call store_error_filename(this%input_fname) - call ustop() end if - if (count_errors() > 0) then - call store_error_filename(this%input_fname) - call ustop() + ! Write message if reduced grid + if (this%nodes < this%nodesuser) then + write (this%iout, fmtnr) this%nodesuser, this%nodes end if ! Array size is now known, so allocate @@ -772,8 +784,6 @@ subroutine grid_finalize(this) if (this%idomain(k) > 0) then this%nodereduced(node) = noder noder = noder + 1 - elseif (this%idomain(k) < 0) then - this%nodereduced(node) = -1 else this%nodereduced(node) = 0 end if @@ -794,15 +804,18 @@ subroutine grid_finalize(this) end do end if - ! Copy bottom into bot + ! Move bottom into bot and put length into disbase%area + ! and set x and y center coordinates do node = 1, this%nodesuser - this%bot(node) = this%bottom(node) + noder = node + if (this%nodes < this%nodesuser) noder = this%nodereduced(node) + if (noder <= 0) cycle + this%bot(noder) = this%bottom(node) + this%area(noder) = this%length(node) end do - ! Assign area in DisBaseType as length - do node = 1, this%nodesuser - this%area(node) = this%length(node) - end do + ! create connectivity using vertices and cell1d + call this%create_connections() ! -- Return return @@ -836,19 +849,19 @@ subroutine allocate_arrays(this) end subroutine allocate_arrays subroutine create_connections(this) - ! -- modules - ! -- dummy + ! modules + ! dummy class(Disv1dType) :: this - ! -- local + ! local integer(I4B) :: nrsize - ! - ! -- create and fill the connections object + + ! create and fill the connections object nrsize = 0 if (this%nodes < this%nodesuser) nrsize = this%nodes - ! - ! -- Allocate connections object + + ! Allocate connections object allocate (this%con) - ! + ! Build connectivity based on vertices call this%con%disv1dconnections_verts(this%name_model, this%nodes, & this%nodesuser, nrsize, this%nvert, & @@ -1117,7 +1130,7 @@ subroutine disv1d_da(this) call mem_deallocate(this%bottom) call mem_deallocate(this%idomain) ! - ! -- cdl hack for arrays for vertices and cell2d blocks + ! -- cdl hack for arrays for vertices and cell1d blocks if (deallocate_vertices) then call mem_deallocate(this%vertices) call mem_deallocate(this%fdc) diff --git a/src/Model/Discretization/Disv2d.f90 b/src/Model/Discretization/Disv2d.f90 index 4102ce881cd..3cfe6f26972 100644 --- a/src/Model/Discretization/Disv2d.f90 +++ b/src/Model/Discretization/Disv2d.f90 @@ -313,10 +313,10 @@ subroutine source_dimensions(this) ! ! -- Allocate vertices array call mem_allocate(this%vertices, 2, this%nvert, 'VERTICES', this%memoryPath) - call mem_allocate(this%cellxy, 2, this%nodes, 'CELLXY', this%memoryPath) + call mem_allocate(this%cellxy, 2, this%nodesuser, 'CELLXY', this%memoryPath) ! ! -- initialize all cells to be active (idomain = 1) - do j = 1, this%nodes + do j = 1, this%nodesuser this%idomain(j) = 1 end do ! @@ -332,7 +332,7 @@ subroutine log_dimensions(this, found) write (this%iout, '(1x,a)') 'Setting Discretization Dimensions' ! if (found%nodes) then - write (this%iout, '(4x,a,i0)') 'NODES = ', this%nodes + write (this%iout, '(4x,a,i0)') 'NODES = ', this%nodesuser end if ! if (found%nvert) then @@ -386,23 +386,24 @@ end subroutine log_griddata !> @brief Finalize grid (check properties, allocate arrays, compute connections) !< subroutine grid_finalize(this) - ! -- dummy + ! dummy class(Disv2dType) :: this - ! -- locals + ! locals integer(I4B) :: node, noder, j, ncell_count - ! -- formats + ! formats character(len=*), parameter :: fmtnr = & "(/1x, 'The specified IDOMAIN results in a reduced number of cells.',& &/1x, 'Number of user nodes: ',I0,& &/1X, 'Number of nodes in solution: ', I0, //)" - ! - ! -- count active cells + + ! count active cells and set nodes to that number ncell_count = 0 - do j = 1, this%nodes + do j = 1, this%nodesuser if (this%idomain(j) > 0) ncell_count = ncell_count + 1 end do - ! - ! -- Check to make sure nodes is a valid number + this%nodes = ncell_count + + ! Check to make sure nodes is a valid number if (ncell_count == 0) then call store_error('Model does not have any active nodes. & &Ensure IDOMAIN array has some values greater & @@ -410,21 +411,22 @@ subroutine grid_finalize(this) call store_error_filename(this%input_fname) end if - if (count_errors() > 0) then - call store_error_filename(this%input_fname) + ! Write message if reduced grid + if (this%nodes < this%nodesuser) then + write (this%iout, fmtnr) this%nodesuser, this%nodes end if - ! - ! -- Array size is now known, so allocate + + ! Array size is now known, so allocate call this%allocate_arrays() - ! - ! -- Fill the nodereduced array with the reduced nodenumber, or - ! a negative number to indicate it is a pass-through cell, or - ! a zero to indicate that the cell is excluded from the - ! solution. + + ! Fill the nodereduced array with the reduced nodenumber, or + ! a negative number to indicate it is a pass-through cell, or + ! a zero to indicate that the cell is excluded from the + ! solution. if (this%nodes < this%nodesuser) then node = 1 noder = 1 - do j = 1, this%nodes + do j = 1, this%nodesuser if (this%idomain(j) > 0) then this%nodereduced(node) = noder noder = noder + 1 @@ -434,12 +436,12 @@ subroutine grid_finalize(this) node = node + 1 end do end if - ! - ! -- allocate and fill nodeuser if a reduced grid + + ! allocate and fill nodeuser if a reduced grid if (this%nodes < this%nodesuser) then node = 1 noder = 1 - do j = 1, this%nodes + do j = 1, this%nodesuser if (this%idomain(j) > 0) then this%nodeuser(noder) = node noder = noder + 1 @@ -448,15 +450,10 @@ subroutine grid_finalize(this) end do end if - ! Copy bottom into bot - do node = 1, this%nodesuser - this%bot(node) = this%bottom(node) - end do - - ! -- Move bottom into bot - ! and set x and y center coordinates + ! Move bottom into bot + ! and set x and y center coordinates node = 0 - do j = 1, this%nodes + do j = 1, this%nodesuser node = node + 1 noder = node if (this%nodes < this%nodesuser) noder = this%nodereduced(node) @@ -465,10 +462,10 @@ subroutine grid_finalize(this) this%xc(noder) = this%cellxy(1, j) this%yc(noder) = this%cellxy(2, j) end do - ! - ! -- Build connections + + ! Build connections call this%connect() - ! + end subroutine grid_finalize !> @brief Load grid vertices from IDM into package @@ -576,7 +573,7 @@ subroutine source_cell2d(this) ! ! -- set cell centers if (associated(cell_x) .and. associated(cell_y)) then - do i = 1, this%nodes + do i = 1, this%nodesuser this%cellxy(1, i) = cell_x(i) this%cellxy(2, i) = cell_y(i) end do @@ -608,7 +605,7 @@ subroutine connect(this) narea_lt_zero = 0 ! ! -- Assign the cell area - do j = 1, this%nodes + do j = 1, this%nodesuser area = this%get_cell2d_area(j) noder = this%get_nodenumber(j, 0) if (noder > 0) this%area(noder) = area @@ -651,7 +648,7 @@ subroutine connect(this) if (this%nodes < this%nodesuser) nrsize = this%nodes allocate (this%con) call this%con%disvconnections(this%name_model, this%nodes, & - this%nodes, 1, nrsize, & + this%nodesuser, 1, nrsize, & this%nvert, this%vertices, this%iavert, & this%javert, this%cellxy, & this%bot, this%bot, & @@ -737,13 +734,13 @@ subroutine write_grb(this, icelltype) write (txt, '(3a, i0)') 'VERTICES ', 'DOUBLE ', 'NDIM 2 2 ', this%nvert txt(lentxt:lentxt) = new_line('a') write (iunit) txt - write (txt, '(3a, i0)') 'CELLX ', 'DOUBLE ', 'NDIM 1 ', this%nodes + write (txt, '(3a, i0)') 'CELLX ', 'DOUBLE ', 'NDIM 1 ', this%nodesuser txt(lentxt:lentxt) = new_line('a') write (iunit) txt - write (txt, '(3a, i0)') 'CELLY ', 'DOUBLE ', 'NDIM 1 ', this%nodes + write (txt, '(3a, i0)') 'CELLY ', 'DOUBLE ', 'NDIM 1 ', this%nodesuser txt(lentxt:lentxt) = new_line('a') write (iunit) txt - write (txt, '(3a, i0)') 'IAVERT ', 'INTEGER ', 'NDIM 1 ', this%nodes + 1 + write (txt, '(3a, i0)') 'IAVERT ', 'INTEGER ', 'NDIM 1 ', this%nodesuser + 1 txt(lentxt:lentxt) = new_line('a') write (iunit) txt write (txt, '(3a, i0)') 'JAVERT ', 'INTEGER ', 'NDIM 1 ', size(this%javert) @@ -773,8 +770,8 @@ subroutine write_grb(this, icelltype) write (iunit) this%angrot ! angrot write (iunit) this%bottom ! botm write (iunit) this%vertices ! vertices - write (iunit) (this%cellxy(1, i), i=1, this%nodes) ! cellx - write (iunit) (this%cellxy(2, i), i=1, this%nodes) ! celly + write (iunit) (this%cellxy(1, i), i=1, this%nodesuser) ! cellx + write (iunit) (this%cellxy(2, i), i=1, this%nodesuser) ! celly write (iunit) this%iavert ! iavert write (iunit) this%javert ! javert write (iunit) this%con%iausr ! iausr @@ -835,18 +832,18 @@ end subroutine nodeu_to_array !> @brief Get reduced node number from user node number !< function get_nodenumber_idx1(this, nodeu, icheck) result(nodenumber) - ! -- return + ! return integer(I4B) :: nodenumber - ! -- dummy + ! dummy class(Disv2dType), intent(in) :: this integer(I4B), intent(in) :: nodeu integer(I4B), intent(in) :: icheck - ! -- local - ! - ! -- check the node number if requested + ! local + + ! check the node number if requested if (icheck /= 0) then - ! - ! -- If within valid range, convert to reduced nodenumber + + ! If within valid range, convert to reduced nodenumber if (nodeu < 1 .or. nodeu > this%nodesuser) then nodenumber = 0 write (errmsg, '(a,i0,a,i0,a)') & @@ -861,7 +858,7 @@ function get_nodenumber_idx1(this, nodeu, icheck) result(nodenumber) nodenumber = nodeu if (this%nodes < this%nodesuser) nodenumber = this%nodereduced(nodeu) end if - ! + end function get_nodenumber_idx1 !> @brief Get normal vector components between the cell and a given neighbor @@ -986,13 +983,13 @@ end subroutine allocate_scalars !> @brief Allocate and initialize arrays !< subroutine allocate_arrays(this) - ! -- dummy + ! dummy class(Disv2dType) :: this - ! - ! -- Allocate arrays in DisBaseType (mshape, top, bot, area) + + ! Allocate arrays in DisBaseType (mshape, top, bot, area) call this%DisBaseType%allocate_arrays() ! - ! -- Allocate arrays for DisvType + ! Allocate arrays for DisvType if (this%nodes < this%nodesuser) then call mem_allocate(this%nodeuser, this%nodes, 'NODEUSER', this%memoryPath) call mem_allocate(this%nodereduced, this%nodesuser, 'NODEREDUCED', & @@ -1001,9 +998,10 @@ subroutine allocate_arrays(this) call mem_allocate(this%nodeuser, 1, 'NODEUSER', this%memoryPath) call mem_allocate(this%nodereduced, 1, 'NODEREDUCED', this%memoryPath) end if - ! -- Initialize - this%mshape(1) = this%nodes - ! + + ! Initialize + this%mshape(1) = this%nodesuser + end subroutine allocate_arrays !> @brief Get the signed area of the cell diff --git a/src/Model/ModelUtilities/Connections.f90 b/src/Model/ModelUtilities/Connections.f90 index 350f209531c..95b09d9ec1f 100644 --- a/src/Model/ModelUtilities/Connections.f90 +++ b/src/Model/ModelUtilities/Connections.f90 @@ -1,8 +1,8 @@ module ConnectionsModule use ArrayReadersModule, only: ReadArray - use KindModule, only: DP, I4B - use ConstantsModule, only: LENMODELNAME, LENMEMPATH, DHALF + use KindModule, only: DP, I4B, LGP + use ConstantsModule, only: LENMODELNAME, LENMEMPATH, DHALF, DONE use MessageModule, only: write_message use SimVariablesModule, only: errmsg use BlockParserModule, only: BlockParserType @@ -954,21 +954,20 @@ end subroutine disuconnections !> @brief Fill the connections object for a disv1d package from vertices !! - !! todo: Still need to handle hwva - !! todo: Only unreduced disv1d grids are allowed at the moment + !! Note that nothing is done for hwva !! !< subroutine disv1dconnections_verts(this, name_model, nodes, nodesuser, & nrsize, nvert, & vertices, iavert, javert, & - cellxyz, cellfdc, nodereduced, nodeuser, & + cellxy, cellfdc, nodereduced, nodeuser, & reach_length) - ! -- modules + ! modules use ConstantsModule, only: DHALF, DZERO, DTHREE, DTWO, DPI use SparseModule, only: sparsematrix use GeomUtilModule, only: get_node use MemoryManagerModule, only: mem_reallocate - ! -- dummy + ! dummy class(ConnectionsType) :: this character(len=*), intent(in) :: name_model integer(I4B), intent(in) :: nodes @@ -978,36 +977,35 @@ subroutine disv1dconnections_verts(this, name_model, nodes, nodesuser, & real(DP), dimension(3, nvert), intent(in) :: vertices integer(I4B), dimension(:), intent(in) :: iavert integer(I4B), dimension(:), intent(in) :: javert - real(DP), dimension(3, nodesuser), intent(in) :: cellxyz - !integer(I4B), dimension(2, nodesuser), intent(in) :: centerverts + real(DP), dimension(2, nodesuser), intent(in) :: cellxy real(DP), dimension(nodesuser), intent(in) :: cellfdc integer(I4B), dimension(:), intent(in) :: nodereduced integer(I4B), dimension(:), intent(in) :: nodeuser real(DP), dimension(:), intent(in) :: reach_length !< length of each reach - ! -- local + ! local integer(I4B), dimension(:), allocatable :: itemp integer(I4B), dimension(:), allocatable :: iavertcells integer(I4B), dimension(:), allocatable :: javertcells type(sparsematrix) :: sparse, vertcellspm - integer(I4B) :: n, m, i, j, ierror - ! - ! -- Allocate scalars + integer(I4B) :: i, j, ierror + + ! Allocate scalars call this%allocate_scalars(name_model) - ! - ! -- Set scalars + + ! Set scalars this%nodes = nodes this%ianglex = 1 - ! - ! -- Create a sparse matrix array with a row for each vertex. The columns - ! in the sparse matrix contains the cells that include that vertex. - ! This array will be used to determine cell connectivity. + + ! Create a sparse matrix array with a row for each vertex. The columns + ! in the sparse matrix contains the cells that include that vertex. + ! This array will be used to determine cell connectivity. allocate (itemp(nvert)) do i = 1, nvert itemp(i) = 4 end do - call vertcellspm%init(nvert, nodes, itemp) + call vertcellspm%init(nvert, nodesuser, itemp) deallocate (itemp) - do j = 1, nodes + do j = 1, nodesuser do i = iavert(j), iavert(j + 1) - 1 call vertcellspm%addconnection(javert(i), j, 1) end do @@ -1017,82 +1015,118 @@ subroutine disv1dconnections_verts(this, name_model, nodes, nodesuser, & allocate (javertcells(vertcellspm%nnz)) call vertcellspm%filliaja(iavertcells, javertcells, ierror) call vertcellspm%destroy() - ! - ! -- Call routine to build a sparse matrix of the connections - call vertexconnectl(this%nodes, nrsize, 6, nodes, sparse, & + + ! Call routine to build a sparse matrix of the connections + call vertexconnectl(this%nodes, nrsize, 6, nodesuser, sparse, & iavertcells, javertcells, nodereduced) - n = sparse%nnz - m = this%nodes this%nja = sparse%nnz this%njas = (this%nja - this%nodes) / 2 - ! - ! -- cleanup memory - deallocate (iavertcells) - deallocate (javertcells) - ! - ! -- Allocate index arrays of size nja and symmetric arrays + + ! Allocate index arrays of size nja and symmetric arrays call this%allocate_arrays() - ! - ! -- Fill the IA and JA arrays from sparse, then destroy sparse + + ! Fill the IA and JA arrays from sparse, then destroy sparse call sparse%sort() call sparse%filliaja(this%ia, this%ja, ierror) call sparse%destroy() - ! - ! -- fill the isym and jas arrays + + ! fill the isym and jas arrays call fillisym(this%nodes, this%nja, this%ia, this%ja, this%isym) call filljas(this%nodes, this%nja, this%ia, this%ja, this%isym, this%jas) - ! Fill disv1d symmetric arrays - ! todo: need to handle cell center shifted from center of reach + ! fill the ihc, cl1, and cl2 arrays call fill_disv1d_symarrays(this%ia, this%ja, this%jas, reach_length, & - this%ihc, this%cl1, this%cl2) + this%ihc, this%cl1, this%cl2, & + nrsize, nodereduced, nodeuser, cellfdc, & + iavert, javert, iavertcells, javertcells) - ! - ! -- Fill symmetric discretization arrays (ihc,cl1,cl2,hwva,anglex) - ! do n = 1, this%nodes - ! do ipos = this%ia(n) + 1, this%ia(n + 1) - 1 - ! m = this%ja(ipos) - ! if(m < n) cycle - ! call geol%cprops(n, m, this%hwva(this%jas(ipos)), & - ! this%cl1(this%jas(ipos)), this%cl2(this%jas(ipos))) - ! enddo - ! enddo + ! cleanup memory + deallocate (iavertcells) + deallocate (javertcells) - ! - ! -- If reduced system, then need to build iausr and jausr, otherwise point - ! them to ia and ja. + ! If reduced system, then need to build iausr and jausr, otherwise point + ! them to ia and ja. call this%iajausr(nrsize, nodesuser, nodereduced, nodeuser) - ! - ! -- Return - return + end subroutine disv1dconnections_verts !> @brief Fill symmetric connection arrays for disv1d !< - subroutine fill_disv1d_symarrays(ia, ja, jas, reach_length, ihc, cl1, cl2) + subroutine fill_disv1d_symarrays(ia, ja, jas, cell_length, ihc, cl1, cl2, & + nrsize, nodereduced, nodeuser, fdc, & + iavert, javert, iavertcells, javertcells) ! dummy integer(I4B), dimension(:), intent(in) :: ia !< csr pointer array integer(I4B), dimension(:), intent(in) :: ja !< csr array integer(I4B), dimension(:), intent(in) :: jas !< csr symmetric array - real(DP), dimension(:), intent(in) :: reach_length !< length of each reach + real(DP), dimension(:), intent(in) :: cell_length !< length of each cell (all cells) integer(I4B), dimension(:), intent(out) :: ihc !< horizontal connection flag real(DP), dimension(:), intent(out) :: cl1 !< distance from n to shared face with m real(DP), dimension(:), intent(out) :: cl2 !< distance from m to shared face with n + integer(I4B), intent(in) :: nrsize !< great than zero indicated reduced nodes present + integer(I4B), dimension(:), intent(in) :: nodereduced !< map user node to reduced node number + integer(I4B), dimension(:), intent(in) :: nodeuser !< map user reduced node to user node number + real(DP), dimension(:), intent(in) :: fdc !< fractional distance along cell to reach cell center + integer(I4B), dimension(:), intent(in) :: iavert ! csr index array of size (nodeuser + 1) for javert + integer(I4B), dimension(:), intent(in) :: javert ! csr array containing vertex numbers that define each cell + integer(I4B), dimension(:), intent(in) :: iavertcells ! csr index array of size (nvert + 1) for javert + integer(I4B), dimension(:), intent(in) :: javertcells ! csr array containing cells numbers that referenced for a vertex + ! local - integer(I4B) :: n - integer(I4B) :: m + integer(I4B) :: nr, nu + integer(I4B) :: mr, mu integer(I4B) :: ipos integer(I4B) :: isympos + real(DP) :: f ! loop through and set array values - do n = 1, size(reach_length) - do ipos = ia(n) + 1, ia(n + 1) - 1 - m = ja(ipos) - if (m < n) cycle + do nu = 1, size(cell_length) + + ! find reduced node number and cycle if cell does not exist + nr = nu + if (nrsize > 0) nr = nodereduced(nu) + if (nr <= 0) cycle + + ! visit each cell connected to reduced cell nr + do ipos = ia(nr) + 1, ia(nr + 1) - 1 + + ! process upper triangle by skipping mr cells less than nr + mr = ja(ipos) + if (mr < nr) cycle + + mu = mr + if (nrsize > 0) mu = nodeuser(mr) + isympos = jas(ipos) ihc(isympos) = 1 - cl1(isympos) = DHALF * reach_length(n) - cl2(isympos) = DHALF * reach_length(m) + + ! if cell m is connected to the downstream end of cell n, then use + ! 1 - fdc times the cell length, otherwise use fdc * length + if (fdc(nu) == DHALF) then + f = DHALF + else + if (connected_down_n(nu, mu, iavert, javert, iavertcells, & + javertcells)) then + f = (DONE - fdc(nu)) + else + f = fdc(nu) + end if + end if + cl1(isympos) = f * cell_length(nu) + + ! do the opposite for the cl2 distance as it is relative to m + if (fdc(mu) == DHALF) then + f = DHALF + else + if (connected_down_n(mu, nu, iavert, javert, iavertcells, & + javertcells)) then + f = (DONE - fdc(mu)) + else + f = fdc(mu) + end if + end if + cl2(isympos) = f * cell_length(mu) + end do end do end subroutine fill_disv1d_symarrays @@ -1370,43 +1404,42 @@ end subroutine vertexconnect subroutine vertexconnectl(nodes, nrsize, maxnnz, nodeuser, sparse, & iavertcells, javertcells, & nodereduced) - ! -- modules + ! modules use SparseModule, only: sparsematrix use GeomUtilModule, only: get_node - ! -- dummy - integer(I4B), intent(in) :: nodes - integer(I4B), intent(in) :: nrsize - integer(I4B), intent(in) :: maxnnz - integer(I4B), intent(in) :: nodeuser - type(SparseMatrix), intent(inout) :: sparse - integer(I4B), dimension(:), intent(in) :: nodereduced - integer(I4B), dimension(:), intent(in) :: iavertcells - integer(I4B), dimension(:), intent(in) :: javertcells - ! -- local + ! dummy + integer(I4B), intent(in) :: nodes !< number of active nodes + integer(I4B), intent(in) :: nrsize !< if > 0 then reduced grid + integer(I4B), intent(in) :: maxnnz !< max number of non zeros + integer(I4B), intent(in) :: nodeuser !< number of user nodes + type(SparseMatrix), intent(inout) :: sparse !< sparse matrix object + integer(I4B), dimension(:), intent(in) :: nodereduced !< map from user to reduced node + integer(I4B), dimension(:), intent(in) :: iavertcells !< csr ia index array for vertices + integer(I4B), dimension(:), intent(in) :: javertcells !< csr list of cells that use each vertex + ! local integer(I4B), dimension(:), allocatable :: rowmaxnnz integer(I4B) :: i, k, nr, mr, nvert integer(I4B) :: con - ! - ! -- Allocate and fill the ia and ja arrays + + ! Setup a sparse object allocate (rowmaxnnz(nodes)) do i = 1, nodes rowmaxnnz(i) = maxnnz end do call sparse%init(nodes, nodes, rowmaxnnz) deallocate (rowmaxnnz) - do nr = 1, nodes - ! - ! -- Process diagonal + + ! Fill the sparse diagonal + do nr = 1, nodeuser mr = nr if (nrsize > 0) mr = nodereduced(mr) if (mr <= 0) cycle call sparse%addconnection(mr, mr, 1) end do - ! - ! -- Go through each vertex and connect up all the cells that use - ! this vertex in their definition. - nvert = size(iavertcells) - 1 + ! Go through each vertex and connect up all the cells that use + ! this vertex in their definition. + nvert = size(iavertcells) - 1 do i = 1, nvert ! loop through cells that share the vertex do k = iavertcells(i), iavertcells(i + 1) - 2 @@ -1423,9 +1456,7 @@ subroutine vertexconnectl(nodes, nrsize, maxnnz, nodeuser, sparse, & end do end do end do - ! - ! -- return - return + end subroutine vertexconnectl !> @brief routine to set a value in the mask array (which has the same shape @@ -1485,4 +1516,37 @@ subroutine iac_to_ia(iac, ia) return end subroutine iac_to_ia + !> @brief Is cell m is connected to the downstream end of cell n + !< + function connected_down_n(nu, mu, iavert, javert, iavertcells, javertcells) & + result(connected_down) + ! dummy + integer(I4B), intent(in) :: nu ! user nodenumber for cell n + integer(I4B), intent(in) :: mu ! user nodenumber for connected cell m + integer(I4B), dimension(:), intent(in) :: iavert ! csr index array of size (nodeuser + 1) for javert + integer(I4B), dimension(:), intent(in) :: javert ! csr array containing vertex numbers that define each cell + integer(I4B), dimension(:), intent(in) :: iavertcells ! csr index array of size (nvert + 1) for javert + integer(I4B), dimension(:), intent(in) :: javertcells ! csr array containing cells numbers that referenced for a vertex + ! return + logical(LGP) :: connected_down + ! - local + integer(I4B) :: ipos + integer(I4B) :: ivert_down + + ! ivert_down is the last vertex for node nu; the last vertex is considered + ! to be the dowstream end of node nu + ivert_down = javert(iavert(nu + 1) - 1) + + ! look through vertex ivert_down, and see if mu is in it; if so, then + ! that means mu is connected to the downstream end of nu + connected_down = .false. + do ipos = iavertcells(ivert_down), iavertcells(ivert_down + 1) - 1 + if (javertcells(ipos) == mu) then + connected_down = .true. + exit + end if + end do + + end function connected_down_n + end module ConnectionsModule