Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <array>
7 : #include <cstddef>
8 : #include <optional>
9 : #include <random>
10 : #include <vector>
11 :
12 : #include "DataStructures/DataBox/DataBox.hpp"
13 : #include "DataStructures/Tensor/EagerMath/RaiseOrLowerIndex.hpp"
14 : #include "DataStructures/Variables.hpp"
15 : #include "Domain/Tags.hpp"
16 : #include "Domain/TagsTimeDependent.hpp"
17 : #include "Evolution/DgSubcell/ActiveGrid.hpp"
18 : #include "Evolution/DgSubcell/Tags/ActiveGrid.hpp"
19 : #include "Evolution/DgSubcell/Tags/Coordinates.hpp"
20 : #include "Evolution/DgSubcell/Tags/Jacobians.hpp"
21 : #include "Evolution/DgSubcell/Tags/Mesh.hpp"
22 : #include "Evolution/Particles/MonteCarlo/InverseJacobianInertialToFluidCompute.hpp"
23 : #include "Evolution/Particles/MonteCarlo/MortarData.hpp"
24 : #include "Evolution/Particles/MonteCarlo/NeutrinoInteractionTable.hpp"
25 : #include "Evolution/Particles/MonteCarlo/Tags.hpp"
26 : #include "Evolution/Particles/MonteCarlo/TemplatedLocalFunctions.hpp"
27 : #include "Evolution/Systems/GeneralizedHarmonic/Tags.hpp"
28 : #include "Parallel/AlgorithmExecution.hpp"
29 : #include "Parallel/GlobalCache.hpp"
30 : #include "PointwiseFunctions/GeneralRelativity/DerivativeSpatialMetric.hpp"
31 : #include "PointwiseFunctions/GeneralRelativity/InverseSpacetimeMetric.hpp"
32 : #include "PointwiseFunctions/GeneralRelativity/SpacetimeNormalVector.hpp"
33 : #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
34 : #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
35 : #include "PointwiseFunctions/Hydro/Tags.hpp"
36 : #include "Time/Tags/TimeStepId.hpp"
37 : #include "Time/Time.hpp"
38 : #include "Time/TimeStepId.hpp"
39 :
40 : namespace Particles::MonteCarlo {
41 :
42 : /// Mutator advancing neutrinos by a single step
43 : template <size_t EnergyBins, size_t NeutrinoSpecies>
44 1 : struct TimeStepMutator {
45 0 : static const size_t Dim = 3;
46 :
47 0 : using return_tags =
48 : tmpl::list<Particles::MonteCarlo::Tags::PacketsOnElement,
49 : Particles::MonteCarlo::Tags::CouplingTildeTau<DataVector>,
50 : Particles::MonteCarlo::Tags::CouplingTildeRhoYe<DataVector>,
51 : Particles::MonteCarlo::Tags::CouplingTildeS<DataVector, Dim>,
52 : Particles::MonteCarlo::Tags::RandomNumberGenerator,
53 : Particles::MonteCarlo::Tags::DesiredPacketEnergyAtEmission<
54 : NeutrinoSpecies>>;
55 : // To do : check carefully DG vs Subcell quantities... everything should
56 : // be on the Subcell grid!
57 0 : using argument_tags = tmpl::list<
58 : ::Tags::TimeStepId, ::Tags::Next<::Tags::TimeStepId>,
59 : hydro::Tags::GrmhdEquationOfState,
60 : Particles::MonteCarlo::Tags::InteractionRatesTable<EnergyBins,
61 : NeutrinoSpecies>,
62 : hydro::Tags::ElectronFraction<DataVector>,
63 : hydro::Tags::RestMassDensity<DataVector>,
64 : hydro::Tags::Temperature<DataVector>,
65 : hydro::Tags::LorentzFactor<DataVector>,
66 : hydro::Tags::SpatialVelocity<DataVector, 3, Frame::Inertial>,
67 : gr::Tags::Lapse<DataVector>,
68 : gr::Tags::Shift<DataVector, Dim, Frame::Inertial>,
69 : gh::Tags::Phi<DataVector, 3, Frame::Inertial>,
70 : gr::Tags::SpatialMetric<DataVector, Dim, Frame::Inertial>,
71 : gr::Tags::InverseSpatialMetric<DataVector, Dim, Frame::Inertial>,
72 : gr::Tags::SqrtDetSpatialMetric<DataVector>,
73 : Particles::MonteCarlo::Tags::CellLightCrossingTime<DataVector>,
74 : evolution::dg::subcell::Tags::Mesh<Dim>,
75 : evolution::dg::subcell::Tags::Coordinates<Dim, Frame::ElementLogical>,
76 : domain::Tags::MeshVelocity<Dim>,
77 : evolution::dg::subcell::fd::Tags::InverseJacobianLogicalToInertial<Dim>,
78 : evolution::dg::subcell::fd::Tags::DetInverseJacobianLogicalToInertial,
79 : domain::Tags::InverseJacobian<Dim + 1, Frame::Inertial, Frame::Fluid>,
80 : domain::Tags::Jacobian<Dim + 1, Frame::Inertial, Frame::Fluid>,
81 : Particles::MonteCarlo::Tags::MortarDataTag<Dim>>;
82 :
83 0 : static void apply(
84 : const gsl::not_null<std::vector<Packet>*> packets,
85 : const gsl::not_null<Scalar<DataVector>*> coupling_tilde_tau,
86 : const gsl::not_null<Scalar<DataVector>*> coupling_tilde_rho_ye,
87 : const gsl::not_null<tnsr::i<DataVector, Dim>*> coupling_tilde_s,
88 : const gsl::not_null<std::mt19937*> random_number_generator,
89 : const gsl::not_null<std::array<DataVector, NeutrinoSpecies>*>
90 : single_packet_energy,
91 : const TimeStepId& current_step_id, const TimeStepId& next_step_id,
92 :
93 : const EquationsOfState::EquationOfState<true, 3>& equation_of_state,
94 : const NeutrinoInteractionTable<EnergyBins, NeutrinoSpecies>&
95 : interaction_table,
96 : const Scalar<DataVector>& electron_fraction,
97 : const Scalar<DataVector>& rest_mass_density,
98 : const Scalar<DataVector>& temperature,
99 : const Scalar<DataVector>& lorentz_factor,
100 : const tnsr::I<DataVector, Dim, Frame::Inertial>& spatial_velocity,
101 : const Scalar<DataVector>& lapse,
102 : const tnsr::I<DataVector, Dim, Frame::Inertial>& shift,
103 : const tnsr::iaa<DataVector, Dim, Frame::Inertial>& phi,
104 :
105 : const tnsr::ii<DataVector, Dim, Frame::Inertial>& spatial_metric,
106 : const tnsr::II<DataVector, Dim, Frame::Inertial>& inv_spatial_metric,
107 : const Scalar<DataVector>& sqrt_determinant_spatial_metric,
108 : const Scalar<DataVector>& cell_light_crossing_time, const Mesh<Dim>& mesh,
109 : const tnsr::I<DataVector, Dim, Frame::ElementLogical>& mesh_coordinates,
110 : const std::optional<tnsr::I<DataVector, Dim, Frame::Inertial>>&
111 : mesh_velocity,
112 : const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
113 : Frame::Inertial>&
114 : inverse_jacobian_logical_to_inertial,
115 : const Scalar<DataVector>& det_inverse_jacobian_logical_to_inertial,
116 : const InverseJacobian<DataVector, Dim + 1, Frame::Inertial, Frame::Fluid>&
117 : inertial_to_fluid_inverse_jacobian,
118 : const Jacobian<DataVector, Dim + 1, Frame::Inertial, Frame::Fluid>&
119 : inertial_to_fluid_jacobian,
120 : const MortarData<Dim>& mortar_data) {
121 : // Number of ghost zones for MC is assumed to be 1 for now.
122 : const size_t num_ghost_zones = 1;
123 : // Get information stored in various databox containers in
124 : // the format expected by take_time_step_on_element
125 : const double start_time = current_step_id.step_time().value();
126 : const double end_time = next_step_id.step_time().value();
127 : Scalar<DataVector> det_jacobian_logical_to_inertial(lapse);
128 : get(det_jacobian_logical_to_inertial) =
129 : 1.0 / get(det_inverse_jacobian_logical_to_inertial);
130 : const DirectionalIdMap<Dim, std::optional<DataVector>>&
131 : electron_fraction_ghost = mortar_data.electron_fraction;
132 : const DirectionalIdMap<Dim, std::optional<DataVector>>&
133 : baryon_density_ghost = mortar_data.rest_mass_density;
134 : const DirectionalIdMap<Dim, std::optional<DataVector>>& temperature_ghost =
135 : mortar_data.temperature;
136 : const DirectionalIdMap<Dim, std::optional<DataVector>>&
137 : cell_light_crossing_time_ghost = mortar_data.cell_light_crossing_time;
138 :
139 : // Calculate temporary tensors needed for MC evolution
140 : using deriv_lapse = ::Tags::deriv<gr::Tags::Lapse<DataVector>,
141 : tmpl::size_t<3>, Frame::Inertial>;
142 : using deriv_shift = ::Tags::deriv<gr::Tags::Shift<DataVector, 3>,
143 : tmpl::size_t<3>, Frame::Inertial>;
144 : using deriv_spatial_metric =
145 : ::Tags::deriv<gr::Tags::SpatialMetric<DataVector, 3>, tmpl::size_t<3>,
146 : Frame::Inertial>;
147 : using deriv_inverse_spatial_metric =
148 : ::Tags::deriv<gr::Tags::InverseSpatialMetric<DataVector, 3>,
149 : tmpl::size_t<3>, Frame::Inertial>;
150 : using temporary_tags = tmpl::list<
151 : hydro::Tags::LowerSpatialFourVelocity<DataVector, Dim, Frame::Inertial>,
152 : gr::Tags::SpacetimeNormalVector<DataVector, 3>,
153 : gr::Tags::InverseSpacetimeMetric<DataVector, 3>, deriv_lapse,
154 : deriv_shift, deriv_spatial_metric, deriv_inverse_spatial_metric>;
155 : Variables<temporary_tags> temp_tags{mesh.number_of_grid_points(), 0.0};
156 :
157 : // u_i = \gamma_{ij} v^j W
158 : auto& lower_spatial_four_velocity =
159 : get<hydro::Tags::LowerSpatialFourVelocity<DataVector, Dim,
160 : Frame::Inertial>>(temp_tags);
161 : raise_or_lower_index(make_not_null(&lower_spatial_four_velocity),
162 : spatial_velocity, spatial_metric);
163 : for (size_t i = 0; i < Dim; i++) {
164 : lower_spatial_four_velocity.get(i) *= get(lorentz_factor);
165 : }
166 : // For the metric, we adapt the calculations performed for the time
167 : // derivative of in GhGrMhd. First get n^a and g^ab
168 : auto& spacetime_normal_vector =
169 : get<gr::Tags::SpacetimeNormalVector<DataVector, 3>>(temp_tags);
170 : auto& inv_spacetime_metric =
171 : get<gr::Tags::InverseSpacetimeMetric<DataVector, 3>>(temp_tags);
172 : gr::spacetime_normal_vector(make_not_null(&spacetime_normal_vector), lapse,
173 : shift);
174 : gr::inverse_spacetime_metric(make_not_null(&inv_spacetime_metric), lapse,
175 : shift, inv_spatial_metric);
176 :
177 : auto& d_lapse = get<deriv_lapse>(temp_tags);
178 : auto& d_shift = get<deriv_shift>(temp_tags);
179 : // Temporary store phi_iab n^a n^b in d_lapse. This is phi_two_normals in GH
180 : for (size_t i = 0; i < Dim; i++) {
181 : for (size_t a = 0; a < Dim + 1; a++) {
182 : for (size_t b = 0; b < Dim + 1; b++) {
183 : d_lapse.get(i) += phi.get(i, a, b) * spacetime_normal_vector.get(a) *
184 : spacetime_normal_vector.get(b);
185 : }
186 : }
187 : }
188 : // Shift derivative using stored quantity in d_lapse
189 : // We use d_i shift^j =
190 : // (g^{j+1 b} phi_{iba} n^a + n^{j+1} phi_{iab} n^a n_b) * lapse
191 : // as in TimeDerivative.hpp in GhGrMhd
192 : for (size_t i = 0; i < Dim; i++) {
193 : for (size_t j = 0; j < Dim; j++) {
194 : d_shift.get(i, j) +=
195 : d_lapse.get(i) * spacetime_normal_vector.get(j + 1);
196 : for (size_t a = 0; a < Dim + 1; a++) {
197 : for (size_t b = 0; b < Dim + 1; b++) {
198 : d_shift.get(i, j) += inv_spacetime_metric.get(j + 1, b) *
199 : phi.get(i, b, a) *
200 : spacetime_normal_vector.get(a);
201 : }
202 : }
203 : d_shift.get(i, j) *= get(lapse);
204 : }
205 : }
206 : // Now use d_i lapse = - lapse * 0.5 * phi_{iab} n^a n^b
207 : // As we already stored phi_{iab} n^a n^b in d_i lapse,
208 : // we just multiply by (-0.5 * lapse)
209 : for (size_t i = 0; i < Dim; i++) {
210 : d_lapse.get(i) *= (-0.5) * get(lapse);
211 : }
212 :
213 : // Extract d_i \gamma_{jk} from phi_{i,j+1,k+1}
214 : auto& d_spatial_metric = get<deriv_spatial_metric>(temp_tags);
215 : for (size_t i = 0; i < Dim; i++) {
216 : for (size_t j = 0; j < Dim; j++) {
217 : for (size_t k = j; k < Dim; k++) {
218 : d_spatial_metric.get(i, j, k) = phi.get(i, j + 1, k + 1);
219 : }
220 : }
221 : }
222 :
223 : auto& d_inv_spatial_metric = get<deriv_inverse_spatial_metric>(temp_tags);
224 : gr::deriv_inverse_spatial_metric(make_not_null(&d_inv_spatial_metric),
225 : inv_spatial_metric, d_spatial_metric);
226 :
227 : TemplatedLocalFunctions<EnergyBins, NeutrinoSpecies> templated_functions;
228 : templated_functions.take_time_step_on_element(
229 : packets, coupling_tilde_tau, coupling_tilde_rho_ye, coupling_tilde_s,
230 : random_number_generator, single_packet_energy, start_time, end_time,
231 : equation_of_state, interaction_table, electron_fraction,
232 : rest_mass_density, temperature, lorentz_factor,
233 : lower_spatial_four_velocity, lapse, shift, d_lapse, d_shift,
234 : d_inv_spatial_metric, spatial_metric, inv_spatial_metric,
235 : sqrt_determinant_spatial_metric, cell_light_crossing_time, mesh,
236 : mesh_coordinates, num_ghost_zones, mesh_velocity,
237 : inverse_jacobian_logical_to_inertial, det_jacobian_logical_to_inertial,
238 : inertial_to_fluid_jacobian, inertial_to_fluid_inverse_jacobian,
239 : electron_fraction_ghost, baryon_density_ghost, temperature_ghost,
240 : cell_light_crossing_time_ghost);
241 : }
242 : };
243 :
244 : namespace Actions {
245 :
246 : /// Action taking a single time step of the Monte-Carlo evolution
247 : /// algorithm, assuming that the fluid and metric data in the ghost
248 : /// zones have been communicated and that packets are on the elements
249 : /// that owns them.
250 : template <size_t EnergyBins, size_t NeutrinoSpecies>
251 1 : struct TakeTimeStep {
252 : template <typename DbTags, typename... InboxTags, typename ArrayIndex,
253 : typename ActionList, typename ParallelComponent,
254 : typename Metavariables>
255 0 : static Parallel::iterable_action_return_t apply(
256 : db::DataBox<DbTags>& box, tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
257 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
258 : const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
259 : const ParallelComponent* const /*meta*/) {
260 : ASSERT(db::get<evolution::dg::subcell::Tags::ActiveGrid>(box) ==
261 : evolution::dg::subcell::ActiveGrid::Subcell,
262 : "MC assumes that we are using the Subcell grid!");
263 :
264 : db::mutate_apply(TimeStepMutator<EnergyBins, NeutrinoSpecies>{},
265 : make_not_null(&box));
266 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
267 : }
268 : };
269 :
270 : } // namespace Actions
271 : } // namespace Particles::MonteCarlo
|
