feat(netcdf): support for a structured export and gridded input arrays (#1934)

* support for a structured export and gridded input arrays

* rebuild makefiles

* fix msvs build

* add netcdf files to make exclude list

* additional netcdf input error handling

* restore meson build

* some cleanup

* few additional doxygen comments

* add xarray import to tests

* run ruff linter

* run ruff formatter

---------

Co-authored-by: mjreno <[email protected]>

Author: mjreno

autotest/test_netcdf_gwe_cnd.py

@@ -7,6 +7,7 @@
 # Imports
 
 import os
+import subprocess
 
 import numpy as np
 import pytest
@@ -20,31 +21,116 @@
     raise Exception(msg)
 
 try:
+    import xarray as xa
     import xugrid as xu
 except ImportError:
-    pytest.skip("xuarray not found", allow_module_level=True)
+    pytest.skip("xarray and xugrid not found", allow_module_level=True)
 
 from framework import TestFramework
+from test_gwe_cnd import cases
 
-cases = ["cnd01"]
 
-
-def build_models(idx, test):
+def build_models(idx, test, export, gridded_input):
     from test_gwe_cnd import build_models as build
 
     sim, dummy = build(idx, test)
     sim.tdis.start_date_time = "2041-01-01T00:00:00-05:00"
     gwe = sim.gwe[0]
-    gwe.name_file.export_netcdf = "ugrid"
+    gwe.name_file.export_netcdf = export
     gwe.dis.export_array_netcdf = True
     gwe.ic.export_array_netcdf = True
     gwe.cnd.export_array_netcdf = True
+
+    name = "gwe-" + cases[idx]
+
+    # netcdf config
+    ncf = flopy.mf6.ModflowUtlncf(
+        gwe.dis,
+        filename=f"{name}.dis.ncf",
+    )
+
     return sim, dummy
 
 
-def check_output(idx, test):
+def check_output(idx, test, export, gridded_input):
     from test_gwe_cnd import check_output as check
 
+    name = "gwe-" + test.name
+
+    if gridded_input == "netcdf":
+        # re-run the simulation with model netcdf input
+        input_fname = f"{name}.nc"
+        nc_fname = f"{name}.{export}.nc"
+        subprocess.run(
+            ["mv", test.workspace / input_fname, test.workspace / nc_fname]
+        )
+
+        with open(test.workspace / f"{name}.nam", "w") as f:
+            f.write("BEGIN options\n")
+            f.write("  SAVE_FLOWS\n")
+            f.write(f"  EXPORT_NETCDF {export}\n")
+            f.write(f"  NETCDF  FILEIN {name}.{export}.nc\n")
+            f.write("END options\n\n")
+            f.write("BEGIN packages\n")
+            f.write(f"  DIS6  {name}.dis  dis\n")
+            f.write(f"  IC6  {name}.ic  ic\n")
+            f.write(f"  ADV6  {name}.adv  adv\n")
+            f.write(f"  CND6  {name}.cnd  cnd\n")
+            f.write(f"  EST6  {name}.est  est\n")
+            f.write(f"  CTP6  {name}.ctp  ctp-1\n")
+            f.write(f"  SSM6  {name}.ssm  ssm\n")
+            f.write(f"  OC6  {name}.oc  oc\n")
+            f.write("END packages\n")
+
+        with open(test.workspace / f"{name}.dis", "w") as f:
+            f.write("BEGIN options\n")
+            f.write("  NOGRB\n")
+            f.write("  EXPORT_ARRAY_NETCDF\n")
+            f.write(f"  NCF6  FILEIN  {name}.dis.ncf\n")
+            f.write("END options\n\n")
+            f.write("BEGIN dimensions\n")
+            f.write("  NLAY  1\n")
+            f.write("  NROW  1\n")
+            f.write("  NCOL  101\n")
+            f.write("END dimensions\n\n")
+            f.write("BEGIN griddata\n")
+            f.write("  delr NETCDF\n")
+            f.write("  delc NETCDF\n")
+            f.write("  top NETCDF\n")
+            f.write("  botm NETCDF\n")
+            f.write("  idomain NETCDF\n")
+            f.write("END griddata\n\n")
+
+        with open(test.workspace / f"{name}.ic", "w") as f:
+            f.write("BEGIN options\n")
+            f.write("  EXPORT_ARRAY_NETCDF\n")
+            f.write("END options\n\n")
+            f.write("BEGIN griddata\n")
+            f.write("  strt NETCDF\n")
+            f.write("END griddata\n")
+
+        with open(test.workspace / f"{name}.cnd", "w") as f:
+            f.write("BEGIN options\n")
+            f.write("  XT3D_OFF\n")
+            f.write("  EXPORT_ARRAY_NETCDF\n")
+            f.write("END options\n\n")
+            f.write("BEGIN griddata\n")
+            f.write("  alh  NETCDF\n")
+            f.write("  ath1  NETCDF\n")
+            f.write("  ktw  NETCDF\n")
+            f.write("  kts  NETCDF\n")
+            f.write("END griddata\n")
+
+        success, buff = flopy.run_model(
+            test.targets["mf6"],
+            test.workspace / "mfsim.nam",
+            model_ws=test.workspace,
+            report=True,
+        )
+
+        assert success
+        test.success = success
+
     check(idx, test)
 
     # read transport results from GWE model
@@ -63,8 +149,11 @@ def check_output(idx, test):
 
     # Check NetCDF output
     nc_fpth = os.path.join(test.workspace, f"{gwename}.nc")
-    ds = xu.open_dataset(nc_fpth)
-    xds = ds.ugrid.to_dataset()
+    if export == "ugrid":
+        ds = xu.open_dataset(nc_fpth)
+        xds = ds.ugrid.to_dataset()
+    elif export == "structured":
+        xds = xa.open_dataset(nc_fpth)
 
     # Compare NetCDF head arrays with binary headfile temperatures
     gwe = test.sims[0].gwe[0]
@@ -78,12 +167,20 @@ def check_output(idx, test):
     for i in range(nper):
         for j in range(int(pd[i][1])):
             rec = cobj.get_data(kstpkper=(j, i))
-            for l in range(nlay):
+            if export == "ugrid":
+                for l in range(nlay):
+                    assert np.allclose(
+                        np.array(rec[l]).flatten(),
+                        xds[f"temperature_l{l+1}"][timestep, :].data,
+                    ), f"NetCDF-temperature comparison failure in timestep {timestep+1}"
+                timestep += 1
+            elif export == "structured":
                 assert np.allclose(
-                    np.array(rec[l]).flatten(),
-                    xds[f"temperature_l{l+1}"][timestep, :].data,
+                    # np.array(rec).flatten(),
+                    np.array(rec),
+                    xds["temperature"][timestep, :].data,
                 ), f"NetCDF-temperature comparison failure in timestep {timestep+1}"
-            timestep += 1
+                timestep += 1
 
     vlist = [
         "dis_delr",
@@ -97,15 +194,6 @@ def check_output(idx, test):
         "cnd_ktw_l",
         "cnd_kts_l",
     ]
-    layer_vlist = [
-        "dis_botm_l",
-        "dis_idomain_l",
-        "ic_strt_l",
-        "cnd_alh_l",
-        "cnd_ath1_l",
-        "cnd_ktw_l",
-        "cnd_kts_l",
-    ]
 
     # Compare NetCDF package input arrays with FloPy arrays
     gwe = test.sims[0].gwe[0]
@@ -115,14 +203,22 @@ def check_output(idx, test):
         array_name = tokens[1].split("_")[0]
         package = getattr(gwe, package_name)
         b = getattr(package, array_name).array
-        if var in layer_vlist:
-            for l in range(nlay):
+        if export == "ugrid":
+            if var.endswith("_l"):
+                for l in range(nlay):
+                    assert np.allclose(
+                        np.array(b[l]).flatten(), xds[f"{var}{l+1}"].data
+                    ), f"NetCDF input array comparison failure, variable={var}{l+1}"
+            else:
                 assert np.allclose(
-                    np.array(b[l]).flatten(), xds[f"{var}{l+1}"].data
-                ), f"NetCDF input array comparison failure, variable={var}{l+1}"
-        else:
+                    np.array(b).flatten(), xds[var].data
+                ), f"NetCDF input array comparison failure, variable={var}"
+        elif export == "structured":
+            var = var.replace("_l", "")
             assert np.allclose(
-                np.array(b).flatten(), xds[var].data
+                # np.array(b).flatten(), xds[var].data
+                np.array(b),
+                xds[var].data,
             ), f"NetCDF input array comparison failure, variable={var}"
 
 
@@ -131,12 +227,14 @@ def check_output(idx, test):
     "idx, name",
     list(enumerate(cases)),
 )
-def test_mf6model(idx, name, function_tmpdir, targets):
+@pytest.mark.parametrize("export", ["ugrid", "structured"])
+@pytest.mark.parametrize("gridded_input", ["ascii", "netcdf"])
+def test_mf6model(idx, name, function_tmpdir, targets, export, gridded_input):
     test = TestFramework(
         name=name,
         workspace=function_tmpdir,
-        build=lambda t: build_models(idx, t),
-        check=lambda t: check_output(idx, t),
+        build=lambda t: build_models(idx, t, export, gridded_input),
+        check=lambda t: check_output(idx, t, export, gridded_input),
         targets=targets,
     )
     test.run()