Skip to content

Commit

Permalink
double shear case
Browse files Browse the repository at this point in the history
  • Loading branch information
@ChrisZYJ
ChrisZYJ committed Sep 26, 2024
1 parent 42cfdf0 commit f0c5f9c
Show file tree
Hide file tree
Showing 3 changed files with 154 additions and 0 deletions.
102 changes: 102 additions & 0 deletions case/double_shear/case.py
View file
Original file line number Diff line number Diff line change
@@ -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))
38 changes: 38 additions & 0 deletions src/pre_process/include/2dHardcodedIC.fpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down
14 changes: 14 additions & 0 deletions src/simulation/m_weno.fpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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) &
Expand Down

0 comments on commit f0c5f9c

Please sign in to comment.