From f0c5f9cfa014d984d0227eadf445e9baa26785f2 Mon Sep 17 00:00:00 2001 From: ChrisZYJ Date: Thu, 26 Sep 2024 02:01:46 -0700 Subject: [PATCH] double shear case --- case/double_shear/case.py | 102 ++++++++++++++++++++++ src/pre_process/include/2dHardcodedIC.fpp | 38 ++++++++ src/simulation/m_weno.fpp | 14 +++ 3 files changed, 154 insertions(+) create mode 100644 case/double_shear/case.py diff --git a/case/double_shear/case.py b/case/double_shear/case.py new file mode 100644 index 0000000000..2d32f68cb3 --- /dev/null +++ b/case/double_shear/case.py @@ -0,0 +1,102 @@ +#!/usr/bin/env python3 + +import math +import json + +# Define non-dimensional parameters based on Table 2 +Ma = 0.1 # Mach number +Re = 10000 # Reynolds number +Pr = 0.73 # Prandtl number +mygamma = 1.4 # Specific heat ratio (renamed to avoid duplication) + +# Dynamic viscosity based on Reynolds number +Mu = 1.0 / Re # Mu = 1/Re = 1/10000 = 0.0001 + +# Pressure based on Mach number +p = 1.0 / (mygamma * Ma**2) # p = 1 / (1.4 * 0.01) = 71.42857142857143 + +# Grid resolution for theta=120 +m = 160 +n = 160 + +# Configuring case dictionary +case_dict = { + # Logistics ================================================ + 'run_time_info' : 'T', + # ========================================================== + + # Computational Domain Parameters ========================== + # Unit square domain for Double Periodic Shear Layer + 'x_domain%beg' : 0.0, + 'x_domain%end' : 1.0, + 'y_domain%beg' : 0.0, + 'y_domain%end' : 1.0, + 'm' : m, # Grid resolution + 'n' : n, + 'p' : 0, + 'dt' : 2.E-04, + 't_step_start' : 0, + 't_step_stop' : 5000, + 't_step_save' : 10, + # ========================================================== + + # Simulation Algorithm Parameters ========================== + 'num_patches' : 1, + 'model_eqns' : 2, # Assuming Euler equations; adjust if different + 'alt_soundspeed' : 'F', + 'num_fluids' : 1, + 'adv_alphan' : 'T', + 'mpp_lim' : 'F', + 'mixture_err' : 'T', + 'time_stepper' : 3, + 'weno_order' : 5, + 'weno_eps' : 1.E-16, + # 'mapped_weno' : 'T', + 'teno' : 'T', + 'teno_CT' : 1.E-6, + 'null_weights' : 'F', + 'mp_weno' : 'F', + 'weno_Re_flux' : 'T', + 'weno_avg' : 'T', + 'riemann_solver' : 2, + 'wave_speeds' : 1, + 'avg_state' : 2, + 'bc_x%beg' : -1, # Periodic boundaries + 'bc_x%end' : -1, + 'bc_y%beg' : -1, + 'bc_y%end' : -1, + # ========================================================== + + # Formatted Database Files Structure Parameters ============ + 'format' : 1, + 'precision' : 2, + 'prim_vars_wrt' : 'T', + 'parallel_io' : 'T', + 'fd_order' : 2, + 'omega_wrt(3)' : 'T', + # ========================================================== + + # Patch 1: Double Periodic Shear Layer ====================== + 'patch_icpp(1)%hcid' : 205, # Set to Case 205 + 'patch_icpp(1)%geometry' : 7, # 2D hard-coded + 'patch_icpp(1)%x_centroid' : 0.5, # Centroid at (0.5, 0.5) for periodic domain + 'patch_icpp(1)%y_centroid' : 0.5, + 'patch_icpp(1)%length_x' : 1.0, # Domain length in x + 'patch_icpp(1)%length_y' : 1.0, # Domain length in y + 'patch_icpp(1)%vel(1)' : 0.0, # Dummy + 'patch_icpp(1)%vel(2)' : 0.0, # Dummy + 'patch_icpp(1)%pres' : p, # Pressure based on Mach number + # Single-phase flow parameters + 'patch_icpp(1)%alpha_rho(1)' : 1.0, # rhoH = rhoL = rho = 1.0 + 'patch_icpp(1)%alpha(1)' : 1.0, # Volume fraction alpha = 1.0 + # ========================================================== + + # Fluids Physical Parameters =============================== + 'fluid_pp(1)%gamma' : mygamma, # Specific heat ratio (renamed) + 'fluid_pp(1)%pi_inf' : 0.0, + 'fluid_pp(1)%Re(1)' : Re, # Shear viscosity: Re + # ========================================================== +} + +# Output the case dictionary as JSON +print(json.dumps(case_dict, indent=4)) diff --git a/src/pre_process/include/2dHardcodedIC.fpp b/src/pre_process/include/2dHardcodedIC.fpp index 7d6ee9602b..5498ec53cd 100644 --- a/src/pre_process/include/2dHardcodedIC.fpp +++ b/src/pre_process/include/2dHardcodedIC.fpp @@ -5,6 +5,12 @@ real(kind(0d0)) :: rhoH, rhoL, pRef, pInt, h, lam, wl, amp, intH, alph + ! Define non-dimensional parameters for case 205 + real(kind(0d0)) :: Ma, Re, Pr, mygamma + real(kind(0d0)) :: theta + real(kind(0d0)) :: p, rho, u, v + real(kind(0d0)) :: E + eps = 1e-9 #:enddef @@ -100,6 +106,38 @@ q_prim_vf(E_idx)%sf(i, j, 0) = pInt + rhoL*9.81*(intH - y_cc(j)) end if + case (205) ! Double Periodic Shear Layer + + Ma = 0.1 + Re = 10000d0 + Pr = 0.73 + mygamma = 1.4d0 + + ! Compute pressure based on Mach number + p = 1.0d0 / (mygamma * Ma**2) + + ! Define theta (only theta=120 is needed) + theta = 120.0d0 + + ! Initialize velocity components based on y-coordinate + if (y_cc(j) <= 0.5d0) then + u = tanh(theta * (y_cc(j) - 0.25d0)) + else + u = tanh(theta * (0.75d0 - y_cc(j))) + end if + + ! Initialize v-component based on x-coordinate (theta=120) + v = 0.05d0 * sin(2.0d0 * pi * x_cc(i)) + + ! Set momentum components + q_prim_vf(momxb)%sf(i, j, 0) = u + q_prim_vf(momxe)%sf(i, j, 0) = v + + ! Compute energy based on pressure and kinetic energy + E = p / (mygamma - 1.0d0) + 0.5d0 * (u**2 + v**2) + q_prim_vf(E_idx)%sf(i, j, 0) = E + + case default if (proc_rank == 0) then call s_int_to_str(patch_id, iStr) diff --git a/src/simulation/m_weno.fpp b/src/simulation/m_weno.fpp index 6a05b9d9c3..974a2aa9e9 100644 --- a/src/simulation/m_weno.fpp +++ b/src/simulation/m_weno.fpp @@ -666,6 +666,20 @@ contains do j = is1_weno%beg, is1_weno%end !$acc loop seq do i = 1, v_size + ! if (i >= mom_idx%beg .and. i <= mom_idx%end .and. i - mom_idx%beg + 1 /= weno_dir) then ! Use higher-order central difference for vL and vR if i != j for du_i/dx_j + ! ! Compute higher-order central differences for vL and vR using 5-point stencil + ! vL_rs_vf_${XYZ}$(j, k, l, i) = ( v_rs_ws_${XYZ}$(j-2, k, l, i) & + ! - 8d0*v_rs_ws_${XYZ}$(j-1, k, l, i) & + ! + 8d0*v_rs_ws_${XYZ}$(j+1, k, l, i) & + ! - v_rs_ws_${XYZ}$(j+2, k, l, i)) / 12d0 + + ! vR_rs_vf_${XYZ}$(j, k, l, i) = ( v_rs_ws_${XYZ}$(j-1, k, l, i) & + ! - 8d0*v_rs_ws_${XYZ}$(j, k, l, i) & + ! + 8d0*v_rs_ws_${XYZ}$(j+2, k, l, i) & + ! - v_rs_ws_${XYZ}$(j+3, k, l, i)) / 12d0 + ! cycle + ! end if + ! reconstruct from left side dvd(1) = v_rs_ws_${XYZ}$ (j + 2, k, l, i) &