|
| 1 | +import math |
| 2 | +import os |
| 3 | +import shutil |
| 4 | +from traceback import format_exc |
| 5 | +from warnings import warn |
| 6 | + |
| 7 | +import numpy as np |
| 8 | +import pytest |
| 9 | +from modflow_devtools.markers import requires_exe, requires_pkg |
| 10 | + |
| 11 | +import flopy |
| 12 | +from flopy.utils import import_optional_dependency |
| 13 | +from flopy.utils.datautil import DatumUtil |
| 14 | + |
| 15 | + |
| 16 | +class TestInfo: |
| 17 | + def __init__( |
| 18 | + self, |
| 19 | + original_simulation_folder, |
| 20 | + netcdf_simulation_folder, |
| 21 | + netcdf_output_file, |
| 22 | + ): |
| 23 | + self.original_simulation_folder = original_simulation_folder |
| 24 | + self.netcdf_simulation_folder = netcdf_simulation_folder |
| 25 | + self.netcdf_output_file = netcdf_output_file |
| 26 | + |
| 27 | + |
| 28 | +@requires_pkg("xarray") |
| 29 | +@requires_exe("mf6") |
| 30 | +@pytest.mark.regression |
| 31 | +def test_load_netcdf_gwfsto01(function_tmpdir, example_data_path): |
| 32 | + xr = import_optional_dependency("xarray") |
| 33 | + data_path_base = example_data_path / "mf6" / "netcdf" |
| 34 | + tests = { |
| 35 | + "test_gwf_sto01_mesh": TestInfo( |
| 36 | + "gwf_sto01", |
| 37 | + "gwf_sto01_write", |
| 38 | + "gwf_sto01.ugrid.nc", |
| 39 | + ), |
| 40 | + "test_gwf_sto01_structured": TestInfo( |
| 41 | + "gwf_sto01", |
| 42 | + "gwf_sto01_write", |
| 43 | + "gwf_sto01.structured.nc", |
| 44 | + ), |
| 45 | + } |
| 46 | + ws = function_tmpdir / "ws" |
| 47 | + for base_folder, test_info in tests.items(): |
| 48 | + print(f"RUNNING TEST: {base_folder}") |
| 49 | + data_path = os.path.join( |
| 50 | + data_path_base, base_folder, test_info.original_simulation_folder |
| 51 | + ) |
| 52 | + # copy example data into working directory |
| 53 | + base_model_folder = os.path.join(ws, f"{base_folder}_base") |
| 54 | + test_model_folder = os.path.join(ws, f"{base_folder}_test") |
| 55 | + shutil.copytree(data_path, base_model_folder) |
| 56 | + # load example |
| 57 | + sim = flopy.mf6.MFSimulation.load(sim_ws=base_model_folder) |
| 58 | + # change simulation path |
| 59 | + sim.set_sim_path(test_model_folder) |
| 60 | + # write example simulation to reset path |
| 61 | + sim.write_simulation() |
| 62 | + |
| 63 | + # compare generated files |
| 64 | + gen_files = [ |
| 65 | + f |
| 66 | + for f in os.listdir(test_model_folder) |
| 67 | + if os.path.isfile(os.path.join(test_model_folder, f)) |
| 68 | + ] |
| 69 | + base_files = [ |
| 70 | + f |
| 71 | + for f in os.listdir(base_model_folder) |
| 72 | + if os.path.isfile(os.path.join(base_model_folder, f)) |
| 73 | + ] |
| 74 | + # assert len(gen_files) == len(base_files) |
| 75 | + for f in base_files: |
| 76 | + print(f"cmp => {f}") |
| 77 | + base = os.path.join(base_model_folder, f) |
| 78 | + gen = os.path.join(test_model_folder, f) |
| 79 | + if f != test_info.netcdf_output_file: |
| 80 | + # "gwf_sto01.dis.ncf": # TODO wkt string missing on write? |
| 81 | + with open(base, "r") as file1, open(gen, "r") as file2: |
| 82 | + # Skip first line |
| 83 | + next(file1) |
| 84 | + next(file2) |
| 85 | + |
| 86 | + for line1, line2 in zip(file1, file2): |
| 87 | + assert line1 == line2 |
| 88 | + else: |
| 89 | + # TODO compare nc files |
| 90 | + assert os.path.exists(gen) |
| 91 | + continue |
| 92 | + |
| 93 | + |
| 94 | +@requires_pkg("xarray") |
| 95 | +@requires_exe("mf6") |
| 96 | +@pytest.mark.regression |
| 97 | +def test_create_netcdf_gwfsto01(function_tmpdir, example_data_path): |
| 98 | + xr = import_optional_dependency("xarray") |
| 99 | + |
| 100 | + cases = ["structured", "mesh2d"] |
| 101 | + |
| 102 | + # static model data |
| 103 | + # temporal discretization |
| 104 | + nper = 31 |
| 105 | + perlen = [1.0] + [365.2500000 for _ in range(nper - 1)] |
| 106 | + nstp = [1] + [6 for _ in range(nper - 1)] |
| 107 | + tsmult = [1.0] + [1.3 for _ in range(nper - 1)] |
| 108 | + tdis_rc = [] |
| 109 | + for i in range(nper): |
| 110 | + tdis_rc.append((perlen[i], nstp[i], tsmult[i])) |
| 111 | + |
| 112 | + # spatial discretization data |
| 113 | + nlay, nrow, ncol = 3, 10, 10 |
| 114 | + delr, delc = 1000.0, 2000.0 |
| 115 | + botm = [-100, -150.0, -350.0] |
| 116 | + strt = 0.0 |
| 117 | + |
| 118 | + # calculate hk |
| 119 | + hk1fact = 1.0 / 50.0 |
| 120 | + hk1 = np.ones((nrow, ncol), dtype=float) * 0.5 * hk1fact |
| 121 | + hk1[0, :] = 1000.0 * hk1fact |
| 122 | + hk1[-1, :] = 1000.0 * hk1fact |
| 123 | + hk1[:, 0] = 1000.0 * hk1fact |
| 124 | + hk1[:, -1] = 1000.0 * hk1fact |
| 125 | + hk = [20.0, hk1, 5.0] |
| 126 | + |
| 127 | + # calculate vka |
| 128 | + vka = [1e6, 7.5e-5, 1e6] |
| 129 | + |
| 130 | + # all cells are active and layer 1 is convertible |
| 131 | + ib = 1 |
| 132 | + laytyp = [1, 0, 0] |
| 133 | + |
| 134 | + # solver options |
| 135 | + nouter, ninner = 500, 300 |
| 136 | + hclose, rclose, relax = 1e-9, 1e-6, 1.0 |
| 137 | + newtonoptions = "NEWTON" |
| 138 | + imsla = "BICGSTAB" |
| 139 | + |
| 140 | + # chd data |
| 141 | + c = [] |
| 142 | + c6 = [] |
| 143 | + ccol = [3, 4, 5, 6] |
| 144 | + for j in ccol: |
| 145 | + c.append([0, nrow - 1, j, strt, strt]) |
| 146 | + c6.append([(0, nrow - 1, j), strt]) |
| 147 | + cd = {0: c} |
| 148 | + cd6 = {0: c6} |
| 149 | + maxchd = len(cd[0]) |
| 150 | + |
| 151 | + # pumping well data |
| 152 | + wr = [0, 0, 0, 0, 1, 1, 2, 2, 3, 3] |
| 153 | + wc = [0, 1, 8, 9, 0, 9, 0, 9, 0, 0] |
| 154 | + wrp = [2, 2, 3, 3] |
| 155 | + wcp = [5, 6, 5, 6] |
| 156 | + wq = [-14000.0, -8000.0, -5000.0, -3000.0] |
| 157 | + d = [] |
| 158 | + d6 = [] |
| 159 | + for r, c, q in zip(wrp, wcp, wq): |
| 160 | + d.append([2, r, c, q]) |
| 161 | + d6.append([(2, r, c), q]) |
| 162 | + wd = {1: d} |
| 163 | + wd6 = {1: d6} |
| 164 | + maxwel = len(wd[1]) |
| 165 | + |
| 166 | + # recharge data |
| 167 | + q = 3000.0 / (delr * delc) |
| 168 | + v = np.zeros((nrow, ncol), dtype=float) |
| 169 | + for r, c in zip(wr, wc): |
| 170 | + v[r, c] = q |
| 171 | + rech = {0: v} |
| 172 | + |
| 173 | + # storage and compaction data |
| 174 | + ske = [6e-4, 3e-4, 6e-4] |
| 175 | + |
| 176 | + # build |
| 177 | + for name in cases: |
| 178 | + # name = cases[0] |
| 179 | + ws = function_tmpdir / name |
| 180 | + |
| 181 | + # build MODFLOW 6 files |
| 182 | + sim = flopy.mf6.MFSimulation( |
| 183 | + sim_name=name, version="mf6", exe_name="mf6", sim_ws=ws |
| 184 | + ) |
| 185 | + # create tdis package |
| 186 | + tdis = flopy.mf6.ModflowTdis( |
| 187 | + sim, time_units="DAYS", nper=nper, perioddata=tdis_rc |
| 188 | + ) |
| 189 | + |
| 190 | + # create gwf model |
| 191 | + top = 0.0 |
| 192 | + zthick = [top - botm[0], botm[0] - botm[1], botm[1] - botm[2]] |
| 193 | + elevs = [top] + botm |
| 194 | + |
| 195 | + gwf = flopy.mf6.ModflowGwf( |
| 196 | + sim, modelname=name, newtonoptions=newtonoptions, save_flows=True |
| 197 | + ) |
| 198 | + |
| 199 | + # create iterative model solution and register the gwf model with it |
| 200 | + ims = flopy.mf6.ModflowIms( |
| 201 | + sim, |
| 202 | + print_option="SUMMARY", |
| 203 | + outer_dvclose=hclose, |
| 204 | + outer_maximum=nouter, |
| 205 | + under_relaxation="NONE", |
| 206 | + inner_maximum=ninner, |
| 207 | + inner_dvclose=hclose, |
| 208 | + rcloserecord=rclose, |
| 209 | + linear_acceleration=imsla, |
| 210 | + scaling_method="NONE", |
| 211 | + reordering_method="NONE", |
| 212 | + relaxation_factor=relax, |
| 213 | + ) |
| 214 | + sim.register_ims_package(ims, [gwf.name]) |
| 215 | + |
| 216 | + dis = flopy.mf6.ModflowGwfdis( |
| 217 | + gwf, |
| 218 | + nlay=nlay, |
| 219 | + nrow=nrow, |
| 220 | + ncol=ncol, |
| 221 | + delr=delr, |
| 222 | + delc=delc, |
| 223 | + top=top, |
| 224 | + botm=botm, |
| 225 | + filename=f"{name}.dis", |
| 226 | + ) |
| 227 | + |
| 228 | + # initial conditions |
| 229 | + ic = flopy.mf6.ModflowGwfic(gwf, strt=strt, filename=f"{name}.ic") |
| 230 | + |
| 231 | + # node property flow |
| 232 | + npf = flopy.mf6.ModflowGwfnpf( |
| 233 | + gwf, |
| 234 | + save_flows=False, |
| 235 | + icelltype=laytyp, |
| 236 | + k=hk, |
| 237 | + k33=vka, |
| 238 | + ) |
| 239 | + # storage |
| 240 | + sto = flopy.mf6.ModflowGwfsto( |
| 241 | + gwf, |
| 242 | + save_flows=False, |
| 243 | + iconvert=laytyp, |
| 244 | + ss=ske, |
| 245 | + sy=0, |
| 246 | + storagecoefficient=None, |
| 247 | + steady_state={0: True}, |
| 248 | + transient={1: True}, |
| 249 | + ) |
| 250 | + |
| 251 | + # recharge |
| 252 | + rch = flopy.mf6.ModflowGwfrcha(gwf, readasarrays=True, recharge=rech) |
| 253 | + |
| 254 | + # wel file |
| 255 | + wel = flopy.mf6.ModflowGwfwel( |
| 256 | + gwf, |
| 257 | + print_input=True, |
| 258 | + print_flows=True, |
| 259 | + maxbound=maxwel, |
| 260 | + stress_period_data=wd6, |
| 261 | + save_flows=False, |
| 262 | + ) |
| 263 | + |
| 264 | + # chd files |
| 265 | + chd = flopy.mf6.modflow.mfgwfchd.ModflowGwfchd( |
| 266 | + gwf, maxbound=maxchd, stress_period_data=cd6, save_flows=False |
| 267 | + ) |
| 268 | + |
| 269 | + # output control |
| 270 | + oc = flopy.mf6.ModflowGwfoc( |
| 271 | + gwf, |
| 272 | + budget_filerecord=f"{name}.cbc", |
| 273 | + head_filerecord=f"{name}.hds", |
| 274 | + headprintrecord=[("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL")], |
| 275 | + saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], |
| 276 | + printrecord=[("HEAD", "LAST"), ("BUDGET", "ALL")], |
| 277 | + ) |
| 278 | + |
| 279 | + sim.write_simulation(netcdf=name) |
| 280 | + |
| 281 | + try: |
| 282 | + success, buff = flopy.run_model( |
| 283 | + "mf6", |
| 284 | + ws / "mfsim.nam", |
| 285 | + model_ws=ws, |
| 286 | + report=True, |
| 287 | + ) |
| 288 | + except Exception: |
| 289 | + warn( |
| 290 | + "MODFLOW 6 serial test", |
| 291 | + name, |
| 292 | + f"failed with error:\n{format_exc()}", |
| 293 | + ) |
| 294 | + success = False |
| 295 | + |
| 296 | + assert success |
0 commit comments