diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 637978aed..e2479f70d 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -28,6 +28,7 @@ jobs: kilonova_1d_1dgrid, kilonova_1d_3dgrid, kilonova_2d_2dgrid, + kilonova_2d_2dgrid_barnesthermalisation, kilonova_2d_2dgrid_expansionopac, kilonova_2d_2dgrid_xcomgammaphotoion, kilonova_2d_3dgrid, diff --git a/artisoptions_christinenonthermal.h b/artisoptions_christinenonthermal.h index 9243a2541..4d5b917af 100644 --- a/artisoptions_christinenonthermal.h +++ b/artisoptions_christinenonthermal.h @@ -133,8 +133,6 @@ constexpr bool USE_CALCULATED_MEANATOMICWEIGHT = false; constexpr bool WRITE_PARTIAL_EMISSIONABSORPTIONSPEC = false; -constexpr bool INSTANT_PARTICLE_DEPOSITION = true; - constexpr enum timestepsizemethods TIMESTEP_SIZE_METHOD = TIMESTEP_SIZES_LOGARITHMIC; constexpr double FIXED_TIMESTEP_WIDTH = -1.; @@ -151,6 +149,8 @@ constexpr float RPKT_BOUNDBOUND_THERMALISATION_PROBABILITY = -1.; constexpr bool USE_XCOM_GAMMAPHOTOION = false; +constexpr auto PARTICLE_THERMALISATION_SCHEME = ThermalisationScheme::INSTANT; + constexpr auto GAMMA_THERMALISATION_SCHEME = ThermalisationScheme::DETAILED; // NOLINTEND(modernize*,misc-unused-parameters) diff --git a/artisoptions_classic.h b/artisoptions_classic.h index e04f1ff3a..1900ec104 100644 --- a/artisoptions_classic.h +++ b/artisoptions_classic.h @@ -131,8 +131,6 @@ constexpr bool USE_CALCULATED_MEANATOMICWEIGHT = false; constexpr bool WRITE_PARTIAL_EMISSIONABSORPTIONSPEC = false; -constexpr bool INSTANT_PARTICLE_DEPOSITION = true; - constexpr enum timestepsizemethods TIMESTEP_SIZE_METHOD = TIMESTEP_SIZES_LOGARITHMIC; constexpr double FIXED_TIMESTEP_WIDTH = -1.; @@ -149,6 +147,8 @@ constexpr float RPKT_BOUNDBOUND_THERMALISATION_PROBABILITY = -1.; constexpr bool USE_XCOM_GAMMAPHOTOION = false; +constexpr auto PARTICLE_THERMALISATION_SCHEME = ThermalisationScheme::INSTANT; + constexpr auto GAMMA_THERMALISATION_SCHEME = ThermalisationScheme::DETAILED; // NOLINTEND(modernize*,misc-unused-parameters) diff --git a/artisoptions_kilonova_lte.h b/artisoptions_kilonova_lte.h index 871236f09..3d38dda29 100644 --- a/artisoptions_kilonova_lte.h +++ b/artisoptions_kilonova_lte.h @@ -130,8 +130,6 @@ constexpr bool USE_CALCULATED_MEANATOMICWEIGHT = true; constexpr bool WRITE_PARTIAL_EMISSIONABSORPTIONSPEC = false; -constexpr bool INSTANT_PARTICLE_DEPOSITION = false; - constexpr enum timestepsizemethods TIMESTEP_SIZE_METHOD = TIMESTEP_SIZES_LOGARITHMIC; constexpr double FIXED_TIMESTEP_WIDTH = -1.; @@ -148,6 +146,8 @@ constexpr float RPKT_BOUNDBOUND_THERMALISATION_PROBABILITY = -1.; constexpr bool USE_XCOM_GAMMAPHOTOION = false; +constexpr auto PARTICLE_THERMALISATION_SCHEME = ThermalisationScheme::DETAILED; + constexpr auto GAMMA_THERMALISATION_SCHEME = ThermalisationScheme::DETAILED; // NOLINTEND(modernize*,misc-unused-parameters) diff --git a/artisoptions_nltenebular.h b/artisoptions_nltenebular.h index ec1f7b4f5..8194ad7db 100644 --- a/artisoptions_nltenebular.h +++ b/artisoptions_nltenebular.h @@ -142,8 +142,6 @@ constexpr bool USE_CALCULATED_MEANATOMICWEIGHT = false; constexpr bool WRITE_PARTIAL_EMISSIONABSORPTIONSPEC = false; -constexpr bool INSTANT_PARTICLE_DEPOSITION = true; - constexpr enum timestepsizemethods TIMESTEP_SIZE_METHOD = TIMESTEP_SIZES_LOGARITHMIC; constexpr double FIXED_TIMESTEP_WIDTH = -1.; @@ -160,6 +158,8 @@ constexpr float RPKT_BOUNDBOUND_THERMALISATION_PROBABILITY = -1.; constexpr bool USE_XCOM_GAMMAPHOTOION = false; +constexpr auto PARTICLE_THERMALISATION_SCHEME = ThermalisationScheme::INSTANT; + constexpr auto GAMMA_THERMALISATION_SCHEME = ThermalisationScheme::DETAILED; // NOLINTEND(modernize*,misc-unused-parameters) diff --git a/artisoptions_nltewithoutnonthermal.h b/artisoptions_nltewithoutnonthermal.h index 5978d1f6d..e144d015c 100644 --- a/artisoptions_nltewithoutnonthermal.h +++ b/artisoptions_nltewithoutnonthermal.h @@ -134,8 +134,6 @@ constexpr bool USE_CALCULATED_MEANATOMICWEIGHT = false; constexpr bool WRITE_PARTIAL_EMISSIONABSORPTIONSPEC = false; -constexpr bool INSTANT_PARTICLE_DEPOSITION = true; - constexpr enum timestepsizemethods TIMESTEP_SIZE_METHOD = TIMESTEP_SIZES_LOGARITHMIC; constexpr double FIXED_TIMESTEP_WIDTH = 0.1; @@ -152,6 +150,8 @@ constexpr float RPKT_BOUNDBOUND_THERMALISATION_PROBABILITY = -1.; constexpr bool USE_XCOM_GAMMAPHOTOION = false; +constexpr auto PARTICLE_THERMALISATION_SCHEME = ThermalisationScheme::INSTANT; + constexpr auto GAMMA_THERMALISATION_SCHEME = ThermalisationScheme::DETAILED; // NOLINTEND(modernize*,misc-unused-parameters) diff --git a/constants.h b/constants.h index a9d9d9c0a..87b94c019 100644 --- a/constants.h +++ b/constants.h @@ -67,6 +67,6 @@ enum timestepsizemethods { TIMESTEP_SIZES_CONSTANT_THEN_LOGARITHMIC = 3, }; -enum class ThermalisationScheme { DETAILED, BARNES_GLOBAL, BARNES_LOCAL }; +enum class ThermalisationScheme { INSTANT, DETAILED, BARNES }; #endif diff --git a/gammapkt.cc b/gammapkt.cc index d084a92e9..99f5009a8 100644 --- a/gammapkt.cc +++ b/gammapkt.cc @@ -969,35 +969,31 @@ void transport_gamma(Packet &pkt, double t2) } } -void barnes_thermalization(Packet &pkt, bool local) +void barnes_thermalisation(Packet &pkt) // Barnes treatment: packet is either getting absorbed immediately and locally // creating a k-packet or it escapes. The absorption probability matches the // Barnes thermalization efficiency, for expressions see the original paper: // https://ui.adsabs.harvard.edu/abs/2016ApJ...829..110B { // compute thermalization efficiency (= absorption probability) - constexpr double mean_gamma_opac = 0.1; - - // determine average initial density - double V_0 = 0.; - for (int nonemptymgi = 0; nonemptymgi < grid::get_nonempty_npts_model(); nonemptymgi++) { - const int mgi = grid::get_mgi_of_nonemptymgi(nonemptymgi); - V_0 += grid::get_modelcell_assocvolume_tmin(mgi); - } - double rho_0 = 0.; - if (!local) { - rho_0 = grid::mtot_input / V_0; - } else { - rho_0 = grid::get_rho_tmin(grid::get_cell_modelgridindex(pkt.where)); - } - - const double R_0 = pow(3 * V_0 / (4 * PI), 1 / 3.); - const double t_0 = grid::get_t_model(); - const double tau_ineff = sqrt(rho_0 * R_0 * pow(t_0, 2) * mean_gamma_opac); + // 0.1 is an average value to fit the analytic approximations from the paper. + // Alternative: Distinguish between low-E (kappa = 1) or high-E (kappa = 0.05) + // packets. + // constexpr double mean_gamma_opac = 0.1; + + // determine average initial density via kinetic energy + const double E_kin = grid::get_ejecta_kinetic_energy(); + const double v_ej = sqrt(E_kin * 2 / grid::mtot_input); + const double t_0 = globals::tmin; + + // const double t_ineff = sqrt(rho_0 * R_0 * pow(t_0, 2) * mean_gamma_opac); + const double t_ineff = 1.4 * 86400. * sqrt(grid::mtot_input / (5.e-3 * 1.989 * 1.e33)) * ((0.2 * 29979200000) / v_ej); // get current time const double t = t_0 + pkt.prop_time; - const double tau = pow(tau_ineff / t, 2.); + const double tau = pow(t_ineff / t, 2.); const double f_gamma = 1. - exp(-tau); + assert_always(f_gamma >= 0.); + assert_always(f_gamma <= 1.); // either absorb packet or let it escape if (rng_uniform() < f_gamma) { @@ -1014,10 +1010,8 @@ void barnes_thermalization(Packet &pkt, bool local) void do_gamma(Packet &pkt, double t2) { if constexpr (GAMMA_THERMALISATION_SCHEME == ThermalisationScheme::DETAILED) { transport_gamma(pkt, t2); - } else if constexpr (GAMMA_THERMALISATION_SCHEME == ThermalisationScheme::BARNES_GLOBAL) { - barnes_thermalization(pkt, false); - } else if constexpr (GAMMA_THERMALISATION_SCHEME == ThermalisationScheme::BARNES_LOCAL) { - barnes_thermalization(pkt, true); + } else if constexpr (GAMMA_THERMALISATION_SCHEME == ThermalisationScheme::BARNES) { + barnes_thermalisation(pkt); } else { __builtin_unreachable(); } diff --git a/grid.cc b/grid.cc index 9a72d9aac..5f3c194cc 100644 --- a/grid.cc +++ b/grid.cc @@ -626,7 +626,7 @@ static void set_elem_stable_abund_from_total(const int mgi, const int element, c modelgrid[mgi].composition[element].abundance = isofracsum + massfracstable; } -static auto get_cellradialposmid(const int cellindex) -> double +auto get_cellradialposmid(const int cellindex) -> double // get the radial distance from the origin to the centre of the cell at time tmin { if (GRID_TYPE == GRID_SPHERICAL1D) { diff --git a/grid.h b/grid.h index a10798e0b..72bba2793 100644 --- a/grid.h +++ b/grid.h @@ -82,6 +82,7 @@ void set_elements_uppermost_ion(int modelgridindex, int element, int newvalue); [[nodiscard]] auto get_cellcoordmax(int cellindex, int axis) -> double; [[nodiscard]] auto get_cellcoordmin(int cellindex, int axis) -> double; [[nodiscard]] auto get_cellcoordpointnum(int cellindex, int axis) -> int; +[[nodiscard]] auto get_cellradialposmid(int cellindex) -> double; [[nodiscard]] auto get_coordcellindexincrement(int axis) -> int; [[nodiscard]] auto get_rho_tmin(int modelgridindex) -> float; [[nodiscard]] auto get_rho(int modelgridindex) -> float; @@ -156,6 +157,19 @@ inline void change_cell(Packet &pkt, const int snext) } } +inline auto get_ejecta_kinetic_energy() { + double E_kin = 0.; + for (int nonemptymgi = 0; nonemptymgi < grid::get_nonempty_npts_model(); nonemptymgi++) { + const int mgi = grid::get_mgi_of_nonemptymgi(nonemptymgi); + const int assoc_cells = grid::get_numassociatedcells(mgi); + double M_cell = grid::get_rho_tmin(mgi) * grid::get_modelcell_assocvolume_tmin(mgi); + const double radial_pos = grid::modelgrid[mgi].initial_radial_pos_sum / assoc_cells; + E_kin += 0.5 * M_cell * std::pow(radial_pos / globals::tmin, 2); + } + + return E_kin; +} + } // namespace grid #endif // GRIDINIT_H diff --git a/tests/kilonova_2d_2dgrid_barnesthermalisation_inputfiles/phixsdata_v2.txt b/tests/kilonova_2d_2dgrid_barnesthermalisation_inputfiles/phixsdata_v2.txt new file mode 100644 index 000000000..249e1c191 --- /dev/null +++ b/tests/kilonova_2d_2dgrid_barnesthermalisation_inputfiles/phixsdata_v2.txt @@ -0,0 +1,2 @@ +100 + 3.0000000e-02 diff --git a/tests/kilonova_2d_2dgrid_barnesthermalisation_inputfiles/results_md5_final.txt b/tests/kilonova_2d_2dgrid_barnesthermalisation_inputfiles/results_md5_final.txt new file mode 100644 index 000000000..0f9fd3ab3 --- /dev/null +++ b/tests/kilonova_2d_2dgrid_barnesthermalisation_inputfiles/results_md5_final.txt @@ -0,0 +1,324 @@ +94090a73b642cd79b4208d34570e5c23 absorption.out +fb7beca851f29fb6d65385671fffdb57 absorption_res_00.out +adaa2826953ac4096930dff9dacf0d5c absorption_res_01.out +7344690a948ddfc603dcced766c4b5b6 absorption_res_02.out +39911c110ee79aa5fa51d48971cefa17 absorption_res_03.out +bc714fe852d7e234590b25a6533e45ac absorption_res_04.out +4b7a4f11a11e57e2fb66d787fe8660bb absorption_res_05.out +106adab262f9234cac82adce9dac3918 absorption_res_06.out +ae48c6f51a67aafd9833d220f97192c5 absorption_res_07.out +97230f16fe902a5aaac5210a625661c5 absorption_res_08.out +66d668b41921cfe581c832cafe89e1db absorption_res_09.out +fbddf2cfd97af2e57aa385fa24d6201c absorption_res_10.out +0fcfda4869243258019a72f48b7867d5 absorption_res_11.out +f2d3979d49991144be2bdbefe836fb69 absorption_res_12.out +c02020563bf66e9ff380188d71a2a866 absorption_res_13.out +f42e4ba4121ea49e2edaa623d4b04b4b absorption_res_14.out +0ce7160a2da38507a7c9047b82afb91c absorption_res_15.out +a5e7ec3545fbfe942ce45004ffbdbe2d absorption_res_16.out +7eda519a04e42d82d63695a85848ede6 absorption_res_17.out +728817befeaedc60520b6904f11ad076 absorption_res_18.out +228c256ae49356fc327b2535b9f3fbff absorption_res_19.out +44e19bbb31dcef1446bc47b51da98947 absorption_res_20.out +112dac2ba4efb244fd7d241ccd4f97b0 absorption_res_21.out +ca5898dd263c4dad666015c7d281428b absorption_res_22.out +6515511bb9956e602a37ce8c70c15896 absorption_res_23.out +117cf9220638ab943df515369a4d8abc absorption_res_24.out +106a896343829f47b949dd43d4c175b6 absorption_res_25.out +0603d8c4009649962a21aa77e1136fc9 absorption_res_26.out +541374ddc59ac621e01f6a3526bfcb3f absorption_res_27.out +eb41c81645f17a6e47333971acb71c0e absorption_res_28.out +639a404bc36ba87355e37bfbb5ac2941 absorption_res_29.out +a2138feb96a829dd99473beed7f58e1d absorption_res_30.out +c1268d1421702dc119b03cfbcb807bec absorption_res_31.out +de34aac06e303e71fff926197a61ed05 absorption_res_32.out +3ede93fec5dbc4975eebd18ee9f7588e absorption_res_33.out +1c79444ab7361d99535b462059a27dbc absorption_res_34.out +143969f3fc05282a4f6b40d6a0780933 absorption_res_35.out +626fa790afc439822038ec24f8df215d absorption_res_36.out +7d30c56f39c63acf9c39c16aa23ae0fe absorption_res_37.out +1aa9a31d3d341a059273107c67f7622d absorption_res_38.out +5827e42c2f836c35630d7d79c9421320 absorption_res_39.out +01820c7ba5a95b29fd799d66cff87203 absorption_res_40.out +308bdab5337138fb4cab9e5983de16ea absorption_res_41.out +b477e21ea7e5403a947a825db27e2f5e absorption_res_42.out +c721d8eaeb6dfc36ee7a513da3446718 absorption_res_43.out +acc3ef7d155705c7f073e42d85fa3490 absorption_res_44.out +c1e0e28812ed0361bfadc6a755085d34 absorption_res_45.out +c31e19aeba478d0c29283f031ab31b7b absorption_res_46.out +41eee3c7220528ddaa420f67598bbf37 absorption_res_47.out +59d3e1724b45e8435900eabdd6ffc81e absorption_res_48.out +7ab8c4b2087685de82373c887f024af8 absorption_res_49.out +075bc2932af9227566eb71315954b917 absorption_res_50.out +1f42e0c35c43e94c2f0e926fc9cf7303 absorption_res_51.out +763bc79eb2c4f5e3c7f0c4592fba14a2 absorption_res_52.out +4055b805831ce3e2ca8e16d93becef32 absorption_res_53.out +416edf1a831f7d301c5f62c2c7fa3537 absorption_res_54.out +c631b917488b0568bdf0f8c281e1b366 absorption_res_55.out +a5e6a181252ad097ffb8fb8ac5398acd absorption_res_56.out +b6071144acd24f7d7776e1babbbd872c absorption_res_57.out +ab81627695045115888bf42af140630c absorption_res_58.out +1c80e46950996e823537121918bfaa42 absorption_res_59.out +559eed3c52f311717ebedb467d910a06 absorption_res_60.out +efe106d8e76799d9e23804ae88e1058e absorption_res_61.out +325d7777e283a1d82e3a9bd659b28626 absorption_res_62.out +0a610b52650ad8427ecdd1e462b71ebc absorption_res_63.out +b4a708d6f51f539e0dd38d3c96278732 absorption_res_64.out +0bb0e1b87091720282a5e4fc044a8eb9 absorption_res_65.out +c544d469c5e9124a1eef1b45348d7d26 absorption_res_66.out +f9d00f24af556ec4265690df01cf96e6 absorption_res_67.out +1324fe23d3e62925883bc529e0013755 absorption_res_68.out +dd9e1eae994a5fcc6e845ceb85b9800d absorption_res_69.out +52f82d4f08ed82ffb07ab84747f26bd5 absorption_res_70.out +a2165f13ec6971a0c6f0cc18244fad2f absorption_res_71.out +c0fbc6862fc513b7eab8eddfc2710e4e absorption_res_72.out +db03e654cc52773398683d95eae0dc12 absorption_res_73.out +69ca416174283f14ba95563020afcd97 absorption_res_74.out +c6b475da35c6064920bf4a810a018395 absorption_res_75.out +3504b37263b345bdeddbf90f5ee55719 absorption_res_76.out +3fcfb6b66759ae77a910509eb12d325e absorption_res_77.out +85bdf8e7d96616263fc41a83f9a60ccf absorption_res_78.out +bf5ae82ea6cea7a4d230851267d09f7a absorption_res_79.out +36633e6926a72dbacb392827450cdca6 absorption_res_80.out +386eb6bdea23245350d722e1f95397da absorption_res_81.out +30e63a6bf62c7ea7747d5487a6003642 absorption_res_82.out +aa1fc13ee5053d107cf548eba48fda8d absorption_res_83.out +3d2029da9f210996dbeb5bb56d4c26d1 absorption_res_84.out +178614cc4780a91ac6392940d35e89eb absorption_res_85.out +31560406504646fd0b1ff4d11832d595 absorption_res_86.out +2b85200ebdafe053deb0aec2f5e9dc63 absorption_res_87.out +177c9bf208a0d74ca486cc49b492ee00 absorption_res_88.out +62a5af963168051fe917b9560cae72fe absorption_res_89.out +d99248eb10a567a7619dceb5eb3dfe29 absorption_res_90.out +9fa4ce77ced9d040cf007494961a28f1 absorption_res_91.out +c77c5d43e9731017dd9afda49b926d94 absorption_res_92.out +2a74c3574ca783275a5915f2f294e9d2 absorption_res_93.out +7bc356661b1b2eaaad505cba1ad1642b absorption_res_94.out +3ed73811b0fb9e65b65018861e05078d absorption_res_95.out +b143b567d2185ac9d423c25663dbbf8f absorption_res_96.out +2c3499c1092c036f15f6732a30dc5548 absorption_res_97.out +2467844c3d885cce6c235057526bdcd3 absorption_res_98.out +e7e78274efcde7d1034d7c56d459ff15 absorption_res_99.out +897316929176464ebc9ad085f31e7284 bflist.out +cd7d9b8340014d5bcc934e71ed6a223b deposition.out +4a5c4759d39721abe5a6ee49d31da3e6 emission.out +cc572d39e872d0104401bb4c4b0fbe92 emission_res_00.out +7bd493fd4e5c762d32d4a140d175300e emission_res_01.out +0e7d16d1a2db4742d78b2eab46c2be58 emission_res_02.out +105ed3f21aa264952c083b060361159b emission_res_03.out +f6864f61f86b1b4acd61147a5f2667df emission_res_04.out +517611008c0c318a2dd2ffde63ed63dc emission_res_05.out +fc42248fa52970762b22ee4face386ae emission_res_06.out +e120802527bb7f9c39af967db50507e0 emission_res_07.out +7fa2a19e9d2fcd7c91891a03e16da841 emission_res_08.out +6b9b541503d4a4bba2696f93862afecf emission_res_09.out +f41b88d0c5747763b0b39263bee9194f emission_res_10.out +5ec7f5ab71a81e1931135e9645e513fb emission_res_11.out +77a5ebc034e53efe24dca03b7527e107 emission_res_12.out +9ff5e8665ba28bfe42aa7f0dbd90e440 emission_res_13.out +327d902a0aebd019a246fe88f137bb51 emission_res_14.out +e4ba77dd58f4bfac1f29cd2c72014da1 emission_res_15.out +e2fc750b6b7161d40ae91ce39300e181 emission_res_16.out +3b1976f718ad58d67278d8b323bede10 emission_res_17.out +ddeb8c804520a45a3a9052683c9a661c emission_res_18.out +e536ec04fe43d80ecf999b73a3b7bf0d emission_res_19.out +67709c86acf269ae150e6fe9e3f207fc emission_res_20.out +cca17d9aec474b7774c95c7155034204 emission_res_21.out +2e7f7cdaee0ad6583bca0133a117a44c emission_res_22.out +e4eb9b7709d2969fff0123a2f7376a5f emission_res_23.out +2e010aa1dde94daf4141ebd091045387 emission_res_24.out +b87219aa7c6440d3f99ddfb120669a41 emission_res_25.out +e5e760be1eea94aa006c10ebcc0c4c6d emission_res_26.out +f31be96f0d7e0e215de61ba079a29c35 emission_res_27.out +b4e210592643bffaa0c5c4ae4e32a2f0 emission_res_28.out +4d66fc69e47b9a66ef6ac1b3e9cc210a emission_res_29.out +b3e3ea9af0ffe16be2c0b091bcb8079f emission_res_30.out +75bd85f2c4c63fb5a52e3cbe33989349 emission_res_31.out +8e0eb10a22811bad8499a782f5b9b3ef emission_res_32.out +fc25207bb79ced180f75a177093ed4a3 emission_res_33.out +d81722d5b10552f34e54cbcc33d70fa4 emission_res_34.out +040fcc8945a22105f9be9c4ebc8280b2 emission_res_35.out +61df5fd7d564441f3f8bf51cded8fd43 emission_res_36.out +1a26f2fc12671af7b0674b6aef129a42 emission_res_37.out +ebb5ed5317273a7ca9f98d1f039c7782 emission_res_38.out +b6b904b46f171b48edd7a7f51ebe4528 emission_res_39.out +8140deeafa0552d7744b04b7ad2dd667 emission_res_40.out +cb19d22e9727ffe492372ab54c25315f emission_res_41.out +a701a6db53451f120de811926dc9fc61 emission_res_42.out +2881ad70bed0d4c959d5ec8c1d51c527 emission_res_43.out +7ebda8c9b9d62ea8c1e50e8914e750ce emission_res_44.out +ff7691998e62c2b551820c01a0437634 emission_res_45.out +7fa565cfaa5031dd381e9f43c486bdc0 emission_res_46.out +e70ca2b14a85818bf17b0d8fe0a25a47 emission_res_47.out +43a51b13050b40da664b211e7b8df13f emission_res_48.out +ebd1ded8c011d5314204028e58c95bc4 emission_res_49.out +91038b64eea91f0b2c93f8ec50ea8697 emission_res_50.out +d2b4262a1a9faf8a32d811589057849f emission_res_51.out +0722ccaa187b4976240224efa365ebd3 emission_res_52.out +7a5a265a228451caed1bac42e41ff303 emission_res_53.out +8ae134a7d3a5ea9ce848c93e7b4fd4e4 emission_res_54.out +cea3be41bbeb139b6d74f32a6ec73fb6 emission_res_55.out +aed3743a86c454c7d672cc9917db222f emission_res_56.out +9c2fac120c5f4a46dff749ef170d5d58 emission_res_57.out +7dfdf248e4297a58dea69d6757d41d40 emission_res_58.out +3152d17723b7da39a56da8f28f67e02f emission_res_59.out +b398b8235c6956c6f9d1071046186c37 emission_res_60.out +ec71ba6aef599b8850c5064ba938810c emission_res_61.out +fdec8ef659eb96f2f03c1eb0717911b6 emission_res_62.out +e2d5d197d5ac590356320471bd977f0e emission_res_63.out +1213f561aca21663891309fc81d2cc30 emission_res_64.out +a9972b67ae0b4ed5217b85289f98d12c emission_res_65.out +714f59f20f2e1e66c197e18b38f7ec37 emission_res_66.out +5701feb4ec2d85c1be73bcb035354a71 emission_res_67.out +ec835a4729b4f2773ac3762a272531dd emission_res_68.out +49bb10b6895de80c3b1e0f1e2fb818b0 emission_res_69.out +3600530580ba5f49da20eeb8fd404b40 emission_res_70.out +0028c79eb1ce191e1eb8aa3ea0741dc7 emission_res_71.out +d263493bbdf347fd0ceb2ebf0759c6b8 emission_res_72.out +bcb5ec7bacd141128b0f59ec10c77c07 emission_res_73.out +cdf0999af79b8b80d422a29f7badfb25 emission_res_74.out +4c382bdafc53e258c29ea2440c3ff811 emission_res_75.out +e5cf051dbbe37c6920c22aae42c1df93 emission_res_76.out +6ee260eb40de74f181c27a7099cf8e20 emission_res_77.out +b2fa3fa154fd575f7fba680654fa6490 emission_res_78.out +3c4ef9b642470e19803e2ac5b5e051a3 emission_res_79.out +e1cc1a7caf77f2a300e3973f42aa3bb8 emission_res_80.out +f51f46d3bca818193cb17320e6bc6e42 emission_res_81.out +14b5601403fecbc8159303bd8d5473a4 emission_res_82.out +14e6d024a8fd42c52888a4ee707e06d3 emission_res_83.out +85f730ac2fc51460e999070c29a64f95 emission_res_84.out +53679f48e5d37e4a3d10365c06317908 emission_res_85.out +506ceb097ba4f3093f4e5fd3a06a58a7 emission_res_86.out +d5267ed66ae0cd646bf8aa5e810980a9 emission_res_87.out +2c5dd3ed2469736bcc4b0e7ee7ed48b7 emission_res_88.out +359c5c6abce22f742e389070465656a2 emission_res_89.out +b6ab5515f84ccb9d88458501d3ac7440 emission_res_90.out +90b25ded9658ad9932a34f42f919a528 emission_res_91.out +c0fd74a87912fe4e4bd64acbf8e668d5 emission_res_92.out +6b299c031abac6dc80ef69096000be5c emission_res_93.out +f33732a20a11af9684a4db3c0eef17b3 emission_res_94.out +664472b4310c3761b8967ea73e865219 emission_res_95.out +49ad0ddc60b06f0bad2f4ee8dce1f729 emission_res_96.out +1e7af5fe6b9e7ff6134b1599299af251 emission_res_97.out +d35808e987c7aef6b22a5416bf17c7b0 emission_res_98.out +8b3f2bd8ee37437d0f7dfdd323c1a6e5 emission_res_99.out +81e206ca6fccc71f5f15af8b85e67ec5 emissiontrue.out +b5fb06ef79ab1001e58006c827162a29 emissiontrue_res_00.out +21e3fb9e41f14362034c3025021c020d emissiontrue_res_01.out +b08e0d7dca103924d22431a79342cd84 emissiontrue_res_02.out +194ac5603dc2a0c41f272780bc161ae9 emissiontrue_res_03.out +8040e5fa9003b6c2955d28db357ceb83 emissiontrue_res_04.out +ed70e0209be4b1c70a4235f48715cd91 emissiontrue_res_05.out +6e4ebc782100a987ea85413f8b7db5a6 emissiontrue_res_06.out +9c39333b2f73f63b378214b6bff0a39b emissiontrue_res_07.out +0213c3aaa1cdf7ae2357fcd48f783fd6 emissiontrue_res_08.out +3a8030378ea116936a291b00115741dd emissiontrue_res_09.out +e3780c7ca484041cc3b59b84e4504248 emissiontrue_res_10.out +a458655ca53641c549cfdca52f83b877 emissiontrue_res_11.out +0ecb29fb03203e805a364149e04523d6 emissiontrue_res_12.out +f46c89692a7453a41a7f7e7a35fd7d06 emissiontrue_res_13.out +a49fc4929d4ebae2cdd22d121b74b061 emissiontrue_res_14.out +26d637af05f9385287ea2f02fa7d8f9b emissiontrue_res_15.out +1338ed4e2fe5447bb3b7585acf57d026 emissiontrue_res_16.out +c0174a2839545a8366ba996957d8c07b emissiontrue_res_17.out +aa919856cb214f0b8add0e01f3e788d0 emissiontrue_res_18.out +b339e492d8862313b95ad4f99a66b89e emissiontrue_res_19.out +59aacbd232b1ceb029c07829a3dd65ac emissiontrue_res_20.out +c1d66a041b4b03c2902bdcd8ac1fe33d emissiontrue_res_21.out +0b17ab8ca6fe0cd88ce8b4d64643e9e4 emissiontrue_res_22.out +abe8de0a2cc6096bfaf23c0769391ccd emissiontrue_res_23.out +6d9274d2ff1fd0b71b9cbbea2084214b emissiontrue_res_24.out +ec3875fbcb728d453ed99be585b22ddf emissiontrue_res_25.out +c5aea4f89bd310772f0d2c681fb40bfd emissiontrue_res_26.out +6f7e9eb20600e6287d4d8fadc71f85ba emissiontrue_res_27.out +754fd1862f2affc4dbd982f69de3707e emissiontrue_res_28.out +c111024f3e1184c5ac30faa959ea4fd6 emissiontrue_res_29.out +6e0dceef8b17d1a0fe1628702517f300 emissiontrue_res_30.out +8b43a167dbabd3cb68e653816fc2e4b7 emissiontrue_res_31.out +a3c3c3ad10570ca68e797feea0c23164 emissiontrue_res_32.out +4fbc6f9bf7502b90cd8360b6b5148d8a emissiontrue_res_33.out +2c7ac8715615eddc8205a08191532340 emissiontrue_res_34.out +753452c80b805d6d88606dfb90d7452f emissiontrue_res_35.out +5fbe80041ad2f8f11191016b3271dfdd emissiontrue_res_36.out +38268b799ab43aa5946498f0a27f7e3e emissiontrue_res_37.out +2b9b3b102e9948eb29d5c4f157cbacf2 emissiontrue_res_38.out +39a8c05e0d385c6c10b794bfb09857d7 emissiontrue_res_39.out +9e3b740705a9e8a9e0ab1783bff00e38 emissiontrue_res_40.out +2d5777c6a9b3718094775352330aeb62 emissiontrue_res_41.out +c264104715625233089b2ac6df291247 emissiontrue_res_42.out +4d9e84788c0bf2aada603025a88956fb emissiontrue_res_43.out +2805238c9d564cbd15247159c130579b emissiontrue_res_44.out +ff30bace38faf2cd3c9b23b6fb31d48a emissiontrue_res_45.out +112d2a753fd0d0b950f3e376bbe1563e emissiontrue_res_46.out +b8e6f601c2834a8dd0a7d1adf2c82144 emissiontrue_res_47.out +0b6c4816459aab0614439437d947d336 emissiontrue_res_48.out +4887b8ce41f4ee7d2f265e2d83218b6f emissiontrue_res_49.out +d216be13765ab0aafe41bd786b26aee2 emissiontrue_res_50.out +744e72ac70721a955522d66ae9eee3bc emissiontrue_res_51.out +4691f93494336dc680191201da3b7e22 emissiontrue_res_52.out +b13dd07587371f5db1761b23efe53940 emissiontrue_res_53.out +9525b7ffcf22e5d040a41491b2872136 emissiontrue_res_54.out +7609a67586253f6fa95cdfeb314565bd emissiontrue_res_55.out +034487f40428c3fb6ccdae85839d9b0b emissiontrue_res_56.out +500c500d5d30b1809ea47107f854a0c5 emissiontrue_res_57.out +0de771068d8c4bebf36e09f0b9881c93 emissiontrue_res_58.out +553ab62fb8f2b57a7742982ade62b043 emissiontrue_res_59.out +99ee5791f3f09c08a15d98168d713859 emissiontrue_res_60.out +dc6abb0e032a4a2cc015a9ad63096b14 emissiontrue_res_61.out +168fab67369e28df4b628602a416dff0 emissiontrue_res_62.out +e0d7d797e264cf2ae0f0d65cdb7af562 emissiontrue_res_63.out +ae573e3a59a960aee5150b72ca779d0c emissiontrue_res_64.out +63c3e3598c109f8005476876efbb2ec1 emissiontrue_res_65.out +0f3fd9d9e337a16dfca92733b53a9704 emissiontrue_res_66.out +db9fca748a0eaad85bee94937f060b11 emissiontrue_res_67.out +669c0988ec271b5cbf34eea76c0d36a4 emissiontrue_res_68.out +8142a70ab2cbe54fe1593be9c5f5265b emissiontrue_res_69.out +d800404628e346b4c604123416d0619d emissiontrue_res_70.out +d77fbfec0fb8eaa44b4ebb10aa583941 emissiontrue_res_71.out +fd396095df75c8057d1f4965787c77fc emissiontrue_res_72.out +56209362fc78524261fbbc9b287e4c66 emissiontrue_res_73.out +6691d9fc1f4125ec29f0b2ee4aacb390 emissiontrue_res_74.out +543d617d5644bcf5eec598431086ea74 emissiontrue_res_75.out +18d07da54360b10b01dcc09ba5a2b2bd emissiontrue_res_76.out +e7fdcd69e1cfb754783ab16fcf000454 emissiontrue_res_77.out +f5f64d841da92ac5cf64799902d62b52 emissiontrue_res_78.out +4b7232550f8a30462fa80ad0206f9a6d emissiontrue_res_79.out +24949075ee7070dd1fe5d0e358beeb5c emissiontrue_res_80.out +13c99407609964cf927928e266babd44 emissiontrue_res_81.out +b05c9252aba7a3257c9c06c81c9990a0 emissiontrue_res_82.out +5fe373eb02284dbedd17d5f41d92bf26 emissiontrue_res_83.out +b65866401ca0201901ea71f8d5c09157 emissiontrue_res_84.out +49176cc16a676fa807403ff1a06ad3e2 emissiontrue_res_85.out +2c8c261fdda28f8f7195cea9958b2740 emissiontrue_res_86.out +6f774234794470bb853ce14ebdcaa44a emissiontrue_res_87.out +d5f9631b947fe103497da7499b8da18f emissiontrue_res_88.out +dbc97ce31d7542505032a0fefc21138c emissiontrue_res_89.out +703832d6bcfba95085ce88a1cec9c5a4 emissiontrue_res_90.out +23d2eea0eb23f0010b7e91065505f2f4 emissiontrue_res_91.out +c050766c9cd2c7ed57fdc02b808e066a emissiontrue_res_92.out +2204a74ff67b7e458d6ee1f55d0d6d49 emissiontrue_res_93.out +e9c98f92ef74f6bc4fe7258e721a4374 emissiontrue_res_94.out +6c0a8e96d005296d35d0780f20775164 emissiontrue_res_95.out +3bb8a8663fa200efb866baaa02841fc1 emissiontrue_res_96.out +0f6a82a362f741a00821f218b6410ff1 emissiontrue_res_97.out +ae37411bb87929ebba53f7d98e028fb2 emissiontrue_res_98.out +6035bce2fcd58b263b879ca0a666437a emissiontrue_res_99.out +7d6e090f3e7cda4461a3965dc38a211e gamma_light_curve.out +484c6a1384f4bf203f31ce715419c6b1 gamma_spec.out +2b769145664780d4f90b07f963588536 gammalinelist.out +29ac1cb06f3139df7cbca0bbdb426e1c grid.out +26e7f83d49072c35020f3e170b91334b light_curve.out +6d8c77ab3451b894939484918b64973d light_curve_res.out +9773becefdfe4afffb041b703210c67c linestat.out +d80043ce6bedca56511b9fd3f7c0a49f modelgridrankassignments.out +824fa63d9d0e8a16b56e60d80ef49a1e packets00_0000.out +ea2bb9d00a595ddea110c5a1c33ba3ca packets00_0001.out +0dea2332a5243b956945abc873b41216 packets00_0002.out +36fc78247a45a5fe182ccc34ac0fc0af packets00_0003.out +6eef5da0bab37e29fd2cb3b9ae04832c spec.out +787ea9b17d6b86d3236c12b89f49c9a3 spec_res.out +a351f1711fecd60c023d0ba7332092db timesteps.out +266a3473fefd978b5ce3fe46bb3c14bd job1/estimators_0000.out +feacc4f2b8eeed8bfadd9aa7dc88242d job1/estimators_0001.out +ce09e41e172a11bfa0a3285dd6fd9580 job1/estimators_0002.out +81dce9f8a1eb0ec858ceccd6fac0df6f job1/estimators_0003.out diff --git a/tests/kilonova_2d_2dgrid_barnesthermalisation_inputfiles/results_md5_job0.txt b/tests/kilonova_2d_2dgrid_barnesthermalisation_inputfiles/results_md5_job0.txt new file mode 100644 index 000000000..7c2b4e448 --- /dev/null +++ b/tests/kilonova_2d_2dgrid_barnesthermalisation_inputfiles/results_md5_job0.txt @@ -0,0 +1,18 @@ +897316929176464ebc9ad085f31e7284 bflist.out +6bd54ee9700306c4f1da537fe8939aeb deposition.out +267fa0f38b2763b29eb4bebb9c3d110c gamma_light_curve.out +2b769145664780d4f90b07f963588536 gammalinelist.out +29ac1cb06f3139df7cbca0bbdb426e1c grid.out +a5f6eed3aa3de1375784bbf175e42494 light_curve.out +9773becefdfe4afffb041b703210c67c linestat.out +d80043ce6bedca56511b9fd3f7c0a49f modelgridrankassignments.out +878862f601178dd48aa4dcf2945f3cd0 packets00_0000.out +7e2db0149c0241f41716846e5390c382 packets00_0001.out +c1977627040085c730a8ce9f41b0daa8 packets00_0002.out +37f35d601bb085b82e3bf9554d4c315f packets00_0003.out +bc999ebbcee373fd633d2f0dcdcc3b67 spec.out +a351f1711fecd60c023d0ba7332092db timesteps.out +424b6e13fe6f6d93ac13890d1a1edf71 job0/estimators_0000.out +64488872c5b47c516c6cc692c668550a job0/estimators_0001.out +a956fa38a510440ba22a60acab3b4859 job0/estimators_0002.out +c5158a986c7b9f377ce1f03c32c46ac2 job0/estimators_0003.out diff --git a/tests/setup_kilonova_2d_2dgrid_barnesthermalisation.sh b/tests/setup_kilonova_2d_2dgrid_barnesthermalisation.sh new file mode 100755 index 000000000..a6811a044 --- /dev/null +++ b/tests/setup_kilonova_2d_2dgrid_barnesthermalisation.sh @@ -0,0 +1,41 @@ +#!/usr/bin/env zsh + +set -x + +runfolder=kilonova_2d_2dgrid_barnesthermalisation_testrun + +mkdir -p $runfolder + +if [ ! -f atomicdata_feconi.tar.xz ]; then curl -O https://theory.gsi.de/~lshingle/artis_http_public/artis/atomicdata_feconi.tar.xz; fi + +tar -xf atomicdata_feconi.tar.xz --directory $runfolder/ + +# same input files as the other test run +rsync -av kilonova_2d_3dgrid_inputfiles/ $runfolder/ + +# for the checksum files +rsync -av --ignore-times kilonova_2d_2dgrid_barnesthermalisation_inputfiles/ $runfolder/ + +cp ../data/* $runfolder/ + +cp ../artisoptions_kilonova_lte.h $runfolder/artisoptions.h + +cd $runfolder + +xz -dv -T0 *.xz + +sed -i'' -e 's/constexpr int MPKTS.*/constexpr int MPKTS = 80000;/g' artisoptions.h + +sed -i'' -e 's/constexpr int GRID_TYPE.*/constexpr int GRID_TYPE = GRID_CYLINDRICAL2D;/g' artisoptions.h + +sed -i'' -e 's/constexpr int TABLESIZE.*/constexpr int TABLESIZE = 20;/g' artisoptions.h +sed -i'' -e 's/constexpr double MINTEMP.*/constexpr double MINTEMP = 1000.;/g' artisoptions.h +sed -i'' -e 's/constexpr double MAXTEMP.*/constexpr double MAXTEMP = 20000.;/g' artisoptions.h + +sed -i'' -e 's/constexpr auto PARTICLE_THERMALISATION_SCHEME.*/constexpr auto PARTICLE_THERMALISATION_SCHEME = ThermalisationScheme::BARNES;/g' artisoptions.h + +sed -i'' -e 's/constexpr auto GAMMA_THERMALISATION_SCHEME.*/constexpr auto GAMMA_THERMALISATION_SCHEME = ThermalisationScheme::BARNES;/g' artisoptions.h + +cd - + +set +x diff --git a/update_grid.cc b/update_grid.cc index 3825716d9..a7e9a2cce 100644 --- a/update_grid.cc +++ b/update_grid.cc @@ -897,6 +897,7 @@ void update_grid_cell(const int mgi, const int nts, const int nts_prev, const in /// Update elemental abundances with radioactive decays decay::update_abundances(mgi, nts, globals::timesteps[nts].mid); + if (globals::opacity_case == 6) { grid::calculate_kappagrey(); } diff --git a/update_packets.cc b/update_packets.cc index f58cd5f79..f30239da0 100644 --- a/update_packets.cc +++ b/update_packets.cc @@ -29,11 +29,30 @@ namespace { void do_nonthermal_predeposit(Packet &pkt, const int nts, const double t2) { double en_deposited = pkt.e_cmf; + const double ts = pkt.prop_time; - if constexpr (INSTANT_PARTICLE_DEPOSITION) { + if constexpr (PARTICLE_THERMALISATION_SCHEME == ThermalisationScheme::INSTANT) { // absorption happens pkt.type = TYPE_NTLEPTON; + } else if constexpr (PARTICLE_THERMALISATION_SCHEME == ThermalisationScheme::BARNES) { + const double E_kin = grid::get_ejecta_kinetic_energy(); + const double v_ej = sqrt(E_kin * 2 / grid::mtot_input); + + const double prefactor = (pkt.pellet_decaytype == decay::DECAYTYPE_ALPHA) ? 7.74 : 7.4; + const double tau_ineff = + prefactor * 86400 * sqrt(grid::mtot_input / (5.e-3 * 1.989 * 1.e33)) * pow((0.2 * 29979200000) / v_ej, 3. / 2.); + const double f_p = log(1 + 2. * ts * ts / tau_ineff / tau_ineff) / (2. * ts * ts / tau_ineff / tau_ineff); + assert_always(f_p >= 0.); + assert_always(f_p <= 1.); + if (rng_uniform() < f_p) { + pkt.type = TYPE_NTLEPTON; + } else { + en_deposited = 0.; + pkt.type = TYPE_ESCAPE; + grid::change_cell(pkt, -99); + } } else { + // local, detailed absorption following Shingles+2023 const double rho = grid::get_rho(grid::get_cell_modelgridindex(pkt.where)); // endot is energy loss rate (positive) in [erg/s]