Information Geometric Regularization

CMakeLists.txt

@@ -497,7 +497,7 @@ function(MFC_SETUP_TARGET)
                 endforeach()
 
                 target_compile_options(${a_target}
-                    PRIVATE -gpu=keep,ptxinfo,lineinfo
+                    PRIVATE -gpu=keep,ptxinfo,lineinfo,fastmath
                 )
 
                 # GH-200 Unified Memory Support

benchmarks/igr/case.py

@@ -0,0 +1,113 @@
+#!/usr/bin/env python3
+# Benchmark igr_T_viscous_T
+# Additional Benchmarked Features
+# - igr : T
+# - viscous : T
+# - igr_order : 5
+
+import json, math, argparse
+
+parser = argparse.ArgumentParser(prog="Benchmarking Case 5", description="This MFC case was created for the purposes of benchmarking MFC.", formatter_class=argparse.ArgumentDefaultsHelpFormatter)
+
+parser.add_argument("--mfc", type=json.loads, default="{}", metavar="DICT", help="MFC's toolchain's internal state.")
+parser.add_argument("--gbpp", type=int, metavar="MEM", default=16, help="Adjusts the problem size per rank to fit into [MEM] GB of GPU memory per GPU.")
+
+ARGS = vars(parser.parse_args())
+DICT = ARGS["mfc"]
+
+size = 1 if DICT["gpu"] else 0
+
+ppg = 8000000 / 16.0
+procs = DICT["nodes"] * DICT["tasks_per_node"]
+ncells = math.floor(ppg * procs * ARGS["gbpp"])
+s = math.floor((ncells) ** (1 / 3))
+Nx, Ny, Nz = s, s, s
+
+Re = 1600
+L = 1
+P0 = 101325
+rho0 = 1
+C0 = math.sqrt(1.4 * P0)
+V0 = 0.1 * C0
+mu = V0 * L / Re
+
+cfl = 0.5
+dx = 2 * math.pi * L / (Nx + 1)
+
+dt = cfl * dx / (C0)
+
+tC = L / V0
+tEnd = 20 * tC
+
+Nt = int(tEnd / dt)
+
+# Configuring case dictionary
+print(
+    json.dumps(
+        {
+            # Logistics
+            "run_time_info": "T",
+            # Computational Domain Parameters
+            "x_domain%beg": -math.pi * L,
+            "x_domain%end": math.pi * L,
+            "y_domain%beg": -math.pi * L,
+            "y_domain%end": math.pi * L,
+            "z_domain%beg": -math.pi * L,
+            "z_domain%end": math.pi * L,
+            "m": Nx,
+            "n": Ny,
+            "p": Nz,
+            "cyl_coord": "F",
+            "dt": dt,
+            "t_step_start": 0,
+            "t_step_stop": int(20 * (95 * size + 5)),
+            "t_step_save": int(20 * (95 * size + 5)),
+            # Simulation Algorithm Parameters
+            "num_patches": 1,
+            "model_eqns": 2,
+            "num_fluids": 1,
+            "time_stepper": 3,
+            "bc_x%beg": -1,
+            "bc_x%end": -1,
+            "bc_y%beg": -1,
+            "bc_y%end": -1,
+            "bc_z%beg": -1,
+            "bc_z%end": -1,
+            "igr": "T",
+            "igr_order": 5,
+            "igr_iter_solver": 1,
+            "num_igr_iters": 3,
+            "num_igr_warm_start_iters": 3,
+            "alf_factor": 10,
+            "viscous": "T",
+            # Formatted Database Files Structure Parameters
+            "format": 1,
+            "precision": 2,
+            "prim_vars_wrt": "T",
+            "omega_wrt(1)": "T",
+            "omega_wrt(2)": "T",
+            "omega_wrt(3)": "T",
+            "qm_wrt": "T",
+            "fd_order": 4,
+            "parallel_io": "T",
+            # Patch 1: Background (AIR - 2)
+            "patch_icpp(1)%geometry": 9,
+            "patch_icpp(1)%x_centroid": 0,
+            "patch_icpp(1)%y_centroid": 0,
+            "patch_icpp(1)%z_centroid": 0,
+            "patch_icpp(1)%length_x": 2 * math.pi * L,
+            "patch_icpp(1)%length_y": 2 * math.pi * L,
+            "patch_icpp(1)%length_z": 2 * math.pi * L,
+            "patch_icpp(1)%vel(1)": f"{V0}*sin(x/{L})*cos(y/{L})*sin(z/{L})",
+            "patch_icpp(1)%vel(2)": f"-{V0}*cos(x/{L})*sin(y/{L})*sin(z/{L})",
+            "patch_icpp(1)%vel(3)": 0,
+            "patch_icpp(1)%pres": f"{P0} + ({rho0}*{V0}**2/16)*(cos(2*x/{L}) + cos(2*y/{L}))*(cos(2*z/{L}) + 2)",
+            "patch_icpp(1)%alpha_rho(1)": 1,
+            "patch_icpp(1)%alpha(1)": 1,
+            # Fluids Physical Parameters
+            "fluid_pp(1)%gamma": 1.0e00 / (1.4 - 1),
+            "fluid_pp(1)%pi_inf": 0,
+            "fluid_pp(1)%Re(1)": 1 / mu,
+        }
+    )
+)

docs/documentation/case.md

@@ -408,6 +408,13 @@ Details of implementation of viscosity in MFC can be found in [Coralic (2015)](r
 | `surface_tension`      | Logical | Activate surface tension |
 | `viscous`              | Logical | Activate viscosity |
 | `hypoelasticity`       | Logical | Activate hypoelasticity* |
+| `igr`                  | Logical | Enable solution via information geometric regularization (IGR) [Cao (2024)](references.md) |
+| `igr_order`            | Integer | Order of reconstruction for IGR [3,5] |
+| `alf_factor`           | Real    | Alpha factor for IGR entropic pressure (default 10) |
+| `igr_pres_lim`         | Logical | Limit IGR pressure to avoid negative values (default F) |
+| `igr_iter_solver`      | Integer | Solution method for IGR elliptic solve [1] Jacobi [2] Gauss-Seidel |
+| `num_igr_iters`        | Integer | Number of iterations for for the IGR elliptic solve (default 2) |
+| `num_igr_warm_start_iters` | Integer | Number of iterations for the IGR elliptic solve at the first time step (default 50) |
 
 - \* Options that work only with `model_eqns = 2`.
 - † Options that work only with ``cyl_coord = 'F'``.
@@ -501,7 +508,7 @@ This option requires `weno_Re_flux` to be true because cell boundary values are
 | `type`*                | Integer | The geometry of the patch. [1]: Line [2]: Circle [3]: Rectangle |
 | `x[y,z]_centroid`*     | Real    | Centroid of the boundary patch in the x[y,z]-direction          |
 | `length_x[y,z]`*       | Real    | Length of the boundary patch in the x[y,z]-direction            |
-| `radius`*              | Real    | Radius of the boundary patch                                    |
+| `radiue`*              | Real    | Radius of the boundary patch                                    |
 *: These parameters should be prepended with `patch_bc(j)%` where $j$ is the patch index.
 
 Boundary condition patches can be used with the following boundary condition types:

docs/documentation/references.md

@@ -10,6 +10,8 @@
 
 - <a id="Bryngelson19">Bryngelson, S. H., Schmidmayer, K., Coralic, V., Meng, J. C., Maeda, K., and Colonius, T. (2019). Mfc: An open-source high-order multi-component, multi-phase, and multi-scale compressible flow solver. arXiv preprint arXiv:1907.10512.</a>
 
+- <a id="Cao24">Cao A. and. Sch&auml;fer F. (2024). Information Geometric Regularization of the Barotropic Euler Equations. arXiv preprint arXiv:2308.14127.</a>
+
 - <a id="Chen22">Chen, S. S., Li, J. P., Li, Z., Yuan, W., & Gao, Z. H. (2022). Anti-dissipation pressure correction under low Mach numbers for Godunov-type schemes. Journal of Computational Physics, 456, 111027. </a>
 
 - <a id="Childs12">Childs, H., Brugger, E., Whitlock, B., Meredith, J., Ahern, S., Pugmire, D., Biagas, K., Miller, M., Harrison, C., Weber, G. H., Krishnan, H., Fogal, T., Sanderson, A., Garth, C., Bethel, E. W., Camp, D., R¨ubel, O., Durant, M., Favre, J. M., and Navr´atil, P. (2012). VisIt: An End-User Tool For Visualizing and Analyzing Very Large Data. In High Performance Visualization–Enabling Extreme-Scale Scientific Insight, pages 357–372.</a>
@@ -62,4 +64,4 @@
 
 - <a id="Powell94">Powell, K. G. (1994). An approximate Riemann solver for magnetohydrodynamics: (That works in more than one dimension). In Upwind and high-resolution schemes (pp. 570-583). Springer.</a>
 
-- <a id="Cao19">Cao, S., Zhang, Y., Liao, D., Zhong, P., and Wang, K. G. (2019). Shock-induced damage and dynamic fracture in cylindrical bodies submerged in liquid. International Journal of Solids and Structures, 169:55–71. Elsevier.</a>
+- <a id="Cao19">Cao, S., Zhang, Y., Liao, D., Zhong, P., and Wang, K. G. (2019). Shock-induced damage and dynamic fracture in cylindrical bodies submerged in liquid. International Journal of Solids and Structures, 169:55–71. Elsevier.</a>

examples/2D_IGR_2fluid/case.py

@@ -0,0 +1,105 @@
+
+#!/usr/bin/env python3
+# This case file demonstrates the Laplace pressure jump of a water droplet in air. The laplace pressure jump
+# in 2D is given by delta = sigma / r where delta is the pressure jump, sigma is the surface tension coefficient,
+# and r is the radius of the droplet. The results of this simulation agree with theory to well within 1%
+# relative error.
+
+import math
+import json
+
+l = 1
+eps = 1e-6
+
+# Numerical setup
+r0 = l / 4
+x0 = -l
+x1 = l
+y0 = -l
+y1 = l
+
+Nx = 200
+Ny = 200
+
+dt = 1e-4
+
+# Configuration case dictionary
+data = {
+    # Logistics
+    "run_time_info": "T",
+    # Computational Domain
+    "x_domain%beg": x0,
+    "x_domain%end": x1,
+    "y_domain%beg": y0,
+    "y_domain%end": y1,
+    "m": Nx,
+    "n": Ny,
+    "p": 0,
+    "cyl_coord": "F",
+    "dt": dt,
+    "t_step_start": 0,
+    "t_step_stop": 5000,
+    "t_step_save": 50,
+    # Simulation Algorithm
+    "model_eqns": 2,
+    "alt_soundspeed": "F",
+    "mixture_err": "T",
+    "mpp_lim": "F",
+    "time_stepper": 3,
+    "bc_x%beg": -1,
+    "bc_x%end": -1,
+    "bc_y%beg": -1,
+    "bc_y%end": -1,
+    "num_patches": 2,
+    "num_fluids": 2,
+    "elliptic_smoothing": "T",
+    "elliptic_smoothing_iters": 10,
+    "igr": "T",
+    "igr_order": 5,
+    "igr_iter_solver": 1,
+    "num_igr_iters": 3,
+    "num_igr_warm_start_iters": 30,
+    "alf_factor": 10,
+    "viscous": "T",
+    # Database Structure Parameters
+    "format": 1,
+    "precision": 2,
+    "prim_vars_wrt": "T",
+    "parallel_io": "T",
+    # Fluid Parameters (Gas)
+    "fluid_pp(1)%gamma": 1.0e00 / (1.4e00 - 1.0e00),
+    "fluid_pp(1)%pi_inf": 0.0e00,
+    "fluid_pp(1)%Re(1)": 1e5,
+    "fluid_pp(2)%gamma": 1.0e00 / (1.4e00 - 1.0e00),
+    "fluid_pp(2)%pi_inf": 0.0e00,
+    "fluid_pp(2)%Re(1)": 1e5,
+    # Ambient pressure
+    "patch_icpp(1)%geometry": 3,
+    "patch_icpp(1)%x_centroid": 0,
+    "patch_icpp(1)%y_centroid": 0,
+    "patch_icpp(1)%length_x": 2,
+    "patch_icpp(1)%length_y": 2,
+    "patch_icpp(1)%vel(1)": 0.0,
+    "patch_icpp(1)%vel(2)": 0.0,
+    "patch_icpp(1)%vel(3)": 0.0,
+    "patch_icpp(1)%pres": 1,
+    "patch_icpp(1)%alpha_rho(1)": 1 - eps,
+    "patch_icpp(1)%alpha(1)": 1 - eps,
+    "patch_icpp(1)%alpha_rho(2)": eps,
+    "patch_icpp(1)%alpha(2)": eps,
+    # High pressure
+    "patch_icpp(2)%alter_patch(1)": "T",
+    "patch_icpp(2)%geometry": 2,
+    "patch_icpp(2)%x_centroid": 0,
+    "patch_icpp(2)%y_centroid": 0,
+    "patch_icpp(2)%radius": r0,
+    "patch_icpp(2)%vel(1)": 0.0,
+    "patch_icpp(2)%vel(2)": 0.0,
+    "patch_icpp(2)%pres": 9.518,
+    "patch_icpp(2)%alpha_rho(1)": eps,
+    "patch_icpp(2)%alpha(1)": eps,
+    "patch_icpp(2)%alpha_rho(2)": (1 - eps) * 5,
+    "patch_icpp(2)%alpha(2)": 1 - eps,
+}
+
+print(json.dumps(data))

examples/2D_IGR_triple_point/case.py

@@ -0,0 +1,102 @@
+#!/usr/bin/env python3
+import math
+import json
+
+eps = 1e-8
+Nx = 699
+Ny = 299
+
+print(
+    json.dumps(
+        {
+            # Logistics
+            "run_time_info": "F",
+            # Computational Domain Parameters
+            "x_domain%beg": 0,
+            "x_domain%end": 7,
+            "y_domain%beg": 0,
+            "y_domain%end": 3,
+            "m": int(Nx),
+            "n": int(Ny),
+            "p": 0,
+            "cfl_adap_dt": "T",
+            "cfl_target": 0.8,
+            "n_start": 0,
+            "t_stop": 4.0,
+            "t_save": 0.04,
+            # Simulation Algorithm Parameters
+            "num_patches": 3,
+            "model_eqns": 2,
+            "alt_soundspeed": "F",
+            "num_fluids": 2,
+            "time_stepper": 3,
+            "elliptic_smoothing": "T",
+            "elliptic_smoothing_iters": 50,
+            "igr": "T",
+            "igr_order": 5,
+            "igr_iter_solver": 1,
+            "num_igr_iters": 3,
+            "num_igr_warm_start_iters": 30,
+            "alf_factor": 10,
+            "bc_x%beg": -3,  # 11,
+            "bc_x%end": -3,  # 12
+            "bc_y%beg": -3,
+            "bc_y%end": -3,
+            # Formatted Database Files Structure Parameters
+            "format": 1,
+            "precision": 2,
+            "prim_vars_wrt": "T",
+            "schlieren_wrt": "T",
+            "fd_order": 4,
+            "schlieren_alpha(1)": 0.5,
+            "schlieren_alpha(2)": 0.5,
+            "parallel_io": "T",
+            # Patch 1: Left state
+            "patch_icpp(1)%geometry": 3,
+            "patch_icpp(1)%x_centroid": 0.5,
+            "patch_icpp(1)%y_centroid": 3,
+            "patch_icpp(1)%length_x": 1,
+            "patch_icpp(1)%length_y": 6,
+            "patch_icpp(1)%vel(1)": 0.0,
+            "patch_icpp(1)%vel(2)": 0.0,
+            "patch_icpp(1)%pres": 1.0,
+            "patch_icpp(1)%alpha_rho(1)": (1 - eps) * 1.0,
+            "patch_icpp(1)%alpha_rho(2)": eps,
+            "patch_icpp(1)%alpha(1)": 1 - eps,
+            "patch_icpp(1)%alpha(2)": eps,
+            # Patch 2: Top right state
+            "patch_icpp(2)%geometry": 3,
+            "patch_icpp(2)%alter_patch(1)": "T",
+            "patch_icpp(2)%x_centroid": 4,
+            "patch_icpp(2)%y_centroid": 2.25,
+            "patch_icpp(2)%length_x": 6,
+            "patch_icpp(2)%length_y": 1.5,
+            "patch_icpp(2)%vel(1)": 0.0,
+            "patch_icpp(2)%vel(2)": 0.0,
+            "patch_icpp(2)%pres": 0.1,
+            "patch_icpp(2)%alpha_rho(1)": (1 - eps) * 0.125,
+            "patch_icpp(2)%alpha_rho(2)": eps,
+            "patch_icpp(2)%alpha(1)": 1 - eps,
+            "patch_icpp(2)%alpha(2)": eps,
+            # Patch 3: Bottom right state
+            "patch_icpp(3)%geometry": 3,
+            "patch_icpp(3)%alter_patch(1)": "T",
+            "patch_icpp(3)%x_centroid": 4,
+            "patch_icpp(3)%y_centroid": 0.75,
+            "patch_icpp(3)%length_x": 6,
+            "patch_icpp(3)%length_y": 1.5,
+            "patch_icpp(3)%vel(1)": 0.0,
+            "patch_icpp(3)%vel(2)": 0.0,
+            "patch_icpp(3)%pres": 0.1,
+            "patch_icpp(3)%alpha_rho(1)": eps,
+            "patch_icpp(3)%alpha_rho(2)": (1 - eps) * 1.0,
+            "patch_icpp(3)%alpha(1)": eps,  # 0.95
+            "patch_icpp(3)%alpha(2)": 1 - eps,  # 0.05,
+            # Fluids Physical Parameters
+            "fluid_pp(1)%gamma": 1.0 / (1.5 - 1.0),
+            "fluid_pp(1)%pi_inf": 0,
+            "fluid_pp(2)%gamma": 1.0 / (1.4 - 1.0),
+            "fluid_pp(2)%pi_inf": 0.0,
+        }
+    )
+)