Skip to content

Commit cadae5d

Browse files
FredJenningsAstroFredJenningsAstro
committed
Apply formatting
1 parent ab22991 commit cadae5d

File tree

5 files changed

+1097
-975
lines changed

5 files changed

+1097
-975
lines changed

‎src/dust/dust.cpp‎

Lines changed: 133 additions & 121 deletions
Original file line numberDiff line numberDiff line change
@@ -268,19 +268,31 @@ Dust::Dust(parthenon::ParameterInput *pin, parthenon::StateDescriptor *hydro_pkg
268268
pin->GetReal("dust", "init_run_stellar_injection_time");
269269
}
270270

271-
if(dust_cooling_mode_str_ == "Dwek_Werner1981" || dust_cooling_mode_str_ == "Dwek_Werner1981_INTEGRATED"){
272-
if(dust_cooling_mode_str_ == "Dwek_Werner1981"){dust_cooling_mode_ = DustCoolingMode::DWEKWERNER1981;}
273-
if(dust_cooling_mode_str_ == "Dwek_Werner1981_INTEGRATED"){dust_cooling_mode_ = DustCoolingMode::DWEKWERNER1981_INTEGRATED;}
274-
// precompute these coefficients here to prevent extra computation in the kernel
275-
dwek_werner_coeff_a_code_units_ = 5.38 * Kokkos::pow(10, -18) * (cm3_to_code_vol_ * erg_to_code_energy_/seconds_to_code_time_);
276-
dwek_werner_coeff_b_code_units_ = 3.37 * Kokkos::pow(10, -13) * (cm3_to_code_vol_ * erg_to_code_energy_/seconds_to_code_time_);
277-
dwek_werner_coeff_c_code_units_ = 6.48 * Kokkos::pow(10, -6) * (cm3_to_code_vol_ * erg_to_code_energy_/seconds_to_code_time_);
278-
dwek_werner_regime_coeff_ = 2.71 * Kokkos::pow(10, 8);
279-
} else if (dust_cooling_mode_str_ == "off"){
280-
printf("Setting dust_cooling_mode_ = DustCoolingMode::OFF; \n");
281-
dust_cooling_mode_ = DustCoolingMode::OFF;
282-
} else {PARTHENON_FAIL("Invalid dust cooling specified")}
283-
271+
if (dust_cooling_mode_str_ == "Dwek_Werner1981" ||
272+
dust_cooling_mode_str_ == "Dwek_Werner1981_INTEGRATED") {
273+
if (dust_cooling_mode_str_ == "Dwek_Werner1981") {
274+
dust_cooling_mode_ = DustCoolingMode::DWEKWERNER1981;
275+
}
276+
if (dust_cooling_mode_str_ == "Dwek_Werner1981_INTEGRATED") {
277+
dust_cooling_mode_ = DustCoolingMode::DWEKWERNER1981_INTEGRATED;
278+
}
279+
// precompute these coefficients here to prevent extra computation in the kernel
280+
dwek_werner_coeff_a_code_units_ =
281+
5.38 * Kokkos::pow(10, -18) *
282+
(cm3_to_code_vol_ * erg_to_code_energy_ / seconds_to_code_time_);
283+
dwek_werner_coeff_b_code_units_ =
284+
3.37 * Kokkos::pow(10, -13) *
285+
(cm3_to_code_vol_ * erg_to_code_energy_ / seconds_to_code_time_);
286+
dwek_werner_coeff_c_code_units_ =
287+
6.48 * Kokkos::pow(10, -6) *
288+
(cm3_to_code_vol_ * erg_to_code_energy_ / seconds_to_code_time_);
289+
dwek_werner_regime_coeff_ = 2.71 * Kokkos::pow(10, 8);
290+
} else if (dust_cooling_mode_str_ == "off") {
291+
printf("Setting dust_cooling_mode_ = DustCoolingMode::OFF; \n");
292+
dust_cooling_mode_ = DustCoolingMode::OFF;
293+
} else {
294+
PARTHENON_FAIL("Invalid dust cooling specified")
295+
}
284296

285297
// Create device views of the required vectors for e.g. the cooling functions
286298
grain_midbin_sizes_microm_ =
@@ -562,74 +574,75 @@ Dust::Dust(parthenon::ParameterInput *pin, parthenon::StateDescriptor *hydro_pkg
562574
if (!file) {
563575
std::cerr << "Error: Could not open the file!" << std::endl;
564576
PARTHENON_FAIL("Error: Could not open the file!");
565-
}
566-
// Write each element on a new line (column format)
567-
for (int i = 0; i < num_grainsize_bins_; ++i) {
568-
file << host_grainsize_bin_edges_microm[i] << " - " << host_grainsize_bin_edges_microm[i+1] << ": " << host_agb_normalised_silicate_number_distribution_array[i] << std::endl;
569-
}
570-
// Close the file
571-
file.close();
572577
}
578+
// Write each element on a new line (column format)
579+
for (int i = 0; i < num_grainsize_bins_; ++i) {
580+
file << host_grainsize_bin_edges_microm[i] << " - "
581+
<< host_grainsize_bin_edges_microm[i + 1] << ": "
582+
<< host_agb_normalised_silicate_number_distribution_array[i] << std::endl;
573583
}
574-
} // if(agb_winds_)
575-
576-
577-
578-
579-
580-
// // Precompute cooling tables. Only do on one rank since we will write to file. just run on host
581-
582-
// if (parthenon::Globals::my_rank == 0){
583-
// std::string column_name;
584-
// std::ostringstream oss;
585-
// // oss << std::setw(3) << std::setfill('0') << dust_i;
586-
// std::string folder_path;
587-
// folder_path = "./dust_cool_table/";
588-
// oss.str(""); // clear the string buffer
589-
// oss.clear(); // reset error/EOF flags
590-
// // Check if folder exists, if not create it
591-
// std::filesystem::create_directories(folder_path);
592-
// std::ofstream cool_file;
593-
// std::string file_tag = "dust_cooling_table.txt";
594-
// cool_file.open(folder_path + file_tag, std::ofstream::app);
595-
// if( !cool_file.is_open() )
596-
// std::cerr << "Error: Unable to open file 'dust_cooling_table.txt' for writing." << std::endl;
597-
598-
// for(int gs_i = 0; gs_i < num_grainsize_bins_; gs_i++){
599-
// if(gs_i == 0){
600-
// cool_file << "T (K)" << "|" ;
601-
// }
602-
// cool_file << std::setprecision(5) << host_grain_midbin_sizes_microm[gs_i] << "|" << "\n";
603-
// }
604-
605-
// printf("num_grainsize_bins_=%d \n",num_grainsize_bins_);
606-
// for(int temp = 0.; temp <= 9.; temp+= 0.1){
607-
// for(int gs_i = 0; gs_i < num_grainsize_bins_; gs_i++){
608-
// const Real dust_de_dt_this_grain_bin = dust::PreComputeDwekWernerGrainCooling(
609-
// Kokkos::pow(10., temp),
610-
// gs_i,
611-
// dwek_werner_regime_coeff_,
612-
// dwek_werner_coeff_a_code_units_,
613-
// dwek_werner_coeff_b_code_units_,
614-
// dwek_werner_coeff_c_code_units_,
615-
// host_grain_midbin_sizes_microm
616-
// );
617-
// printf("Kokkos::pow(10., temp)=%e gs_i = %d dust_de_dt_this_grain_bin=%e \n", Kokkos::pow(10., temp), gs_i, dust_de_dt_this_grain_bin);
618-
// if(temp == 0){
619-
// cool_file << std::setprecision(5) << temp << " " ;
620-
// }
621-
// cool_file << std::setprecision(5) << Kokkos::log10(Kokkos::abs(dust_de_dt_this_grain_bin)) << " " ;
622-
// }
623-
// cool_file << std::endl;
624-
// }
625-
// cool_file.close();
626-
627-
// }
628-
629-
630-
631-
hydro_pkg->AddParam<>("dust", *this); // So we can access the instance everywhere we can access hydro_pkg
632-
} //Dust::Dust
584+
// Close the file
585+
file.close();
586+
}
587+
}
588+
} // if(agb_winds_)
589+
590+
// // Precompute cooling tables. Only do on one rank since we will write to file. just
591+
// run on host
592+
593+
// if (parthenon::Globals::my_rank == 0){
594+
// std::string column_name;
595+
// std::ostringstream oss;
596+
// // oss << std::setw(3) << std::setfill('0') << dust_i;
597+
// std::string folder_path;
598+
// folder_path = "./dust_cool_table/";
599+
// oss.str(""); // clear the string buffer
600+
// oss.clear(); // reset error/EOF flags
601+
// // Check if folder exists, if not create it
602+
// std::filesystem::create_directories(folder_path);
603+
// std::ofstream cool_file;
604+
// std::string file_tag = "dust_cooling_table.txt";
605+
// cool_file.open(folder_path + file_tag, std::ofstream::app);
606+
// if( !cool_file.is_open() )
607+
// std::cerr << "Error: Unable to open file 'dust_cooling_table.txt' for writing."
608+
// << std::endl;
609+
610+
// for(int gs_i = 0; gs_i < num_grainsize_bins_; gs_i++){
611+
// if(gs_i == 0){
612+
// cool_file << "T (K)" << "|" ;
613+
// }
614+
// cool_file << std::setprecision(5) << host_grain_midbin_sizes_microm[gs_i] << "|"
615+
// << "\n";
616+
// }
617+
618+
// printf("num_grainsize_bins_=%d \n",num_grainsize_bins_);
619+
// for(int temp = 0.; temp <= 9.; temp+= 0.1){
620+
// for(int gs_i = 0; gs_i < num_grainsize_bins_; gs_i++){
621+
// const Real dust_de_dt_this_grain_bin = dust::PreComputeDwekWernerGrainCooling(
622+
// Kokkos::pow(10., temp),
623+
// gs_i,
624+
// dwek_werner_regime_coeff_,
625+
// dwek_werner_coeff_a_code_units_,
626+
// dwek_werner_coeff_b_code_units_,
627+
// dwek_werner_coeff_c_code_units_,
628+
// host_grain_midbin_sizes_microm
629+
// );
630+
// printf("Kokkos::pow(10., temp)=%e gs_i = %d dust_de_dt_this_grain_bin=%e \n",
631+
// Kokkos::pow(10., temp), gs_i, dust_de_dt_this_grain_bin); if(temp == 0){
632+
// cool_file << std::setprecision(5) << temp << " " ;
633+
// }
634+
// cool_file << std::setprecision(5) <<
635+
// Kokkos::log10(Kokkos::abs(dust_de_dt_this_grain_bin)) << " " ;
636+
// }
637+
// cool_file << std::endl;
638+
// }
639+
// cool_file.close();
640+
641+
// }
642+
643+
hydro_pkg->AddParam<>(
644+
"dust", *this); // So we can access the instance everywhere we can access hydro_pkg
645+
} // Dust::Dust
633646

634647
// Function to measure dust masses and cooling rates in radii and gas T bins and write to
635648
// file
@@ -689,25 +702,25 @@ void Dust::MeasureAndRecordHistory(parthenon::MeshData<parthenon::Real> *md,
689702
Real cm3_to_code_vol_ = units.cm() * units.cm() * units.cm();
690703
Real msun_to_code_mass = units.msun();
691704

692-
int we_have_dust_cooling;
693-
auto dust_cooling_mode_ = this->dust_cooling_mode_;
694-
// std::optional<Real> nH_to_ne;
695-
if (hydro_pkg->Param<bool>("dust_on")){
696-
we_have_dust_cooling = 1;
697-
// nH_to_ne = hydro_pkg->Param<Real>("nH_to_ne");
698-
699-
switch(dust_cooling_mode_) {
700-
case dust::DustCoolingMode::OFF:
701-
we_have_dust_cooling = 0;
702-
break;
703-
case dust::DustCoolingMode::DWEKWERNER1981:
704-
break;
705-
case dust::DustCoolingMode::DWEKWERNER1981_INTEGRATED:
706-
break;
707-
}
708-
}else{
709-
we_have_dust_cooling = 0;
710-
}
705+
int we_have_dust_cooling;
706+
auto dust_cooling_mode_ = this->dust_cooling_mode_;
707+
// std::optional<Real> nH_to_ne;
708+
if (hydro_pkg->Param<bool>("dust_on")) {
709+
we_have_dust_cooling = 1;
710+
// nH_to_ne = hydro_pkg->Param<Real>("nH_to_ne");
711+
712+
switch (dust_cooling_mode_) {
713+
case dust::DustCoolingMode::OFF:
714+
we_have_dust_cooling = 0;
715+
break;
716+
case dust::DustCoolingMode::DWEKWERNER1981:
717+
break;
718+
case dust::DustCoolingMode::DWEKWERNER1981_INTEGRATED:
719+
break;
720+
}
721+
} else {
722+
we_have_dust_cooling = 0;
723+
}
711724

712725
int dust_scalar_idx_start = hydro_pkg->Param<int>("dust_scalar_idx_start");
713726
int dust_scalar_idx_end = hydro_pkg->Param<int>("dust_scalar_idx_end");
@@ -762,25 +775,24 @@ void Dust::MeasureAndRecordHistory(parthenon::MeshData<parthenon::Real> *md,
762775
const int disable_all_gas_cooling_for_testing =
763776
hydro_pkg->Param<int>("disable_all_gas_cooling_for_testing");
764777

765-
int dust_piecewise_mode_int = 0;
766-
if(this->piecewise_mode_ == dust::DustPiecewiseMode::LINEAR){
767-
dust_piecewise_mode_int = 1;
768-
} else if(this->piecewise_mode_ == dust::DustPiecewiseMode::LOGLINEAR){
769-
dust_piecewise_mode_int = 2;
770-
}
771-
772-
773-
Kokkos::parallel_for(
774-
"DustHst",
775-
Kokkos::MDRangePolicy<Kokkos::Rank<5>>(
776-
parthenon::DevExecSpace(), {0, 0, kb.s, jb.s, ib.s},
777-
{cons_pack.GetDim(5), num_dust_bins, kb.e + 1, jb.e + 1, ib.e + 1}),
778-
// {1, 1, 1, 1, ib.e + 1 - ib.s}),
779-
KOKKOS_LAMBDA(const int b, const int dust_i, const int k, const int j, const int i) {
780-
// dust_i runs from 0 to number of dust types * number of size bins. It will be
781-
// half the size of the total number of dust variables in the cons pack, which
782-
// stores both mass and number density for each bin
778+
int dust_piecewise_mode_int = 0;
779+
if (this->piecewise_mode_ == dust::DustPiecewiseMode::LINEAR) {
780+
dust_piecewise_mode_int = 1;
781+
} else if (this->piecewise_mode_ == dust::DustPiecewiseMode::LOGLINEAR) {
782+
dust_piecewise_mode_int = 2;
783+
}
783784

785+
Kokkos::parallel_for(
786+
"DustHst",
787+
Kokkos::MDRangePolicy<Kokkos::Rank<5>>(
788+
parthenon::DevExecSpace(), {0, 0, kb.s, jb.s, ib.s},
789+
{cons_pack.GetDim(5), num_dust_bins, kb.e + 1, jb.e + 1, ib.e + 1}),
790+
// {1, 1, 1, 1, ib.e + 1 - ib.s}),
791+
KOKKOS_LAMBDA(const int b, const int dust_i, const int k, const int j,
792+
const int i) {
793+
// dust_i runs from 0 to number of dust types * number of size bins. It will be
794+
// half the size of the total number of dust variables in the cons pack, which
795+
// stores both mass and number density for each bin
784796

785797
auto f_a_gaseous_cool_rate = scatter_f_gaseous_cool_rate.access();
786798
auto f_a_dust_cool_rate = scatter_f_dust_cool_rate.access();
@@ -1108,11 +1120,11 @@ void Dust::MeasureAndRecordHistory(parthenon::MeshData<parthenon::Real> *md,
11081120
}
11091121
}
11101122
}
1111-
} else{
1112-
// Currently can only do dust history with Tabular Cooling
1113-
;
1114-
}
1115-
}; // Dust::MeasureAndRecordHistory
1123+
} else {
1124+
// Currently can only do dust history with Tabular Cooling
1125+
;
1126+
}
1127+
}; // Dust::MeasureAndRecordHistory
11161128

11171129
std::vector<double> Dust::get_r_bin_edges() const { return this->r_bin_edges_; };
11181130

0 commit comments

Comments
 (0)