Use AST module to help post-process assignments in scripts

notebooks/ex-gwf-drn-p01.ipynb

@@ -1393,7 +1393,7 @@
     "        sfr = flopy.mf6.ModflowGwfsfr(\n",
     "            gwf,\n",
     "            pname=\"SFR-1\",\n",
-    "            unit_conversion=1.486,\n",
+    "            length_conversion=3.28081,\n",
     "            mover=True,\n",
     "            nreaches=len(sfr_pakdata),\n",
     "            packagedata=sfr_pakdata,\n",

notebooks/ex-gwf-lak-p02.ipynb

@@ -422,16 +422,6 @@
     "SFR package"
    ]
   },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "e827da47",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "sfr_conv = 128390.4"
-   ]
-  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -1024,7 +1014,8 @@
     "        flopy.mf6.ModflowGwfsfr(\n",
     "            gwf,\n",
     "            pname=\"SFR-1\",\n",
-    "            unit_conversion=sfr_conv,\n",
+    "            time_conversion=86400.000,\n",
+    "            length_conversion=3.28081,\n",
     "            mover=True,\n",
     "            print_stage=True,\n",
     "            print_flows=True,\n",

notebooks/ex-gwf-lgr.ipynb

@@ -918,7 +918,7 @@
     "            save_flows=True,\n",
     "            mover=True,\n",
     "            pname=\"SFR-parent\",\n",
-    "            unit_conversion=86400.0,\n",
+    "            time_conversion=86400.0,\n",
     "            boundnames=False,\n",
     "            nreaches=len(connsp),\n",
     "            packagedata=pkdat,\n",
@@ -1028,7 +1028,7 @@
     "            save_flows=True,\n",
     "            mover=True,\n",
     "            pname=\"SFR-child\",\n",
-    "            unit_conversion=86400.00,\n",
+    "            time_conversion=86400.00,\n",
     "            boundnames=False,\n",
     "            nreaches=len(connsc),\n",
     "            packagedata=pkdatc,\n",

notebooks/ex-gwf-sagehen.ipynb

@@ -763,7 +763,7 @@
     "            save_flows=True,\n",
     "            mover=True,\n",
     "            pname=\"SFR-1\",\n",
-    "            unit_conversion=86400.0,\n",
+    "            time_conversion=86400.0,\n",
     "            boundnames=True,\n",
     "            nreaches=len(conns),\n",
     "            packagedata=pkdat,\n",

notebooks/ex-gwf-sfr-p01.ipynb

@@ -1102,7 +1102,7 @@
     "        )\n",
     "        sfr = flopy.mf6.ModflowGwfsfr(\n",
     "            gwf,\n",
-    "            unit_conversion=1.486,\n",
+    "            length_conversion=3.28081,\n",
     "            nreaches=len(sfr_pakdata),\n",
     "            packagedata=sfr_pakdata,\n",
     "            connectiondata=sfr_conn,\n",

notebooks/ex-gwf-sfr-p01b.ipynb

@@ -4191,7 +4191,7 @@
     "            mover=True,\n",
     "            pname=\"SFR-1\",\n",
     "            maximum_depth_change=0.1e-05,\n",
-    "            unit_conversion=1.486,\n",
+    "            length_conversion=3.28081,\n",
     "            nreaches=len(sfr_pakdata),\n",
     "            packagedata=sfr_pakdata,\n",
     "            connectiondata=sfr_conn,\n",

notebooks/ex-gwt-mt3dms-p01.ipynb

@@ -672,25 +672,29 @@
     "        # Instantiating MODFLOW 6 transport mass storage package (formerly \"reaction\" package in MT3DMS)\n",
     "        if retardation != 1.0:\n",
     "            sorption = \"linear\"\n",
+    "            bulk_density = rhob\n",
     "            kd = (\n",
     "                (retardation - 1.0) * prsity / rhob\n",
     "            )  # prsity & rhob defined in\n",
     "        else:  # global variable section\n",
     "            sorption = None\n",
-    "            kd = 1.0\n",
+    "            bulk_density = None\n",
+    "            kd = None\n",
     "        if decay != 0.0:\n",
     "            first_order_decay = True\n",
+    "            decay_arg = decay\n",
     "        else:\n",
     "            first_order_decay = False\n",
+    "            decay_arg = None\n",
     "        flopy.mf6.ModflowGwtmst(\n",
     "            gwt,\n",
     "            porosity=prsity,\n",
     "            sorption=sorption,\n",
-    "            bulk_density=rhob,\n",
+    "            bulk_density=bulk_density,\n",
     "            distcoef=kd,\n",
     "            first_order_decay=first_order_decay,\n",
-    "            decay=decay,\n",
-    "            decay_sorbed=decay,\n",
+    "            decay=decay_arg,\n",
+    "            decay_sorbed=decay_arg,\n",
     "            filename=\"{}.mst\".format(gwtname),\n",
     "        )\n",
     "\n",

notebooks/ex-gwt-mt3dms-p02.ipynb

@@ -188,10 +188,10 @@
     "parameter_units = {\n",
     "    \"beta\": \"$s^{-1}$\",\n",
     "    \"sorption\": \"text string\",\n",
-    "    \"Kf\": \"$\\mu g L g^{} mg^{-1}$\",\n",
+    "    \"Kf\": r\"$\\mu g L g^{} mg^{-1}$\",\n",
     "    \"a\": \"unitless\",\n",
     "    \"Kl\": \"$L mg^{-1}$\",\n",
-    "    \"S\": \"$\\mu g g^{-1}$\",\n",
+    "    \"S\": r\"$\\mu g g^{-1}$\",\n",
     "}"
    ]
   },
@@ -419,10 +419,15 @@
     "        sp2 = a\n",
     "    if S is not None:\n",
     "        sp2 = S\n",
+    "    volfracim = 0.\n",
+    "    if beta is not None:\n",
+    "        if beta > 0:\n",
+    "            volfracim = bulk_density / (bulk_density + porosity)\n",
+    "\n",
     "    flopy.mf6.ModflowGwtmst(\n",
     "        gwt,\n",
-    "        porosity=porosity,\n",
-    "        bulk_density=bulk_density,\n",
+    "        porosity=porosity / (1.0 - volfracim),\n",
+    "        bulk_density=bulk_density if sorption else None,\n",
     "        sorption=sorption,\n",
     "        distcoef=distcoef,\n",
     "        sp2=sp2,\n",
@@ -433,8 +438,11 @@
     "    )\n",
     "    if beta is not None:\n",
     "        if beta > 0:\n",
-    "            thetaim = bulk_density\n",
-    "            flopy.mf6.ModflowGwtist(gwt, thetaim=thetaim, zetaim=beta)\n",
+    "            porosity_im = bulk_density / volfracim\n",
+    "            flopy.mf6.ModflowGwtist(gwt,\n",
+    "                                    volfrac=volfracim,\n",
+    "                                    porosity=porosity_im,\n",
+    "                                    zetaim=beta)\n",
     "    pd = [\n",
     "        (\"GWFHEAD\", \"../mf6gwf/flow.hds\".format(), None),\n",
     "        (\"GWFBUDGET\", \"../mf6gwf/flow.bud\", None),\n",

notebooks/ex-gwt-mt3dsupp632.ipynb

@@ -216,6 +216,7 @@
     "botm = 0  # Layer bottom elevation ($m$)\n",
     "specific_discharge = 0.06  # Specific discharge ($md^{-1}$)\n",
     "longitudinal_dispersivity = 10  # Longitudinal dispersivity ($m$)\n",
+    "volfrac = 0.2  # volume fraction that is immobile domain (unitless)\n",
     "porosity = 0.2  # Porosity of mobile domain (unitless)\n",
     "porosity_immobile = 0.05  # Porosity of immobile domain (unitless)\n",
     "bulk_density = 4.0  # Bulk density ($gL^{-1})$\n",
@@ -355,19 +356,23 @@
     "        first_order_decay = True\n",
     "    if distribution_coefficient > 0:\n",
     "        sorption = \"linear\"\n",
+    "        bd = bulk_density\n",
+    "        kd = distribution_coefficient\n",
     "    else:\n",
     "        sorption = None\n",
+    "        bd = None\n",
+    "        kd = None\n",
     "    flopy.mf6.ModflowGwtic(gwt, strt=0)\n",
     "    flopy.mf6.ModflowGwtmst(\n",
     "        gwt,\n",
     "        zero_order_decay=zero_order_decay,\n",
     "        first_order_decay=first_order_decay,\n",
     "        sorption=sorption,\n",
-    "        porosity=porosity,\n",
+    "        porosity=porosity / (1.0 - volfrac),\n",
     "        decay=decay,\n",
     "        decay_sorbed=decay_sorbed,\n",
-    "        bulk_density=bulk_density,\n",
-    "        distcoef=distribution_coefficient,\n",
+    "        bulk_density=bd,\n",
+    "        distcoef=kd,\n",
     "    )\n",
     "    istsorption = sorption is not None\n",
     "    if dual_domain:\n",
@@ -376,11 +381,12 @@
     "            zero_order_decay=zero_order_decay,\n",
     "            first_order_decay=first_order_decay,\n",
     "            sorption=istsorption,\n",
-    "            thetaim=porosity_immobile,\n",
+    "            porosity=porosity_immobile / volfrac,\n",
+    "            volfrac=volfrac,\n",
     "            zetaim=zeta_im,\n",
     "            decay=decay,\n",
     "            decay_sorbed=decay_sorbed,\n",
-    "            bulk_density=bulk_density,\n",
+    "            bulk_density=bd,\n",
     "            distcoef=distribution_coefficient,\n",
     "        )\n",
     "    flopy.mf6.ModflowGwtadv(gwt)\n",

notebooks/ex-gwt-prudic2004t2.ipynb

@@ -459,7 +459,8 @@
     "        budget_filerecord=name + \".sfr.bud\",\n",
     "        mover=True,\n",
     "        pname=\"SFR-1\",\n",
-    "        unit_conversion=128390.00,\n",
+    "        time_conversion=86400.000,\n",
+    "        length_conversion=3.28081,\n",
     "        boundnames=True,\n",
     "        nreaches=len(sfrconnectiondata),\n",
     "        packagedata=sfrpackagedata,\n",

scripts/ex-gwt-mt3dms-p02.py

@@ -69,10 +69,10 @@
 parameter_units = {
     "beta": "$s^{-1}$",
     "sorption": "text string",
-    "Kf": "$\mu g L g^{} mg^{-1}$",
+    "Kf": r"$\mu g L g^{} mg^{-1}$",
     "a": "unitless",
     "Kl": "$L mg^{-1}$",
-    "S": "$\mu g g^{-1}$",
+    "S": r"$\mu g g^{-1}$",
 }
 
 # Model units

scripts/process-scripts.py

@@ -1,3 +1,4 @@
+import ast
 import os
 import sys
 import re
@@ -66,6 +67,7 @@ def make_notebooks():
             "jupytext",
             "--from py",
             "--to ipynb",
+            "--update",
             "-o",
             opth,
             tpth,
@@ -126,39 +128,33 @@ def make_tables():
         )
         # do a little processing
         with open(file) as f:
-            lines = f.read().splitlines()
+            txt = f.read()
+        doc = ast.parse(txt, file)
+        lines = txt.splitlines()
+
+        # gather all assignments, organize by (starting) lineno and name
+        assign_bylineno = {}
+        assign_byname = {}
+        for node in doc.body:
+            if isinstance(node, ast.Assign):
+                target = node.targets[0]
+                if isinstance(target, ast.Name):
+                    assign_bylineno[node.lineno] = node
+                    assign_byname[target.id] = node.value
 
         # parse optional parameters dictionary
         parameters = None
-        tag = "parameters ="
-        for idx, line in enumerate(lines):
-            if line.lower().startswith(tag):
-                dict_str = line[len(tag) + 1 :].strip()
-                for jdx in range(idx + 1, len(lines)):
-                    if len(lines[jdx].strip()) < 1:
-                        break
-                    dict_str += " " + lines[jdx].strip()
-
-                # parse the dictionary string into a dictionary
-                parameters = eval(dict_str)
-                break
+        if "parameters" in assign_byname:
+            obj = assign_byname["parameters"]
+            parameters = eval(ast.unparse(obj))
 
         # parse optional parameter units dictionary, if parameters are
         # specified in the script
         if parameters is not None:
             parameter_units = None
-            tag = "parameter_units ="
-            for idx, line in enumerate(lines):
-                if line.lower().startswith(tag):
-                    dict_str = line[len(tag) + 1 :].strip()
-                    for jdx in range(idx + 1, len(lines)):
-                        if len(lines[jdx].strip()) < 1:
-                            break
-                        dict_str += " " + lines[jdx].strip()
-
-                    # parse the dictionary string into a dictionary
-                    parameter_units = eval(dict_str)
-                    break
+            if "parameter_units" in assign_byname:
+                obj = assign_byname["parameter_units"]
+                parameter_units = eval(ast.unparse(obj))
 
         # create scenario table if parameters are specified in the script
         if parameters is not None:
@@ -208,25 +204,38 @@ def make_tables():
 
         # tables
         table_number = 0
-        for idx, line in enumerate(lines):
-            table_text = []
-            table_value = []
+        scanning_table = False
+        table_text = []
+        table_value = []
+        for lineno, line in enumerate(lines, 1):
             if line.lower().startswith("# table"):
-                for table_line in lines[idx + 1 :]:
-                    # skip empty lines
-                    if len(table_line.strip()) < 1:
-                        continue
-                    if table_line.startswith("# "):
+                scanning_table = True
+                continue
+            if scanning_table:
+                if line.startswith("# "):
+                    scanning_table = False
+                if lineno in assign_bylineno:
+                    obj = assign_bylineno[lineno]
+                    # get comment from at least one line
+                    for idx in range(obj.lineno - 1, obj.end_lineno):
+                        table_line = lines[idx]
+                        if (ipos := table_line.find("# ")) > 0:
+                            table_text.append(table_line[ipos + 1 :].strip())
+                            break
+                    else:
+                        scanning_table = False
                         break
-                    ipos = table_line.find("# ")
-                    if ipos > 0:
-                        table_text.append(table_line[ipos + 1 :].strip())
-                        table_value.append(
-                            _replace_quotes(
-                                table_line[0:ipos].split("=")[1].strip()
-                            )
-                        )
-            if len(table_text) > 0:
+                    # get content on right side of "="
+                    val = obj.value
+                    if val.lineno == val.end_lineno:
+                        value = line[val.col_offset:val.end_col_offset]
+                    else:
+                        # string may span more than one line
+                        value = obj.value.value
+                    if value.startswith("'") or value.startswith('"'):
+                        value = _replace_quotes(value)
+                    table_value.append(value)
+            if not scanning_table and len(table_text) > 0:
                 table_number += 1
                 tab_name = "{}-{:02d}".format(basename, table_number)
                 caption = "Model parameters for example {}.".format(basename)
@@ -250,6 +259,10 @@ def make_tables():
                 # finalize table
                 f.close()
 
+                # reset these if there is more than one table
+                table_text = []
+                table_value = []
+
 
 def get_ordered_examples():
     print("creating a ordered list of examples from the LaTeX document")