Merge MUSCL_fix branch + upstream master (4 commits)

Author: ChrisZYJ

CMakeLists.txt

@@ -636,7 +636,7 @@ exit 0
                     target_compile_options(${a_target} PRIVATE -fopenmp)
                     target_link_options(${a_target} PRIVATE -fopenmp)
                 elseif(CMAKE_Fortran_COMPILER_ID STREQUAL "LLVMFlang")
-                    target_compile_options(${a_target} PRIVATE -fopenmp --offload-arch=gfx90a)
+                    target_compile_options(${a_target} PRIVATE -fopenmp --offload-arch=gfx90a -fopenmp-target-fast -fopenmp-assume-threads-oversubscription -fopenmp-assume-teams-oversubscription)
                     target_link_options(${a_target} PRIVATE -fopenmp --offload-arch=gfx90a)
                 endif()
             endif()
@@ -1026,4 +1026,4 @@ site_name(SITE_NAME)
 
 configure_file(
     "${CMAKE_CURRENT_SOURCE_DIR}/toolchain/cmake/configuration.cmake.in"
-    "${CMAKE_CURRENT_BINARY_DIR}/configuration.txt")
+    "${CMAKE_CURRENT_BINARY_DIR}/configuration.txt")

docs/documentation/case.md

@@ -445,6 +445,7 @@ See @ref equations "Equations" for the mathematical models these parameters cont
 | `mp_weno`                  | Logical | Monotonicity preserving WENO |
 | `muscl_order`              | Integer | MUSCL order [1,2] |
 | `muscl_lim`                | Integer | MUSCL Slope Limiter: [1] minmod; [2] monotonized central; [3] Van Albada; [4] Van Leer; [5] SUPERBEE |
+| `muscl_eps`                | Real    | MUSCL limiter slope-product threshold (default: hard-coded thresholds; set to 0 for textbook behavior) |
 | `flux_lim`                 | Integer | Flux limiter for post-process: [1] minmod; [2] MUSCL; [3] OSPRE; [4] SUPERBEE |
 | `int_comp`                 | Logical | THINC Interface Compression |
 | `ic_eps`                   | Real    | Interface compression threshold (default: 1e-4) |
@@ -546,6 +547,10 @@ It is recommended to set `weno_eps` to $10^{-6}$ for WENO-JS, and to $10^{-40}$
 - `muscl_lim` specifies the slope limiter that is used in 2nd order MUSCL Reconstruction by an integer from 1 through 5.
 `muscl_lim = 1`, `2`, `3`, `4`, and `5` correspond to minmod, monotonized central, Van Albada, Van Leer, and SUPERBEE, respectively.
 
+- `muscl_eps` controls the slope-product activation threshold for all MUSCL limiters.
+When not set (default), the threshold is 1e-9 for minmod/MC, and 1e-6 for others.
+Setting `muscl_eps = 0` gives textbook limiter behavior where limiters activate whenever both slopes have the same sign.
+
 - `int_comp` activates interface compression using THINC used in MUSCL Reconstruction, with control parameters (`ic_eps`, and `ic_beta`).
 
 - `riemann_solver` specifies the choice of the Riemann solver that is used in simulation by an integer from 1 through 4.
@@ -1081,8 +1086,20 @@ When ``cyl_coord = 'T'`` is set in 2D the following constraints must be met:
 
 - `cantera_file` specifies the chemical mechanism file. If the file is part of the standard Cantera library, only the filename is required. Otherwise, the file must be located in the same directory as your `case.py` file
 
+### 18. Chemistry-Specific Boundary Conditions
+
+| Parameter          | Type    | Description                                                                 |
+| ---:               | :----:  | :---                                                                        |
+| `bc_[x,y,z]%%isothermal_in`    | Logical | Enable isothermal wall at the domain entrance (minimum coordinate).         |
+| `bc_[x,y,z]%%isothermal_out`   | Logical | Enable isothermal wall at the domain exit (maximum coordinate).             |
+| `bc_[x,y,z]%%Twall_in`         | Real    | Temperature [K] of the entrance isothermal wall.                            |
+| `bc_[x,y,z]%%Twall_out`        | Real    | Temperature [K] of the exit isothermal wall.                                |
+
+This boundary condition can be used for fixed-temperature (isothermal) walls at the domain extremities. It is exclusively available for reacting flows and requires chemistry to be enabled. It properly evaluates heat and species fluxes at the interface when ``chemistry = 'T'``, ``chem_params%%diffusion = 'T'``, and the corresponding domain boundary is set to a slip wall (`bc_[x,y,z]%%[beg,end]` = -15) or a no-slip wall (`bc_[x,y,z]%%[beg,end]` = -16).
+
+
 
-### 18. GPU Performance (NVIDIA UVM)
+### 19. GPU Performance (NVIDIA UVM)
 
 | Parameter                  | Type    | Description                                              |
 | ---:                       | :---:   | :---                                                     |

examples/2D_Thermal_Flatplate/case.py

@@ -0,0 +1,76 @@
+#!/usr/bin/env python3
+import json
+import math
+
+import cantera as ct
+
+Lx = 0.05
+Ly = 0.05
+
+ctfile = "h2o2.yaml"
+sol_L = ct.Solution(ctfile)
+sol_L.TPY = 1125, ct.one_atm, "O2:0.21,N2:0.79"
+# Configuring case dictionary
+case = {
+    "run_time_info": "T",
+    "x_domain%beg": 0.0,
+    "x_domain%end": Lx,
+    "y_domain%beg": 0.0,
+    "y_domain%end": Ly,
+    "m": 699,
+    "n": 699,
+    "p": 0,
+    "dt": 4.0e-08,
+    "t_step_start": 0,
+    "t_step_stop": 75000,
+    "t_step_save": 4500,
+    "num_patches": 1,
+    "model_eqns": 2,
+    "alt_soundspeed": "F",
+    "num_fluids": 1,
+    "mpp_lim": "F",
+    "mixture_err": "T",
+    "time_stepper": 3,
+    "mp_weno": "F",
+    "weno_order": 5,
+    "weno_eps": 1e-16,
+    "riemann_solver": 2,
+    "wave_speeds": 1,
+    "avg_state": 2,
+    "bc_x%beg": -7,
+    "bc_x%end": -3,
+    "bc_y%beg": -16,
+    "bc_y%end": -3,
+    "bc_y%isothermal_in": "T",
+    "bc_y%Twall_in": 600.0,
+    "format": 1,
+    "precision": 2,
+    "prim_vars_wrt": "T",
+    "parallel_io": "T",
+    "chemistry": "T",
+    "chem_params%diffusion": "T",
+    "chem_params%reactions": "F",
+    "cantera_file": ctfile,
+    "chem_wrt_T": "T",
+    "patch_icpp(1)%geometry": 3,
+    "patch_icpp(1)%hcid": 291,
+    "patch_icpp(1)%x_centroid": Lx / 2,
+    "patch_icpp(1)%y_centroid": Ly / 2,
+    "patch_icpp(1)%length_x": Lx,
+    "patch_icpp(1)%length_y": Ly,
+    "patch_icpp(1)%vel(1)": 0,
+    "patch_icpp(1)%vel(2)": 0,
+    "patch_icpp(1)%pres": 101325,
+    "patch_icpp(1)%alpha_rho(1)": 1.00,
+    "patch_icpp(1)%alpha(1)": 1,
+    "fluid_pp(1)%gamma": 1.0e00 / (1.4e00 - 1.0e00),
+    "fluid_pp(1)%pi_inf": 0.0e00,
+    "viscous": "T",
+    "fluid_pp(1)%Re(1)": 100000,
+}
+for i in range(len(sol_L.Y)):
+    case[f"chem_wrt_Y({i + 1})"] = "T"
+    case[f"patch_icpp(1)%Y({i + 1})"] = sol_L.Y[i]
+
+if __name__ == "__main__":
+    print(json.dumps(case))

src/common/include/1dHardcodedIC.fpp

@@ -51,6 +51,21 @@
         molar_mass_inv = y1/31.998_wp + y2/18.01508_wp + y3/16.04256_wp + y4/28.0134_wp
 
         q_prim_vf(eqn_idx%cont%beg)%sf(i, 0, 0) = 1.01325_wp*(10.0_wp)**5/(temp*8.3144626_wp*1000.0_wp*molar_mass_inv)
+
+    case(191)  ! 1D Dual Isothermal case
+
+        q_prim_vf(eqn_idx%E)%sf(i, 0, 0) = 101325.0_wp
+        q_prim_vf(eqn_idx%mom%beg)%sf(i, 0, 0) = 0.0_wp
+        q_prim_vf(eqn_idx%species%beg)%sf(i, 0, 0) = 1.0_wp
+
+        if (x_cc(i) <= 0.025_wp) then
+            temp = 700.0_wp + ((1000.0_wp - 700.0_wp)/0.025_wp)*x_cc(i)
+        else
+            temp = 1200.0_wp + ((900.0_wp - 1000.0_wp)/0.025_wp)*(x_cc(i) - 0.025_wp)
+        end if
+
+        molar_mass_inv = 1.0_wp/2.01588_wp
+        q_prim_vf(eqn_idx%cont%beg)%sf(i, 0, 0) = 101325.0_wp/(temp*8.3144626_wp*1000.0_wp*molar_mass_inv)
     case default
         call s_int_to_str(patch_id, iStr)
         call s_mpi_abort("Invalid hcid specified for patch " // trim(iStr))

src/common/include/2dHardcodedIC.fpp

@@ -8,8 +8,16 @@
     real(wp) :: sinA, cosA
     real(wp) :: r_sq
 
+    ! # 291 - Shear/Thermal Layer Case
+    real(wp) :: delta_shear, u_max, u_mean
+    real(wp) :: T_wall, T_inf, P_atm, T_loc
+    real(wp) :: delta_th, R_mix
+    real(wp) :: Y_N2, Y_O2, MW_N2, MW_O2
+    real(wp) :: bottom_blend_u, bottom_blend_T
+
     ! # 207
     real(wp) :: sigma, gauss1, gauss2
+
     ! # 208
     real(wp) :: ei, d, fsm, alpha_air, alpha_sf6
 
@@ -305,6 +313,32 @@
             q_prim_vf(eqn_idx%mom%beg + 1)%sf(i, j, &
                       & 0) = 112.99092883944267*((0.1/0.3))*x_cc(i)*exp(0.5*(1 - sqrt(x_cc(i)**2 + y_cc(j)**2)))
         end if
+    case (291)  ! Isothermal Flat Plate
+        T_inf = 1125.0_wp
+        T_wall = 600.0_wp
+        P_atm = 101325.0_wp
+
+        ! Boundary/Shear Layer thicknesses
+        delta_th = 0.0003_wp  ! Thermal BL thickness
+        delta_shear = 8e-3_wp  ! Velocity BL thickness
+
+        u_max = 50.0_wp  ! Freestream Velocity (m/s)
+
+        MW_N2 = 28.0134e-3_wp
+        MW_O2 = 31.999e-3_wp
+        Y_N2 = 0.767_wp
+        Y_O2 = 0.233_wp
+        R_mix = 8.314462618_wp*((Y_N2/MW_N2) + (Y_O2/MW_O2))
+        bottom_blend_u = tanh(y_cc(j)/delta_shear)
+        bottom_blend_T = tanh(y_cc(j)/delta_th)
+        u_mean = u_max*bottom_blend_u
+        T_loc = T_wall + (T_inf - T_wall)*bottom_blend_T
+        q_prim_vf(eqn_idx%cont%beg)%sf(i, j, 0) = P_atm/(R_mix*T_loc)
+        q_prim_vf(eqn_idx%mom%beg)%sf(i, j, 0) = u_mean
+        q_prim_vf(eqn_idx%mom%end)%sf(i, j, 0) = 0.0_wp
+        q_prim_vf(eqn_idx%E)%sf(i, j, 0) = P_atm
+        q_prim_vf(eqn_idx%species%beg)%sf(i, j, 0) = Y_O2
+        q_prim_vf(eqn_idx%species%end)%sf(i, j, 0) = Y_N2
     case default
         if (proc_rank == 0) then
             call s_int_to_str(patch_id, iStr)