Add Herschel-Bulkley non-Newtonian viscosity model

.github/workflows/frontier_amd/bench.sh

@@ -0,0 +1 @@
+../frontier/bench.sh

.github/workflows/frontier_amd/build.sh

@@ -0,0 +1 @@
+../frontier/build.sh

.github/workflows/frontier_amd/submit.sh

@@ -0,0 +1 @@
+../frontier/submit.sh

.github/workflows/frontier_amd/test.sh

@@ -0,0 +1 @@
+../frontier/test.sh

examples/2D_lid_driven_cavity_nn/case.py

@@ -0,0 +1,116 @@
+#!/usr/bin/env python3
+"""
+2D Lid-Driven Cavity Flow with Herschel-Bulkley Non-Newtonian Fluid
+
+Re = 5000, n = 1.5 (shear-thickening) with mesh stretching near walls.
+
+HB model: mu = (tau0/gdot)*(1 - exp(-m*gdot)) + K * gdot^(n-1)
+Re_gen = rho * U^(2-n) * L^n / K = 1 * 1^0.5 * 1^1.5 / 0.0002 = 5000
+
+Mesh stretching: cosh-based clustering near all 4 walls (x_a, x_b, y_a, y_b).
+"""
+import json
+
+eps = 1e-6
+
+# HB model parameters
+tau0 = 0.0          # Yield stress (set to 0 for power-law fluid)
+K = 0.0002          # Consistency index (Re=5000: K = 1/5000)
+nn = 1.5            # Flow behavior index (shear-thickening)
+mu_min = 0.00002    # K * gdot_min^(n-1) = 0.0002 * (0.01)^0.5
+mu_max = 0.0632     # K * gdot_max^(n-1) = 0.0002 * (1e5)^0.5
+hb_m = 1000.0       # Papanastasiou regularization parameter
+mu_bulk = 0.0
+
+lid_velocity = 1.0  # Lid velocity (m/s)
+
+# Configuring case dictionary
+print(
+    json.dumps(
+        {
+            # Logistics
+            "run_time_info": "T",
+            # Computational Domain Parameters
+            "x_domain%beg": 0.0,
+            "x_domain%end": 1.0,
+            "y_domain%beg": 0.0,
+            "y_domain%end": 1.0,
+            "m": 255,
+            "n": 255,
+            "p": 0,
+            "cfl_adap_dt": "T",
+            "cfl_target": 0.5,
+            "n_start": 0,
+            "t_stop": 50.0,
+            "t_save": 25.0,
+            # Simulation Algorithm Parameters
+            "num_patches": 1,
+            "model_eqns": 2,
+            "alt_soundspeed": "F",
+            "num_fluids": 2,
+            "mpp_lim": "F",
+            "mixture_err": "T",
+            "time_stepper": 3,
+            "weno_order": 5,
+            "weno_eps": 1e-16,
+            "mapped_weno": "T",
+            "weno_Re_flux": "T",
+            "mp_weno": "T",
+            "weno_avg": "T",
+            "riemann_solver": 2,
+            "wave_speeds": 1,
+            "avg_state": 2,
+            "bc_x%beg": -16,
+            "bc_x%end": -16,
+            "bc_y%beg": -16,
+            "bc_y%end": -16,
+            "bc_y%ve1": lid_velocity,
+            "viscous": "T",
+            # Formatted Database Files Structure Parameters
+            "format": 1,
+            "precision": 2,
+            "prim_vars_wrt": "T",
+            "omega_wrt(3)": "T",
+            "fd_order": 4,
+            "parallel_io": "T",
+            # Patch 1: Base
+            "patch_icpp(1)%geometry": 3,
+            "patch_icpp(1)%x_centroid": 0.5,
+            "patch_icpp(1)%y_centroid": 0.5,
+            "patch_icpp(1)%length_x": 1.0,
+            "patch_icpp(1)%length_y": 1.0,
+            "patch_icpp(1)%vel(1)": 0,
+            "patch_icpp(1)%vel(2)": 0.0,
+            "patch_icpp(1)%pres": 1e3,
+            "patch_icpp(1)%alpha_rho(1)": 0.5,
+            "patch_icpp(1)%alpha(1)": 0.5,
+            "patch_icpp(1)%alpha_rho(2)": 0.5,
+            "patch_icpp(1)%alpha(2)": 0.5,
+            # Fluids Physical Parameters
+            # Fluid 1:
+            "fluid_pp(1)%gamma": 1.0 / (1.4 - 1.0),
+            "fluid_pp(1)%pi_inf": 0.0,
+            "fluid_pp(1)%Re(1)": 1.0 / K,
+            "fluid_pp(1)%non_newtonian": "T",
+            "fluid_pp(1)%tau0": tau0,
+            "fluid_pp(1)%K": K,
+            "fluid_pp(1)%nn": nn,
+            "fluid_pp(1)%mu_max": mu_max,
+            "fluid_pp(1)%mu_min": mu_min,
+            "fluid_pp(1)%mu_bulk": mu_bulk,
+            "fluid_pp(1)%hb_m": hb_m,
+            # Fluid 2:
+            "fluid_pp(2)%gamma": 1.0 / (1.4 - 1.0),
+            "fluid_pp(2)%pi_inf": 0.0,
+            "fluid_pp(2)%Re(1)": 1.0 / K,
+            "fluid_pp(2)%non_newtonian": "T",
+            "fluid_pp(2)%tau0": tau0,
+            "fluid_pp(2)%K": K,
+            "fluid_pp(2)%nn": nn,
+            "fluid_pp(2)%mu_max": mu_max,
+            "fluid_pp(2)%mu_min": mu_min,
+            "fluid_pp(2)%mu_bulk": mu_bulk,
+            "fluid_pp(2)%hb_m": hb_m,
+        }
+    )
+)

examples/2D_poiseuille_nn/case.py

@@ -0,0 +1,158 @@
+#!/usr/bin/env python3
+"""
+2D Poiseuille Flow with Herschel-Bulkley Non-Newtonian Fluid
+
+Pressure-driven channel flow between two no-slip walls, driven by a constant
+body force in the streamwise (x) direction. Periodic BCs in x, no-slip in y.
+
+HB model: mu = (tau0/gdot)*(1 - exp(-m*gdot)) + K * gdot^(n-1)
+  - tau0: yield stress
+  - K:    consistency index
+  - n:    flow behavior index (< 1 shear-thinning, > 1 shear-thickening)
+  - m:    Papanastasiou regularization parameter
+
+For Newtonian Poiseuille validation, set tau0=0, nn=1, K=mu.
+The analytical solution is: u(y) = (G/(2*mu)) * y * (H - y)
+where G = rho * g_x is the effective pressure gradient.
+"""
+import json
+import math
+
+# === Channel geometry (square domain) ===
+L = 1.0         # Channel length (streamwise, x)
+H = 1.0         # Channel height (wall-normal, y)
+
+# === Grid resolution ===
+Nx = 24         # Cells in x (streamwise, minimal — periodic)
+Ny = 81         # Cells in y (wall-normal)
+
+# === Fluid properties ===
+rho = 1.0       # Density
+p0 = 1e5        # Reference pressure (high for low Mach)
+gamma = 1.4     # Ratio of specific heats
+
+# Sound speed and CFL
+c = math.sqrt(gamma * p0 / rho)
+dx = L / (Nx + 1)
+cfl = 0.3
+dt = cfl * dx / c
+
+# === Body force (pressure gradient substitute) ===
+# G = rho * g_x acts as dp/dx driving force
+g_x = 0.5
+
+# === HB non-Newtonian model parameters ===
+tau0 = 0.0          # Yield stress (set 0 for power-law)
+K = 0.1            # Consistency index
+nn = 2.0          # Flow behavior index (< 1 = shear-thinning)
+hb_m = 1000.0       # Papanastasiou regularization parameter
+mu_min = 1e-4       # Minimum viscosity bound
+mu_max = 10.0       # Maximum viscosity bound
+mu_bulk = 0.0       # Bulk viscosity
+
+# Reference Re based on consistency index (used as baseline)
+Re_ref = 1.0 / K    # = 100
+
+# === Time control ===
+t_end = 10.0        # End time (allow flow to reach steady state)
+t_save = 5.0        # Save interval
+
+eps = 1e-6
+
+# Configuring case dictionary
+print(
+    json.dumps(
+        {
+            # Logistics
+            "run_time_info": "T",
+            # Computational Domain Parameters
+            "x_domain%beg": 0.0,
+            "x_domain%end": L,
+            "y_domain%beg": 0.0,
+            "y_domain%end": H,
+            "m": Nx,
+            "n": Ny,
+            "p": 0,
+            "cfl_adap_dt": "T",
+            "cfl_target": cfl,
+            "n_start": 0,
+            "t_stop": t_end,
+            "t_save": t_save,
+            # Simulation Algorithm Parameters
+            "num_patches": 1,
+            "model_eqns": 2,
+            "alt_soundspeed": "F",
+            "num_fluids": 2,
+            "mpp_lim": "F",
+            "mixture_err": "T",
+            "time_stepper": 3,
+            "weno_order": 5,
+            "weno_eps": 1e-16,
+            "mapped_weno": "T",
+            "weno_Re_flux": "T",
+            "mp_weno": "T",
+            "weno_avg": "T",
+            "riemann_solver": 2,
+            "wave_speeds": 1,
+            "avg_state": 2,
+            # Boundary Conditions
+            # Periodic in x (streamwise), no-slip walls in y
+            "bc_x%beg": -1,
+            "bc_x%end": -1,
+            "bc_y%beg": -16,
+            "bc_y%end": -16,
+            # Viscous
+            "viscous": "T",
+            # Body Force (drives the flow like a pressure gradient)
+            "bf_x": "T",
+            "g_x": g_x,
+            "k_x": 0.0,
+            "w_x": 0.0,
+            "p_x": 0.0,
+            # Formatted Database Files Structure Parameters
+            "format": 1,
+            "precision": 2,
+            "prim_vars_wrt": "T",
+            "omega_wrt(3)": "T",
+            "fd_order": 4,
+            "parallel_io": "T",
+            # Patch 1: Entire channel domain (initially at rest)
+            "patch_icpp(1)%geometry": 3,
+            "patch_icpp(1)%x_centroid": L / 2.0,
+            "patch_icpp(1)%y_centroid": H / 2.0,
+            "patch_icpp(1)%length_x": L,
+            "patch_icpp(1)%length_y": H,
+            "patch_icpp(1)%vel(1)": 0.0,
+            "patch_icpp(1)%vel(2)": 0.0,
+            "patch_icpp(1)%pres": p0,
+            "patch_icpp(1)%alpha_rho(1)": rho * 0.5,
+            "patch_icpp(1)%alpha(1)": 0.5,
+            "patch_icpp(1)%alpha_rho(2)": rho * 0.5,
+            "patch_icpp(1)%alpha(2)": 0.5,
+            # Fluid 1: HB non-Newtonian fluid
+            "fluid_pp(1)%gamma": 1.0 / (gamma - 1.0),
+            "fluid_pp(1)%pi_inf": 0.0,
+            "fluid_pp(1)%Re(1)": Re_ref,
+            "fluid_pp(1)%non_newtonian": "T",
+            "fluid_pp(1)%tau0": tau0,
+            "fluid_pp(1)%K": K,
+            "fluid_pp(1)%nn": nn,
+            "fluid_pp(1)%hb_m": hb_m,
+            "fluid_pp(1)%mu_min": mu_min,
+            "fluid_pp(1)%mu_max": mu_max,
+            "fluid_pp(1)%mu_bulk": mu_bulk,
+            # Fluid 2: same properties (single-phase effectively)
+            "fluid_pp(2)%gamma": 1.0 / (gamma - 1.0),
+            "fluid_pp(2)%pi_inf": 0.0,
+            "fluid_pp(2)%Re(1)": Re_ref,
+            "fluid_pp(2)%non_newtonian": "T",
+            "fluid_pp(2)%tau0": tau0,
+            "fluid_pp(2)%K": K,
+            "fluid_pp(2)%nn": nn,
+            "fluid_pp(2)%hb_m": hb_m,
+            "fluid_pp(2)%mu_min": mu_min,
+            "fluid_pp(2)%mu_max": mu_max,
+            "fluid_pp(2)%mu_bulk": mu_bulk,
+        }
+    )
+)

src/common/m_derived_types.fpp

@@ -367,6 +367,14 @@ module m_derived_types
         real(wp) :: qv      !< reference energy per unit mass for SGEOS, q (see Le Metayer (2004))
         real(wp) :: qvp     !< reference entropy per unit mass for SGEOS, q' (see Le Metayer (2004))
         real(wp) :: G
+        logical :: non_newtonian  !< Non-Newtonian fluid flag
+        real(wp) :: tau0          !< Yield stress (Herschel-Bulkley model)
+        real(wp) :: K             !< Consistency index (Herschel-Bulkley model)
+        real(wp) :: nn            !< Flow behavior index (Herschel-Bulkley model)
+        real(wp) :: mu_max        !< Maximum viscosity limit (shear)
+        real(wp) :: mu_min        !< Minimum viscosity limit (shear)
+        real(wp) :: mu_bulk       !< Bulk viscosity for non-Newtonian fluids
+        real(wp) :: hb_m          !< Papanastasiou regularization parameter
     end type physical_parameters
 
     !> Derived type annexing the physical parameters required for sub-grid bubble models

src/post_process/m_global_parameters.fpp

@@ -434,6 +434,14 @@ contains
             fluid_pp(i)%qv = 0._wp
             fluid_pp(i)%qvp = 0._wp
             fluid_pp(i)%G = dflt_real
+            fluid_pp(i)%non_newtonian = .false.
+            fluid_pp(i)%tau0 = 0._wp
+            fluid_pp(i)%K = 0._wp
+            fluid_pp(i)%nn = 1._wp
+            fluid_pp(i)%mu_max = dflt_real
+            fluid_pp(i)%mu_min = 0._wp
+            fluid_pp(i)%mu_bulk = 0._wp
+            fluid_pp(i)%hb_m = 1000._wp
         end do
 
         ! Subgrid bubble parameters

src/post_process/m_mpi_proxy.fpp

@@ -130,6 +130,10 @@ contains
             call MPI_BCAST(fluid_pp(i)%qv, 1, mpi_p, 0, MPI_COMM_WORLD, ierr)
             call MPI_BCAST(fluid_pp(i)%qvp, 1, mpi_p, 0, MPI_COMM_WORLD, ierr)
             call MPI_BCAST(fluid_pp(i)%G, 1, mpi_p, 0, MPI_COMM_WORLD, ierr)
+            call MPI_BCAST(fluid_pp(i)%non_newtonian, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
+            #:for VAR in ['tau0', 'K', 'nn', 'mu_max', 'mu_min', 'mu_bulk', 'hb_m']
+                call MPI_BCAST(fluid_pp(i)%${VAR}$, 1, mpi_p, 0, MPI_COMM_WORLD, ierr)
+            #:endfor
         end do
 
         ! Subgrid bubble parameters

src/pre_process/m_global_parameters.fpp

@@ -592,6 +592,14 @@ contains
             fluid_pp(i)%qv = 0._wp
             fluid_pp(i)%qvp = 0._wp
             fluid_pp(i)%G = 0._wp
+            fluid_pp(i)%non_newtonian = .false.
+            fluid_pp(i)%tau0 = 0._wp
+            fluid_pp(i)%K = 0._wp
+            fluid_pp(i)%nn = 1._wp
+            fluid_pp(i)%mu_max = dflt_real
+            fluid_pp(i)%mu_min = 0._wp
+            fluid_pp(i)%mu_bulk = 0._wp
+            fluid_pp(i)%hb_m = 1000._wp
         end do
 
         Bx0 = dflt_real

src/pre_process/m_mpi_proxy.fpp

@@ -143,6 +143,11 @@ contains
                 call MPI_BCAST(fluid_pp(i)%${VAR}$, 1, mpi_p, 0, MPI_COMM_WORLD, ierr)
             #:endfor
 
+            call MPI_BCAST(fluid_pp(i)%non_newtonian, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
+            #:for VAR in ['tau0', 'K', 'nn', 'mu_max', 'mu_min', 'mu_bulk', 'hb_m']
+                call MPI_BCAST(fluid_pp(i)%${VAR}$, 1, mpi_p, 0, MPI_COMM_WORLD, ierr)
+            #:endfor
+
             call MPI_BCAST(simplex_params%perturb_dens(i), 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
             call MPI_BCAST(simplex_params%perturb_dens_freq(i), 1, mpi_p, 0, MPI_COMM_WORLD, ierr)
             call MPI_BCAST(simplex_params%perturb_dens_scale(i), 1, mpi_p, 0, MPI_COMM_WORLD, ierr)