Skip to content

Commit

Permalink
add Metadata::Restart
Browse files Browse the repository at this point in the history
  • Loading branch information
@BenWibking
BenWibking committed Mar 5, 2024
1 parent 96ff9be commit 3c4f7e1
Show file tree
Hide file tree
Showing 2 changed files with 34 additions and 34 deletions.
50 changes: 27 additions & 23 deletions src/pgen/precipitator.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -570,15 +570,19 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *pkg
auto &hydro_pkg = pkg;

/// add gravitational potential field
auto m = Metadata({Metadata::Cell, Metadata::OneCopy}, std::vector<int>({1}));
auto m = Metadata({Metadata::Cell, Metadata::OneCopy, Metadata::Restart},
std::vector<int>({1}));
pkg->AddField("grav_phi", m);
m = Metadata({Metadata::Cell, Metadata::OneCopy}, std::vector<int>({1}));
m = Metadata({Metadata::Cell, Metadata::OneCopy, Metadata::Restart},
std::vector<int>({1}));
pkg->AddField("grav_phi_zface", m);

// add hydrostatic pressure, density fields
m = Metadata({Metadata::Cell, Metadata::OneCopy}, std::vector<int>({1}));
m = Metadata({Metadata::Cell, Metadata::OneCopy, Metadata::Restart},
std::vector<int>({1}));
pkg->AddField("pressure_hse", m);
m = Metadata({Metadata::Cell, Metadata::OneCopy}, std::vector<int>({1}));
m = Metadata({Metadata::Cell, Metadata::OneCopy, Metadata::Restart},
std::vector<int>({1}));
pkg->AddField("density_hse", m);

/// add derived fields
Expand Down Expand Up @@ -634,7 +638,7 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *pkg
const Real L = std::min({x1max - x1min, x2max - x2min, x3max - x3min});

const Real gam = pin->GetReal("hydro", "gamma");
hydro_pkg->AddParam("gamma", gam); // adiabatic index
hydro_pkg->AddParam("gamma", gam, parthenon::Params::Mutability::Restart); // adiabatic index

/************************************************************
* Initialize magic heating
Expand All @@ -643,44 +647,44 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *pkg
// store PI error integral over time
parthenon::ParArray1D<Real> PI_error_integral("PI_error_integral",
REDUCTION_ARRAY_SIZE);
pkg->AddParam("PI_error_integral", PI_error_integral, true);
pkg->AddParam("PI_error_integral", PI_error_integral, parthenon::Params::Mutability::Restart);

// feedback loop target temperature
const Real PI_target_T = pin->GetReal("precipitator", "thermostat_temperature");
pkg->AddParam("PI_controller_temperature", PI_target_T);
pkg->AddParam("PI_controller_temperature", PI_target_T, parthenon::Params::Mutability::Restart);

// K_p feedback loop constant [dimensionless]
const Real PI_Kp = pin->GetReal("precipitator", "thermostat_Kp");
pkg->AddParam("PI_controller_Kp", PI_Kp);
pkg->AddParam("PI_controller_Kp", PI_Kp, parthenon::Params::Mutability::Restart);

// K_i feedback loop constant [dimensionless)
const Real PI_Ki = pin->GetReal("precipitator", "thermostat_Ki");
pkg->AddParam("PI_controller_Ki", PI_Ki);
pkg->AddParam("PI_controller_Ki", PI_Ki, parthenon::Params::Mutability::Restart);

/************************************************************
* Initialize the hydrostatic profile
************************************************************/
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);
hydro_pkg->AddParam<>("uniform_init", uniform_init, parthenon::Params::Mutability::Restart);
if (uniform_init == 1) {
const auto &uniform_init_height = pin->GetReal("precipitator", "uniform_init_height");
hydro_pkg->AddParam<>("uniform_init_height", uniform_init_height);
hydro_pkg->AddParam<>("uniform_init_height", uniform_init_height, parthenon::Params::Mutability::Restart);
}

const PrecipitatorProfile P_rho_profile(filename);
hydro_pkg->AddParam<>("precipitator_profile", P_rho_profile);
hydro_pkg->AddParam<>("precipitator_profile", P_rho_profile, parthenon::Params::Mutability::Restart);

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

// 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); // smoothing scale (code units)
h_smooth_heatcool, parthenon::Params::Mutability::Restart); // smoothing scale (code units)

// read perturbation parameters
const int kx_max = pin->GetInteger("precipitator", "perturb_kx");
Expand All @@ -689,11 +693,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);
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);
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);

// density perturbations
if (amp > 0.) {
Expand Down Expand Up @@ -765,12 +769,12 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *pkg
});

// save modes in hydro_pkg
hydro_pkg->AddParam("drho_hat", drho_hat);
hydro_pkg->AddParam("drho_hat", drho_hat, parthenon::Params::Mutability::Restart);
}

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

if (sigma_v > 0) {
auto k_peak_v = pin->GetReal("precipitator/driving", "k_peak");
Expand All @@ -787,7 +791,7 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *pkg

auto vertical_driving_only =
pin->GetOrAddBoolean("precipitator/driving", "vertical_driving_only", false);
hydro_pkg->AddParam("vertical_driving_only", vertical_driving_only);
hydro_pkg->AddParam("vertical_driving_only", vertical_driving_only, parthenon::Params::Mutability::Restart);

// when only v_z is desired, ensure that always k_z == 0
const bool xy_modes_only = vertical_driving_only;
Expand All @@ -796,7 +800,7 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *pkg

auto few_modes_ft = FewModesFT(pin, hydro_pkg, "precipitator_perturb_v", num_modes_v,
k_vec_v, k_peak_v, sol_weight_v, t_corr, rseed_v);
hydro_pkg->AddParam<>("precipitator/few_modes_ft_v", few_modes_ft);
hydro_pkg->AddParam<>("precipitator/few_modes_ft_v", few_modes_ft, parthenon::Params::Mutability::Restart);
}

if (parthenon::Globals::my_rank == 0) {
Expand Down
18 changes: 7 additions & 11 deletions src/utils/few_modes_ft.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -89,25 +89,24 @@ FewModesFT::FewModesFT(parthenon::ParameterInput *pin, parthenon::StateDescripto

// read var_hat
{
std::cout << "reading var_hat...\n";
//std::cout << "reading var_hat...\n";
auto var_hat_host = Kokkos::create_mirror_view(var_hat_);
for (int i = 0; i < 3; i++) {
for (int m = 0; m < num_modes; m++) {
auto real = pin->GetReal("few_modes_ft", "var_hat_" + std::to_string(i) + "_" +
std::to_string(m) + "_r");
auto imag = pin->GetReal("few_modes_ft", "var_hat_" + std::to_string(i) + "_" +
std::to_string(m) + "_i");
std::cout << "(" << i << "," << m << "): " << real << "\t" << imag << "\n";
//std::cout << "(" << i << "," << m << "): " << real << "\t" << imag << "\n";
var_hat_host(i, m) = Complex(real, imag);
}
}
Kokkos::deep_copy(var_hat_, var_hat_host);
}

#if 0
// read var_hat_new
{
std::cout << "reading var_hat_new...\n";
//std::cout << "reading var_hat_new...\n";
auto var_hat_new_host = Kokkos::create_mirror_view(var_hat_new_);
for (int i = 0; i < 3; i++) {
for (int m = 0; m < num_modes; m++) {
Expand All @@ -117,27 +116,26 @@ FewModesFT::FewModesFT(parthenon::ParameterInput *pin, parthenon::StateDescripto
auto imag_new =
pin->GetReal("few_modes_ft", "var_hat_new_" + std::to_string(i) + "_" +
std::to_string(m) + "_i");
std::cout << "(" << i << "," << m << "): " << real_new << "\t" << imag_new
<< "\n";
//std::cout << "(" << i << "," << m << "): " << real_new << "\t" << imag_new
// << "\n";
var_hat_new_host(i, m) = Complex(real_new, imag_new);
}
}
Kokkos::deep_copy(var_hat_new_, var_hat_new_host);
}
#endif

// read rng state
{
std::istringstream iss(pin->GetString("few_modes_ft", "state_rng"));
iss >> rng_;
std::cout << "rng state: " << rng_ << "\n";
//std::cout << "rng state: " << rng_ << "\n";
}

// read dist
{
std::istringstream iss(pin->GetString("few_modes_ft", "state_dist"));
iss >> dist_;
std::cout << "dist state: " << dist_ << "\n";
//std::cout << "dist state: " << dist_ << "\n";
}

} else {
Expand Down Expand Up @@ -175,7 +173,6 @@ void FewModesFT::SaveStateBeforeOutput(Mesh *mesh, ParameterInput *pin) {
}
}

#if 0
// var_hat_new_
{
auto var_hat_new_host = Kokkos::create_mirror_view_and_copy(parthenon::HostMemSpace(), var_hat_new_);
Expand All @@ -196,7 +193,6 @@ void FewModesFT::SaveStateBeforeOutput(Mesh *mesh, ParameterInput *pin) {
}
}
}
#endif

// save rng
{
Expand Down

0 comments on commit 3c4f7e1

Please sign in to comment.