Skip to content

Commit

Permalink
Fix Deardorff overwrite of theta_v component.
Browse files Browse the repository at this point in the history
  • Loading branch information
@AMLattanzi
AMLattanzi committed Jun 11, 2024
1 parent 79bd2a0 commit a333067
Show file tree
Hide file tree
Showing 3 changed files with 180 additions and 44 deletions.
100 changes: 100 additions & 0 deletions Exec/DevTests/Bomex/input_Kessler
View file
Original file line number Diff line number Diff line change
@@ -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
32 changes: 18 additions & 14 deletions Exec/DevTests/Bomex/input_SAM
View file
Original file line number Diff line number Diff line change
@@ -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

Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
Expand Down
92 changes: 62 additions & 30 deletions Source/Diffusion/ComputeTurbulentViscosity.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -281,31 +313,31 @@ 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;
case EddyDiff::KE_h:
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);
// NOTE: Vertical component already set in the valid and ghost regions
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];
});
}
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];
mu_turb(i,j,k,indx_v) = mu_turb(i,j,k,indx);
});
break;
}
Expand Down Expand Up @@ -342,37 +374,37 @@ 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;
case EddyDiff::KE_h:
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;
}
Expand Down

0 comments on commit a333067

Please sign in to comment.