From 33417f70e5c975a81b729e7899b5b8375daf4b3c Mon Sep 17 00:00:00 2001 From: "Aaron M. Lattanzi" <103702284+AMLattanzi@users.noreply.github.com> Date: Tue, 11 Jun 2024 15:22:10 -0700 Subject: [PATCH] Fix Deardorff overwrite of theta_v component. (#1641) * Fix Deardorff overwrite of theta_v component. * fix indentation. * fix logic error. --- Exec/DevTests/Bomex/input_Kessler | 100 +++++++ Exec/DevTests/Bomex/input_SAM | 32 ++- .../Diffusion/ComputeTurbulentViscosity.cpp | 263 ++++++++++-------- 3 files changed, 266 insertions(+), 129 deletions(-) create mode 100644 Exec/DevTests/Bomex/input_Kessler diff --git a/Exec/DevTests/Bomex/input_Kessler b/Exec/DevTests/Bomex/input_Kessler new file mode 100644 index 000000000..00ec1ee2f --- /dev/null +++ b/Exec/DevTests/Bomex/input_Kessler @@ -0,0 +1,100 @@ +# ------------------ INPUTS TO MAIN PROGRAM ------------------- +stop_time = 43200 + +amrex.fpe_trap_invalid = 1 + +fabarray.mfiter_tile_size = 1024 1024 1024 + +# PROBLEM SIZE & GEOMETRY +geometry.prob_extent = 3200 3200 4000 +amr.n_cell = 32 32 100 + +geometry.is_periodic = 1 1 0 + +# MOST BOUNDARY (DEFAULT IS ADIABATIC FOR THETA) +zlo.type = "Most" +erf.most.flux_type = "custom" +erf.most.ustar = 0.28 # ustar +erf.most.tstar = 8.0e-3 # theta flux +erf.most.qstar = 5.2e-5 # qv flux +erf.most.z0 = 0.1 +erf.most.zref = 20.0 + +// NOTE: This should have a qv grad too (use hoextrapcc?!) +zhi.type = "SlipWall" +zhi.theta_grad = 0.00365 + +# TIME STEP CONTROL +erf.fixed_dt = 0.5 # fixed time step depending on grid resolution + +# DIAGNOSTICS & VERBOSITY +erf.sum_interval = 1 # timesteps between computing mass +erf.v = 1 # verbosity in ERF.cpp +amr.v = 1 # verbosity in Amr.cpp +erf.data_log = "surf" "mean" "flux" "subgrid" +erf.profile_int = 120 + +# REFINEMENT / REGRIDDING +amr.max_level = 0 # maximum level number allowed + +# CHECKPOINT FILES +erf.check_file = chk # root name of checkpoint file +erf.check_int = 600 # number of timesteps between checkpoints + +# PLOTFILES +erf.plot_file_1 = plt # prefix of plotfile name +erf.plot_int_1 = 120 # number of timesteps between plotfiles +erf.plot_vars_1 = density rhotheta x_velocity y_velocity z_velocity pressure temp theta qt qp qv qc qi + +# SOLVER CHOICE +erf.alpha_T = 0.0 +erf.alpha_C = 0.0 +erf.use_gravity = true + +erf.dycore_horiz_adv_type = Upwind_5th +erf.dycore_vert_adv_type = Upwind_5th +erf.dryscal_horiz_adv_type = WENOZ5 +erf.dryscal_vert_adv_type = WENOZ5 +erf.moistscal_horiz_adv_type = WENOZ5 +erf.moistscal_vert_adv_type = WENOZ5 + +erf.moisture_model = "Kessler" + +erf.molec_diff_type = "None" +erf.les_type = "Smagorinsky" +erf.Cs = 0.17 +erf.Pr_t = 0.333333 +erf.Sc_t = 0.333333 + +erf.init_type = "input_sounding" +erf.init_sounding_ideal = true + +erf.add_custom_rhotheta_forcing = true +erf.add_custom_moisture_forcing = true +erf.add_custom_geostrophic_profile = true +erf.add_custom_w_subsidence = true +erf.custom_forcing_uses_primitive_vars = false + +# Higher values of perturbations lead to instability +# Instability seems to be coming from BC +prob.U_0_Pert_Mag = 0.01 +prob.V_0_Pert_Mag = 0.01 +prob.W_0_Pert_Mag = 0.0 + +prob.pert_ref_height = 1600.0 +prob.T_0_Pert_Mag = 0.1 +prob.qv_0_Pert_Mag = 0.000025 + +prob.advection_heating_rate = -2.3148E-5 +prob.source_cutoff = 1500.0 +prob.source_cutoff_transition = 1500.0 + +prob.advection_moisture_rate = -1.2E-8 +prob.moisture_source_cutoff = 300.0 +prob.moisture_source_cutoff_transition = 200.0 + +prob.wbar_sub_max = -0.0065 +prob.wbar_cutoff_max = 1500.0 +prob.wbar_cutoff_min = 2100.0 + +prob.custom_TKE = true diff --git a/Exec/DevTests/Bomex/input_SAM b/Exec/DevTests/Bomex/input_SAM index a108e3de5..9d30122eb 100644 --- a/Exec/DevTests/Bomex/input_SAM +++ b/Exec/DevTests/Bomex/input_SAM @@ -1,13 +1,13 @@ # ------------------ INPUTS TO MAIN PROGRAM ------------------- -stop_time = 21600 +stop_time = 43200 amrex.fpe_trap_invalid = 1 fabarray.mfiter_tile_size = 1024 1024 1024 # PROBLEM SIZE & GEOMETRY -geometry.prob_extent = 6400 6400 4000 -amr.n_cell = 64 64 100 +geometry.prob_extent = 3200 3200 4000 +amr.n_cell = 32 32 100 geometry.is_periodic = 1 1 0 @@ -17,20 +17,22 @@ erf.most.flux_type = "custom" erf.most.ustar = 0.28 # ustar erf.most.tstar = 8.0e-3 # theta flux erf.most.qstar = 5.2e-5 # qv flux -erf.most.z0 = 0.1 +erf.most.z0 = 0.1 erf.most.zref = 20.0 +// NOTE: This should have a qv grad too (use hoextrapcc?!) zhi.type = "SlipWall" - +zhi.theta_grad = 0.00365 + # TIME STEP CONTROL -erf.fixed_dt = 1.0 # fixed time step depending on grid resolution +erf.fixed_dt = 0.5 # fixed time step depending on grid resolution # DIAGNOSTICS & VERBOSITY erf.sum_interval = 1 # timesteps between computing mass erf.v = 1 # verbosity in ERF.cpp amr.v = 1 # verbosity in Amr.cpp erf.data_log = "surf" "mean" "flux" "subgrid" -erf.profile_int = 60 +erf.profile_int = 120 # REFINEMENT / REGRIDDING amr.max_level = 0 # maximum level number allowed @@ -41,12 +43,12 @@ erf.check_int = 600 # number of timesteps between checkpoints # PLOTFILES erf.plot_file_1 = plt # prefix of plotfile name -erf.plot_int_1 = 60 # number of timesteps between plotfiles +erf.plot_int_1 = 120 # number of timesteps between plotfiles erf.plot_vars_1 = density rhotheta x_velocity y_velocity z_velocity pressure temp theta qt qp qv qc qi # SOLVER CHOICE erf.alpha_T = 0.0 -erf.alpha_C = 1.0 +erf.alpha_C = 0.0 erf.use_gravity = true erf.dycore_horiz_adv_type = Upwind_5th @@ -60,15 +62,17 @@ erf.moisture_model = "SAM" erf.molec_diff_type = "None" erf.les_type = "Smagorinsky" -erf.Cs = 0.1 +erf.Cs = 0.17 +erf.Pr_t = 0.333333 +erf.Sc_t = 0.333333 erf.init_type = "input_sounding" erf.init_sounding_ideal = true -erf.add_custom_rhotheta_forcing = true -erf.add_custom_moisture_forcing = true -erf.add_custom_geostrophic_profile = true -erf.add_custom_w_subsidence = true +erf.add_custom_rhotheta_forcing = true +erf.add_custom_moisture_forcing = true +erf.add_custom_geostrophic_profile = true +erf.add_custom_w_subsidence = true erf.custom_forcing_uses_primitive_vars = false # Higher values of perturbations lead to instability diff --git a/Source/Diffusion/ComputeTurbulentViscosity.cpp b/Source/Diffusion/ComputeTurbulentViscosity.cpp index db7742162..cf681c12a 100644 --- a/Source/Diffusion/ComputeTurbulentViscosity.cpp +++ b/Source/Diffusion/ComputeTurbulentViscosity.cpp @@ -110,101 +110,101 @@ void ComputeTurbulentViscosityLES (const MultiFab& Tau11, const MultiFab& Tau22, //*********************************************************************************** else if (turbChoice.les_type == LESType::Deardorff) { - const Real l_C_k = turbChoice.Ck; - const Real l_C_e = turbChoice.Ce; - const Real l_C_e_wall = turbChoice.Ce_wall; - const Real Ce_lcoeff = amrex::max(0.0, l_C_e - 1.9*l_C_k); - const Real l_abs_g = const_grav; - const Real l_inv_theta0 = 1.0 / turbChoice.theta_ref; + const Real l_C_k = turbChoice.Ck; + const Real l_C_e = turbChoice.Ce; + const Real l_C_e_wall = turbChoice.Ce_wall; + const Real Ce_lcoeff = amrex::max(0.0, l_C_e - 1.9*l_C_k); + const Real l_abs_g = const_grav; + const Real l_inv_theta0 = 1.0 / turbChoice.theta_ref; #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif - for ( MFIter mfi(eddyViscosity,TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - Box bxcc = mfi.tilebox(); + for ( MFIter mfi(eddyViscosity,TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + Box bxcc = mfi.tilebox(); - const Array4& mu_turb = eddyViscosity.array(mfi); - const Array4& hfx_x = Hfx1.array(mfi); - const Array4& hfx_y = Hfx2.array(mfi); - const Array4& hfx_z = Hfx3.array(mfi); - const Array4& diss = Diss.array(mfi); + const Array4& mu_turb = eddyViscosity.array(mfi); + const Array4& hfx_x = Hfx1.array(mfi); + const Array4& hfx_y = Hfx2.array(mfi); + const Array4& hfx_z = Hfx3.array(mfi); + const Array4& diss = Diss.array(mfi); - const Array4 &cell_data = cons_in.array(mfi); + const Array4 &cell_data = cons_in.array(mfi); - Array4 mf_u = mapfac_u.array(mfi); - Array4 mf_v = mapfac_v.array(mfi); + Array4 mf_u = mapfac_u.array(mfi); + Array4 mf_v = mapfac_v.array(mfi); - Array4 z_nd_arr = (use_terrain) ? z_phys_nd->const_array(mfi) : Array4{}; + Array4 z_nd_arr = (use_terrain) ? z_phys_nd->const_array(mfi) : Array4{}; - ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Real dxInv = cellSizeInv[0]; - Real dyInv = cellSizeInv[1]; - Real dzInv = cellSizeInv[2]; - if (use_terrain) { - // the terrain grid is only deformed in z for now - dzInv /= Compute_h_zeta_AtCellCenter(i,j,k, cellSizeInv, z_nd_arr); - } - Real cellVolMsf = 1.0 / (dxInv * mf_u(i,j,0) * dyInv * mf_v(i,j,0) * dzInv); - Real DeltaMsf = std::pow(cellVolMsf,1.0/3.0); - - // Calculate stratification-dependent mixing length (Deardorff 1980) - Real eps = std::numeric_limits::epsilon(); - Real dtheta_dz; - if (use_most && k==klo) { - if (exp_most) { - dtheta_dz = ( cell_data(i,j,k+1,RhoTheta_comp)/cell_data(i,j,k+1,Rho_comp) - - cell_data(i,j,k ,RhoTheta_comp)/cell_data(i,j,k ,Rho_comp) )*dzInv; - } else { - dtheta_dz = 0.5 * (-3 * cell_data(i,j,k ,RhoTheta_comp) - / cell_data(i,j,k ,Rho_comp) - + 4 * cell_data(i,j,k+1,RhoTheta_comp) - / cell_data(i,j,k+1,Rho_comp) - - cell_data(i,j,k+2,RhoTheta_comp) - / cell_data(i,j,k+2,Rho_comp) ) * dzInv; - } - } else { - dtheta_dz = 0.5 * ( cell_data(i,j,k+1,RhoTheta_comp)/cell_data(i,j,k+1,Rho_comp) - - cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp) )*dzInv; - } - Real E = cell_data(i,j,k,RhoKE_comp) / cell_data(i,j,k,Rho_comp); - Real strat = l_abs_g * dtheta_dz * l_inv_theta0; // stratification - Real length; - if (strat <= eps) { - length = DeltaMsf; - } else { - length = 0.76 * std::sqrt(E / strat); - // mixing length should be _reduced_ for stable stratification - length = amrex::min(length, DeltaMsf); - // following WRF, make sure the mixing length isn't too small - length = amrex::max(length, 0.001 * DeltaMsf); - } + ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + Real dxInv = cellSizeInv[0]; + Real dyInv = cellSizeInv[1]; + Real dzInv = cellSizeInv[2]; + if (use_terrain) { + // the terrain grid is only deformed in z for now + dzInv /= Compute_h_zeta_AtCellCenter(i,j,k, cellSizeInv, z_nd_arr); + } + Real cellVolMsf = 1.0 / (dxInv * mf_u(i,j,0) * dyInv * mf_v(i,j,0) * dzInv); + Real DeltaMsf = std::pow(cellVolMsf,1.0/3.0); + + // Calculate stratification-dependent mixing length (Deardorff 1980) + Real eps = std::numeric_limits::epsilon(); + Real dtheta_dz; + if (use_most && k==klo) { + if (exp_most) { + dtheta_dz = ( cell_data(i,j,k+1,RhoTheta_comp)/cell_data(i,j,k+1,Rho_comp) + - cell_data(i,j,k ,RhoTheta_comp)/cell_data(i,j,k ,Rho_comp) )*dzInv; + } else { + dtheta_dz = 0.5 * (-3 * cell_data(i,j,k ,RhoTheta_comp) + / cell_data(i,j,k ,Rho_comp) + + 4 * cell_data(i,j,k+1,RhoTheta_comp) + / cell_data(i,j,k+1,Rho_comp) + - cell_data(i,j,k+2,RhoTheta_comp) + / cell_data(i,j,k+2,Rho_comp) ) * dzInv; + } + } else { + dtheta_dz = 0.5 * ( cell_data(i,j,k+1,RhoTheta_comp)/cell_data(i,j,k+1,Rho_comp) + - cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp) )*dzInv; + } + Real E = cell_data(i,j,k,RhoKE_comp) / cell_data(i,j,k,Rho_comp); + Real strat = l_abs_g * dtheta_dz * l_inv_theta0; // stratification + Real length; + if (strat <= eps) { + length = DeltaMsf; + } else { + length = 0.76 * std::sqrt(E / strat); + // mixing length should be _reduced_ for stable stratification + length = amrex::min(length, DeltaMsf); + // following WRF, make sure the mixing length isn't too small + length = amrex::max(length, 0.001 * DeltaMsf); + } - // Calculate eddy diffusivities - // K = rho * C_k * l * KE^(1/2) - mu_turb(i,j,k,EddyDiff::Mom_h) = cell_data(i,j,k,Rho_comp) * l_C_k * length * std::sqrt(E); - mu_turb(i,j,k,EddyDiff::Mom_v) = mu_turb(i,j,k,EddyDiff::Mom_h); - // KH = (1 + 2*l/delta) * mu_turb - mu_turb(i,j,k,EddyDiff::Theta_v) = (1.+2.*length/DeltaMsf) * mu_turb(i,j,k,EddyDiff::Mom_v); - - // Calculate SFS quantities - // - dissipation - Real Ce; - if ((l_C_e_wall > 0) && (k==0)) { - Ce = l_C_e_wall; - } else { - Ce = 1.9*l_C_k + Ce_lcoeff*length / DeltaMsf; - } - diss(i,j,k) = cell_data(i,j,k,Rho_comp) * Ce * std::pow(E,1.5) / length; - // - heat flux - // (Note: If using ERF_EXPLICIT_MOST_STRESS, the value at k=0 will - // be overwritten when BCs are applied) - hfx_x(i,j,k) = 0.0; - hfx_y(i,j,k) = 0.0; - hfx_z(i,j,k) = -mu_turb(i,j,k,EddyDiff::Theta_v) * dtheta_dz; // (rho*w)' theta' [kg m^-2 s^-1 K] - }); - } + // Calculate eddy diffusivities + // K = rho * C_k * l * KE^(1/2) + mu_turb(i,j,k,EddyDiff::Mom_h) = cell_data(i,j,k,Rho_comp) * l_C_k * length * std::sqrt(E); + mu_turb(i,j,k,EddyDiff::Mom_v) = mu_turb(i,j,k,EddyDiff::Mom_h); + // KH = (1 + 2*l/delta) * mu_turb + mu_turb(i,j,k,EddyDiff::Theta_v) = (1.+2.*length/DeltaMsf) * mu_turb(i,j,k,EddyDiff::Mom_v); + + // Calculate SFS quantities + // - dissipation + Real Ce; + if ((l_C_e_wall > 0) && (k==0)) { + Ce = l_C_e_wall; + } else { + Ce = 1.9*l_C_k + Ce_lcoeff*length / DeltaMsf; + } + diss(i,j,k) = cell_data(i,j,k,Rho_comp) * Ce * std::pow(E,1.5) / length; + // - heat flux + // (Note: If using ERF_EXPLICIT_MOST_STRESS, the value at k=0 will + // be overwritten when BCs are applied) + hfx_x(i,j,k) = 0.0; + hfx_y(i,j,k) = 0.0; + hfx_z(i,j,k) = -mu_turb(i,j,k,EddyDiff::Theta_v) * dtheta_dz; // (rho*w)' theta' [kg m^-2 s^-1 K] + }); + } } // Extrapolate Kturb in x/y, fill remaining elements (relevent to lev==0) @@ -271,6 +271,38 @@ void ComputeTurbulentViscosityLES (const MultiFab& Tau11, const MultiFab& Tau22, }); } + // Copy Theta_v component into lateral ghost cells if using Deardorff (populated above) + if (use_KE) { + if (i_lo == domain.smallEnd(0)) { + ParallelFor(planex, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + int lj = amrex::min(amrex::max(j, domain.smallEnd(1)), domain.bigEnd(1)); + mu_turb(i_lo-i, j, k, EddyDiff::Theta_v) = mu_turb(i_lo, lj, k, EddyDiff::Theta_v); + }); + } + if (i_hi == domain.bigEnd(0)) { + ParallelFor(planex, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + int lj = amrex::min(amrex::max(j, domain.smallEnd(1)), domain.bigEnd(1)); + mu_turb(i_hi+i, j, k, EddyDiff::Theta_v) = mu_turb(i_hi, lj, k, EddyDiff::Theta_v); + }); + } + if (j_lo == domain.smallEnd(1)) { + ParallelFor(planey, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + int li = amrex::min(amrex::max(i, domain.smallEnd(0)), domain.bigEnd(0)); + mu_turb(i, j_lo-j, k, EddyDiff::Theta_v) = mu_turb(li, j_lo, k, EddyDiff::Theta_v); + }); + } + if (j_hi == domain.bigEnd(1)) { + ParallelFor(planey, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + int li = amrex::min(amrex::max(i, domain.smallEnd(0)), domain.bigEnd(0)); + mu_turb(i, j_hi+j, k, EddyDiff::Theta_v) = mu_turb(li, j_hi, k, EddyDiff::Theta_v); + }); + } + } + // refactor the code to eliminate the need for ifdef's for (auto n = 1; n < (EddyDiff::NumDiffs-1)/2; ++n) { int offset = (EddyDiff::NumDiffs-1)/2; @@ -281,10 +313,10 @@ void ComputeTurbulentViscosityLES (const MultiFab& Tau11, const MultiFab& Tau22, if(use_QKE) { ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - int indx = n; - int indx_v = indx + offset; - mu_turb(i,j,k,indx) = mu_turb(i,j,k,EddyDiff::Mom_h) * fac_ptr[indx-1]; - mu_turb(i,j,k,indx_v) = mu_turb(i,j,k,indx); + int indx = n; + int indx_v = indx + offset; + mu_turb(i,j,k,indx) = mu_turb(i,j,k,EddyDiff::Mom_h) * fac_ptr[indx-1]; + mu_turb(i,j,k,indx_v) = mu_turb(i,j,k,indx); }); } break; @@ -292,20 +324,21 @@ void ComputeTurbulentViscosityLES (const MultiFab& Tau11, const MultiFab& Tau22, if (use_KE) { ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - int indx = n; - int indx_v = indx + offset; - mu_turb(i,j,k,indx) = mu_turb(i,j,k,EddyDiff::Mom_h) * fac_ptr[indx-1]; - mu_turb(i,j,k,indx_v) = mu_turb(i,j,k,indx); + int indx = n; + int indx_v = indx + offset; + mu_turb(i,j,k,indx) = mu_turb(i,j,k,EddyDiff::Mom_h) * fac_ptr[indx-1]; + mu_turb(i,j,k,indx_v) = mu_turb(i,j,k,indx); }); } break; default: ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - int indx = n; - int indx_v = indx + offset; - mu_turb(i,j,k,indx) = mu_turb(i,j,k,EddyDiff::Mom_h) * fac_ptr[indx-1]; - mu_turb(i,j,k,indx_v) = mu_turb(i,j,k,indx); + int indx = n; + int indx_v = indx + offset; + mu_turb(i,j,k,indx) = mu_turb(i,j,k,EddyDiff::Mom_h) * fac_ptr[indx-1]; + // NOTE: Theta_v has already been set for Deardorff + if (!(indx_v == EddyDiff::Theta_v && use_KE)) mu_turb(i,j,k,indx_v) = mu_turb(i,j,k,indx); }); break; } @@ -342,12 +375,12 @@ void ComputeTurbulentViscosityLES (const MultiFab& Tau11, const MultiFab& Tau22, if(use_QKE) { ParallelFor(planez, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - int indx = n; - int indx_v = indx + offset; - mu_turb(i, j, k_lo-k, indx ) = mu_turb(i, j, k_lo, indx ); - mu_turb(i, j, k_hi+k, indx ) = mu_turb(i, j, k_hi, indx ); - mu_turb(i, j, k_lo-k, indx_v) = mu_turb(i, j, k_lo, indx_v); - mu_turb(i, j, k_hi+k, indx_v) = mu_turb(i, j, k_hi, indx_v); + int indx = n; + int indx_v = indx + offset; + mu_turb(i, j, k_lo-k, indx ) = mu_turb(i, j, k_lo, indx ); + mu_turb(i, j, k_hi+k, indx ) = mu_turb(i, j, k_hi, indx ); + mu_turb(i, j, k_lo-k, indx_v) = mu_turb(i, j, k_lo, indx_v); + mu_turb(i, j, k_hi+k, indx_v) = mu_turb(i, j, k_hi, indx_v); }); } break; @@ -355,24 +388,24 @@ void ComputeTurbulentViscosityLES (const MultiFab& Tau11, const MultiFab& Tau22, if (use_KE) { ParallelFor(planez, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - int indx = n; - int indx_v = indx + offset; - mu_turb(i, j, k_lo-k, indx ) = mu_turb(i, j, k_lo, indx ); - mu_turb(i, j, k_hi+k, indx ) = mu_turb(i, j, k_hi, indx ); - mu_turb(i, j, k_lo-k, indx_v) = mu_turb(i, j, k_lo, indx_v); - mu_turb(i, j, k_hi+k, indx_v) = mu_turb(i, j, k_hi, indx_v); + int indx = n; + int indx_v = indx + offset; + mu_turb(i, j, k_lo-k, indx ) = mu_turb(i, j, k_lo, indx ); + mu_turb(i, j, k_hi+k, indx ) = mu_turb(i, j, k_hi, indx ); + mu_turb(i, j, k_lo-k, indx_v) = mu_turb(i, j, k_lo, indx_v); + mu_turb(i, j, k_hi+k, indx_v) = mu_turb(i, j, k_hi, indx_v); }); } break; default: ParallelFor(planez, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - int indx = n ; - int indx_v = indx + offset; - mu_turb(i, j, k_lo-k, indx ) = mu_turb(i, j, k_lo, indx ); - mu_turb(i, j, k_hi+k, indx ) = mu_turb(i, j, k_hi, indx ); - mu_turb(i, j, k_lo-k, indx_v) = mu_turb(i, j, k_lo, indx_v); - mu_turb(i, j, k_hi+k, indx_v) = mu_turb(i, j, k_hi, indx_v); + int indx = n; + int indx_v = indx + offset; + mu_turb(i, j, k_lo-k, indx ) = mu_turb(i, j, k_lo, indx ); + mu_turb(i, j, k_hi+k, indx ) = mu_turb(i, j, k_hi, indx ); + mu_turb(i, j, k_lo-k, indx_v) = mu_turb(i, j, k_lo, indx_v); + mu_turb(i, j, k_hi+k, indx_v) = mu_turb(i, j, k_hi, indx_v); }); break; }