Skip to content

Commit

Permalink
update driving params
Browse files Browse the repository at this point in the history
  • Loading branch information
@BenWibking
BenWibking committed Jun 13, 2024
1 parent 673d89b commit 0a2a063
Show file tree
Hide file tree
Showing 2 changed files with 27 additions and 33 deletions.
6 changes: 3 additions & 3 deletions inputs/precipitator_lowres.in
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,12 +11,12 @@ dt = 50.0 # time increment between outputs
<parthenon/output2>
file_type = hdf5 # HDF5 data dump
variables = prim, drho_over_rho, dP_over_P, dK_over_K, dT_over_T, entropy, temperature, tcool_over_tff, dv_x, dv_y, dv_z, accel_x, accel_y, accel_z
dt = 100.0 # Time increment between outputs
dt = 50.0 # Time increment between outputs
id = prim # Name to append to output

<parthenon/output3>
file_type = rst # Binary data dump
dt = 2000.0 # time increment between outputs
dt = 250.0 # time increment between outputs
id = restart

<parthenon/time>
Expand Down Expand Up @@ -107,5 +107,5 @@ num_modes = 12
vertical_driving_only = true
sol_weight = 1.0
rseed = 42
t_corr = 50.0 # Myr
t_corr = 500.0 # Myr
max_height = 80.0 # kpc
54 changes: 24 additions & 30 deletions src/pgen/precipitator.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -296,8 +296,8 @@ void TurbSrcTerm(MeshData<Real> *md, const parthenon::SimTime /*time*/, const Re

if (sigma_v > 0) {
// generate perturbations
auto few_modes_ft = hydro_pkg->Param<FewModesFT>("precipitator/few_modes_ft_v");
few_modes_ft.Generate(md, dt, "tmp_perturb");
auto *few_modes_ft = hydro_pkg->MutableParam<FewModesFT>("precipitator/few_modes_ft_v");
few_modes_ft->Generate(md, dt, "tmp_perturb");

// normalize perturbations
Real v2_sum{};
Expand Down Expand Up @@ -704,12 +704,10 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *pkg
************************************************************/
const auto &filename = pin->GetString("precipitator", "hse_profile_filename");
const auto &uniform_init = pin->GetInteger("precipitator", "uniform_init");
hydro_pkg->AddParam<>("uniform_init", uniform_init,
parthenon::Params::Mutability::Restart);
hydro_pkg->AddParam<>("uniform_init", uniform_init);
if (uniform_init == 1) {
const auto &uniform_init_height = pin->GetReal("precipitator", "uniform_init_height");
hydro_pkg->AddParam<>("uniform_init_height", uniform_init_height,
parthenon::Params::Mutability::Restart);
hydro_pkg->AddParam<>("uniform_init_height", uniform_init_height);
}

const PrecipitatorProfile P_rho_profile(filename);
Expand All @@ -718,16 +716,14 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *pkg

const auto enable_heating_str =
pin->GetOrAddString("precipitator", "enable_heating", "none");
hydro_pkg->AddParam<>("enable_heating", enable_heating_str,
parthenon::Params::Mutability::Restart);
hydro_pkg->AddParam<>("enable_heating", enable_heating_str);

// read smoothing height for heating/cooling
const Real h_smooth_heatcool =
pin->GetReal("precipitator", "h_smooth_heatcool_kpc") * units.kpc();

hydro_pkg->AddParam<Real>(
"h_smooth_heatcool", h_smooth_heatcool,
parthenon::Params::Mutability::Restart); // smoothing scale (code units)
"h_smooth_heatcool", h_smooth_heatcool); // smoothing scale (code units)

// read perturbation parameters
const int kx_max = pin->GetInteger("precipitator", "perturb_kx");
Expand All @@ -736,13 +732,11 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *pkg
const int expo = pin->GetInteger("precipitator", "perturb_exponent");
const Real amp = pin->GetReal("precipitator", "perturb_sin_drho_over_rho");

hydro_pkg->AddParam<int>("perturb_kx", kx_max, parthenon::Params::Mutability::Restart);
hydro_pkg->AddParam<int>("perturb_ky", ky_max, parthenon::Params::Mutability::Restart);
hydro_pkg->AddParam<int>("perturb_kz", kz_max, parthenon::Params::Mutability::Restart);
hydro_pkg->AddParam<int>("perturb_exponent", expo,
parthenon::Params::Mutability::Restart);
hydro_pkg->AddParam<Real>("perturb_amplitude", amp,
parthenon::Params::Mutability::Restart);
hydro_pkg->AddParam<int>("perturb_kx", kx_max);
hydro_pkg->AddParam<int>("perturb_ky", ky_max);
hydro_pkg->AddParam<int>("perturb_kz", kz_max);
hydro_pkg->AddParam<int>("perturb_exponent", expo);
hydro_pkg->AddParam<Real>("perturb_amplitude", amp);

// density perturbations
if (amp > 0.) {
Expand Down Expand Up @@ -819,7 +813,7 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *pkg

// setup velocity perturbations
const auto sigma_v = pin->GetOrAddReal("precipitator/driving", "sigma_v", 0.0);
hydro_pkg->AddParam<>("sigma_v", sigma_v, parthenon::Params::Mutability::Restart);
hydro_pkg->AddParam<>("sigma_v", sigma_v);

if (sigma_v > 0) {
// the maximum height at which to drive turbulence
Expand Down Expand Up @@ -897,14 +891,14 @@ void ProblemGenerator(MeshBlock *pmb, parthenon::ParameterInput *pin) {
}

// Get HSE profile and parameters
const auto &P_rho_profile =
hydro_pkg->Param<PrecipitatorProfile>("precipitator_profile");
const auto *P_rho_profile =
hydro_pkg->MutableParam<PrecipitatorProfile>("precipitator_profile");

// Initialize phase factors for velocity perturbations
const auto sigma_v = hydro_pkg->Param<Real>("sigma_v");
if (sigma_v > 0) {
auto few_modes_ft = hydro_pkg->Param<FewModesFT>("precipitator/few_modes_ft_v");
few_modes_ft.SetPhases(pmb, pin);
auto *few_modes_ft = hydro_pkg->MutableParam<FewModesFT>("precipitator/few_modes_ft_v");
few_modes_ft->SetPhases(pmb, pin);
}

// Get (density) perturbation parameters
Expand Down Expand Up @@ -970,7 +964,7 @@ void ProblemGenerator(MeshBlock *pmb, parthenon::ParameterInput *pin) {
// Calculate height
const Real zcen = coords.Xc<3>(k);
const Real zcen_cgs = std::abs(zcen) * code_length_cgs;
const Real phi_i = P_rho_profile.phi(zcen_cgs);
const Real phi_i = P_rho_profile->phi(zcen_cgs);
grav_phi(0, k, j, i) = phi_i / code_potential_cgs;
});
}
Expand Down Expand Up @@ -999,7 +993,7 @@ void ProblemGenerator(MeshBlock *pmb, parthenon::ParameterInput *pin) {
KOKKOS_LAMBDA(const int, const int k, const int j, const int i) {
// Calculate height
const Real zmin_cgs = std::abs(coords.Xf<3>(k)) * code_length_cgs;
const Real phi_iminus = P_rho_profile.phi(zmin_cgs) / code_potential_cgs;
const Real phi_iminus = P_rho_profile->phi(zmin_cgs) / code_potential_cgs;
grav_phi_zface(0, k, j, i) = phi_iminus;
});
}
Expand All @@ -1013,8 +1007,8 @@ void ProblemGenerator(MeshBlock *pmb, parthenon::ParameterInput *pin) {
DEFAULT_LOOP_PATTERN, "SetHydrostaticProfileCells", parthenon::DevExecSpace(), 0,
0, kbp.s, kbp.e, jbp.s, jbp.e, ibp.s, ibp.e,
KOKKOS_LAMBDA(const int, const int k, const int j, const int i) {
auto p_hse = [=](Real z) { return P_rho_profile.P(z); };
auto rho_hse = [=](Real z) { return P_rho_profile.rho(z); };
auto p_hse = [=](Real z) { return P_rho_profile->P(z); };
auto rho_hse = [=](Real z) { return P_rho_profile->rho(z); };
pressure_hse(0, k, j, i) = p_hse(uniform_init_height) / code_pressure_cgs;
density_hse(0, k, j, i) = rho_hse(uniform_init_height) / code_density_cgs;
});
Expand All @@ -1032,8 +1026,8 @@ void ProblemGenerator(MeshBlock *pmb, parthenon::ParameterInput *pin) {
const Real dz_cgs = dx3 * code_length_cgs;

parthenon::math::quadrature::gauss<Real, 7> quad;
auto p_hse = [=](Real z) { return P_rho_profile.P(z); };
auto rho_hse = [=](Real z) { return P_rho_profile.rho(z); };
auto p_hse = [=](Real z) { return P_rho_profile->P(z); };
auto rho_hse = [=](Real z) { return P_rho_profile->rho(z); };
const Real P_hse_avg = quad.integrate(p_hse, zmin_cgs, zmax_cgs) / dz_cgs;
const Real rho_hse_avg = quad.integrate(rho_hse, zmin_cgs, zmax_cgs) / dz_cgs;

Expand Down Expand Up @@ -1077,8 +1071,8 @@ void UserMeshWorkBeforeOutput(Mesh *mesh, ParameterInput *pin,

const auto sigma_v = pkg->Param<Real>("sigma_v");
if (sigma_v > 0) {
auto few_modes_ft = pkg->Param<FewModesFT>("precipitator/few_modes_ft_v");
few_modes_ft.SaveStateBeforeOutput(mesh, pin);
auto *few_modes_ft = pkg->MutableParam<FewModesFT>("precipitator/few_modes_ft_v");
few_modes_ft->SaveStateBeforeOutput(mesh, pin);
}

const Units units(pin);
Expand Down

0 comments on commit 0a2a063

Please sign in to comment.