Named values for enumerated case parameters and AST-based analytic IC codegen (#1550)

Author: sbryngelson

docs/documentation/case.md

@@ -98,6 +98,11 @@ There are multiple sets of parameters that must be specified in the python input
 Items 8, 9, 10, 11 and 12 are optional sets of parameters that activate the acoustic source model, ensemble-averaged bubble model, initial velocity field setup, phase change, artificial Mach number respectively.
 Definition of the parameters is described in the following subsections.
 
+Enumerated parameters accept named values as well as integer codes: `"riemann_solver": "hllc"`
+is equivalent to `"riemann_solver": 2`. Defined names appear in each parameter's table entry
+(e.g. `1` (`hll`), `2` (`hllc`)). Existing case files can be rewritten to named syntax with
+`./mfc.sh validate case.py --migrate`.
+
 ### 1. Runtime {#sec-runtime}
 
 | Parameter        | Type           | Description                               |
@@ -213,24 +218,34 @@ Some parameters, as described above, can be defined by analytical functions in t
 
 where `alpha_rho` is defined with the `1 + 0.1*sin(20*x*pi)` function.
 
-The expressions are interpreted as Fortran code, in the execution context and scope of the function that defines the patch.
-Additionally, the following variables are made available as shorthand:
+Expressions use **Python syntax** and are parsed at case load time.
+Syntax errors and unknown variable or function names are immediate, named errors — they are reported before any Fortran is compiled, with a message identifying the offending expression and listing the available names.
+
+The following variables are available in IC patch expressions:
 
 | Shorthand | Expands To               | Shorthand | Expands To                | Shorthand | Expands To               |
 | --------- | ------------------------ | --------- | ------------------------- | --------- | ------------------------ |
 | `x`       | `x_cc(i)`                | `lx`      | The patch's `length_x`    | `xc`      | The patch's `x_centroid` |
 | `y`       | `y_cc(j)`                | `ly`      | The patch's `length_y`    | `yc`      | The patch's `y_centroid` |
 | `z`       | `z_cc(k)`                | `lz`      | The patch's `length_z`    | `zc`      | The patch's `z_centroid` |
 | `eps`     | The patch's `epsilon`    | `beta`    | The patch's `beta`        | `radii`   | The patch's `radii`      |
-| `tau_e`   | The patch's `tau_e`      | `r`       | The patch's `radius`      | `pi` and `e` | \f$\pi\f$ and \f$e\f$ |
+| `tau_e`   | The patch's `tau_e`      | `r`       | The patch's `radius`      | `pi`      | \f$\pi\f$ (Fortran constant from `m_constants`) |
 
 where $(i,j,k)$ are the grid-indices of the current cell in each coordinate direction.
 
+The allowed functions are the standard Fortran intrinsics:
+`sin`, `cos`, `tan`, `asin`, `acos`, `atan`, `atan2`,
+`sinh`, `cosh`, `tanh`, `exp`, `log`, `log10`, `sqrt`,
+`abs`, `min`, `max`, `mod`, `sign`.
+
+**Euler's number:** bare `e` is **not** a variable.
+Write `exp(1.0)` or a numeric literal instead (e.g. `2.718281828`).
+
 In the example above, the following code is generated:
 
 ```f90
 if (patch_id == 2) then
-    q_prim_vf(eqn_idx%cont%beg)%sf(i, 0, 0) = 1 + 0.1*sin(20*x_cc(i)*3.141592653589793)
+    q_prim_vf(eqn_idx%cont%beg)%sf(i, 0, 0) = 1 + 0.1 * sin(20 * x_cc(i) * pi)
 end if
 ```
 
@@ -355,9 +370,13 @@ Additional details on this specification can be found in [NACA airfoil](https://
 
 - `moving_ibm` sets the method by which movement will be applied to the immersed boundary. Using 0 will result in no movement. Using 1 will result 1-way coupling where the boundary moves at a constant rate and applied forces to the fluid based upon it's own motion. In 1-way coupling, the fluid does not apply forces back onto the IB. Using 2 will result in 2-way coupling, where the boundary pushes on the fluid and the fluid pushes back on the boundary via pressure and viscous forces. If external forces are applied, the boundary will also experience those forces.
 
-- `vel(i)` is the initial linear velocity of the IB in the x, y, z direction for i=1, 2, 3. When `moving_ibm` equals 2, this velocity is just the starting speed of the object, which will then accelerate due to external forces. If `moving_ibm` equals 1, then this is constant if it is a number, or can be described analytically with an expression. 
+- `vel(i)` is the initial linear velocity of the IB in the x, y, z direction for i=1, 2, 3. When `moving_ibm` equals 2, this velocity is just the starting speed of the object, which will then accelerate due to external forces. If `moving_ibm` equals 1, then this is constant if it is a number, or can be described analytically with an expression.
+
+- `angular_vel(i)` is the initial angular velocity of the IB about the x, y, z axes for i=1, 2, 3 in radians per second. When `moving_ibm` equals 2, this rotation rate is just the starting rate of the object, which will then change due to external torques. If `moving_ibm` equals 1, then this is constant if it is a number, or can be described analytically with an expression.
 
-- `angular_vel(i)` is the initial angular velocity of the IB about the x, y, z axes for i=1, 2, 3 in radians per second. When `moving_ibm` equals 2, this rotation rate is just the starting rate of the object, which will then change due to external torques. If `moving_ibm` equals 1, then this is constant if it is a number, or can be described analytically with an expression. 
+  Moving-IB analytic expressions use the same Python syntax and error-reporting as IC patch expressions (see the "Analytical Definition of Primitive Variables" section above).
+  Available variables: `x` (`x_cc(i)`), `y` (`y_cc(j)`), `z` (`z_cc(k)`), `t` (current simulation time), and `r` (the IB patch radius).
+  The same intrinsic functions and `pi` constant apply; bare `e` is not available.
 
 - `coefficient_of_restitution` is a number from 0 (exclusive) to 1 (inclusive) describing how elastic IB collisions are. 0 is for perfectly inellastic collisions while 1 is for perfectly ellastic collisions.
 

examples/0D_bubblecollapse_adap/case.py

@@ -69,25 +69,25 @@
             "t_step_save": t_save,
             # Simulation Algorithm Parameters
             "num_patches": 1,
-            "model_eqns": 2,
+            "model_eqns": "5eq",
             "alt_soundspeed": "F",
             "num_fluids": 1,
             "mpp_lim": "F",
             "mixture_err": "F",
-            "time_stepper": 3,
+            "time_stepper": "rk3",
             "weno_order": 5,
             "weno_eps": 1.0e-16,
             "mapped_weno": "T",
             "null_weights": "F",
             "mp_weno": "T",
-            "riemann_solver": 2,
-            "wave_speeds": 1,
-            "avg_state": 2,
+            "riemann_solver": "hllc",
+            "wave_speeds": "direct",
+            "avg_state": "arithmetic",
             "bc_x%beg": -1,
             "bc_x%end": -1,
             # Formatted Database Files Structure Parameters
-            "format": 1,
-            "precision": 2,
+            "format": "silo",
+            "precision": "double",
             "prim_vars_wrt": "T",
             "parallel_io": "T",
             "fd_order": 1,
@@ -103,7 +103,7 @@
             "patch_icpp(1)%v0": 0.0,
             # Bubbles
             "bubbles_euler": "T",
-            "bubble_model": 2,
+            "bubble_model": "keller_miksis",
             # Nondimensional numbers
             "bub_pp%R0ref": R0ref / x0,
             "bub_pp%p0ref": p0ref / p0,

examples/1D_acoustic_dipole/case.py

@@ -19,26 +19,26 @@
             "t_step_save": 1,
             # Simulation Algorithm Parameters
             "num_patches": 1,
-            "model_eqns": 2,
+            "model_eqns": "5eq",
             "alt_soundspeed": "F",
             "num_fluids": 1,
             "mpp_lim": "F",
             "mixture_err": "F",
-            "time_stepper": 3,
+            "time_stepper": "rk3",
             "weno_order": 5,
             "weno_eps": 1.0e-16,
             "teno": "T",
             "teno_CT": 1e-8,
             "null_weights": "F",
             "mp_weno": "F",
-            "riemann_solver": 2,
-            "wave_speeds": 1,
-            "avg_state": 2,
+            "riemann_solver": "hllc",
+            "wave_speeds": "direct",
+            "avg_state": "arithmetic",
             "bc_x%beg": -6,
             "bc_x%end": -6,
             # Formatted Database Files Structure Parameters
-            "format": 1,
-            "precision": 2,
+            "format": "silo",
+            "precision": "double",
             "prim_vars_wrt": "T",
             "rho_wrt": "T",
             "parallel_io": "T",

examples/1D_acoustic_gauss_sigmadist/case.py

@@ -19,26 +19,26 @@
             "t_step_save": 1,
             # Simulation Algorithm Parameters
             "num_patches": 1,
-            "model_eqns": 2,
+            "model_eqns": "5eq",
             "alt_soundspeed": "F",
             "num_fluids": 1,
             "mpp_lim": "F",
             "mixture_err": "F",
-            "time_stepper": 3,
+            "time_stepper": "rk3",
             "weno_order": 5,
             "weno_eps": 1.0e-16,
             "teno": "T",
             "teno_CT": 1e-8,
             "null_weights": "F",
             "mp_weno": "F",
-            "riemann_solver": 2,
-            "wave_speeds": 1,
-            "avg_state": 2,
+            "riemann_solver": "hllc",
+            "wave_speeds": "direct",
+            "avg_state": "arithmetic",
             "bc_x%beg": -6,
             "bc_x%end": -6,
             # Formatted Database Files Structure Parameters
-            "format": 1,
-            "precision": 2,
+            "format": "silo",
+            "precision": "double",
             "prim_vars_wrt": "T",
             "rho_wrt": "T",
             "parallel_io": "T",

examples/1D_acoustic_gauss_sigmatime/case.py

@@ -19,26 +19,26 @@
             "t_step_save": 1,
             # Simulation Algorithm Parameters
             "num_patches": 1,
-            "model_eqns": 2,
+            "model_eqns": "5eq",
             "alt_soundspeed": "F",
             "num_fluids": 1,
             "mpp_lim": "F",
             "mixture_err": "F",
-            "time_stepper": 3,
+            "time_stepper": "rk3",
             "weno_order": 5,
             "weno_eps": 1.0e-16,
             "teno": "T",
             "teno_CT": 1e-8,
             "null_weights": "F",
             "mp_weno": "F",
-            "riemann_solver": 2,
-            "wave_speeds": 1,
-            "avg_state": 2,
+            "riemann_solver": "hllc",
+            "wave_speeds": "direct",
+            "avg_state": "arithmetic",
             "bc_x%beg": -6,
             "bc_x%end": -6,
             # Formatted Database Files Structure Parameters
-            "format": 1,
-            "precision": 2,
+            "format": "silo",
+            "precision": "double",
             "prim_vars_wrt": "T",
             "rho_wrt": "T",
             "parallel_io": "T",

examples/1D_acoustic_sine_frequency/case.py

@@ -19,26 +19,26 @@
             "t_step_save": 1,
             # Simulation Algorithm Parameters
             "num_patches": 1,
-            "model_eqns": 2,
+            "model_eqns": "5eq",
             "alt_soundspeed": "F",
             "num_fluids": 1,
             "mpp_lim": "F",
             "mixture_err": "F",
-            "time_stepper": 3,
+            "time_stepper": "rk3",
             "weno_order": 5,
             "weno_eps": 1.0e-16,
             "teno": "T",
             "teno_CT": 1e-8,
             "null_weights": "F",
             "mp_weno": "F",
-            "riemann_solver": 2,
-            "wave_speeds": 1,
-            "avg_state": 2,
+            "riemann_solver": "hllc",
+            "wave_speeds": "direct",
+            "avg_state": "arithmetic",
             "bc_x%beg": -6,
             "bc_x%end": -6,
             # Formatted Database Files Structure Parameters
-            "format": 1,
-            "precision": 2,
+            "format": "silo",
+            "precision": "double",
             "prim_vars_wrt": "T",
             "rho_wrt": "T",
             "parallel_io": "T",

examples/1D_acoustic_sine_wavelength/case.py

@@ -19,26 +19,26 @@
             "t_step_save": 5,
             # Simulation Algorithm Parameters
             "num_patches": 1,
-            "model_eqns": 2,
+            "model_eqns": "5eq",
             "alt_soundspeed": "F",
             "num_fluids": 1,
             "mpp_lim": "F",
             "mixture_err": "F",
-            "time_stepper": 3,
+            "time_stepper": "rk3",
             "weno_order": 5,
             "weno_eps": 1.0e-16,
             "teno": "T",
             "teno_CT": 1e-8,
             "null_weights": "F",
             "mp_weno": "F",
-            "riemann_solver": 2,
-            "wave_speeds": 1,
-            "avg_state": 2,
+            "riemann_solver": "hllc",
+            "wave_speeds": "direct",
+            "avg_state": "arithmetic",
             "bc_x%beg": -6,
             "bc_x%end": -6,
             # Formatted Database Files Structure Parameters
-            "format": 1,
-            "precision": 2,
+            "format": "silo",
+            "precision": "double",
             "prim_vars_wrt": "T",
             "rho_wrt": "T",
             "parallel_io": "T",

examples/1D_advection_convergence/case.py

@@ -39,13 +39,13 @@
 
 if args.muscl:
     scheme_params = {
-        "recon_type": 2,
-        "muscl_order": 2,
+        "recon_type": "muscl",
+        "muscl_order": "second_order",
         "muscl_lim": args.muscl_lim,
     }
 else:
     scheme_params = {
-        "recon_type": 1,
+        "recon_type": "weno",
         "weno_order": args.order,
         "weno_eps": 1.0e-16,
         "weno_Re_flux": "F",
@@ -69,19 +69,19 @@
             "t_step_stop": Nt,
             "t_step_save": Nt,
             "num_patches": 1,
-            "model_eqns": 2,
+            "model_eqns": "5eq",
             "alt_soundspeed": "F",
             "num_fluids": 2,
             "mpp_lim": "F",
             "mixture_err": "F",
-            "time_stepper": 3,
-            "riemann_solver": 2,
-            "wave_speeds": 1,
-            "avg_state": 2,
+            "time_stepper": "rk3",
+            "riemann_solver": "hllc",
+            "wave_speeds": "direct",
+            "avg_state": "arithmetic",
             "bc_x%beg": -1,
             "bc_x%end": -1,
-            "format": 1,
-            "precision": 2,
+            "format": "silo",
+            "precision": "double",
             "prim_vars_wrt": "T",
             "parallel_io": "F",
             "patch_icpp(1)%geometry": 1,

examples/1D_brio_wu/case.py

@@ -22,25 +22,25 @@
             "t_step_save": 100,
             # Simulation Algorithm Parameters
             "num_patches": 2,
-            "model_eqns": 2,
+            "model_eqns": "5eq",
             "alt_soundspeed": "F",
             "num_fluids": 1,
             "mpp_lim": "F",
             "mixture_err": "F",
-            "time_stepper": 3,
+            "time_stepper": "rk3",
             "weno_order": 5,
             "mapped_weno": "T",
             "weno_eps": 1.0e-16,
             "null_weights": "F",
             "mp_weno": "F",
-            "riemann_solver": 1,
-            "wave_speeds": 1,
-            "avg_state": 2,
+            "riemann_solver": "hll",
+            "wave_speeds": "direct",
+            "avg_state": "arithmetic",
             "bc_x%beg": -3,
             "bc_x%end": -3,
             # Formatted Database Files Structure Parameters
-            "format": 1,
-            "precision": 2,
+            "format": "silo",
+            "precision": "double",
             "prim_vars_wrt": "T",
             "rho_wrt": "T",
             "parallel_io": "T",

examples/1D_brio_wu_hlld/case.py

@@ -23,25 +23,25 @@
             "t_step_save": 100,
             # Simulation Algorithm Parameters
             "num_patches": 2,
-            "model_eqns": 2,
+            "model_eqns": "5eq",
             "alt_soundspeed": "F",
             "num_fluids": 1,
             "mpp_lim": "F",
             "mixture_err": "F",
-            "time_stepper": 3,
+            "time_stepper": "rk3",
             "weno_order": 5,
             "mapped_weno": "T",
             "weno_eps": 1.0e-16,
             "null_weights": "F",
             "mp_weno": "F",
-            "riemann_solver": 4,
-            "wave_speeds": 1,
-            "avg_state": 2,
+            "riemann_solver": "hlld",
+            "wave_speeds": "direct",
+            "avg_state": "arithmetic",
             "bc_x%beg": -3,
             "bc_x%end": -3,
             # Formatted Database Files Structure Parameters
-            "format": 1,
-            "precision": 2,
+            "format": "silo",
+            "precision": "double",
             "prim_vars_wrt": "T",
             "rho_wrt": "T",
             "parallel_io": "T",