Skip to content

Commit d640fa2

Browse files
FredJenningsAstroFredJenningsAstro
committed
Formatting
1 parent 2f65823 commit d640fa2

File tree

8 files changed

+4792
-4158
lines changed

8 files changed

+4792
-4158
lines changed

‎src/dust/dust.cpp‎

Lines changed: 1252 additions & 1116 deletions
Large diffs are not rendered by default.

‎src/dust/dust.hpp‎

Lines changed: 2796 additions & 2346 deletions
Large diffs are not rendered by default.

‎src/hydro/hydro.cpp‎

Lines changed: 83 additions & 91 deletions
Original file line numberDiff line numberDiff line change
@@ -38,7 +38,6 @@
3838
#include "srcterms/tabular_cooling.hpp"
3939
#include "utils/error_checking.hpp"
4040

41-
4241
// Dust headers
4342
#include "../dust/dust.hpp"
4443

@@ -239,11 +238,11 @@ TaskStatus AddUnsplitSources(MeshData<Real> *md, const SimTime &tm, const Real b
239238
if (enable_cooling == Cooling::tabular) {
240239
const TabularCooling &tabular_cooling =
241240
hydro_pkg->Param<TabularCooling>("tabular_cooling");
242-
const Real current_time = tm.time; // FJJ needed for writing dust history files from within subcycling
241+
const Real current_time =
242+
tm.time; // FJJ needed for writing dust history files from within subcycling
243243
tabular_cooling.SrcTerm(md, beta_dt, current_time);
244244
}
245245

246-
247246
if (ProblemSourceUnsplit != nullptr) {
248247
ProblemSourceUnsplit(md, tm, beta_dt);
249248
}
@@ -722,21 +721,17 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
722721
PARTHENON_FAIL("AthenaPK hydro: Unknown EOS");
723722
}
724723

725-
726724
/************************************************************
727725
* Read Dust 1 before tabular cooling to ensure correct scheme
728726
************************************************************/
729727
auto dust_on = pin->GetOrAddBoolean("dust", "active", false);
730728
pkg->AddParam<bool>("dust_on", dust_on);
731-
if(!dust_on){
729+
if (!dust_on) {
732730
pkg->AddParam<std::string>("dust_time_integrator", "none");
733731
pkg->AddParam<bool>("dust_subcycle_with_cooling", "false");
734732
pkg->AddParam<int>("dust_num_grains_sizes", 0);
735-
736733
}
737734

738-
739-
740735
/************************************************************
741736
* Read Tabular Cooling
742737
************************************************************/
@@ -797,100 +792,97 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
797792
prim_labels.emplace_back("scalar_" + std::to_string(i));
798793
running_scalar_idx += 1;
799794
}
800-
if(nscalars >=1){
801-
pkg->AddParam<int>("jet_scalar_idx", nhydro);
802-
} else{
795+
if (nscalars >= 1) {
796+
pkg->AddParam<int>("jet_scalar_idx", nhydro);
797+
} else {
803798
pkg->AddParam<int>("jet_scalar_idx", -1);
804799
}
805800

806-
807-
808-
809-
810801
/************************************************************
811802
* Read Dust 2
812803
************************************************************/
813804
dust::Dust dust(pin, pkg.get());
814-
auto disable_all_gas_cooling_for_testing = pin->GetOrAddBoolean("dust", "disable_all_gas_cooling_for_testing", false);
815-
pkg->AddParam<int>("disable_all_gas_cooling_for_testing", disable_all_gas_cooling_for_testing ? 1 : 0);
816-
817-
if(dust_on){
818-
819-
820-
821-
822-
auto num_grain_size_bins = pin->GetInteger("dust", "num_grainsize_bins");
823-
auto carbonaceous_grains = pin->GetOrAddBoolean("dust", "carbonaceous_grains", false);
824-
auto silicate_grains = pin->GetOrAddBoolean("dust", "silicate_grains", false);
825-
auto dust_subcycle_with_cooling = pin->GetOrAddBoolean("dust", "subcycle", true);
826-
auto dust_time_integrator = pin->GetString("dust", "time_integrator");
827-
PARTHENON_REQUIRE(dust_time_integrator == "euler" || dust_time_integrator == "heun", "Invalid dust time integrator!");
828-
auto max_frac_change_in_dust_bin = pin->GetOrAddReal("dust", "max_dM_in_bin", -1.);
829-
pkg->AddParam<bool>("dust_subcycle_with_cooling", dust_subcycle_with_cooling);
830-
pkg->AddParam<std::string>("dust_time_integrator", dust_time_integrator);
831-
std::vector<std::string> dust_var_names = {"density"};
832-
833-
int num_dust_vars = dust_var_names.size();
834-
int num_grain_compositions = (carbonaceous_grains + silicate_grains);
835-
int num_dust_bins = num_grain_compositions*num_grain_size_bins;
836-
int num_dust_scalars = 2 * num_dust_bins * num_dust_vars; // 2 since we must store both mass and density for the piecewise linear/loglinear schemes
837-
838-
pkg->AddParam<Real>("max_frac_change_in_dust_bin", max_frac_change_in_dust_bin);
839-
840-
// pkg->AddParam<int>("dust_num_grains_bins", num_dust_scalars);
841-
pkg->AddParam<int>("dust_num_grain_compositions", num_grain_compositions);
842-
pkg->AddParam<int>("dust_num_grains_sizes", num_grain_size_bins);
843-
844-
pkg->AddParam<int>("dust_scalar_idx_start", running_scalar_idx); // taking into account the jet_scalar
845-
pkg->AddParam<int>("dust_scalar_idx_end", running_scalar_idx+num_dust_scalars-1);
846-
running_scalar_idx += num_dust_scalars;
847-
848-
// Athena++ code paper: "**Ci (primitive variable) is the specific density of each
849-
// scalar and (ρCi) is the mass of each scalar species (conserved variable).**"
850-
// in ConstoPrim for tracer: prim(n, k, j, i) = cons(n, k, j, i) * di. with di = 1/cons(IDN, k, j, i);
851-
// So cons = prim * density
852-
853-
854-
855-
std::vector<std::string> dust_grain_compositions_names = {};
856-
// Order is carbonaceous_grains, silicate_grains
857-
if(carbonaceous_grains){
858-
dust_grain_compositions_names.emplace_back("carbonaceous");
859-
for (auto j = 0; j < num_grain_size_bins; j++) {
860-
cons_labels.emplace_back("dust_carbonaceous_grain_number_density_" + std::to_string(j));
861-
prim_labels.emplace_back("dust_carbonaceous_grain_number_density_" + std::to_string(j));
862-
cons_labels.emplace_back("dust_carbonaceous_grain_mass_density_" + std::to_string(j));
863-
prim_labels.emplace_back("dust_carbonaceous_grain_mass_density_" + std::to_string(j));
864-
}
865-
}
866-
if(silicate_grains){
867-
dust_grain_compositions_names.emplace_back("silicate");
868-
for (auto j = 0; j < num_grain_size_bins; j++) {
869-
cons_labels.emplace_back("dust_silicate_grain_number_density_" + std::to_string(j));
870-
prim_labels.emplace_back("dust_silicate_grain_number_density_" + std::to_string(j));
871-
cons_labels.emplace_back("dust_silicate_grain_mass_density_" + std::to_string(j));
872-
prim_labels.emplace_back("dust_silicate_grain_mass_density_" + std::to_string(j));
805+
auto disable_all_gas_cooling_for_testing =
806+
pin->GetOrAddBoolean("dust", "disable_all_gas_cooling_for_testing", false);
807+
pkg->AddParam<int>("disable_all_gas_cooling_for_testing",
808+
disable_all_gas_cooling_for_testing ? 1 : 0);
809+
810+
if (dust_on) {
811+
812+
auto num_grain_size_bins = pin->GetInteger("dust", "num_grainsize_bins");
813+
auto carbonaceous_grains = pin->GetOrAddBoolean("dust", "carbonaceous_grains", false);
814+
auto silicate_grains = pin->GetOrAddBoolean("dust", "silicate_grains", false);
815+
auto dust_subcycle_with_cooling = pin->GetOrAddBoolean("dust", "subcycle", true);
816+
auto dust_time_integrator = pin->GetString("dust", "time_integrator");
817+
PARTHENON_REQUIRE(dust_time_integrator == "euler" || dust_time_integrator == "heun",
818+
"Invalid dust time integrator!");
819+
auto max_frac_change_in_dust_bin = pin->GetOrAddReal("dust", "max_dM_in_bin", -1.);
820+
pkg->AddParam<bool>("dust_subcycle_with_cooling", dust_subcycle_with_cooling);
821+
pkg->AddParam<std::string>("dust_time_integrator", dust_time_integrator);
822+
std::vector<std::string> dust_var_names = {"density"};
823+
824+
int num_dust_vars = dust_var_names.size();
825+
int num_grain_compositions = (carbonaceous_grains + silicate_grains);
826+
int num_dust_bins = num_grain_compositions * num_grain_size_bins;
827+
int num_dust_scalars =
828+
2 * num_dust_bins * num_dust_vars; // 2 since we must store both mass and density
829+
// for the piecewise linear/loglinear schemes
830+
831+
pkg->AddParam<Real>("max_frac_change_in_dust_bin", max_frac_change_in_dust_bin);
832+
833+
// pkg->AddParam<int>("dust_num_grains_bins", num_dust_scalars);
834+
pkg->AddParam<int>("dust_num_grain_compositions", num_grain_compositions);
835+
pkg->AddParam<int>("dust_num_grains_sizes", num_grain_size_bins);
836+
837+
pkg->AddParam<int>("dust_scalar_idx_start",
838+
running_scalar_idx); // taking into account the jet_scalar
839+
pkg->AddParam<int>("dust_scalar_idx_end", running_scalar_idx + num_dust_scalars - 1);
840+
running_scalar_idx += num_dust_scalars;
841+
842+
// Athena++ code paper: "**Ci (primitive variable) is the specific density of each
843+
// scalar and (ρCi) is the mass of each scalar species (conserved variable).**"
844+
// in ConstoPrim for tracer: prim(n, k, j, i) = cons(n, k, j, i) * di. with di =
845+
// 1/cons(IDN, k, j, i); So cons = prim * density
846+
847+
std::vector<std::string> dust_grain_compositions_names = {};
848+
// Order is carbonaceous_grains, silicate_grains
849+
if (carbonaceous_grains) {
850+
dust_grain_compositions_names.emplace_back("carbonaceous");
851+
for (auto j = 0; j < num_grain_size_bins; j++) {
852+
cons_labels.emplace_back("dust_carbonaceous_grain_number_density_" +
853+
std::to_string(j));
854+
prim_labels.emplace_back("dust_carbonaceous_grain_number_density_" +
855+
std::to_string(j));
856+
cons_labels.emplace_back("dust_carbonaceous_grain_mass_density_" +
857+
std::to_string(j));
858+
prim_labels.emplace_back("dust_carbonaceous_grain_mass_density_" +
859+
std::to_string(j));
860+
}
861+
}
862+
if (silicate_grains) {
863+
dust_grain_compositions_names.emplace_back("silicate");
864+
for (auto j = 0; j < num_grain_size_bins; j++) {
865+
cons_labels.emplace_back("dust_silicate_grain_number_density_" +
866+
std::to_string(j));
867+
prim_labels.emplace_back("dust_silicate_grain_number_density_" +
868+
std::to_string(j));
869+
cons_labels.emplace_back("dust_silicate_grain_mass_density_" + std::to_string(j));
870+
prim_labels.emplace_back("dust_silicate_grain_mass_density_" + std::to_string(j));
871+
}
873872
}
874-
}
875-
876-
// Use ints for easy capture into Kokkos lambdas
877-
pkg->AddParam("dust_carbonaceous_grains", carbonaceous_grains ? 1 : 0 );
878-
pkg->AddParam("dust_silicate_grains", silicate_grains ? 1 : 0 );
879-
pkg->AddParam("dust_grain_compositions_names", dust_grain_compositions_names);
880-
881-
// std::cout << "num_dust_scalars" << num_dust_scalars << std::endl;
882-
} // if dust_on
883-
884-
885-
886-
887-
888873

874+
// Use ints for easy capture into Kokkos lambdas
875+
pkg->AddParam("dust_carbonaceous_grains", carbonaceous_grains ? 1 : 0);
876+
pkg->AddParam("dust_silicate_grains", silicate_grains ? 1 : 0);
877+
pkg->AddParam("dust_grain_compositions_names", dust_grain_compositions_names);
889878

890-
pkg->AddParam("user_nscalars", nscalars); // the user-input nscalars (e.g. without dust...)
891-
nscalars = running_scalar_idx - nhydro;
892-
pkg->AddParam("nscalars", nscalars);
879+
// std::cout << "num_dust_scalars" << num_dust_scalars << std::endl;
880+
} // if dust_on
893881

882+
pkg->AddParam("user_nscalars",
883+
nscalars); // the user-input nscalars (e.g. without dust...)
884+
nscalars = running_scalar_idx - nhydro;
885+
pkg->AddParam("nscalars", nscalars);
894886

895887
Metadata m(
896888
{Metadata::Cell, Metadata::Independent, Metadata::FillGhost, Metadata::WithFluxes},

0 commit comments

Comments
 (0)