diff --git a/.typos.toml b/.typos.toml
index d32ea7e534..6216993cb7 100644
--- a/.typos.toml
+++ b/.typos.toml
@@ -17,6 +17,7 @@ gam = "gam"
Gam = "Gam"
strang = "strang"
Strang = "Strang"
+TKE = "TKE"
[files]
extend-exclude = ["docs/documentation/references*", "tests/", "toolchain/cce_simulation_workgroup_256.sh"]
diff --git a/docs/documentation/case.md b/docs/documentation/case.md
index d000aeb467..7d55e301b3 100644
--- a/docs/documentation/case.md
+++ b/docs/documentation/case.md
@@ -825,7 +825,6 @@ When ``polytropic = 'F'``, the gas compression is modeled as non-polytropic due
| `mixlayer_vel_profile` | Logical | Set the mean streamwise velocity to hyperbolic tangent profile |
| `mixlayer_vel_coef` | Real | Coefficient for the hyperbolic tangent profile of a mixing layer |
| `mixlayer_perturb` | Logical | Perturb the initial velocity field by instability waves |
-| `mixlayer_domain` | Real | Domain size of a mixing layer for the linear stability analysis |
The table lists velocity field parameters.
The parameters are optionally used to define initial velocity profiles and perturbations.
@@ -840,15 +839,13 @@ The parameters are optionally used to define initial velocity profiles and pertu
- `perturb_sph_fluid` specifies the fluid component whose partial density is to be perturbed.
-- `mixlayer_vel_profile` activates setting the mean streamwise velocity to a hyperbolic tangent profile. This option works only for 2D and 3D cases.
+- `mixlayer_vel_profile` activates setting the mean streamwise velocity to a hyperbolic tangent profile. This option works only for `n > 0`.
- `mixlayer_vel_coef` is a parameter for the hyperbolic tangent profile of a mixing layer when `mixlayer_vel_profile = 'T'`. The mean streamwise velocity profile is given as:
-$$ u = patch\_icpp(1)\%vel(1) * tanh(y\_cc * mixlayer\_vel\_profile) $$
+$$ u = \mbox{patch\_icpp(1)\%vel(1)} * \tanh( y_{cc} * \mbox{mixlayer\_vel\_coef}) $$
-- `mixlayer_perturb` activates the perturbation of initial velocity by instability waves obtained from linear stability analysis for a mixing layer with hyperbolic tangent mean streamwise velocity profile. This option only works for `n > 0`, `bc_y%[beg,end] = -6`, `num_fluids = 1`, `model_eqns = 2` and `mixlayer_vel_profile = 'T'`.
-
-- `mixlayer_domain` defines the domain size to compute spatial eigenvalues of the linear instability analysis when `mixlayer_perturb = 'T'`. For example, the spatial eigenvalue in `x` direction in 2D problem will be $2 \pi \alpha / (mixlayer\_domain*patch\_icpp(1)\%length\_y)$ for $\alpha = 1$, $2$ and $4$.
+- `mixlayer_perturb` activates the velocity perturbation for a temporal mixing layer with hyperbolic tangent mean streamwise velocity profile, using an inverter version of the spectrum-based synthetic turbulence generation method proposed by Guo et al. (2023, JFM). This option only works for `p > 0` and `mixlayer_vel_profile = 'T'`.
### 11. Phase Change Model
| Parameter | Type | Description |
diff --git a/docs/documentation/references.md b/docs/documentation/references.md
index 7f2479048a..f97337269e 100644
--- a/docs/documentation/references.md
+++ b/docs/documentation/references.md
@@ -20,6 +20,8 @@
- Gottlieb, S. and Shu, C.-W. (1998). Total variation diminishing runge-kutta schemes. Mathematics of computation of the American Mathematical Society, 67(221):73–85.
+- Guo, H., Jiang, P., Ye, L., & Zhu, Y. (2023). An efficient and low-divergence method for generating inhomogeneous and anisotropic turbulence with arbitrary spectra. Journal of Fluid Mechanics, 970, A2.
+
- Henrick, A. K., Aslam, T. D., and Powers, J. M. (2005). Mapped weighted essentially nonoscillatory schemes: achieving optimal order near critical points. Journal of Computational Physics, 207(2):542–567.
- Borges, R., Carmona, M., Costa, B., and Don, W. S. (2008). An improved weighted essentially non-oscillatory scheme for hyperbolic conservation laws. Journal of computational physics, 227(6):3191–3211.
@@ -34,6 +36,8 @@
- Meng, J. C. C. (2016). Numerical simulations of droplet aerobreakup. PhD thesis, California Institute of Technology.
+- Michalke, A. (1964). On the inviscid instability of the hyperbolictangent velocity profile. Journal of Fluid Mechanics, 19(4), 543-556.
+
- Pirozzoli, S., and Colonius, T. (2013). Generalized characteristic relaxation boundary conditions for unsteady compressible flow simulations. Journal of Computational Physics, 248:109-126.
- Preston, A., Colonius, T., and Brennen, C. (2007). A reduced-order model of diffusive effects on the dynamics of bubbles. Physics of Fluids, 19(12):123302.
diff --git a/examples/2D_mixing_artificial_Ma/case.py b/examples/2D_mixing_artificial_Ma/case.py
index 825ce9a24a..e9d926df71 100644
--- a/examples/2D_mixing_artificial_Ma/case.py
+++ b/examples/2D_mixing_artificial_Ma/case.py
@@ -108,8 +108,6 @@
# Mixing layer
"mixlayer_vel_profile": "T",
"mixlayer_vel_coef": 1.0,
- "mixlayer_domain": 1.0,
- "mixlayer_perturb": "T",
# Artificial Mach number
"pi_fac": pi_fac,
# Fluids Physical Parameters
diff --git a/examples/3D_turb_mixing/case.py b/examples/3D_turb_mixing/case.py
index 386ed621cf..4055e57b28 100644
--- a/examples/3D_turb_mixing/case.py
+++ b/examples/3D_turb_mixing/case.py
@@ -4,26 +4,25 @@
# SURROUNDING FLOW
# Nondimensional parameters
-Re0 = 50.0 # Reynolds number
+Re0 = 320.0 # Reynolds number
M0 = 0.2 # Mach number
# Fluid properties
gamma = 1.4
-pi_inf = 0.0
# Free stream velocity & pressure
u0 = 1.0
pres0 = 1.0 / (gamma * M0**2)
# Domain size
-Lx = 59.0
-Ly = 59.0
-Lz = 59.0
+Lx = 160.0
+Ly = 160.0
+Lz = 80.0
# Number of grid cells
-Nx = 191
-Ny = 191
-Nz = 191
+Nx = 1023
+Ny = 1023
+Nz = 511
# Grid spacing
dx = Lx / float(Nx)
@@ -31,12 +30,12 @@
dz = Lz / float(Nz)
# Time advancement
-cfl = 0.1
-T = 100.0
-dt = cfl * dx / u0
+cfl = 0.5
+T = 150.0
+dt = cfl * dx / (u0 / M0 + 1)
Ntfinal = int(T / dt)
Ntstart = 0
-Nfiles = 100
+Nfiles = 30
t_save = int(math.ceil((Ntfinal - 0) / float(Nfiles)))
Nt = t_save * Nfiles
t_step_start = Ntstart
@@ -53,6 +52,11 @@
"x_domain%end": Lx,
"y_domain%beg": -Ly / 2.0,
"y_domain%end": Ly / 2.0,
+ "stretch_y": "T",
+ "a_y": 2,
+ "y_a": -0.3 * Ly,
+ "y_b": 0.3 * Ly,
+ "loops_y": 2,
"z_domain%beg": 0.0,
"z_domain%end": Lz,
"m": Nx,
@@ -68,10 +72,9 @@
"num_fluids": 1,
"time_stepper": 3,
"weno_order": 5,
- "weno_eps": 1.0e-16,
- "weno_Re_flux": "T",
- "weno_avg": "T",
- "mapped_weno": "T",
+ "weno_eps": 1.0e-40,
+ "weno_Re_flux": "F",
+ "wenoz": "T",
"riemann_solver": 2,
"wave_speeds": 1,
"avg_state": 2,
@@ -85,7 +88,7 @@
# Formatted Database Files Structure Parameters
"format": 1,
"precision": 2,
- "cons_vars_wrt": "T",
+ "cons_vars_wrt": "F",
"prim_vars_wrt": "T",
"parallel_io": "T",
"fd_order": 1,
@@ -99,7 +102,7 @@
"patch_icpp(1)%y_centroid": 0.0,
"patch_icpp(1)%z_centroid": Lz / 2.0,
"patch_icpp(1)%length_x": Lx,
- "patch_icpp(1)%length_y": Ly,
+ "patch_icpp(1)%length_y": 1e6,
"patch_icpp(1)%length_z": Lz,
"patch_icpp(1)%alpha_rho(1)": 1.0,
"patch_icpp(1)%alpha(1)": 1.0,
@@ -109,8 +112,6 @@
"patch_icpp(1)%pres": pres0,
# Mixing layer
"mixlayer_vel_profile": "T",
- "mixlayer_vel_coef": 1.0,
- "mixlayer_domain": 1.0,
"mixlayer_perturb": "T",
# Fluids Physical Parameters
# Surrounding liquid
diff --git a/examples/3D_turb_mixing/turbulence_stat/README.md b/examples/3D_turb_mixing/turbulence_stat/README.md
new file mode 100644
index 0000000000..490da952bf
--- /dev/null
+++ b/examples/3D_turb_mixing/turbulence_stat/README.md
@@ -0,0 +1,23 @@
+# Temporally evolving turbulent mixing layer (3D)
+
+Reference:
+> Bell, J. H., & Mehta, R. D. (1990). Development of a two-stream mixing layer from tripped and untripped boundary layers. AIAA Journal, 28(12), 2034-2042.
+> Rogers, M. M., & Moser, R. D. (1994). Direct simulation of a self‐similar turbulent mixing layer. Physics of Fluids, 6(2), 903-923.
+> Pantano, C., & Sarkar, S. (2002). A study of compressibility effects in the high-speed turbulent shear layer using direct simulation. Journal of Fluid Mechanics, 451, 329-371.
+> Vaghefi, S. N. S. (2014). Simulation and modeling of compressible turbulent mixing layer. State University of New York at Buffalo.
+> Wang, X., Wang, J., & Chen, S. (2022). Compressibility effects on statistics and coherent structures of compressible turbulent mixing layers. Journal of Fluid Mechanics, 947, A38.
+
+## Description
+This directory (`turbulence_stat`) contains sub-directories and Matlab scripts for computing turbulence statistics of 3D temporally evolving turbulent mixing layer as described below:
+
+Files:
+- `set_user_inputs.m`: User input parameters are defined in this file.
+- `run_turbulence.m`: Turbulence statistics (growth rate, Reynolds stress, TKE budget) are computed and plots are generated in this file.
+- `average_tke_over_self_similar.m`: Averaging over self-similar period of TKE budget is performed in this file.
+- `submit_batch_job_*.sh`: Example scripts for batch job submission on DoD Carpenter (PBS system) and NCSA Delta (Slurm system)
+
+Directories:
+- `log`: Log files from batch job are stored in this directory.
+- `reference_data`: Data from reference papers are provided in this directory.
+- `variables`: User inputs from `set_user_inputs.m` are saved as Matlab data (`user_inputs.mat`) in this directory.
+- `results`: Outputs are stored in sub-directories in the directory.
diff --git a/examples/3D_turb_mixing/turbulence_stat/average_tke_over_self_similar.m b/examples/3D_turb_mixing/turbulence_stat/average_tke_over_self_similar.m
new file mode 100644
index 0000000000..ed3b2e9547
--- /dev/null
+++ b/examples/3D_turb_mixing/turbulence_stat/average_tke_over_self_similar.m
@@ -0,0 +1,118 @@
+close all; clear all;
+
+
+% Setup
+disp("Start average_tke_over_self_similar ..."); tic;
+set_user_inputs(); load variables/user_inputs.mat;
+
+ybeg = -5; yend = 5; ny = 101;
+y = linspace(ybeg,yend,ny);
+
+% Array
+T0_averaged = zeros(ny,1);
+P_averaged = zeros(ny,1);
+D_averaged = zeros(ny,1);
+
+% Compute averaged TKE budget
+for q = 1:Nfiles
+ load("results/tke_budget_data/tstep_"+string(timesteps(q))+".mat");
+
+ % Normalization
+ T0 = T0 / (8/mth); % T / (Delta U^3 / mth)
+ P = P / (8/mth); % P / (Delta U^3 / mth)
+ D = D / (8/mth); % D / (Delta U^3 / mth)
+
+ % Interpolation
+ i_start = 1;
+ for j = 1:ny
+ for i = i_start:length(y_norm_mth) - 1
+ if (y_norm_mth(i) <= y(j) && y_norm_mth(i+1) > y(j))
+ T0_averaged(j) = T0_averaged(j) + ((T0(i+1) - T0(i))/(y_norm_mth(i+1) - y_norm_mth(i))*(y(j) - y_norm_mth(i)) + T0(i))/Nfiles;
+ P_averaged(j) = P_averaged(j) + ((P(i+1) - P(i))/(y_norm_mth(i+1) - y_norm_mth(i))*(y(j) - y_norm_mth(i)) + P(i))/Nfiles;
+ D_averaged(j) = D_averaged(j) + ((D(i+1) - D(i))/(y_norm_mth(i+1) - y_norm_mth(i))*(y(j) - y_norm_mth(i)) + D(i))/Nfiles;
+ i_start = i;
+ break;
+ end
+ end
+ end
+end
+
+T_averaged = f_compute_derivative_1d(T0_averaged,y*mth);
+
+% Plot
+plot_tke_budget(T_averaged, P_averaged, D_averaged, y, mth);
+
+disp("End of program"); toc;
+
+%% FUNCTIONS
+% Compute the wall-normal derivative of a discretized function, fun(y)
+function dfunds = f_compute_derivative_1d(fun,s)
+
+ dfunds = zeros(size(fun)); % initialize discrete derivative vector
+
+ % Compute one-sided derivative at the bottom boundary
+ dfunds(1) = (fun(2) - fun(1)) / (s(2) - s(1));
+
+ % Compute one-sided derivative at the top boundary
+ dfunds(end) = (fun(end) - fun(end-1)) / (s(end) - s(end-1));
+
+ % Compute two-sided derivatives for interior points
+ for i = 2:length(s)-1
+ dfunds(i) = (fun(i+1) - fun(i-1)) / (s(i+1) - s(i-1));
+ end
+end
+
+% Plot TKE budget
+function plot_tke_budget(T, P, D, y_norm_mth, mth)
+
+ load variables/user_inputs.mat;
+ load reference_data/reference.mat;
+
+ % Plot
+ f1 = figure("DefaultAxesFontSize",18);
+ set(f1,"Position",[200 200 1000 700]);
+
+ % Present
+ h1 = plot([-100 -100],[-100 -100],'-k','LineWidth',2); hold on; grid on;
+ plot(y_norm_mth,T,'-b','LineWidth',2);
+ plot(y_norm_mth,P,'-g','LineWidth',2);
+ plot(y_norm_mth,D,'-r','LineWidth',2);
+ xlim([-5 5]); xticks([-5:1:5]);
+ ylim([-0.002 0.003]);
+ xlabel('$y/\delta_\theta$','interpreter','latex');
+ set(gca,'TickLabelInterpreter','latex');
+
+ % Pantano & Sarkar (2002)
+ h2 = plot([-100 -100],[-100 -100],'ko','LineWidth',2,'MarkerSize',8);
+ plot(p2002_tke_transport(:,1),p2002_tke_transport(:,2),'bo','LineWidth',2,'MarkerSize',8);
+ plot(p2002_tke_production(:,1),p2002_tke_production(:,2),'go','LineWidth',2,'MarkerSize',8);
+ plot(p2002_tke_dissipation(:,1),p2002_tke_dissipation(:,2),'ro','LineWidth',2,'MarkerSize',8);
+
+ % Rogers & Moser (1994)
+ h3 = plot([-100 -100],[-100 -100],'k^','LineWidth',2,'MarkerSize',8);
+ plot(r1994_tke_transport(:,1),r1994_tke_transport(:,2),'b^','LineWidth',2,'MarkerSize',8);
+ plot(r1994_tke_production(:,1),r1994_tke_production(:,2),'g^','LineWidth',2,'MarkerSize',8);
+ plot(r1994_tke_dissipation(:,1),r1994_tke_dissipation(:,2),'r^','LineWidth',2,'MarkerSize',8);
+
+ % Vaghefi (2014)
+ h4 = plot([-100 -100],[-100 -100],'k+','LineWidth',2,'MarkerSize',8);
+ plot(v2014_tke_transport(:,1),v2014_tke_transport(:,2),'b+','LineWidth',2,'MarkerSize',8);
+ plot(v2014_tke_production(:,1),v2014_tke_production(:,2),'g+','LineWidth',2,'MarkerSize',8);
+ plot(v2014_tke_dissipation(:,1),v2014_tke_dissipation(:,2),'r+','LineWidth',2,'MarkerSize',8);
+
+ % Wang et al. (2022)
+ h5 = plot([-100 -100],[-100 -100],'k*','LineWidth',2,'MarkerSize',8);
+ plot(w2022_tke_transport(:,1),w2022_tke_transport(:,2),'b*','LineWidth',2,'MarkerSize',8);
+ plot(w2022_tke_production(:,1),w2022_tke_production(:,2),'g*','LineWidth',2,'MarkerSize',8);
+ plot(w2022_tke_dissipation(:,1),w2022_tke_dissipation(:,2),'r*','LineWidth',2,'MarkerSize',8);
+
+ legend([h1,h2,h3,h4,h5], {"$\mbox{Present}$", ...
+ "$\mbox{Pantano \& Sarkar (2002)}$", ...
+ "$\mbox{Rogers \& Moser (1994)}$", ...
+ "$\mbox{Vaghefi (2014)}$", ...
+ "$\mbox{Wang et al. (2022)}$"}, ...
+ 'interpreter','latex','location','northeast');
+
+ saveas(f1, "results/tke_budget/avg_self_similar","png");
+ close(f1);
+end
diff --git a/examples/3D_turb_mixing/turbulence_stat/reference_data/reference.mat b/examples/3D_turb_mixing/turbulence_stat/reference_data/reference.mat
new file mode 100644
index 0000000000..8a72ad405c
Binary files /dev/null and b/examples/3D_turb_mixing/turbulence_stat/reference_data/reference.mat differ
diff --git a/examples/3D_turb_mixing/turbulence_stat/results/Reynolds_stress/tstep_11520.png b/examples/3D_turb_mixing/turbulence_stat/results/Reynolds_stress/tstep_11520.png
new file mode 100644
index 0000000000..15340f5f58
Binary files /dev/null and b/examples/3D_turb_mixing/turbulence_stat/results/Reynolds_stress/tstep_11520.png differ
diff --git a/examples/3D_turb_mixing/turbulence_stat/results/growth_rate.png b/examples/3D_turb_mixing/turbulence_stat/results/growth_rate.png
new file mode 100644
index 0000000000..37d9f37676
Binary files /dev/null and b/examples/3D_turb_mixing/turbulence_stat/results/growth_rate.png differ
diff --git a/examples/3D_turb_mixing/turbulence_stat/results/tke_budget/avg_self_similar.png b/examples/3D_turb_mixing/turbulence_stat/results/tke_budget/avg_self_similar.png
new file mode 100644
index 0000000000..87e387565f
Binary files /dev/null and b/examples/3D_turb_mixing/turbulence_stat/results/tke_budget/avg_self_similar.png differ
diff --git a/examples/3D_turb_mixing/turbulence_stat/results/tke_budget/tstep_11520.png b/examples/3D_turb_mixing/turbulence_stat/results/tke_budget/tstep_11520.png
new file mode 100644
index 0000000000..099c12e8b0
Binary files /dev/null and b/examples/3D_turb_mixing/turbulence_stat/results/tke_budget/tstep_11520.png differ
diff --git a/examples/3D_turb_mixing/turbulence_stat/run_turbulence.m b/examples/3D_turb_mixing/turbulence_stat/run_turbulence.m
new file mode 100644
index 0000000000..2a6473a5ca
--- /dev/null
+++ b/examples/3D_turb_mixing/turbulence_stat/run_turbulence.m
@@ -0,0 +1,601 @@
+close all; clear all;
+
+%% Initialization
+disp("Initialize run_turbulence ..."); tic;
+% Create directories
+create_directory();
+% Read user inputs
+set_user_inputs(); load variables/user_inputs.mat;
+% Set grids
+create_grid(); load variables/user_inputs.mat;
+% Allocate array for momentum thickness
+mth = zeros(1,Nfiles);
+toc;
+
+%% Loop over timesteps
+for i = 1:Nfiles
+
+ % Read data in conservative form
+ disp(" "); disp("[1/6] Read data ..."); tic;
+ qp = f_read_data(strcat("../restart_data/lustre_",int2str(timesteps(i)),".dat")); toc;
+
+ % Compute mean quantities
+ disp("[2/6] Compute mean and fluctuating quantities ..."); tic;
+ [qp_mean qp_fluc] = f_compute_qp_mean_fluc(qp); toc;
+
+ % Compute derivatives
+ disp("[3/6] Compute velocity derivatives ..."); tic;
+ [dvel_ds dvelmean_dy dvelfluc_ds] = ...
+ f_compute_vel_derivatives(qp(momxb:momxe,:,:,:), qp_mean(momxb:momxe,:), qp_fluc(momxb:momxe,:,:,:)); toc;
+
+ % Compute mixing layer thickness
+ disp("[4/6] Compute mixing layer thickness ..."); tic;
+ [vth mth(i) y_norm_mth y_norm_vth] = f_compute_mixing_layer_thickness(qp_mean,dvelmean_dy); toc;
+
+ % Compute Reynolds stress
+ disp("[5/6] Compute Reynolds stress ..."); tic;
+ if (Reynolds_stress)
+ f_compute_Reynolds_stress(qp(1,:,:,:), qp_fluc(momxb:momxe,:,:,:), y_norm_vth, timesteps(i));
+ else
+ disp("skipped");
+ end
+ toc;
+
+ % Compute TKE budget
+ disp("[6/6] Compute TKE budget ..."); tic;
+ if (tke_budget)
+ f_compute_tke_budget(dvel_ds, dvelmean_dy, dvelfluc_ds, ...
+ squeeze(qp(1,:,:,:)), qp_fluc(momxb:momxe,:,:,:), squeeze(qp(E_idx,:,:,:)), ...
+ y_norm_mth, mth(i), timesteps(i));
+ else
+ disp("skipped");
+ end
+ toc;
+end
+disp("End of loop");
+
+% Compute mixing layer growth rate
+disp("Compute growth rate ..."); tic;
+plot_growth_rate(time,mth); toc;
+
+disp("End of program");
+
+
+%% FUNCTIONS
+% Create directory for saving results
+function create_directory()
+ if exist("variables", 'dir')
+ delete variables/*.mat;
+ else
+ mkdir(strcat("variables"));
+ end
+ if ~exist(strcat("results"), "dir")
+ mkdir(strcat("results"));
+ end
+ if ~exist(strcat("results/tke_budget"), "dir")
+ mkdir(strcat("results/tke_budget"));
+ end
+ if ~exist(strcat("results/tke_budget_data"), "dir")
+ mkdir(strcat("results/tke_budget_data"));
+ end
+ if ~exist(strcat("results/Reynolds_stress"), "dir")
+ mkdir(strcat("results/Reynolds_stress"));
+ end
+ if ~exist(strcat("results/Reynolds_stress_data"), "dir")
+ mkdir(strcat("results/Reynolds_stress_data"));
+ end
+end
+
+% Create grids
+function create_grid()
+
+ load variables/user_inputs.mat;
+
+ % Grid-related parameters
+ mp = m + 1; mpp = mp + 1;
+ np = n + 1; npp = np + 1;
+ pp = p + 1; ppp = pp + 1;
+ r = mp*np*pp;
+
+ % Allocate memory
+ x_cc = zeros(1,mp); x_cb = zeros(1,mpp);
+ y_cc = zeros(1,np); y_cb = zeros(1,npp);
+ z_cc = zeros(1,pp); z_cb = zeros(1,ppp);
+
+ % x-dir
+ for i = 1:mpp
+ x_cb(i) = (i - 1)*Lx/mp;
+ end
+ dx = x_cb(2:mpp) - x_cb(1:mp);
+ x_cc = (x_cb(2:mpp) + x_cb(1:mp))/2;
+
+ % y-dir
+ for i = 1:npp
+ y_cb(i) = (i - 1)*Ly/np - 0.5*Ly;
+ end
+ if (stretch_y == "T")
+ y_cb = y_cb / Ly;
+ y_a = y_a / Ly;
+ y_b = y_b / Ly;
+ for j = 1:loops_y
+ for i = 1:npp
+ y_cb(i) = y_cb(i)/a_y* ...
+ (a_y + log(cosh(a_y*(y_cb(i) - y_a))) ...
+ + log(cosh(a_y*(y_cb(i) - y_b))) ...
+ - 2*log(cosh(a_y*(y_b - y_a)/2)));
+ end
+ end
+ y_cb = y_cb*Ly;
+ end
+ dy = y_cb(2:npp) - y_cb(1:np);
+ y_cc = (y_cb(2:npp) + y_cb(1:np))/2;
+
+ % z-dir
+ for i = 1:ppp
+ z_cb(i) = (i - 1)*Lz/pp;
+ end
+ dz = z_cb(2:ppp) - z_cb(1:pp);
+ z_cc = (z_cb(2:ppp) + z_cb(1:pp))/2;
+
+ % Save variables
+ save variables/user_inputs.mat
+end
+
+% Read lustre files
+function qp = f_read_data(filename)
+
+ load variables/user_inputs.mat;
+
+ % Read data
+ disp(filename);
+ fileID = fopen(filename,'r');
+ A = fread(fileID,'double');
+ fclose(fileID);
+
+ % Reassign density & velocity components
+ qc = permute(reshape(A, mp, np, pp, sys_size),[4 1 2 3]);
+
+ % Convert conservative variables to primitive variables
+ qp = f_convert_cons_to_prim(qc);
+end
+
+% Convert conservative variables to primitive variables
+function qp = f_convert_cons_to_prim(qc)
+
+ load variables/user_inputs.mat;
+
+ qp = zeros(size(qc));
+
+ % density
+ qp(contxb:contxe,:,:,:) = qc(contxb:contxe,:,:,:);
+
+ dyn_p = 0;
+ for i = momxb:momxe
+ qp(i,:,:,:) = qc(i,:,:,:) ./ qc(1,:,:,:); % Velocity
+ dyn_p = dyn_p + 0.5*qc(i,:,:,:).*qp(i,:,:,:); % Dynamic pressure
+ end
+
+ % pressure
+ qp(E_idx,:,:,:) = (1/Gamma)*(qc(E_idx,:,:,:) - dyn_p - Pi_inf);
+
+end
+
+% Compute mean and fluctuating flow field data
+function [qp_mean qp_fluc] = f_compute_qp_mean_fluc(qp)
+
+ load variables/user_inputs.mat;
+
+ % Compute mean quantities
+ qp_mean = zeros(E_idx,np);
+
+ % mean rho
+ qp_mean(contxb:contxe,:) = squeeze(mean(qp(contxb:contxe,:,:,:),[2 4]));
+
+ % favre-averaged velocity
+ for i = momxb:momxe
+ qp_mean(i,:) = squeeze(mean(qp(1,:,:,:).*qp(i,:,:,:),[2 4])./mean(qp(1,:,:,:),[2 4]));
+ end
+
+ % mean pressure
+ qp_mean(E_idx,:) = squeeze(mean(qp(E_idx,:,:,:),[2 4]));
+
+ % fluctuation
+ qp_fluc = qp(1:E_idx,:,:,:) - permute(repmat(qp_mean, [1 1 mp pp]), [1 3 2 4]);
+end
+
+% Compute velocity derivatives
+function [dvel_ds dvelmean_dy dvelfluc_ds] = f_compute_vel_derivatives(vel, vel_mean, vel_fluc)
+
+ load variables/user_inputs.mat;
+
+ dvelmean_dy = zeros(3,np);
+ dvel_ds = zeros(3,3,mp,np,pp);
+ dvelfluc_ds = zeros(3,3,mp,np,pp);
+
+ % favre-averaged velocity derivatives
+ dvelmean_dy(1,:) = f_compute_derivative_1d(vel_mean(1,:),y_cc);
+ dvelmean_dy(2,:) = f_compute_derivative_1d(vel_mean(2,:),y_cc);
+ dvelmean_dy(3,:) = f_compute_derivative_1d(vel_mean(3,:),y_cc);
+
+ if (tke_budget)
+ % Compute velocity derivatives
+ vel1 = squeeze(vel(1,:,:,:));
+ vel2 = squeeze(vel(2,:,:,:));
+ vel3 = squeeze(vel(3,:,:,:));
+
+ dvel_ds(1,1,:,:,:) = f_compute_derivative_3d(vel1,x_cc,1);
+ dvel_ds(2,1,:,:,:) = f_compute_derivative_3d(vel2,x_cc,1);
+ dvel_ds(3,1,:,:,:) = f_compute_derivative_3d(vel3,x_cc,1);
+
+ dvel_ds(1,2,:,:,:) = f_compute_derivative_3d(vel1,y_cc,2);
+ dvel_ds(2,2,:,:,:) = f_compute_derivative_3d(vel2,y_cc,2);
+ dvel_ds(3,2,:,:,:) = f_compute_derivative_3d(vel3,y_cc,2);
+
+ dvel_ds(1,3,:,:,:) = f_compute_derivative_3d(vel1,z_cc,3);
+ dvel_ds(2,3,:,:,:) = f_compute_derivative_3d(vel2,z_cc,3);
+ dvel_ds(3,3,:,:,:) = f_compute_derivative_3d(vel3,z_cc,3);
+
+ % fluctuating velocity derivatives
+ vel_fluc1 = squeeze(vel_fluc(1,:,:,:));
+ vel_fluc2 = squeeze(vel_fluc(2,:,:,:));
+ vel_fluc3 = squeeze(vel_fluc(3,:,:,:));
+
+ dvelfluc_ds(1,1,:,:,:) = f_compute_derivative_3d(vel_fluc1,x_cc,1);
+ dvelfluc_ds(2,1,:,:,:) = f_compute_derivative_3d(vel_fluc2,x_cc,1);
+ dvelfluc_ds(3,1,:,:,:) = f_compute_derivative_3d(vel_fluc3,x_cc,1);
+
+ dvelfluc_ds(1,2,:,:,:) = f_compute_derivative_3d(vel_fluc1,y_cc,2);
+ dvelfluc_ds(2,2,:,:,:) = f_compute_derivative_3d(vel_fluc2,y_cc,2);
+ dvelfluc_ds(3,2,:,:,:) = f_compute_derivative_3d(vel_fluc3,y_cc,2);
+
+ dvelfluc_ds(1,3,:,:,:) = f_compute_derivative_3d(vel_fluc1,z_cc,3);
+ dvelfluc_ds(2,3,:,:,:) = f_compute_derivative_3d(vel_fluc2,z_cc,3);
+ dvelfluc_ds(3,3,:,:,:) = f_compute_derivative_3d(vel_fluc3,z_cc,3);
+ end
+end
+
+% Compute mixing layer thickness
+function [vth mth y_norm_mth y_norm_vth] = f_compute_mixing_layer_thickness(qp_mean, dvelmean_dy)
+
+ load variables/user_inputs.mat;
+
+ % Compute vorticity thickness
+ vth = 2 / max(abs(dvelmean_dy(1,:)),[],"all");
+ disp("vth: "+num2str(vth));
+
+ % Compute momentum thickness
+ f = qp_mean(1,:) .* (1 - qp_mean(momxb, :)) .* (1 + qp_mean(momxb, :)) / 4;
+ mth = trapz(y_cc,f);
+ disp("mth: "+num2str(mth));
+
+ % Compute mth-normalized y value
+ y_norm_mth = y_cc/mth;
+ y_norm_vth = y_cc/vth;
+
+end
+
+% Compute Reynolds stress
+function f_compute_Reynolds_stress(rho, vel_fluc, y_norm_vth, timestep)
+
+ rho = squeeze(rho);
+ vel_fluc1 = squeeze(vel_fluc(1, :, :, :));
+ vel_fluc2 = squeeze(vel_fluc(2, :, :, :));
+ vel_fluc3 = squeeze(vel_fluc(3, :, :, :));
+
+ % Compute Reynolds stress
+ ruu = mean(rho .* vel_fluc1.^2, [1 3]);
+ rvv = mean(rho .* vel_fluc2.^2, [1 3]);
+ rww = mean(rho .* vel_fluc3.^2, [1 3]);
+ ruv = mean(rho .* vel_fluc1.*vel_fluc2, [1 3]);
+
+ % Plot Reynolds stress
+ plot_Reynolds_stress(ruu, rvv, rww, ruv, y_norm_vth, timestep);
+
+ % Save data
+ save("results/Reynolds_stress_data/tstep_"+string(timestep)+".mat","y_norm_vth","ruu","rvv","rww","ruv");
+
+end
+
+% Compute strain rate tensor
+function [SS S11 S12 S13 S22 S23 S33] = f_compute_strain_rate_tensor(dvel_ds);
+ S11 = squeeze(dvel_ds(1,1,:,:,:));
+ S12 = squeeze(0.5*(dvel_ds(1,2,:,:,:) + dvel_ds(2,1,:,:,:)));
+ S13 = squeeze(0.5*(dvel_ds(1,3,:,:,:) + dvel_ds(3,1,:,:,:)));
+ S22 = squeeze(dvel_ds(2,2,:,:,:));
+ S23 = squeeze(0.5*(dvel_ds(2,3,:,:,:) + dvel_ds(3,2,:,:,:)));
+ S33 = squeeze(dvel_ds(3,3,:,:,:));
+
+ SS = S11.^2 + 2*S12.^2 + 2*S13.^2 + S22.^2 + 2*S23.^2 + S33.^2;
+end
+
+% Compute TKE budget
+function f_compute_tke_budget(dvel_ds, dvelmean_dy, dvelfluc_ds, rho, vel_fluc, pres_fluc, y_norm_mth, mth, timestep)
+
+ load variables/user_inputs.mat;
+
+ % Compute strain-rate tensor
+ [SS S11 S12 S13 S22 S23 S33] = f_compute_strain_rate_tensor(dvel_ds);
+ S21 = S12; S31 = S13; S32 = S23;
+
+ % Viscous stress tensor
+ divU = squeeze(dvel_ds(1,1,:,:,:) + dvel_ds(2,2,:,:,:) + dvel_ds(3,3,:,:,:));
+ sig = zeros(3,3,mp,np,pp);
+ sig(1,1,:,:,:) = 2/Re*(S11 - divU/3);
+ sig(1,2,:,:,:) = 2/Re*(S12);
+ sig(1,3,:,:,:) = 2/Re*(S13);
+ sig(2,1,:,:,:) = 2/Re*(S21);
+ sig(2,2,:,:,:) = 2/Re*(S22 - divU/3);
+ sig(2,3,:,:,:) = 2/Re*(S23);
+ sig(3,1,:,:,:) = 2/Re*(S31);
+ sig(3,2,:,:,:) = 2/Re*(S32);
+ sig(3,3,:,:,:) = 2/Re*(S33 - divU/3);
+
+ sigmean = zeros(3,3,np);
+ for i = 1:3
+ for j = 1:3
+ sigmean(i,j,:) = squeeze(mean(sig(i,j,:,:,:),[3 5]));
+ end
+ end
+ sigfluc = sig - permute(repmat(sigmean, [1 1 1 mp pp]), [1 2 4 3 5]);
+
+ % Transport (T)
+ T1 = -0.5*mean(...
+ rho.*squeeze(vel_fluc(1,:,:,:).*vel_fluc(1,:,:,:).*vel_fluc(2,:,:,:)) ...
+ + rho.*squeeze(vel_fluc(2,:,:,:).*vel_fluc(2,:,:,:).*vel_fluc(2,:,:,:)) ...
+ + rho.*squeeze(vel_fluc(3,:,:,:).*vel_fluc(3,:,:,:).*vel_fluc(2,:,:,:)),[1 3]);
+ T2 = -mean(pres_fluc.*squeeze(vel_fluc(2,:,:,:)),[1 3]);
+ T3 = mean(...
+ squeeze(sigfluc(1,2,:,:,:)).*squeeze(vel_fluc(1,:,:,:)) ...
+ + squeeze(sigfluc(2,2,:,:,:)).*squeeze(vel_fluc(2,:,:,:)) ...
+ + squeeze(sigfluc(3,2,:,:,:)).*squeeze(vel_fluc(3,:,:,:)),[1 3]);
+ T0 = T1 + T2 + T3;
+ T = f_compute_derivative_1d(T0,y_cc);
+
+ % Turbulent production (P)
+ P1 = -mean(rho.*squeeze(vel_fluc(1,:,:,:).*vel_fluc(2,:,:,:)),[1 3]).*squeeze(dvelmean_dy(1,:));
+ P2 = -mean(rho.*squeeze(vel_fluc(2,:,:,:).*vel_fluc(2,:,:,:)),[1 3]).*squeeze(dvelmean_dy(2,:));
+ P3 = -mean(rho.*squeeze(vel_fluc(3,:,:,:).*vel_fluc(2,:,:,:)),[1 3]).*squeeze(dvelmean_dy(3,:));
+ P = P1 + P2 + P3;
+
+ % Dissipation (D)
+ D = -mean(squeeze(sigfluc(1,1,:,:,:).*dvelfluc_ds(1,1,:,:,:) ...
+ + sigfluc(2,1,:,:,:).*dvelfluc_ds(2,1,:,:,:) ...
+ + sigfluc(3,1,:,:,:).*dvelfluc_ds(3,1,:,:,:) ...
+ + sigfluc(1,2,:,:,:).*dvelfluc_ds(1,2,:,:,:) ...
+ + sigfluc(2,2,:,:,:).*dvelfluc_ds(2,2,:,:,:) ...
+ + sigfluc(3,2,:,:,:).*dvelfluc_ds(3,2,:,:,:) ...
+ + sigfluc(1,3,:,:,:).*dvelfluc_ds(1,3,:,:,:) ...
+ + sigfluc(2,3,:,:,:).*dvelfluc_ds(2,3,:,:,:) ...
+ + sigfluc(3,3,:,:,:).*dvelfluc_ds(3,3,:,:,:)), [1 3]);
+
+ % Plot TKE budget
+ plot_tke_budget(T, P, D, y_norm_mth, mth, timestep);
+
+ % Save data
+ save("results/tke_budget_data/tstep_"+string(timestep)+".mat","y_norm_mth","mth","T0","P","D");
+
+end
+
+% Compute the wall-normal derivative of a discretized function, fun(y)
+function dfunds = f_compute_derivative_1d(fun,s)
+
+ dfunds = zeros(size(fun)); % initialize discrete derivative vector
+
+ for j = 1:length(s) % for each index in s
+ % if at bottom boundary, use 1-sided derivative
+ dfunds(1) = (fun(2) - fun(1)) / (s(2) - s(1));
+ % if at top boundary, use 1-sided derivative
+ dfunds(end) = (fun(end) - fun(end-1)) / (s(end) - s(end-1));
+ % otherwise, if in the interior of the domain, 2-sided derivative
+ for i = 2:length(s)-1
+ dfunds(i) = (fun(i+1) - fun(i-1)) / (s(i+1) - s(i-1));
+ end
+ end
+end
+
+% Compute the wall-normal derivative of a discretized function, fun(y)
+function dfunds = f_compute_derivative_3d(fun,s,dir)
+
+ dfunds = zeros(size(fun)); % initialize discrete derivative vector
+
+ if (dir == 1)
+ % Forward difference at start
+ dfunds(1,:,:) = (fun(2,:,:) - fun(1,:,:)) / (s(2) - s(1));
+ % Backward difference at end
+ dfunds(end,:,:) = (fun(end,:,:) - fun(end-1,:,:)) / (s(end) - s(end-1));
+ % Central difference for interior
+ for i = 2:length(s)-1
+ dfunds(i,:,:) = (fun(i+1,:,:) - fun(i-1,:,:)) / (s(i+1) - s(i-1));
+ end
+ elseif (dir == 2)
+ dfunds(:,1,:) = (fun(:,2,:) - fun(:,1,:)) / (s(2) - s(1));
+ dfunds(:,end,:) = (fun(:,end,:) - fun(:,end-1,:)) / (s(end) - s(end-1));
+ for i = 2:length(s)-1
+ dfunds(:,i,:) = (fun(:,i+1,:) - fun(:,i-1,:)) / (s(i+1) - s(i-1));
+ end
+ elseif (dir == 3)
+ dfunds(:,:,1) = (fun(:,:,2) - fun(:,:,1)) / (s(2) - s(1));
+ dfunds(:,:,end) = (fun(:,:,end) - fun(:,:,end-1)) / (s(end) - s(end-1));
+ for i = 2:length(s)-1
+ dfunds(:,:,i) = (fun(:,:,i+1) - fun(:,:,i-1)) / (s(i+1) - s(i-1));
+ end
+ else
+ disp("Wrong dir input");
+ end
+end
+
+% Plot Reynolds stress
+function plot_Reynolds_stress(ruu, rvv, rww, ruv, y_norm_vth, timestep)
+
+ load variables/user_inputs.mat;
+ load reference_data/reference.mat;
+
+ % Sqrt & Normalization
+ ruu = sqrt(ruu) / 2; % sqrt(ruu) / Delta U
+ rvv = sqrt(rvv) / 2; % sqrt(rvv) / Delta U
+ rww = sqrt(rww) / 2; % sqrt(rww) / Delta U
+ ruv(ruv > 0) = 0;
+ ruv = sqrt(-ruv) / 2; % sqrt(-ruv) / Delta U
+
+ % Plot
+ f1 = figure("DefaultAxesFontSize", 18);
+ set(f1,"Position",[200 200 1000 700]);
+
+ % sqrt(ruu)/Delta U
+ subplot(2,2,1);
+ plot(y_norm_vth,ruu,'k-','LineWidth',1); hold on; grid on;
+ plot(b1990_rs_ruu(:,1),b1990_rs_ruu(:,2),'g+','LineWidth',1,'MarkerSize',8);
+ plot(v2014_rs_ruu(:,1),v2014_rs_ruu(:,2),'bo','LineWidth',1,'MarkerSize',8);
+ plot(w2022_rs_ruu(:,1),w2022_rs_ruu(:,2),'r^','LineWidth',1,'MarkerSize',8);
+ axis([-1.5 1.5 0 0.25]);
+ xticks([-1.5:0.5:1.5]);
+ yticks([0:0.05:0.25]);
+ xlabel('$y/\delta_\omega$','Interpreter','latex');
+ ylabel('$\sqrt{\left< \overline{\rho} u^{\prime 2} \right>} / \Delta U$','Interpreter','latex');
+ set(gca,'TickLabelInterpreter','latex');
+
+ % sqrt(rvv)/Delta U
+ subplot(2,2,2);
+ plot(y_norm_vth,rvv,'k-','LineWidth',1); hold on; grid on;
+ plot(b1990_rs_rvv(:,1),b1990_rs_rvv(:,2),'g+','LineWidth',1,'MarkerSize',8);
+ plot(v2014_rs_rvv(:,1),v2014_rs_rvv(:,2),'bo','LineWidth',1,'MarkerSize',8);
+ plot(w2022_rs_rvv(:,1),w2022_rs_rvv(:,2),'r^','LineWidth',1,'MarkerSize',8);
+ axis([-1.5 1.5 0 0.25]);
+ xticks([-1.5:0.5:1.5]);
+ yticks([0:0.05:0.25]);
+ xlabel('$y/\delta_\omega$','Interpreter','latex');
+ ylabel('$\sqrt{\left< \overline{\rho} v^{\prime 2} \right>} / \Delta U$','Interpreter','latex');
+ set(gca,'TickLabelInterpreter','latex');
+
+ % sqrt(rww)/Delta U
+ subplot(2,2,3);
+ plot(y_norm_vth,rww,'k-','LineWidth',1); hold on; grid on;
+ plot(b1990_rs_rww(:,1),b1990_rs_rww(:,2),'g+','LineWidth',1,'MarkerSize',8);
+ plot(v2014_rs_rww(:,1),v2014_rs_rww(:,2),'bo','LineWidth',1,'MarkerSize',8);
+ plot(w2022_rs_rww(:,1),w2022_rs_rww(:,2),'r^','LineWidth',1,'MarkerSize',8);
+ axis([-1.5 1.5 0 0.25]);
+ xticks([-1.5:0.5:1.5]);
+ yticks([0:0.05:0.25]);
+ xlabel('$y/\delta_\omega$','Interpreter','latex');
+ ylabel('$\sqrt{\left< \overline{\rho} w^{\prime 2} \right>} / \Delta U$','Interpreter','latex');
+ set(gca,'TickLabelInterpreter','latex');
+
+ % sqrt(-rvu)/Delta U
+ subplot(2,2,4);
+ plot(y_norm_vth,ruv,'k-','LineWidth',1); hold on; grid on;
+ plot(b1990_rs_ruv(:,1),b1990_rs_ruv(:,2),'g+','LineWidth',1,'MarkerSize',8);
+ plot(v2014_rs_ruv(:,1),v2014_rs_ruv(:,2),'bo','LineWidth',1,'MarkerSize',8);
+ plot(w2022_rs_ruv(:,1),w2022_rs_ruv(:,2),'r^','LineWidth',1,'MarkerSize',8);
+ axis([-1.5 1.5 0 0.25]);
+ xticks([-1.5:0.5:1.5]);
+ yticks([0:0.05:0.25]);
+ xlabel('$y/\delta_\omega$','Interpreter','latex');
+ ylabel('$\sqrt{-\left< \overline{\rho} u^{\prime} v^{\prime}\right>} / \Delta U$','Interpreter','latex');
+ set(gca,'TickLabelInterpreter','latex');
+
+ subplot(2,2,2);
+ legend("$\mbox{Present}$",...
+ "$\mbox{Bell \& Mehta (1990)}$",...
+ "$\mbox{Vaghefi (2014)}$",...
+ "$\mbox{Wang et al. (2022)}$",...
+ 'Interpreter','latex','location','northeast');
+
+ saveas(f1, "results/Reynolds_stress/tstep_"+string(timestep),"png");
+ close(f1);
+
+end
+
+% Plot TKE budget
+function plot_tke_budget(T, P, D, y_norm_mth, mth, timestep)
+
+ % Normalization
+ T = T / (8/mth); % T / (Delta U^3 / mth)
+ P = P / (8/mth); % P / (Delta U^3 / mth)
+ D = D / (8/mth); % D / (Delta U^3 / mth)
+
+ load variables/user_inputs.mat;
+ load reference_data/reference.mat;
+
+ % Plot
+ f1 = figure("DefaultAxesFontSize",18);
+ set(f1,"Position",[200 200 1000 700]);
+
+ % Present
+ h1 = plot([-100 -100],[-100 -100],'-k','LineWidth',2); hold on; grid on;
+ plot(y_norm_mth,T,'-b','LineWidth',2);
+ plot(y_norm_mth,P,'-g','LineWidth',2);
+ plot(y_norm_mth,D,'-r','LineWidth',2);
+ xlim([-5 5]); xticks([-5:1:5]);
+ ylim([-0.002 0.003]);
+ xlabel('$y/\delta_\theta$','interpreter','latex');
+ set(gca,'TickLabelInterpreter','latex');
+
+ % Pantano & Sarkar (2002)
+ h2 = plot([-100 -100],[-100 -100],'ko','LineWidth',2,'MarkerSize',8);
+ plot(p2002_tke_transport(:,1),p2002_tke_transport(:,2),'bo','LineWidth',2,'MarkerSize',8);
+ plot(p2002_tke_production(:,1),p2002_tke_production(:,2),'go','LineWidth',2,'MarkerSize',8);
+ plot(p2002_tke_dissipation(:,1),p2002_tke_dissipation(:,2),'ro','LineWidth',2,'MarkerSize',8);
+
+ % Rogers & Moser (1994)
+ h3 = plot([-100 -100],[-100 -100],'k^','LineWidth',2,'MarkerSize',8);
+ plot(r1994_tke_transport(:,1),r1994_tke_transport(:,2),'b^','LineWidth',2,'MarkerSize',8);
+ plot(r1994_tke_production(:,1),r1994_tke_production(:,2),'g^','LineWidth',2,'MarkerSize',8);
+ plot(r1994_tke_dissipation(:,1),r1994_tke_dissipation(:,2),'r^','LineWidth',2,'MarkerSize',8);
+
+ % Vaghefi (2014)
+ h4 = plot([-100 -100],[-100 -100],'k+','LineWidth',2,'MarkerSize',8);
+ plot(v2014_tke_transport(:,1),v2014_tke_transport(:,2),'b+','LineWidth',2,'MarkerSize',8);
+ plot(v2014_tke_production(:,1),v2014_tke_production(:,2),'g+','LineWidth',2,'MarkerSize',8);
+ plot(v2014_tke_dissipation(:,1),v2014_tke_dissipation(:,2),'r+','LineWidth',2,'MarkerSize',8);
+
+ % Wang et al. (2022)
+ h5 = plot([-100 -100],[-100 -100],'k*','LineWidth',2,'MarkerSize',8);
+ plot(w2022_tke_transport(:,1),w2022_tke_transport(:,2),'b*','LineWidth',2,'MarkerSize',8);
+ plot(w2022_tke_production(:,1),w2022_tke_production(:,2),'g*','LineWidth',2,'MarkerSize',8);
+ plot(w2022_tke_dissipation(:,1),w2022_tke_dissipation(:,2),'r*','LineWidth',2,'MarkerSize',8);
+
+ legend([h1,h2,h3,h4,h5], {"$\mbox{Present}$", ...
+ "$\mbox{Pantano \& Sarkar (2002)}$", ...
+ "$\mbox{Rogers \& Moser (1994)}$", ...
+ "$\mbox{Vaghefi (2014)}$", ...
+ "$\mbox{Wang et al. (2022)}$"}, ...
+ 'interpreter','latex','location','northeast');
+
+ saveas(f1, "results/tke_budget/tstep_"+string(timestep),"png");
+ close(f1);
+end
+
+% Compute and plot growth rate
+function plot_growth_rate(time, mth)
+
+ % Growth rate
+ dmth = (mth(2:end) - mth(1:end-1)) ./ (time(2:end) - time(1:end-1));
+
+ % Normalization
+ dmth = dmth / 2; % (dmth/dt) * (1/Delta U)
+
+ % Integrated growth rate
+ f_growth_rate = figure("DefaultAxesFontSize", 18);
+
+ h1 = plot(time(1:end-1),dmth,'-ko','LineWidth',2); hold on; grid on;
+ h2 = plot([0 max(time)],[0.012 0.012], 'r--','LineWidth',1);
+ h3 = plot([0 max(time)],[0.0135 0.0135], 'b--','LineWidth',1);
+ h4 = plot([0 max(time)],[0.014 0.014], 'g--','LineWidth',1);
+ h5 = plot([0 max(time)],[0.015 0.015], 'c--','LineWidth',1);
+ h6 = plot([0 max(time)],[0.017 0.017], 'm--.','LineWidth',1);
+ axis([0 max(time) 0 0.04]);
+ xlabel('$t U_1 / \delta_{\omega h}^0$','interpreter','latex');
+ ylabel('$\dot{\delta}_{\theta} / \Delta U$','interpreter','latex');
+ legend([h1,h2,h3,h4,h5,h6],{"$\mbox{Present}$",...
+ "$0.012^{[1]}$",... % [1] Baltzer & Livescu (2020, JFM)
+ "$0.0135^{[2]}$",... % [2] Vaghefi (2014, Thesis)
+ "$0.014^{[3,4]}$",... % [3] Rogers & Moser (1993, PoF) [4] Blakeley et al. (2023, JFM)
+ "$0.015^{[5]}$",... % [5] Sharan, Matheou & Dimotakis (2019, JFM)
+ "$0.017^{[6]}$"},... % [6] Almagro et al. (2017, JFM)
+ 'Interpreter','latex','location','northeast');
+ set(gca,'TickLabelInterpreter','latex');
+
+ saveas(f_growth_rate, "results/growth_rate.png");
+ close(f_growth_rate);
+
+end
diff --git a/examples/3D_turb_mixing/turbulence_stat/set_user_inputs.m b/examples/3D_turb_mixing/turbulence_stat/set_user_inputs.m
new file mode 100644
index 0000000000..979fbddc08
--- /dev/null
+++ b/examples/3D_turb_mixing/turbulence_stat/set_user_inputs.m
@@ -0,0 +1,57 @@
+function set_user_inputs()
+
+ %% Options
+ Reynolds_stress = false;
+ tke_budget = false;
+
+ %% Fluid property (Air)
+ gamma_a = 1.4;
+ pi_inf_a = 0;
+ Gamma = 1/(gamma_a - 1);
+ Pi_inf = gamma_a*pi_inf_a/(gamma_a - 1);
+
+ %% Reynolds number
+ Re = 320; % Reynolds number
+
+ %% Grids
+ Lx = 160.0; % Domain size in x-direction
+ Ly = 160.0; % Domain size in y-direction
+ Lz = 80.0; % Domain size in z-direction
+
+ m = 1023; % Number of grids in x-direction (0, ..., m)
+ n = 1023; % Number of grids in y-direction (0, ..., n)
+ p = 511; % Number of grids in z-direction (0, ..., p)
+
+ % Stretched grid in y-direction
+ stretch_y = "T";
+ a_y = 2;
+ y_a = -0.3*Ly;
+ y_b = 0.3*Ly;
+ loops_y = 2;
+
+ %% Time
+ Nt_beg = 384*0;
+ Nt_end = 384*30;
+ Nt_save = 384;
+ dt = 0.013033561420658195;
+ Nfiles = (Nt_end - Nt_beg)/Nt_save + 1;
+ timesteps = Nt_beg:Nt_save:Nt_end;
+ time = timesteps*dt;
+
+ %% Configuration
+ num_fluids = 1; % Number of fluids
+ num_dims = 3; % Number of dimensions
+
+ %% Index
+ contxb = 1;
+ contxe = num_fluids;
+ momxb = contxe + 1;
+ momxe = contxe + num_dims;
+ E_idx = momxe + 1;
+ advxb = E_idx + 1;
+ advxe = E_idx + num_fluids;
+ sys_size = advxe;
+
+ save variables/user_inputs.mat;
+
+end
diff --git a/src/common/m_constants.fpp b/src/common/m_constants.fpp
index 1a77f83f28..32a68ec882 100644
--- a/src/common/m_constants.fpp
+++ b/src/common/m_constants.fpp
@@ -64,6 +64,13 @@ module m_constants
! Relativity
integer, parameter :: relativity_cons_to_prim_max_iter = 100
+ ! Pseudo-random number generator
+ integer, parameter :: modulus = 2**30 - 1
+ integer, parameter :: multiplier = 1664525
+ integer, parameter :: increment = 1013904223
+ integer, parameter :: amplifier = 3**13
+ real(wp), parameter :: decimal_trim = 1.e5_wp
+
! System constants
integer, parameter :: CASE_FILE_ERROR_CODE = 22
diff --git a/src/pre_process/m_checker.fpp b/src/pre_process/m_checker.fpp
index 972f3bcb70..7423d8ed64 100644
--- a/src/pre_process/m_checker.fpp
+++ b/src/pre_process/m_checker.fpp
@@ -188,11 +188,7 @@ contains
"mixlayer_vel_profile requires n > 0")
! Instability wave
- @:PROHIBIT(mixlayer_perturb .and. n == 0, "mixlayer_perturb requires n > 0")
- @:PROHIBIT(mixlayer_perturb .and. model_eqns /= 2, "mixlayer_perturb requires model_eqns = 2")
- @:PROHIBIT(mixlayer_perturb .and. num_fluids > 1, "mixlayer_perturb requires num_fluids = 1")
- @:PROHIBIT(mixlayer_perturb .and. any((/bc_y%beg, bc_y%end/) /= BC_CHAR_NR_SUB_BUFFER), &
- "mixlayer_perturb requires both bc_y%beg and bc_y%end to be 6")
+ @:PROHIBIT(mixlayer_perturb .and. p == 0, "mixlayer_perturb requires p > 0")
@:PROHIBIT(elliptic_smoothing .and. elliptic_smoothing_iters < 1, &
"elliptic_smoothing_iters must be positive")
diff --git a/src/pre_process/m_global_parameters.fpp b/src/pre_process/m_global_parameters.fpp
index bb49f35d0d..2358cd71f4 100644
--- a/src/pre_process/m_global_parameters.fpp
+++ b/src/pre_process/m_global_parameters.fpp
@@ -141,8 +141,11 @@ module m_global_parameters
logical :: mixlayer_vel_profile !< Set hyperbolic tangent streamwise velocity profile
real(wp) :: mixlayer_vel_coef !< Coefficient for the hyperbolic tangent streamwise velocity profile
- real(wp) :: mixlayer_domain !< Domain for the hyperbolic tangent streamwise velocity profile
logical :: mixlayer_perturb !< Superimpose instability waves to surrounding fluid flow
+ integer :: mixlayer_perturb_nk !< Number of Fourier modes for perturbation with mixlayer_perturb flag
+ real(wp) :: mixlayer_perturb_k0 !< Peak wavenumber of prescribed energy spectra with mixlayer_perturb flag
+ !! Default value (k0 = 0.4446) is most unstable mode obtained from linear stability analysis
+ !! See Michalke (1964, JFM) for details
real(wp) :: pi_fac !< Factor for artificial pi_inf
@@ -367,8 +370,9 @@ contains
bubbles_lagrange = .false.
mixlayer_vel_profile = .false.
mixlayer_vel_coef = 1._wp
- mixlayer_domain = 1._wp
mixlayer_perturb = .false.
+ mixlayer_perturb_nk = 100
+ mixlayer_perturb_k0 = 0.4446_wp
perturb_flow = .false.
perturb_flow_fluid = dflt_int
perturb_flow_mag = dflt_real
diff --git a/src/pre_process/m_initial_condition.fpp b/src/pre_process/m_initial_condition.fpp
index 8b7e001ac7..4761d110a6 100644
--- a/src/pre_process/m_initial_condition.fpp
+++ b/src/pre_process/m_initial_condition.fpp
@@ -189,7 +189,7 @@ contains
if (perturb_flow) call s_perturb_surrounding_flow(q_prim_vf)
if (perturb_sph) call s_perturb_sphere(q_prim_vf)
- if (mixlayer_perturb) call s_superposition_instability_wave(q_prim_vf)
+ if (mixlayer_perturb) call s_perturb_mixlayer(q_prim_vf)
if (elliptic_smoothing) call s_elliptic_smoothing(q_prim_vf, bc_type)
! Converting the primitive variables to the conservative ones
diff --git a/src/pre_process/m_mpi_proxy.fpp b/src/pre_process/m_mpi_proxy.fpp
index b9d60e2ac1..b0aabd3658 100644
--- a/src/pre_process/m_mpi_proxy.fpp
+++ b/src/pre_process/m_mpi_proxy.fpp
@@ -47,7 +47,7 @@ contains
& 'weno_order', 'precision', 'perturb_flow_fluid', &
& 'perturb_sph_fluid', 'num_patches', 'thermal', 'nb', 'dist_type',&
& 'R0_type', 'relax_model', 'num_ibs', 'n_start', 'elliptic_smoothing_iters', &
- & 'num_bc_patches' ]
+ & 'num_bc_patches', 'mixlayer_perturb_nk' ]
call MPI_BCAST(${VAR}$, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
#:endfor
@@ -58,7 +58,7 @@ contains
& 'qbmm', 'file_per_process', 'adv_n', 'ib' , 'cfl_adap_dt', &
& 'cfl_const_dt', 'cfl_dt', 'surface_tension', &
& 'hyperelasticity', 'pre_stress', 'elliptic_smoothing', 'viscous',&
- & 'bubbles_lagrange', 'bc_io', 'mhd', 'relativity', 'cont_damage' ]
+ & 'bubbles_lagrange', 'bc_io', 'mhd', 'relativity', 'cont_damage']
call MPI_BCAST(${VAR}$, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
#:endfor
call MPI_BCAST(fluid_rho(1), num_fluids_max, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
@@ -69,8 +69,8 @@ contains
& 'bc_x%end', 'bc_y%beg', 'bc_y%end', 'bc_z%beg', 'bc_z%end', &
& 'perturb_flow_mag', 'pref', 'rhoref', 'poly_sigma', 'R0ref', &
& 'Web', 'Ca', 'Re_inv', 'sigR', 'sigV', 'rhoRV', 'palpha_eps', &
- & 'ptgalpha_eps', 'sigma', 'pi_fac', 'mixlayer_vel_coef', &
- & 'mixlayer_domain', 'Bx0' ]
+ & 'ptgalpha_eps', 'sigma', 'pi_fac', 'mixlayer_vel_coef', 'Bx0', &
+ & 'mixlayer_perturb_k0']
call MPI_BCAST(${VAR}$, 1, mpi_p, 0, MPI_COMM_WORLD, ierr)
#:endfor
diff --git a/src/pre_process/m_perturbation.fpp b/src/pre_process/m_perturbation.fpp
index 0bfef5da83..32b3eaa2b7 100644
--- a/src/pre_process/m_perturbation.fpp
+++ b/src/pre_process/m_perturbation.fpp
@@ -12,44 +12,20 @@ module m_perturbation
use m_mpi_proxy !< Message passing interface (MPI) module proxy
- use m_eigen_solver ! Subroutines to solve eigenvalue problem for
-
use m_boundary_common ! Boundary conditions module
- use m_helper_basic !< Functions to compare floating point numbers
+ use m_helper
use ieee_arithmetic
implicit none
- integer :: mixlayer_nvar ! Number of variables in linear stability analysis solver for mixing layer
- integer, allocatable, dimension(:) :: mixlayer_var ! Index of variables in linear stability analysis solver
- integer :: nbp ! Number of grid cell boundary points in y-direction
- integer :: mixlayer_bc_fd ! Order of finite difference applied at the boundaries of mixing layer
- integer :: n_bc_skip ! Number of points skipped in the linear stability analysis due to the boundary condition
-
real(wp), allocatable, dimension(:, :, :, :) :: q_prim_temp
contains
impure subroutine s_initialize_perturbation_module()
- if (mixlayer_perturb) then
- mixlayer_bc_fd = 2
- nbp = n + 2
- if (model_eqns == 2 .and. num_fluids == 1) then
- n_bc_skip = mixlayer_bc_fd*2
- mixlayer_nvar = 5 ! 1 continuity + 3 momentum + 1 energy
- allocate (mixlayer_var(mixlayer_nvar))
-
- mixlayer_var(1) = contxb
- mixlayer_var(2) = momxb
- mixlayer_var(3) = momxb + 1
- mixlayer_var(4) = momxb + 2
- mixlayer_var(5) = momxb + 3
- end if
- end if
-
if (elliptic_smoothing) then
allocate (q_prim_temp(0:m, 0:n, 0:p, 1:sys_size))
end if
@@ -111,510 +87,8 @@ contains
end do
end do
end do
-
end subroutine s_perturb_surrounding_flow
- !> This subroutine computes velocity perturbations for a temporal mixing
- !! layer with hypertangent mean streamwise velocity profile
- !! obtained from linear stability analysis. For a 2D case,
- !! instability waves with spatial wavenumbers, (4,0), (2,0),
- !! and (1,0) are superposed. For a 3D waves, (4,4), (4,-4),
- !! (2,2), (2,-2), (1,1), (1,-1) areadded on top of 2D waves.
- impure subroutine s_superposition_instability_wave(q_prim_vf)
- type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf
- real(wp), dimension(mixlayer_nvar, 0:m, 0:n, 0:p) :: wave, wave1, wave2, wave_tmp
- real(wp) :: uratio, Ldomain
- integer :: i, j, k, q
-
- uratio = 1._wp/patch_icpp(1)%vel(1)
- Ldomain = mixlayer_domain*patch_icpp(1)%length_y
-
- wave = 0._wp
- wave1 = 0._wp
- wave2 = 0._wp
-
- ! Compute 2D waves
- call s_instability_wave(2*pi*4.0_wp/Ldomain, 0._wp, wave_tmp, 0._wp)
- wave1 = wave1 + wave_tmp
- call s_instability_wave(2*pi*2.0_wp/Ldomain, 0._wp, wave_tmp, 0._wp)
- wave1 = wave1 + wave_tmp
- call s_instability_wave(2*pi*1.0_wp/Ldomain, 0._wp, wave_tmp, 0._wp)
- wave1 = wave1 + wave_tmp
- wave = wave1*0.05_wp
-
- if (p > 0) then
- ! Compute 3D waves with phase shifts.
- call s_instability_wave(2*pi*4.0_wp/Ldomain, 2*pi*4.0_wp/Ldomain, wave_tmp, 2*pi*11._wp/31._wp)
- wave2 = wave2 + wave_tmp
- call s_instability_wave(2*pi*2.0_wp/Ldomain, 2*pi*2.0_wp/Ldomain, wave_tmp, 2*pi*13._wp/31._wp)
- wave2 = wave2 + wave_tmp
- call s_instability_wave(2*pi*1.0_wp/Ldomain, 2*pi*1.0_wp/Ldomain, wave_tmp, 2*pi*17._wp/31._wp)
- wave2 = wave2 + wave_tmp
- call s_instability_wave(2*pi*4.0_wp/Ldomain, -2*pi*4.0_wp/Ldomain, wave_tmp, 2*pi*19._wp/31._wp)
- wave2 = wave2 + wave_tmp
- call s_instability_wave(2*pi*2.0_wp/Ldomain, -2*pi*2.0_wp/Ldomain, wave_tmp, 2*pi*23._wp/31._wp)
- wave2 = wave2 + wave_tmp
- call s_instability_wave(2*pi*1.0_wp/Ldomain, -2*pi*1.0_wp/Ldomain, wave_tmp, 2*pi*29._wp/31._wp)
- wave2 = wave2 + wave_tmp
- wave = wave + 0.15_wp*wave2
- end if
-
- ! Superpose velocity perturbuations (instability waves) to the velocity field
- do k = 0, p
- do j = 0, n
- do i = 0, m
- q_prim_vf(contxb)%sf(i, j, k) = q_prim_vf(contxb)%sf(i, j, k) + wave(mixlayer_var(1), i, j, k) ! rho
- q_prim_vf(momxb)%sf(i, j, k) = q_prim_vf(momxb)%sf(i, j, k) + wave(mixlayer_var(2), i, j, k)/uratio ! u
- q_prim_vf(momxb + 1)%sf(i, j, k) = q_prim_vf(momxb + 1)%sf(i, j, k) + wave(mixlayer_var(3), i, j, k)/uratio ! v
- if (p > 0) then
- q_prim_vf(momxb + 2)%sf(i, j, k) = q_prim_vf(momxb + 2)%sf(i, j, k) + wave(mixlayer_var(4), i, j, k)/uratio ! w
- end if
- q_prim_vf(E_idx)%sf(i, j, k) = q_prim_vf(E_idx)%sf(i, j, k) + wave(mixlayer_var(5), i, j, k)/uratio**2 ! p
-
- if (bubbles_euler .and. (.not. qbmm)) then
- do q = 1, nb
- call s_compute_equilibrium_state(q_prim_vf(E_idx)%sf(i, j, k), R0(q), q_prim_vf(bub_idx%rs(q))%sf(i, j, k))
- end do
- end if
- end do
- end do
- end do
-
- end subroutine s_superposition_instability_wave
-
- !> This subroutine computes equilibrium bubble radius of the perturbed pressure field
- impure subroutine s_compute_equilibrium_state(fP, fR0, fR)
- real(wp), intent(in) :: fP, fR0
- real(wp), intent(inout) :: fR
- real(wp) :: f0, f1
- real(wp) :: gam_b
- integer :: ii
-
- gam_b = 1._wp + 1._wp/fluid_pp(num_fluids + 1)%gamma
-
- ! Loop
- ii = 1
- do while (.true.)
-
- f0 = (Ca + 2._wp/Web)*(fR0/fR)**(3._wp*gam_b) - 2._wp/(Web*fR) + 1._wp - Ca - fP
- f1 = -3._wp*gam_b*(Ca + 2._wp/Web)*(fR0/fR)**(3._wp*gam_b + 1._wp) + 2._wp/(Web*fR**2._wp)
-
- if (abs(f0) <= 1.e-10_wp) then
- ! Converged
- exit
- else
- ! Update radius
- fR = fR - f0/f1
- end if
-
- ! Failed case
- if (ieee_is_nan(f0) .or. &
- ieee_is_nan(f1) .or. &
- ii > 1000 .or. &
- fR < 0._wp) then
-
- print *, "Failed to compute equilibrium radius"
-
- fR = fR0
- exit
- end if
-
- ii = ii + 1
- end do
-
- end subroutine s_compute_equilibrium_state
-
- !> This subroutine computes instability waves for a given set of spatial
- !! wavenumbers (alpha, beta) in x and z directions.
- !! The eigenvalue problem is derived from the linearized
- !! Euler equations with parallel mean flow assumption
- !! (See Sandham 1989 PhD thesis for details).
- pure subroutine s_instability_wave(alpha, beta, wave, shift)
- real(wp), intent(in) :: alpha, beta !< spatial wavenumbers
- real(wp), dimension(mixlayer_nvar, 0:m, 0:n, 0:p), intent(inout) :: wave !< instability wave
- real(wp), intent(in) :: shift !< phase shift
- real(wp), dimension(0:nbp - 1) :: u_mean !< mean density and velocity profiles
- real(wp) :: rho_mean, p_mean !< mean density and pressure
- real(wp), dimension(0:nbp - 1, 0:nbp - 1) :: d !< differential operator in y dir
- real(wp) :: gam, pi_inf, mach, c1, adv
- real(wp) :: xratio, uratio
- integer :: j !< generic loop iterators
-
- xratio = mixlayer_vel_coef
- uratio = 1._wp/patch_icpp(1)%vel(1)
-
- ! Set fluid flow properties
- if (bubbles_euler) then
- adv = patch_icpp(1)%alpha(num_fluids)
- else
- adv = 0._wp
- end if
- gam = 1._wp + 1._wp/fluid_pp(1)%gamma
- pi_inf = fluid_pp(1)%pi_inf*(gam - 1._wp)/gam*uratio**2
- rho_mean = patch_icpp(1)%alpha_rho(1)
- p_mean = patch_icpp(1)%pres*uratio**2
- c1 = sqrt((gam*(p_mean + pi_inf))/(rho_mean*(1._wp - adv)))
- mach = 1._wp/c1
-
- ! Assign mean profiles
- do j = 0, n + 1
- u_mean(j) = tanh(y_cb(j - 1)*xratio)
- end do
-
- ! Compute differential operator in y-dir
- ! based on 2nd order central difference
- d = 0._wp
- d(0, 0) = -1._wp/((y_cb(0) - y_cb(-1))*xratio)
- d(0, 1) = 1._wp/((y_cb(0) - y_cb(-1))*xratio)
- do j = 1, n
- d(j, j - 1) = -1._wp/((y_cb(j) - y_cb(j - 2))*xratio)
- d(j, j + 1) = 1._wp/((y_cb(j) - y_cb(j - 2))*xratio)
- end do
- d(n + 1, n) = -1._wp/((y_cb(n) - y_cb(n - 1))*xratio)
- d(n + 1, n + 1) = 1._wp/((y_cb(n) - y_cb(n - 1))*xratio)
-
- ! Compute
- call s_solve_linear_system(alpha, beta, u_mean, rho_mean, p_mean, d, gam, pi_inf, mach, wave, shift)
-
- end subroutine s_instability_wave
-
- !> This subroutine solves linear system from linear stability analysis and
- !! generate instability waves for the given set of spatial
- !! wave numbers and phase shift.
- pure subroutine s_solve_linear_system(alpha, beta, u_mean, rho_mean, p_mean, d, gam, pi_inf, mach, wave, shift)
- real(wp), intent(in) :: alpha, beta !< spatial wavenumbers
- real(wp), dimension(0:nbp - 1), intent(in) :: u_mean !< mean velocity profiles
- real(wp), intent(in) :: rho_mean, p_mean !< mean density and pressure
- real(wp), dimension(0:nbp - 1, 0:nbp - 1), intent(in) :: d !< differential operator in y dir
- real(wp), intent(in) :: gam, pi_inf, mach, shift
- real(wp), dimension(mixlayer_nvar, 0:m, 0:n, 0:p), intent(inout) :: wave
-
- real(wp), dimension(0:nbp - 1) :: drho_mean, du_mean !< y-derivatives of mean profiles
- real(wp), dimension(0:mixlayer_nvar*nbp - 1, 0:mixlayer_nvar*nbp - 1) :: ar, ai !< matrices for eigenvalue problem
- real(wp), dimension(0:mixlayer_nvar*nbp - 1, 0:mixlayer_nvar*nbp - 1) :: br, bi, ci !< matrices for eigenvalue problem
- real(wp), dimension(0:mixlayer_nvar*n - n_bc_skip - 1, 0:mixlayer_nvar*n - n_bc_skip - 1) :: hr, hi !< matrices for eigenvalue problem
-
- real(wp), dimension(0:mixlayer_nvar*n - n_bc_skip - 1, 0:mixlayer_nvar*n - n_bc_skip - 1) :: zr, zi !< eigenvectors
- real(wp), dimension(0:mixlayer_nvar*n - n_bc_skip - 1) :: wr, wi !< eigenvalues
- real(wp), dimension(0:mixlayer_nvar*n - n_bc_skip - 1) :: fv1, fv2, fv3 !< temporary memory
-
- integer :: ierr
- integer :: j, k !< generic loop iterators
- integer :: ii, jj !< block matrix indices
-
- ! Compute y-derivatives of rho and u
- do j = 0, nbp - 1
- drho_mean(j) = 0
- du_mean(j) = 0
- do k = 0, nbp - 1
- drho_mean(j) = 0._wp
- du_mean(j) = du_mean(j) + d(j, k)*u_mean(k)
- end do
- end do
-
- ! Compute B and C, then A = B + C. Here, A is the matrix for the linear
- ! systems of equation (i.e. we are going to solve x for Ax = lambda x).
- ! Here, B includes components of A without differential operator, and
- ! C includes components of A with differential operator.
- br = 0._wp
- bi = 0._wp
- ci = 0._wp
- do j = 0, nbp - 1
- ii = mixlayer_var(1); jj = mixlayer_var(1); br((ii - 1)*nbp + j, (jj - 1)*nbp + j) = alpha*u_mean(j);
- ii = mixlayer_var(1); jj = mixlayer_var(2); br((ii - 1)*nbp + j, (jj - 1)*nbp + j) = alpha*rho_mean;
- ii = mixlayer_var(1); jj = mixlayer_var(3); bi((ii - 1)*nbp + j, (jj - 1)*nbp + j) = -drho_mean(j);
- ii = mixlayer_var(1); jj = mixlayer_var(4); br((ii - 1)*nbp + j, (jj - 1)*nbp + j) = beta*rho_mean;
- ii = mixlayer_var(2); jj = mixlayer_var(2); br((ii - 1)*nbp + j, (jj - 1)*nbp + j) = alpha*u_mean(j);
- ii = mixlayer_var(2); jj = mixlayer_var(3); bi((ii - 1)*nbp + j, (jj - 1)*nbp + j) = -du_mean(j);
- ii = mixlayer_var(2); jj = mixlayer_var(5); br((ii - 1)*nbp + j, (jj - 1)*nbp + j) = alpha/rho_mean;
- ii = mixlayer_var(3); jj = mixlayer_var(3); br((ii - 1)*nbp + j, (jj - 1)*nbp + j) = alpha*u_mean(j);
- ii = mixlayer_var(4); jj = mixlayer_var(4); br((ii - 1)*nbp + j, (jj - 1)*nbp + j) = alpha*u_mean(j);
- ii = mixlayer_var(4); jj = mixlayer_var(5); br((ii - 1)*nbp + j, (jj - 1)*nbp + j) = beta/rho_mean;
- ii = mixlayer_var(5); jj = mixlayer_var(2); br((ii - 1)*nbp + j, (jj - 1)*nbp + j) = gam*(p_mean + pi_inf)*alpha;
- ii = mixlayer_var(5); jj = mixlayer_var(4); br((ii - 1)*nbp + j, (jj - 1)*nbp + j) = gam*(p_mean + pi_inf)*beta;
- ii = mixlayer_var(5); jj = mixlayer_var(5); br((ii - 1)*nbp + j, (jj - 1)*nbp + j) = alpha*u_mean(j);
- do k = 0, n + 1
- ii = mixlayer_var(1); jj = mixlayer_var(3); ci((ii - 1)*nbp + j, (jj - 1)*nbp + k) = -rho_mean*d(j, k);
- ii = mixlayer_var(3); jj = mixlayer_var(5); ci((ii - 1)*nbp + j, (jj - 1)*nbp + k) = -d(j, k)/rho_mean;
- ii = mixlayer_var(5); jj = mixlayer_var(3); ci((ii - 1)*nbp + j, (jj - 1)*nbp + k) = -gam*(p_mean + pi_inf)*d(j, k);
- end do
- end do
- ar = br
- ai = bi + ci
-
- ! Apply BC to ar and ai matrices
- if (bc_y%beg == BC_CHAR_NR_SUB_BUFFER .and. bc_y%end == BC_CHAR_NR_SUB_BUFFER) then
- ! Nonreflecting subsonic buffer BC
- call s_instability_nonreflecting_subsonic_buffer_bc(ar, ai, hr, hi, rho_mean, mach)
- end if
-
- ! Compute eigenvalues and eigenvectors
- call cg(mixlayer_nvar*n - n_bc_skip, mixlayer_nvar*n - n_bc_skip, hr, hi, wr, wi, zr, zi, fv1, fv2, fv3, ierr)
-
- ! Generate instability wave
- call s_generate_wave(wi, zr, zi, rho_mean, mach, alpha, beta, wave, shift)
-
- end subroutine s_solve_linear_system
-
- !> This subroutine applies non-reflecting subsonic buffer boundary condition
- !! to the linear system of equations (i.e. matrix A).
- pure subroutine s_instability_nonreflecting_subsonic_buffer_bc(ar, ai, hr, hi, rho_mean, mach)
- real(wp), dimension(0:mixlayer_nvar*nbp - 1, 0:mixlayer_nvar*nbp - 1), intent(inout) :: ar, ai !< matrices for eigenvalue problem
- real(wp), dimension(0:mixlayer_nvar*n - n_bc_skip - 1, 0:mixlayer_nvar*n - n_bc_skip - 1), intent(out) :: hr, hi !< matrices for eigenvalue problem
- real(wp), intent(in) :: rho_mean !< mean density profiles
- real(wp), intent(in) :: mach
- real(wp), dimension(0:mixlayer_nvar*n - 1, 0:mixlayer_nvar*n - 1) :: fr, fi !< matrices for eigenvalue problem
- real(wp), dimension(0:mixlayer_nvar*n - n_bc_skip - 1, 0:mixlayer_nvar*n - 1) :: gr, gi !< matrices for eigenvalue problem
- integer :: i, j, k, l, ii, jj
-
- ! Condition 1: v = 0 at BC - no action required here
-
- ! Condition 2: du/dy = 0 at BC
- do j = 0, mixlayer_nvar*nbp - 1
- ! at y_domain%beg
- ii = mixlayer_var(1)*nbp
- ar(j, ii + 1) = ar(j, ii + 1) + ar(j, ii)
- ai(j, ii + 1) = ai(j, ii + 1) + ai(j, ii)
- ! at y_domain%end
- ii = mixlayer_var(1)*nbp + nbp - 1
- ar(j, ii - 1) = ar(j, ii - 1) + ar(j, ii)
- ai(j, ii - 1) = ai(j, ii - 1) + ai(j, ii)
- end do
-
- ! Condition 3: dw/dy = 0 at BC
- do j = 0, mixlayer_nvar*nbp - 1
- ! at y_domain%beg
- ii = (mixlayer_var(3))*nbp
- ar(j, ii + 1) = ar(j, ii + 1) + ar(j, ii)
- ai(j, ii + 1) = ai(j, ii + 1) + ai(j, ii)
- ! at y_domain%end
- ii = (mixlayer_var(3))*nbp + nbp - 1
- ar(j, ii - 1) = ar(j, ii - 1) + ar(j, ii)
- ai(j, ii - 1) = ai(j, ii - 1) + ai(j, ii)
- end do
-
- ! Condition 4: dp/dy +- rho c dv/dy = 0 at BC
- do j = 0, mixlayer_nvar*nbp - 1
- ! at y_domain%beg
- ii = mixlayer_var(4)*nbp
- ar(j, ii + 1) = ar(j, ii + 1) + ar(j, ii)
- ai(j, ii + 1) = ai(j, ii + 1) + ai(j, ii)
- jj = mixlayer_var(2)*nbp
- ar(j, jj + 1) = ar(j, jj + 1) + ar(j, ii)*rho_mean/mach
- ai(j, jj + 1) = ai(j, jj + 1) + ai(j, ii)*rho_mean/mach
- ! at y_domain%end
- ii = mixlayer_var(4)*nbp + nbp - 1
- ar(j, ii - 1) = ar(j, ii - 1) + ar(j, ii)
- ai(j, ii - 1) = ai(j, ii - 1) + ai(j, ii)
- jj = mixlayer_var(2)*nbp + nbp - 1
- ar(j, jj - 1) = ar(j, jj - 1) - ar(j, ii)*rho_mean/mach
- ai(j, jj - 1) = ai(j, jj - 1) - ai(j, ii)*rho_mean/mach
- end do
-
- ! Condition 5: c^2 drho/dy +- dp/dy = 0 at BC
- do j = 0, mixlayer_nvar*nbp - 1
- ! at y_domain%beg
- ii = 0
- ar(j, ii + 1) = ar(j, ii + 1) + ar(j, ii)
- ai(j, ii + 1) = ai(j, ii + 1) + ai(j, ii)
- jj = mixlayer_var(2)*nbp
- ar(j, jj + 1) = ar(j, jj + 1) + ar(j, ii)*rho_mean*mach
- ai(j, jj + 1) = ai(j, jj + 1) + ai(j, ii)*rho_mean*mach
- ! at y_domain%end
- ii = nbp - 1
- ar(j, ii - 1) = ar(j, ii - 1) + ar(j, ii)
- ai(j, ii - 1) = ai(j, ii - 1) + ai(j, ii)
- jj = mixlayer_var(2)*nbp + nbp - 1
- ar(j, jj - 1) = ar(j, jj - 1) - ar(j, ii)*rho_mean*mach
- ai(j, jj - 1) = ai(j, jj - 1) - ai(j, ii)*rho_mean*mach
- end do
-
- ! Remove unnecessary rows of the matrix A (rho, u, v, w, p at the boundaries)
- fr = 0._wp
- fi = 0._wp
- do ii = 1, mixlayer_nvar
- do jj = 1, mixlayer_nvar
- do k = 0, n - 1
- do l = 0, n - 1
- fr((ii - 1)*n + k, (jj - 1)*n + l) = ar((ii - 1)*nbp + k + 1, (jj - 1)*nbp + l + 1)
- fi((ii - 1)*n + k, (jj - 1)*n + l) = ai((ii - 1)*nbp + k + 1, (jj - 1)*nbp + l + 1)
- end do
- end do
- end do
- end do
-
- gr = 0._wp
- gi = 0._wp
- do ii = 1, mixlayer_nvar
- do j = 0, mixlayer_nvar*n - 1
- if (ii <= mixlayer_var(2)) then
- do k = 0, n - 1
- gr((ii - 1)*n + k, j) = fr((ii - 1)*n + k, j)
- gi((ii - 1)*n + k, j) = fi((ii - 1)*n + k, j)
- end do
- elseif (ii == mixlayer_var(3)) then
- do k = 0, n - n_bc_skip - 1
- gr((ii - 1)*n + k, j) = fr((ii - 1)*n + k + mixlayer_bc_fd, j)
- gi((ii - 1)*n + k, j) = fi((ii - 1)*n + k + mixlayer_bc_fd, j)
- end do
- else
- do k = 0, n - 1
- gr((ii - 1)*n - n_bc_skip + k, j) = fr((ii - 1)*n + k, j)
- gi((ii - 1)*n - n_bc_skip + k, j) = fi((ii - 1)*n + k, j)
- end do
- end if
- end do
- end do
-
- hr = 0._wp
- hi = 0._wp
- do i = 0, mixlayer_nvar*n - n_bc_skip - 1
- do jj = 1, mixlayer_nvar
- if (jj <= mixlayer_var(2)) then
- do k = 0, n - 1
- hr(i, (jj - 1)*n + k) = gr(i, (jj - 1)*n + k)
- hi(i, (jj - 1)*n + k) = gi(i, (jj - 1)*n + k)
- end do
- elseif (jj == mixlayer_var(3)) then
- do k = 0, n - n_bc_skip - 1
- hr(i, (jj - 1)*n + k) = gr(i, (jj - 1)*n + k + mixlayer_bc_fd)
- hi(i, (jj - 1)*n + k) = gi(i, (jj - 1)*n + k + mixlayer_bc_fd)
- end do
- else
- do k = 0, n - 1
- hr(i, (jj - 1)*n - n_bc_skip + k) = gr(i, (jj - 1)*n + k)
- hi(i, (jj - 1)*n - n_bc_skip + k) = gi(i, (jj - 1)*n + k)
- end do
- end if
- end do
- end do
-
- end subroutine s_instability_nonreflecting_subsonic_buffer_bc
-
- !> This subroutine generates an instability wave using the most unstable
- !! eigenvalue and corresponding eigenvector among the
- !! given set of eigenvalues and eigenvectors.
- pure subroutine s_generate_wave(wi, zr, zi, rho_mean, mach, alpha, beta, wave, shift)
- real(wp), dimension(0:mixlayer_nvar*n - n_bc_skip - 1), intent(in) :: wi !< eigenvalues
- real(wp), dimension(0:mixlayer_nvar*n - n_bc_skip - 1, 0:mixlayer_nvar*n - n_bc_skip - 1), intent(in) :: zr, zi !< eigenvectors
- real(wp), intent(in) :: rho_mean
- real(wp), dimension(mixlayer_nvar, 0:m, 0:n, 0:p), intent(inout) :: wave
- real(wp), intent(in) :: alpha, beta, mach, shift
- real(wp), dimension(0:mixlayer_nvar*n - n_bc_skip - 1) :: vr, vi, vnr, vni !< most unstable eigenvector
- real(wp), dimension(0:mixlayer_nvar*nbp - 1) :: xbr, xbi !< eigenvectors
- real(wp), dimension(0:mixlayer_nvar*(nbp - 1) - 1) :: xcr, xci !< eigenvectors
- real(wp) :: ang, norm
- real(wp) :: tr, ti, cr, ci !< temporary memory
- real(wp) :: xratio
- integer :: idx
- integer :: i, j, k
-
- xratio = mixlayer_vel_coef
-
- ! Find the most unstable eigenvalue and corresponding eigenvector
- k = 0
- do i = 1, mixlayer_nvar*n - n_bc_skip - 1
- if (wi(i) > wi(k)) then
- k = i
- end if
- end do
- vr = zr(:, k)
- vi = zi(:, k)
-
- ! Normalize the eigenvector by its component with the largest modulus.
- norm = 0._wp
- do i = 0, mixlayer_nvar*n - n_bc_skip - 1
- if (sqrt(vr(i)**2 + vi(i)**2) > norm) then
- idx = i
- norm = sqrt(vr(i)**2 + vi(i)**2)
- end if
- end do
-
- tr = vr(idx)
- ti = vi(idx)
- do i = 0, mixlayer_nvar*n - n_bc_skip - 1
- call cdiv(vr(i), vi(i), tr, ti, cr, ci)
- vnr(i) = cr
- vni(i) = ci
- end do
-
- ! Reassign missing values at boundaries based on the boundary condition
- xbr = 0._wp
- xbi = 0._wp
- do i = 1, mixlayer_nvar
- if (i <= mixlayer_var(2)) then
- do k = 0, n - 1
- xbr((i - 1)*nbp + k + 1) = vnr((i - 1)*n + k)
- xbi((i - 1)*nbp + k + 1) = vni((i - 1)*n + k)
- end do
- elseif (i == mixlayer_var(3)) then
- do k = 0, n - n_bc_skip - 1
- xbr((i - 1)*nbp + mixlayer_bc_fd + k + 1) = vnr((i - 1)*n + k)
- xbi((i - 1)*nbp + mixlayer_bc_fd + k + 1) = vni((i - 1)*n + k)
- end do
- else
- do k = 0, n - 1
- xbr((i - 1)*nbp + k + 1) = vnr((i - 1)*n - n_bc_skip + k)
- xbi((i - 1)*nbp + k + 1) = vni((i - 1)*n - n_bc_skip + k)
- end do
- end if
- end do
-
- ! rho at boundaries
- xbr(0) = xbr(1) + xbr(mixlayer_var(2)*nbp + 1)*rho_mean*mach
- xbi(0) = xbi(1) + xbi(mixlayer_var(2)*nbp + 1)*rho_mean*mach
- xbr(nbp - 1) = xbr(n) - xbr(mixlayer_var(2)*nbp + n)*rho_mean*mach
- xbi(nbp - 1) = xbi(n) - xbi(mixlayer_var(2)*nbp + n)*rho_mean*mach
-
- ! u at boundaries
- xbr(mixlayer_var(1)*nbp) = xbr(mixlayer_var(1)*nbp + 1)
- xbi(mixlayer_var(1)*nbp) = xbi(mixlayer_var(1)*nbp + 1)
- xbr(mixlayer_var(1)*nbp + nbp - 1) = xbr(mixlayer_var(1)*nbp + n)
- xbi(mixlayer_var(1)*nbp + nbp - 1) = xbi(mixlayer_var(1)*nbp + n)
-
- ! w at boundaries
- xbr((mixlayer_var(3))*nbp + 0) = xbr((mixlayer_var(3))*nbp + 1)
- xbi((mixlayer_var(3))*nbp + 0) = xbi((mixlayer_var(3))*nbp + 1)
- xbr((mixlayer_var(3))*nbp + nbp - 1) = xbr((mixlayer_var(3))*nbp + n)
- xbi((mixlayer_var(3))*nbp + nbp - 1) = xbi((mixlayer_var(3))*nbp + n)
-
- ! p at boundaries
- xbr(mixlayer_var(4)*nbp + 0) = xbr(mixlayer_var(4)*nbp + 1) + xbr(mixlayer_var(2)*nbp + 1)*rho_mean/mach
- xbi(mixlayer_var(4)*nbp + 0) = xbi(mixlayer_var(4)*nbp + 1) + xbi(mixlayer_var(2)*nbp + 1)*rho_mean/mach
- xbr(mixlayer_var(4)*nbp + nbp - 1) = xbr(mixlayer_var(4)*nbp + n) - xbr(mixlayer_var(2)*nbp + n)*rho_mean/mach
- xbi(mixlayer_var(4)*nbp + nbp - 1) = xbi(mixlayer_var(4)*nbp + n) - xbi(mixlayer_var(2)*nbp + n)*rho_mean/mach
-
- ! Compute average to get cell-centered values
- xcr = 0._wp
- xci = 0._wp
- do i = 1, mixlayer_nvar
- do k = 0, n
- xcr((i - 1)*(nbp - 1) + k) = 5.e-1_wp*(xbr((i - 1)*nbp + k) + xbr((i - 1)*nbp + k + 1))
- xci((i - 1)*(nbp - 1) + k) = 5.e-1_wp*(xbi((i - 1)*nbp + k) + xbi((i - 1)*nbp + k + 1))
- end do
- end do
-
- ! Generate instability waves in x- and z-directions with phase shifts
- ! wave = Re(eigfunc * exp(i*(alpha*x + beta*z)))
- do i = 0, m
- do j = 0, n
- do k = 0, p
- if (f_approx_equal(beta, 0._wp)) then
- ang = alpha*(x_cc(i)*xratio)
- else
- ang = alpha*(x_cc(i)*xratio) + beta*(z_cc(k)*xratio) + shift
- end if
- wave(mixlayer_var(1), i, j, k) = xcr(j)*cos(ang) - xci(j)*sin(ang) ! rho
- wave(mixlayer_var(2), i, j, k) = xcr(mixlayer_var(1)*(nbp - 1) + j)*cos(ang) - xci(mixlayer_var(1)*(nbp - 1) + j)*sin(ang) ! u
- wave(mixlayer_var(3), i, j, k) = xcr(mixlayer_var(2)*(nbp - 1) + j)*cos(ang) - xci(mixlayer_var(2)*(nbp - 1) + j)*sin(ang) ! v
- wave(mixlayer_var(4), i, j, k) = xcr(mixlayer_var(3)*(nbp - 1) + j)*cos(ang) - xci(mixlayer_var(3)*(nbp - 1) + j)*sin(ang) ! w
- wave(mixlayer_var(5), i, j, k) = xcr(mixlayer_var(4)*(nbp - 1) + j)*cos(ang) - xci(mixlayer_var(4)*(nbp - 1) + j)*sin(ang) ! p
- end do
- end do
- end do
-
- end subroutine s_generate_wave
-
impure subroutine s_elliptic_smoothing(q_prim_vf, bc_type)
type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf
@@ -676,6 +150,151 @@ contains
end subroutine s_elliptic_smoothing
+ !> This subroutine computes velocity perturbations for a temporal mixing
+ !! layer with a hyperbolic tangent mean streamwise velocity
+ !! profile, using an inverter version of the spectrum-based
+ !! synthetic turbulence generation method proposed by
+ !! Guo et al. (2023, JFM).
+ subroutine s_perturb_mixlayer(q_prim_vf)
+ type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf
+ real(wp), dimension(mixlayer_perturb_nk) :: k, Ek
+ real(wp), dimension(3, 3) :: Rij, Lmat
+ real(wp), dimension(3) :: velfluc, sig_tmp, sig, khat, xi
+ real(wp) :: dk, alpha, Eksum, q, uu0, phi
+ integer :: i, j, l, r, ierr
+
+ ! Initialize parameters
+ dk = 1._wp/mixlayer_perturb_nk
+
+ ! Compute prescribed energy spectra
+ Eksum = 0._wp
+ do i = 1, mixlayer_perturb_nk
+ k(i) = dk*i
+ Ek(i) = (k(i)/mixlayer_perturb_k0)**4._wp*exp(-2._wp*(k(i)/mixlayer_perturb_k0)**2._wp)
+ Eksum = Eksum + Ek(i)
+ end do
+
+ ! Main loop
+ do r = 0, n
+ ! Compute prescribed Reynolds stress tensor with about half
+ ! magnitude of its self-similar value
+ Rij(:, :) = 0._wp
+ uu0 = patch_icpp(1)%vel(1)**2._wp &
+ *(1._wp - tanh(y_cc(r)*mixlayer_vel_coef)**2._wp)
+ Rij(1, 1) = 0.05_wp*uu0
+ Rij(2, 2) = 0.03_wp*uu0
+ Rij(3, 3) = 0.03_wp*uu0
+ Rij(1, 2) = -0.02_wp*uu0
+ Rij(2, 1) = Rij(1, 2)
+
+ ! Cholesky decomposition for inhomogeneity and anisotropy
+ Lmat = 0._wp
+ Lmat(1, 1) = sqrt(Rij(1, 1))
+ if (abs(Lmat(1, 1)) < sgm_eps) Lmat(1, 1) = sgm_eps
+ Lmat(2, 1) = Rij(2, 1)/Lmat(1, 1)
+ Lmat(2, 2) = sqrt(Rij(2, 2) - Lmat(2, 1)**2._wp)
+ if (abs(Lmat(2, 2)) <