hehe

Author: Dimitrios Adam

CMakeLists.txt

@@ -240,7 +240,11 @@ if (CMAKE_BUILD_TYPE STREQUAL "Release")
             message(STATUS "LTO/IPO is not available with NVHPC using Unified Memory")
         elseif(CMAKE_Fortran_COMPILER_VERSION VERSION_LESS "23.11")
             message(STATUS "LTO/IPO is not supported in NVHPC Version < 23.11. Use a newer version of NVHPC for best performance.")
-        else()
+        elseif (CMAKE_Fortran_COMPILER_VERSION VERSION_GREATER "24.11" AND CMAKE_Fortran_COMPILER_VERSION VERSION_LESS "25.9") 
+            message(STATUS "LTO/IPO is not supported in NVHPC Version 24.11 to 25.9. Use >=25.9 or (<=24.11 && > 23.11) Performance will be degraded.")
+            set(NVHPC_USE_TWO_PASS_IPO FALSE)
+          else()
+ 
             message(STATUS "Performing IPO using -Mextract followed by -Minline")
             set(NVHPC_USE_TWO_PASS_IPO TRUE)
         endif()

examples/1D_reactive_shocktube/case.py

@@ -110,8 +110,6 @@
 if args.chemistry:
     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]
-        case[f"patch_icpp(2)%Y({i+1})"] = sol_R.Y[i]
 
 if __name__ == "__main__":
     print(json.dumps(case))

examples/2D_riemann_test/case.py

@@ -1,6 +1,8 @@
 #!/usr/bin/env python3
 import json
 import math
+import cantera as ct
+
 Lx = 0.04347
 Ly = 0.01466+3.048/1000
 Cav_Lx = 43.47/1000-21.99/1000-14.19/1000
@@ -9,10 +11,12 @@
 Cav_CenterY = 3.048/2000 
 Cav_Center2X = 2*Cav_CenterX+14.15/1000+21.99/2000
 
+
+ctfile = "h2.yaml"
+sol_L = ct.Solution(ctfile)
+sol_L.TPY = 1125, ct.one_atm, "O2:0.21,N2:0.79"
 # Configuring case dictionary
-print(
-    json.dumps(
-        {
+case =     {    
             # Logistics
             "run_time_info": "T",
             # Computational Domain Parameters
@@ -23,10 +27,10 @@
             "m": 499,
             "n": 499,
             "p": 0,
-            "dt": 8.0e-08,
+            "dt": 4.0e-08,
             "t_step_start": 0,
-            "t_step_stop": 2000,
-            "t_step_save": 50,
+            "t_step_stop": 10000,
+            "t_step_save": 1000,
             # Simulation Algorithm Parameters
             "num_patches": 1,
             "model_eqns": 2,
@@ -45,7 +49,7 @@
             "wave_speeds": 1,
             "avg_state": 2,
             "bc_x%beg": -7,
-            "bc_x%end": -8,
+            "bc_x%end": -3,
             "bc_y%beg": -16,
             "bc_y%end": -16,
             # Formatted Database Files Structure Parameters
@@ -54,7 +58,14 @@
             "prim_vars_wrt": "T",
             "parallel_io": "T",
                  "ib": "T",
-            "num_ibs": 2,
+            "num_ibs": 1,
+            "omega_wrt(3)": "T",
+          "fd_order": 2,
+              "chemistry": "T" ,
+      "chem_params%diffusion": "T",
+      "chem_params%reactions": "F",
+      "cantera_file": ctfile,
+       "chem_wrt_T": "T",
             # Patch 1: Base
             "patch_icpp(1)%geometry": 3,
             "patch_icpp(1)%hcid": 291,
@@ -71,21 +82,26 @@
              "patch_ib(1)%geometry": 3,
             "patch_ib(1)%x_centroid": 0,
             "patch_ib(1)%y_centroid": 0,
-            "patch_ib(1)%length_x": 2*Cav_Lx,
+            "patch_ib(1)%length_x": 0.1,
             "patch_ib(1)%length_y": 2*Cav_Ly,
             "patch_ib(1)%slip": "F",
-            # Patch 2: IBM
-            "patch_ib(2)%geometry": 3,
-            "patch_ib(2)%x_centroid": Lx,
-            "patch_ib(2)%y_centroid": 0,
-            "patch_ib(2)%length_x": 2*21.99/1000,
-            "patch_ib(2)%length_y": 2*Cav_Ly,
-            "patch_ib(2)%slip": "F",
+
+     #       # Patch 2: IBM
+     #       "patch_ib(2)%geometry": 3,
+     #       "patch_ib(2)%x_centroid": Lx,
+     #       "patch_ib(2)%y_centroid": 0,
+     #       "patch_ib(2)%length_x": 2*21.99/1000,
+     #       "patch_ib(2)%length_y": 2*Cav_Ly,
+     #       "patch_ib(2)%slip": "F",
             # Fluids Physical Parameters
             "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))

examples/2D_riemann_test/h2.yaml

@@ -0,0 +1,84 @@
+description: |-
+  Hydrogen-Oxygen submechanism extracted from GRI-Mech 3.0.
+  Modified from the original to include N2.
+  
+  Redlich-Kwong coefficients are based on tabulated critical properties or
+  estimated according to the method of Joback and Reid, "Estimation of pure-
+  component properties from group-contributions," Chem. Eng. Comm. 57 (1987)
+  233-243
+
+generator: ck2yaml
+input-files: [h2o2.inp, gri30_tran.dat]
+cantera-version: 2.5.0
+date: Wed, 11 Dec 2019 16:59:04 -0500
+
+units: {length: cm, time: s, quantity: mol, activation-energy: cal/mol}
+
+phases:
+- name: ohmech
+  thermo: ideal-gas
+  elements: [O, N]
+  species: [ O2,  N2]
+  kinetics: gas
+  transport: mixture-averaged
+  state: {T: 300.0, P: 1 atm}
+
+- name: ohmech-RK
+  thermo: Redlich-Kwong
+  elements: [O,  N]
+  species: [ O2,  N2]
+  kinetics: gas
+  transport: mixture-averaged
+  state: {T: 300.0, P: 1 atm}
+
+species:
+- name: O2
+  composition: {O: 2}
+  thermo:
+    model: NASA7
+    temperature-ranges: [200.0, 1000.0, 3500.0]
+    data:
+    - [3.78245636, -2.99673416e-03, 9.84730201e-06, -9.68129509e-09, 3.24372837e-12,
+      -1063.94356, 3.65767573]
+    - [3.28253784, 1.48308754e-03, -7.57966669e-07, 2.09470555e-10, -2.16717794e-14,
+      -1088.45772, 5.45323129]
+    note: TPIS89
+  transport:
+    model: gas
+    geometry: linear
+    well-depth: 107.4
+    diameter: 3.458
+    polarizability: 1.6
+    rotational-relaxation: 3.8
+  equation-of-state:
+    model: Redlich-Kwong
+    a: 1.74102e+12
+    b: 22.08100907
+
+- name: N2
+  composition: {N: 2}
+  thermo:
+    model: NASA7
+    temperature-ranges: [300.0, 1000.0, 5000.0]
+    data:
+    - [3.298677, 1.4082404e-03, -3.963222e-06, 5.641515e-09, -2.444854e-12,
+      -1020.8999, 3.950372]
+    - [2.92664, 1.4879768e-03, -5.68476e-07, 1.0097038e-10, -6.753351e-15,
+      -922.7977, 5.980528]
+    note: '121286'
+  transport:
+    model: gas
+    geometry: linear
+    well-depth: 97.53
+    diameter: 3.621
+    polarizability: 1.76
+    rotational-relaxation: 4.0
+  equation-of-state:
+    model: Redlich-Kwong
+    a: 1.55976e+12
+    b: 26.81724983
+
+reactions:
+- equation: O2 + N2 + M <=> O2 + N2 + M  # Reaction 1
+  type: three-body
+  rate-constant: {A: 1.2e+17, b: -1.0, Ea: 0.0}

src/common/include/2dHardcodedIC.fpp

@@ -18,6 +18,10 @@
 real(wp) :: y_shear_layer, delta_shear, tanh_arg, mix_factor
 real(wp) :: u_max, p_ref, rho_ref
 real(wp) :: u_mean, v_kick, perturbation_amp, perturbation_k
+real(wp) :: T_wall, T_inf, P_atm, y_height, T_loc,y_dist
+real(wp) :: delta_th,  T_target, rho_target, R_mix
+real(wp) :: Y_N2, Y_O2, MW_N2, MW_O2
+real(wp) :: x_lim_left, x_lim_right
 
     eps = 1.e-9_wp
 #:enddef
@@ -338,53 +342,82 @@ real(wp) :: u_mean, v_kick, perturbation_amp, perturbation_k
 
     case (291)
 
-y_shear_layer = 3.048e-3_wp  
+    ! Setup Values
+T_wall  = 600.0_wp
+T_inf   = 1125.0_wp
+P_atm   = 101325.0_wp  ! 1 atm in Pa
+delta_th = 0.0001_wp    ! Thermal BL thickness (e.g., 2mm - ADJUST THIS)
+
+! Geometry Limits (converting mm to meters if your grid is in meters)
+x_lim_left  = (43.48e-3_wp - 21.99e-3_wp - 14.15e-3_wp)
+x_lim_right = (43.48e-3_wp - 21.99e-3_wp)
+
+! Air Properties (Approximate)
+MW_N2 = 28.0134e-3_wp ! kg/mol
+MW_O2 = 31.999e-3_wp  ! kg/mol
+Y_N2  = 0.767_wp      ! Mass fraction
+Y_O2  = 0.233_wp
+
+! Calculate Mixture Gas Constant R_mix = R_u * sum(Y_i / MW_i)
+R_mix = 8.314462618_wp * ( (Y_N2 / MW_N2) + (Y_O2 / MW_O2) )
+
+y_shear_layer = 3.048e-3_wp
 
 ! 2. Shear Layer Thickness (Smoothing)
-! This spreads the jump over approx 1.0 mm to prevent acoustic shock.
-delta_shear = 0.5e-3_wp      
+delta_shear = 0.8e-3_wp
 
-! 3. Flow Reference Conditions
+! Flow Reference Conditions
 u_max   = 50.0_wp       ! Freestream Velocity (m/s)
-p_ref   = 101325.0_wp   ! Reference Pressure (Pa) - KEEP CONSTANT
-rho_ref = 1.0_wp        ! Reference Density (kg/m3) - KEEP CONSTANT
+p_ref   = 101325.0_wp   ! Reference Pressure (Pa)
+rho_ref = 1.0_wp        ! Reference Density (kg/m3)
+
+! Perturbation Parameters (To trigger vortices/Rossiter modes)
+perturbation_amp = 0.05_wp * u_max   ! Amplitude: 5% of freestream
+perturbation_k   = 200.0_wp          ! Wavenumber: Oscillations along X
+
+
+
+  
+        IF (y_cc(j) > y_shear_layer) THEN
 
-! 4. Perturbation Parameters (To trigger vortices)
-! Amplitude: 5% of freestream velocity
-! Wavenumber: Makes the kick oscillate along X
-perturbation_amp = 0.05_wp * u_max 
-perturbation_k   = 200.0_wp ! Adjust to fit
+            tanh_arg = (y_cc(j) - y_shear_layer) / delta_shear
+            mix_factor = tanh(tanh_arg)
+            u_mean     = u_max * mix_factor
 
-tanh_arg = (y_cc(j) - y_shear_layer) / delta_shear
-        mix_factor = 0.5_wp * (1.0_wp + tanh(tanh_arg))
+            v_kick = 0.0_wp
 
-        ! --- B. Mean Velocity ---
-        u_mean = u_max * mix_factor
+        ELSE
+
+
+            u_mean = 0.0_wp
+            v_kick = 0.0_wp
+
+        END IF
+
+!IF (y_cc(j) > y_shear_layer-0.01) THEN
 
-        ! --- C. Vertical Perturbation (The "Kick") ---
-        ! This adds a small vertical wiggle ONLY near the shear layer.
-        ! exp(...) ensures the kick is zero far away from the interface.
-        ! sin(...) makes it oscillate along the streamwise direction.
 
-        v_kick = perturbation_amp * &
-                 sin(perturbation_k * x_cc(i)) * &
-                 exp( -1.0_wp * (tanh_arg**2) )
+        ! We are in the Left or Right regions: Apply Tanh Profile
+ !       y_dist = y_cc(j) ! Distance from bottom wall
+        
+        ! Tanh Temperature Profile: T = Tw + (Tinf - Tw) * tanh(y/delta)
+  !      T_loc = T_wall + (T_inf - T_wall) * tanh((y_dist-y_shear_layer) / delta_th)
+
 
-        ! --- D. Assign to State Vector ---
-        ! 1. Density (Constant)
-        q_prim_vf(contxb)%sf(i,j,0) = rho_ref
+   ! else
+        ! We are in the middle gap: Uniform Freestream Temperature
+        T_loc = T_inf
+   ! end if
 
-        ! 2. U-Velocity (Smooth Profile)
+        q_prim_vf(contxb)%sf(i,j,0) = P_atm / (R_mix * T_loc)
         q_prim_vf(momxb)%sf(i,j,0) = u_mean
 
-        ! 3. V-Velocity (The Perturbation)
         q_prim_vf(momxe)%sf(i,j,0) = v_kick
 
-        ! 4. W-Velocity (assuming 2D or 0 initially)
-        ! q_prim_vf(momxz)%sf(i,j,0) = 0.0_wp 
-
-        ! 5. Pressure (Constant - Crucial for stability!)
         q_prim_vf(E_idx)%sf(i,j,0)  = p_ref
+
+       q_prim_vf(chemxb)%sf(i,j,0) =   0.21_wp
+        q_prim_vf(chemxe)%sf(i,j,0) =   0.79_wp
     case default
         if (proc_rank == 0) then
             call s_int_to_str(patch_id, iStr)

src/common/m_chemistry.fpp

@@ -248,12 +248,16 @@ contains
                             rho_L = q_prim_qp(1)%sf(x, y, z)
                             rho_R = q_prim_qp(1)%sf(x + offsets(1), y + offsets(2), z + offsets(3))
 
-                            T_L = P_L/rho_L/Rgas_L
+                            T_L = P_L/rho_L/Rgas_L      
                             T_R = P_R/rho_R/Rgas_R
 
                             rho_cell = 0.5_wp*(rho_L + rho_R)
                             dT_dxi = (T_R - T_L)/grid_spacing
 
+                            if (x .eq. 200 .and. idir .eq. 2) then
+                            print *, dT_dxi, y , T_L, T_R
+                            end if
+
                             ! Get transport properties
                             call get_species_mass_diffusivities_mixavg(P_L, T_L, Ys_L, mass_diffusivities_mixavg1)
                             call get_species_mass_diffusivities_mixavg(P_R, T_R, Ys_R, mass_diffusivities_mixavg2)