17 #include "Utilities/ContainerHelpers.hpp"
64 template <
typename... FluxTags,
size_t Dim>
68 const Variables<tmpl::list<FluxTags...>>& fluxes,
const Mesh<Dim>& mesh,
69 const InverseJacobian<DataVector, Dim, Frame::Logical, Frame::Inertial>&
70 det_jac_times_inverse_jacobian) noexcept {
71 if (
UNLIKELY(divergence_of_fluxes->number_of_grid_points() !=
72 fluxes.number_of_grid_points())) {
73 divergence_of_fluxes->initialize(fluxes.number_of_grid_points());
76 const auto apply_matrix_in_first_dim =
77 [](
double* result,
const double*
const input,
const Matrix& matrix,
78 const size_t size,
const bool add_to_result) noexcept {
82 size / matrix.columns(),
96 if constexpr (Dim == 1) {
97 (void)det_jac_times_inverse_jacobian;
99 apply_matrix_in_first_dim(divergence_of_fluxes->data(), fluxes.data(),
101 fluxes.size(),
false);
104 const auto transform_to_logical_frame =
105 [&det_jac_times_inverse_jacobian, &fluxes](
109 auto logical_index_of_jacobian) noexcept {
110 using flux_tag = tmpl::type_from<decltype(flux_tag_v)>;
113 auto& result = get<div_tag>(*result_buffer);
114 const auto& flux = get<flux_tag>(fluxes);
116 for (
size_t result_storage_index = 0;
117 result_storage_index < result.size(); ++result_storage_index) {
118 const auto result_tensor_index =
119 result.get_tensor_index(result_storage_index);
120 const auto flux_x_tensor_index =
prepend(result_tensor_index, 0_st);
121 const auto flux_y_tensor_index =
prepend(result_tensor_index, 1_st);
122 if constexpr (Dim == 2) {
123 result[result_storage_index] =
125 0>(det_jac_times_inverse_jacobian) *
126 flux.get(flux_x_tensor_index) +
128 1>(det_jac_times_inverse_jacobian) *
129 flux.get(flux_y_tensor_index);
131 const auto flux_z_tensor_index =
132 prepend(result_tensor_index, 2_st);
133 result[result_storage_index] =
135 0>(det_jac_times_inverse_jacobian) *
136 flux.get(flux_x_tensor_index) +
138 1>(det_jac_times_inverse_jacobian) *
139 flux.get(flux_y_tensor_index) +
141 2>(det_jac_times_inverse_jacobian) *
142 flux.get(flux_z_tensor_index);
147 if constexpr (Dim == 2) {
148 Variables<tmpl::list<Tags::div<FluxTags>...>> data_buffer{
149 divergence_of_fluxes->number_of_grid_points()};
151 const Matrix& eta_weak_div_matrix =
153 const Matrix& xi_weak_div_matrix =
161 double* div_ptr = divergence_of_fluxes->data();
163 xi_weak_div_matrix.rows(),
164 divergence_of_fluxes->size() / xi_weak_div_matrix.rows());
165 apply_matrix_in_first_dim(data_buffer.data(),
166 divergence_of_fluxes->data(),
167 eta_weak_div_matrix, data_buffer.size(),
false);
169 const size_t chunk_size = Variables<tmpl::list<Tags::div<FluxTags>...>>::
170 number_of_independent_components *
171 eta_weak_div_matrix.rows();
173 data_buffer.size() / chunk_size);
179 apply_matrix_in_first_dim(divergence_of_fluxes->data(),
180 data_buffer.data(), xi_weak_div_matrix,
181 data_buffer.size(),
true);
182 }
else if constexpr (Dim == 3) {
183 Variables<tmpl::list<Tags::div<FluxTags>...>> data_buffer0{
184 divergence_of_fluxes->number_of_grid_points()};
185 Variables<tmpl::list<Tags::div<FluxTags>...>> data_buffer1{
186 divergence_of_fluxes->number_of_grid_points()};
187 constexpr
size_t number_of_independent_components =
188 decltype(data_buffer1)::number_of_independent_components;
190 const Matrix& zeta_weak_div_matrix =
192 const Matrix& eta_weak_div_matrix =
194 const Matrix& xi_weak_div_matrix =
203 xi_weak_div_matrix.rows() * eta_weak_div_matrix.rows();
204 double* result_ptr = data_buffer1.data();
206 data_buffer0.size() / chunk_size);
207 apply_matrix_in_first_dim(data_buffer0.data(), data_buffer1.data(),
208 zeta_weak_div_matrix, data_buffer0.size(),
211 number_of_independent_components * zeta_weak_div_matrix.rows();
212 result_ptr = divergence_of_fluxes->data();
214 data_buffer1.size() / chunk_size);
221 chunk_size = xi_weak_div_matrix.rows();
222 result_ptr = data_buffer1.data();
224 data_buffer1.size() / chunk_size);
225 apply_matrix_in_first_dim(data_buffer0.data(), data_buffer1.data(),
226 eta_weak_div_matrix, data_buffer0.size(),
228 chunk_size = number_of_independent_components *
229 eta_weak_div_matrix.rows() * zeta_weak_div_matrix.rows();
230 result_ptr = data_buffer1.data();
232 data_buffer0.size() / chunk_size);
233 *divergence_of_fluxes += data_buffer1;
239 apply_matrix_in_first_dim(divergence_of_fluxes->data(),
240 data_buffer0.data(), xi_weak_div_matrix,
241 data_buffer0.size(),
true);
243 static_assert(Dim == 1 or Dim == 2 or Dim == 3,
244 "Weak divergence only implemented in 1d, 2d, and 3d.");