Skip to content

Commit a30110d

Browse files
FredJenningsAstroFredJenningsAstro
committed
Minor changes and checks
1 parent 46a6258 commit a30110d

File tree

4 files changed

+98
-129
lines changed

4 files changed

+98
-129
lines changed

‎docs/dust.md‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -15,7 +15,7 @@ active = [true, false]
1515
subcycle = [true, false] -- default = true <br>
1616
time_integrator = [heun, euler] <br>
1717
max_dM_in_bin = 0.1 -- default = -1. <br>
18-
piecewise_method = [loglinear, linear] <br>
18+
piecewise_method = [hybrid_loglinear,loglinear, linear] <br>
1919
cooling = [Dwek_Werner1981_INTEGRATED, Dwek_Werner1981, off] <br>
2020
dust_cool_table_N_Tbins = [-1, or +/ve int] <br>
2121
disable_all_gas_cooling_for_testing = [true, false] -- default = false <br>

‎src/dust/dust.cpp‎

Lines changed: 31 additions & 41 deletions
Original file line numberDiff line numberDiff line change
@@ -268,30 +268,19 @@ 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" ||
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-
dust_cooling_mode_ = DustCoolingMode::OFF;
292-
} else {
293-
PARTHENON_FAIL("Invalid dust cooling specified")
294-
}
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+
295284

296285
// Create device views of the required vectors for e.g. the cooling functions
297286
grain_midbin_sizes_microm_ =
@@ -719,24 +708,25 @@ void Dust::MeasureAndRecordHistory(parthenon::MeshData<parthenon::Real> *md,
719708
const int disable_all_gas_cooling_for_testing =
720709
hydro_pkg->Param<int>("disable_all_gas_cooling_for_testing");
721710

722-
int dust_piecewise_mode_int = 0;
723-
if (this->piecewise_mode_ == dust::DustPiecewiseMode::LINEAR) {
724-
dust_piecewise_mode_int = 1;
725-
} else if (this->piecewise_mode_ == dust::DustPiecewiseMode::LOGLINEAR) {
726-
dust_piecewise_mode_int = 2;
727-
}
711+
int dust_piecewise_mode_int = 0;
712+
if(this->piecewise_mode_ == dust::DustPiecewiseMode::LINEAR){
713+
dust_piecewise_mode_int = 1;
714+
} else if(this->piecewise_mode_ == dust::DustPiecewiseMode::LOGLINEAR){
715+
dust_piecewise_mode_int = 2;
716+
}
717+
718+
719+
Kokkos::parallel_for(
720+
"DustHst",
721+
Kokkos::MDRangePolicy<Kokkos::Rank<5>>(
722+
parthenon::DevExecSpace(), {0, 0, kb.s, jb.s, ib.s},
723+
{cons_pack.GetDim(5), num_dust_bins, kb.e + 1, jb.e + 1, ib.e + 1}),
724+
// {1, 1, 1, 1, ib.e + 1 - ib.s}),
725+
KOKKOS_LAMBDA(const int b, const int dust_i, const int k, const int j, const int i) {
726+
// dust_i runs from 0 to number of dust types * number of size bins. It will be
727+
// half the size of the total number of dust variables in the cons pack, which
728+
// stores both mass and number density for each bin
728729

729-
Kokkos::parallel_for(
730-
"DustHst",
731-
Kokkos::MDRangePolicy<Kokkos::Rank<5>>(
732-
parthenon::DevExecSpace(), {0, 0, kb.s, jb.s, ib.s},
733-
{cons_pack.GetDim(5), num_dust_bins, kb.e + 1, jb.e + 1, ib.e + 1}),
734-
// {1, 1, 1, 1, ib.e + 1 - ib.s}),
735-
KOKKOS_LAMBDA(const int b, const int dust_i, const int k, const int j,
736-
const int i) {
737-
// dust_i runs from 0 to number of dust types * number of size bins. It will be
738-
// half the size of the total number of dust variables in the cons pack, which
739-
// stores both mass and number density for each bin
740730

741731
auto f_a_gaseous_cool_rate = scatter_f_gaseous_cool_rate.access();
742732
auto f_a_dust_cool_rate = scatter_f_dust_cool_rate.access();

‎src/dust/dust.hpp‎

Lines changed: 55 additions & 70 deletions
Original file line numberDiff line numberDiff line change
@@ -829,26 +829,27 @@ struct DustDevice {
829829
int only_sputtering_for_debug;
830830
Real whole_box_extent;
831831

832-
// AGB vars
833-
Real agb_max_radius;
834-
Real gamma_star;
835-
Real sersic_n;
836-
Real sersic_Re;
837-
Real stellar_profile_norm;
838-
Real stellar_mass_cent;
839-
Real stellar_density_profile_r_low;
840-
Real stellar_density_profile_r_up;
841-
Real dust_return_silicates_mass_fraction_per_megayear;
842-
Real dust_return_carbon_mass_fraction_per_megayear;
843-
Real code_to_megayear;
844-
int carbonaceous_grains;
845-
int silicate_grains;
846-
ParArray1D<Real> agb_normalised_carbonaceous_mass_distribution_array;
847-
ParArray1D<Real> agb_normalised_carbonaceous_number_distribution_array;
848-
ParArray1D<Real> agb_normalised_silicate_mass_distribution_array;
849-
ParArray1D<Real> agb_normalised_silicate_number_distribution_array;
850-
StellarRadialProfile stellar_radial_profile;
851-
// fjjcurrent
832+
// AGB vars
833+
Real agb_max_radius;
834+
Real gamma_star;
835+
Real sersic_n;
836+
Real sersic_Re;
837+
Real stellar_profile_norm;
838+
Real stellar_mass_cent;
839+
Real stellar_density_profile_r_low;
840+
Real stellar_density_profile_r_up;
841+
Real dust_return_silicates_mass_fraction_per_megayear;
842+
Real dust_return_carbon_mass_fraction_per_megayear;
843+
Real code_to_megayear;
844+
int carbonaceous_grains;
845+
int silicate_grains;
846+
ParArray1D<Real> agb_normalised_carbonaceous_mass_distribution_array;
847+
ParArray1D<Real> agb_normalised_carbonaceous_number_distribution_array;
848+
ParArray1D<Real> agb_normalised_silicate_mass_distribution_array;
849+
ParArray1D<Real> agb_normalised_silicate_number_distribution_array;
850+
StellarRadialProfile stellar_radial_profile;
851+
852+
//fjjcurrent
852853

853854
// Helper function for calculating the cooling rates for the Linear reconstruction
854855
// method
@@ -1207,7 +1208,6 @@ struct DustDevice {
12071208

12081209

12091210

1210-
12111211

12121212
// top-level function to call for DW dust cooling. Can decide here to do for single size bin or for all bins
12131213
KOKKOS_INLINE_FUNCTION
@@ -1237,6 +1237,9 @@ struct DustDevice {
12371237
integrated_rates
12381238
);
12391239
} else { // over all dust bins
1240+
1241+
// printf("A single_grain_densities.extent(0)=%d grain_midbin_sizes_microm.extent(0)=%d \n ", single_grain_densities.extent(0), grain_midbin_sizes_microm.extent(0));
1242+
KOKKOS_ASSERT(single_grain_densities.extent(0)*grain_midbin_sizes_microm.extent(0)>=1);
12401243
for(int gc_i = 0; gc_i < single_grain_densities.extent(0); gc_i++){ // loop over grain compositions
12411244
for(int gs_i = 0; gs_i < grain_midbin_sizes_microm.extent(0); gs_i++){ // loop over grain sizes
12421245
//get correct index into cons_pack subview for this grain type and size bin
@@ -2567,55 +2570,37 @@ void DustFilladotView(const Real temperature, const int gc_i, const int gs_i, co
25672570
DustCalculateAdotPerBin(temperature, rho, DustDevObj, adot_sputter, adot_accretion,
25682571
adot);
25692572

2570-
const Real whole_box_extent = DustDevObj.whole_box_extent;
2571-
const auto coords = cons_pack.GetCoords(b);
2572-
const auto x = coords.Xc<1>(i);
2573-
const auto y = coords.Xc<2>(j);
2574-
const auto z = coords.Xc<3>(k);
2575-
const auto r = Kokkos::sqrt(x * x + y * y + z * z);
2576-
// Check correct signs. Don;t worry too much if very near a boundary, where densities
2577-
// might go weird
2578-
if (adot_sputter > 0 || adot_sputter != adot_sputter) {
2579-
if (std::abs(x) < 0.9 * whole_box_extent && std::abs(y) < 0.9 * whole_box_extent &&
2580-
std::abs(z) < 0.9 * whole_box_extent) {
2581-
// printf("[FJJ DEBUG] Sputtering is growing grains! x=%e y =%e z=%e r=%e
2582-
// whole_box_extent=%e f_sput=%e rho =%e adot_sputter=%e sput_prefac=%e
2583-
// sput_dens=%e sput_T=%e \n", x,y,z,r, whole_box_extent, f_sput, rho,
2584-
// adot_sputter, sput_prefac, sput_dens, sput_T);
2585-
printf("[FJJ DEBUG] Sputtering is growing grains! x=%e y =%e z=%e "
2586-
"whole_box_extent=%e rho =%e\n",
2587-
x, y, z, whole_box_extent, rho);
2588-
}
2589-
if (std::abs(x) < 0.9 * whole_box_extent && std::abs(y) < 0.9 * whole_box_extent &&
2590-
std::abs(z) < 0.9 * whole_box_extent) { // ignore weird things at box boundary
2591-
// e.g. negative densities
2592-
printf("[FJJ DEBUG] Sputtering is growing grains inside of the boundary! x=%e y "
2593-
"=%e z=%e whole_box_extent=%e rho =%e temperature=%e\n",
2594-
x, y, z, whole_box_extent, rho, temperature);
2595-
PARTHENON_REQUIRE(adot_sputter <= 0, "Sputtering is growing grains!");
2596-
PARTHENON_REQUIRE(adot_sputter == adot_sputter, "adot_sputter is nan!");
2597-
}
2598-
adot_sputter = 0.;
2599-
adot_accretion = 0; // don;t do any dust updates if sputtering already is bad
2600-
adot = 0.;
2601-
}
2602-
if (adot_accretion < 0 || adot_accretion != adot_accretion) {
2603-
if (std::abs(x) < 0.9 * whole_box_extent && std::abs(y) < 0.9 * whole_box_extent &&
2604-
std::abs(z) < 0.9 * whole_box_extent) {
2605-
printf("[FJJ DEBUG] Accretion is shrinking grains! x=%e y =%e z=%e r=%e "
2606-
"whole_box_extent=%e rho =%e adot_accretion=%e \n",
2607-
x, y, z, r, whole_box_extent, rho, adot_accretion);
2608-
}
2609-
2610-
if (std::abs(x) < 0.9 * whole_box_extent && std::abs(y) < 0.9 * whole_box_extent &&
2611-
std::abs(z) < 0.9 * whole_box_extent) { // ignore weird things at box boundary
2612-
// e.g. negative densities
2613-
printf("[FJJ DEBUG] Accretion is shrinking grains inside of the boundary! x=%e y "
2614-
"=%e z=%e r=%e whole_box_extent=%e rho =%e adot_accretion=%e \n",
2615-
x, y, z, r, whole_box_extent, rho, adot_accretion);
2616-
PARTHENON_REQUIRE(adot_accretion >= 0, "Accretion is shrinking grains!");
2617-
PARTHENON_REQUIRE(adot_accretion == adot_accretion, "adot_sputter is nan!");
2618-
}
2573+
const Real whole_box_extent = DustDevObj.whole_box_extent;
2574+
const auto coords = cons_pack.GetCoords(b);
2575+
const auto x = coords.Xc<1>(i);
2576+
const auto y = coords.Xc<2>(j);
2577+
const auto z = coords.Xc<3>(k);
2578+
const auto r = Kokkos::sqrt(x * x + y * y + z * z);
2579+
// Check correct signs. Don;t worry too much if very near a boundary, where densities might go weird
2580+
if(adot_sputter > 0 || adot_sputter != adot_sputter){
2581+
if(std::abs(x) < 0.9*whole_box_extent && std::abs(y) < 0.9*whole_box_extent && std::abs(z) < 0.9*whole_box_extent){
2582+
// printf("[FJJ DEBUG] Sputtering is growing grains! x=%e y =%e z=%e r=%e whole_box_extent=%e f_sput=%e rho =%e adot_sputter=%e sput_prefac=%e sput_dens=%e sput_T=%e \n", x,y,z,r, whole_box_extent, f_sput, rho, adot_sputter, sput_prefac, sput_dens, sput_T);
2583+
// printf("[FJJ DEBUG] Sputtering is growing grains! x=%e y =%e z=%e whole_box_extent=%e rho =%e\n", x,y,z, whole_box_extent, rho);
2584+
}
2585+
if(std::abs(x) < 0.9*whole_box_extent && std::abs(y) < 0.9*whole_box_extent && std::abs(z) < 0.9*whole_box_extent){ // ignore weird things at box boundary e.g. negative densities
2586+
// printf("[FJJ DEBUG] Sputtering is growing grains inside of the boundary! x=%e y =%e z=%e whole_box_extent=%e rho =%e temperature=%e\n", x,y,z, whole_box_extent, rho, temperature);
2587+
// PARTHENON_REQUIRE(adot_sputter <= 0 , "Sputtering is growing grains!");
2588+
// PARTHENON_REQUIRE(adot_sputter == adot_sputter , "adot_sputter is nan!");
2589+
}
2590+
adot_sputter = 0.;
2591+
adot_accretion = 0; // don;t do any dust updates if sputtering already is bad
2592+
adot = 0.;
2593+
}
2594+
if(adot_accretion < 0 || adot_accretion != adot_accretion){
2595+
if(std::abs(x) < 0.9*whole_box_extent && std::abs(y) < 0.9*whole_box_extent && std::abs(z) < 0.9*whole_box_extent){
2596+
// printf("[FJJ DEBUG] Accretion is shrinking grains! x=%e y =%e z=%e r=%e whole_box_extent=%e rho =%e adot_accretion=%e \n", x,y,z,r, whole_box_extent, rho, adot_accretion);
2597+
}
2598+
2599+
if(std::abs(x) < 0.9*whole_box_extent && std::abs(y) < 0.9*whole_box_extent && std::abs(z) < 0.9*whole_box_extent){ // ignore weird things at box boundary e.g. negative densities
2600+
// printf("[FJJ DEBUG] Accretion is shrinking grains inside of the boundary! x=%e y =%e z=%e r=%e whole_box_extent=%e rho =%e adot_accretion=%e \n", x,y,z,r, whole_box_extent, rho, adot_accretion);
2601+
// PARTHENON_REQUIRE(adot_accretion >= 0, "Accretion is shrinking grains!");
2602+
// PARTHENON_REQUIRE(adot_accretion == adot_accretion , "adot_sputter is nan!");
2603+
}
26192604

26202605
adot_sputter = 0.;
26212606
adot_accretion = 0; // don;t do any dust updates if sputtering already is bad

‎src/hydro/srcterms/tabular_cooling.cpp‎

Lines changed: 11 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -409,7 +409,7 @@ void TabularCooling::SubcyclingFixedIntSrcTerm(MeshData<Real> *md, const Real dt
409409
IndexRange jb = md->GetBlockData(0)->GetBoundsJ(IndexDomain::entire);
410410
IndexRange kb = md->GetBlockData(0)->GetBoundsK(IndexDomain::entire);
411411

412-
par_for(
412+
par_for(
413413
DEFAULT_LOOP_PATTERN, "TabularCooling::SubcyclingSplitSrcTerm", DevExecSpace(), 0,
414414
cons_pack.GetDim(5) - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
415415
KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i) {
@@ -734,24 +734,18 @@ void TabularCooling::SubcyclingFixedIntSrcTerm(MeshData<Real> *md, const Real dt
734734
// Latter technically not required if no other tasks follows before
735735
// ConservedToPrim conversion, but keeping it for now (better safe than sorry).
736736
prim(IPR, k, j, i) = rho * internal_e * gm1;
737+
738+
737739
});
738740

739-
if (DustDevObj.agb_winds_on == 1 && DustObj.write_dust_history_to_file_ &&
740-
DustDevObj.dust_subcycle_with_cooling == 1) {
741-
Kokkos::Experimental::contribute(reduction_view_agb_injected_mass_c,
742-
scatter_f_agb_injected_mass_c);
743-
Kokkos::View<double *, Kokkos::LayoutRight, Kokkos::HostSpace>
744-
host_reduction_view_agb_injected_mass_c =
745-
Kokkos::create_mirror_view(reduction_view_agb_injected_mass_c);
746-
Kokkos::deep_copy(host_reduction_view_agb_injected_mass_c,
747-
reduction_view_agb_injected_mass_c);
748-
Kokkos::Experimental::contribute(reduction_view_agb_injected_mass_s,
749-
scatter_f_agb_injected_mass_s);
750-
Kokkos::View<double *, Kokkos::LayoutRight, Kokkos::HostSpace>
751-
host_reduction_view_agb_injected_mass_s =
752-
Kokkos::create_mirror_view(reduction_view_agb_injected_mass_s);
753-
Kokkos::deep_copy(host_reduction_view_agb_injected_mass_s,
754-
reduction_view_agb_injected_mass_s);
741+
742+
if(DustDevObj.agb_winds_on == 1 && DustObj.write_dust_history_to_file_ && DustDevObj.dust_subcycle_with_cooling==1){
743+
Kokkos::Experimental::contribute(reduction_view_agb_injected_mass_c, scatter_f_agb_injected_mass_c);
744+
Kokkos::View<double*, Kokkos::LayoutRight, Kokkos::HostSpace> host_reduction_view_agb_injected_mass_c = Kokkos::create_mirror_view(reduction_view_agb_injected_mass_c);
745+
Kokkos::deep_copy(host_reduction_view_agb_injected_mass_c, reduction_view_agb_injected_mass_c);
746+
Kokkos::Experimental::contribute(reduction_view_agb_injected_mass_s, scatter_f_agb_injected_mass_s);
747+
Kokkos::View<double*, Kokkos::LayoutRight, Kokkos::HostSpace> host_reduction_view_agb_injected_mass_s = Kokkos::create_mirror_view(reduction_view_agb_injected_mass_s);
748+
Kokkos::deep_copy(host_reduction_view_agb_injected_mass_s, reduction_view_agb_injected_mass_s);
755749
Kokkos::Experimental::contribute(reduction_view_stellar_mass, scatter_f_stellar_mass);
756750
Kokkos::View<double *, Kokkos::LayoutRight, Kokkos::HostSpace>
757751
host_reduction_view_stellar_mass =

0 commit comments

Comments
 (0)