SpECTRE Documentation Coverage Report
Current view: top level - Evolution/Systems/NewtonianEuler/Limiters - CharacteristicHelpers.hpp Hit Total Coverage
Commit: 9a905b0737f373631c1b8e8389b8f26e67fa5313 Lines: 6 8 75.0 %
Date: 2024-03-28 09:03:18
Legend: Lines: hit not hit

          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 <limits>
       9             : #include <string>
      10             : #include <type_traits>
      11             : #include <unordered_map>
      12             : #include <utility>
      13             : 
      14             : #include "DataStructures/Tags.hpp"
      15             : #include "DataStructures/Tags/TempTensor.hpp"
      16             : #include "DataStructures/Tensor/Tensor.hpp"
      17             : #include "DataStructures/Variables.hpp"
      18             : #include "Evolution/Systems/NewtonianEuler/Characteristics.hpp"
      19             : #include "Evolution/Systems/NewtonianEuler/SoundSpeedSquared.hpp"
      20             : #include "Evolution/Systems/NewtonianEuler/Tags.hpp"
      21             : #include "NumericalAlgorithms/LinearOperators/MeanValue.hpp"
      22             : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
      23             : #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
      24             : #include "Utilities/Gsl.hpp"
      25             : #include "Utilities/TMPL.hpp"
      26             : #include "Utilities/TaggedTuple.hpp"
      27             : 
      28           0 : namespace NewtonianEuler::Limiters {
      29             : 
      30             : /*!
      31             :  * \brief Compute the transform matrices between the conserved variables and
      32             :  * the characteristic variables of the NewtonianEuler system.
      33             :  *
      34             :  * Wraps calls to `NewtonianEuler::right_eigenvectors` and
      35             :  * `NewtonianEuler::left_eigenvectors`.
      36             :  */
      37             : template <size_t VolumeDim>
      38           1 : std::pair<Matrix, Matrix> right_and_left_eigenvectors(
      39             :     const Scalar<double>& mean_density,
      40             :     const tnsr::I<double, VolumeDim>& mean_momentum,
      41             :     const Scalar<double>& mean_energy,
      42             :     const EquationsOfState::EquationOfState<false, 2>& equation_of_state,
      43             :     const tnsr::i<double, VolumeDim>& unit_vector,
      44             :     bool compute_char_transformation_numerically = false);
      45             : 
      46             : /// @{
      47             : /// \brief Compute characteristic fields from conserved fields
      48             : ///
      49             : /// Note that these functions apply the same transformation to every grid point
      50             : /// in the element, using the same matrix from `left_eigenvectors`.
      51             : /// This is in contrast to the characteristic transformation used in the
      52             : /// GeneralizedHarmonic upwind flux, which is computed pointwise.
      53             : ///
      54             : /// By using a fixed matrix of eigenvectors computed from the cell-averaged
      55             : /// fields, we ensure we can consistently transform back from the characteristic
      56             : /// variables even after applying the limiter (the limiter changes the pointwise
      57             : /// values but preserves the cell averages).
      58             : template <size_t VolumeDim>
      59           1 : void characteristic_fields(
      60             :     gsl::not_null<tuples::TaggedTuple<
      61             :         ::Tags::Mean<NewtonianEuler::Tags::VMinus>,
      62             :         ::Tags::Mean<NewtonianEuler::Tags::VMomentum<VolumeDim>>,
      63             :         ::Tags::Mean<NewtonianEuler::Tags::VPlus>>*>
      64             :         char_means,
      65             :     const tuples::TaggedTuple<
      66             :         ::Tags::Mean<NewtonianEuler::Tags::MassDensityCons>,
      67             :         ::Tags::Mean<NewtonianEuler::Tags::MomentumDensity<VolumeDim>>,
      68             :         ::Tags::Mean<NewtonianEuler::Tags::EnergyDensity>>& cons_means,
      69             :     const Matrix& left);
      70             : 
      71             : template <size_t VolumeDim>
      72           1 : void characteristic_fields(
      73             :     gsl::not_null<Scalar<DataVector>*> char_v_minus,
      74             :     gsl::not_null<tnsr::I<DataVector, VolumeDim>*> char_v_momentum,
      75             :     gsl::not_null<Scalar<DataVector>*> char_v_plus,
      76             :     const Scalar<DataVector>& cons_mass_density,
      77             :     const tnsr::I<DataVector, VolumeDim>& cons_momentum_density,
      78             :     const Scalar<DataVector>& cons_energy_density, const Matrix& left);
      79             : 
      80             : template <size_t VolumeDim>
      81           1 : void characteristic_fields(
      82             :     gsl::not_null<
      83             :         Variables<tmpl::list<NewtonianEuler::Tags::VMinus,
      84             :                              NewtonianEuler::Tags::VMomentum<VolumeDim>,
      85             :                              NewtonianEuler::Tags::VPlus>>*>
      86             :         char_vars,
      87             :     const Variables<tmpl::list<NewtonianEuler::Tags::MassDensityCons,
      88             :                                NewtonianEuler::Tags::MomentumDensity<VolumeDim>,
      89             :                                NewtonianEuler::Tags::EnergyDensity>>& cons_vars,
      90             :     const Matrix& left);
      91             : /// @}
      92             : 
      93             : /// \brief Compute conserved fields from characteristic fields
      94             : ///
      95             : /// Note that this function applies the same transformation to every grid point
      96             : /// in the element, using the same matrix from `right_eigenvectors`.
      97             : /// This is in contrast to the characteristic transformation used in the
      98             : /// GeneralizedHarmonic upwind flux, which is computed pointwise.
      99             : ///
     100             : /// By using a fixed matrix of eigenvectors computed from the cell-averaged
     101             : /// fields, we ensure we can consistently transform back from the characteristic
     102             : /// variables even after applying the limiter (the limiter changes the pointwise
     103             : /// values but preserves the cell averages).
     104             : template <size_t VolumeDim>
     105           1 : void conserved_fields_from_characteristic_fields(
     106             :     gsl::not_null<Scalar<DataVector>*> cons_mass_density,
     107             :     gsl::not_null<tnsr::I<DataVector, VolumeDim>*> cons_momentum_density,
     108             :     gsl::not_null<Scalar<DataVector>*> cons_energy_density,
     109             :     const Scalar<DataVector>& char_v_minus,
     110             :     const tnsr::I<DataVector, VolumeDim>& char_v_momentum,
     111             :     const Scalar<DataVector>& char_v_plus, const Matrix& right);
     112             : 
     113             : /// \brief Apply a limiter to the characteristic fields computed with respect
     114             : /// to each direction in the volume, then take average of results
     115             : ///
     116             : /// When computing the characteristic fields in the volume to pass as inputs to
     117             : /// the limiter, it is not necessarily clear (in more than one dimension) which
     118             : /// unit vector should be used in the characteristic decomposition. This is in
     119             : /// contrast to uses of characteristic fields for, e.g., computing boundary
     120             : /// conditions where the boundary normal provides a clear choice of unit vector.
     121             : //
     122             : /// The common solution to this challenge is to average the results of applying
     123             : /// the limiter with different choices of characteristic decomposition:
     124             : /// - For each direction (e.g. in 3D for each of
     125             : ///   \f$(\hat{x}, \hat{y}, \hat{z})\f$), compute the characteristic fields with
     126             : ///   respect to this direction. This gives `VolumeDim` sets of characteristic
     127             : ///   fields.
     128             : /// - Apply the limiter to each set of characteristic fields, then convert back
     129             : ///   to conserved fields. This gives `VolumeDim` different sets of limited
     130             : ///   conserved fields.
     131             : /// - Average the new conserved fields; this is the result of the limiting
     132             : ///   process.
     133             : ///
     134             : /// This function handles the logic of computing the different characteristic
     135             : /// fields, limiting them, converting back to conserved fields, and averaging.
     136             : ///
     137             : /// The limiter to apply is passed in via a lambda which is responsible for
     138             : /// converting any neighbor data to the characteristic representation, and then
     139             : /// applying the limiter to all NewtonianEuler characteristic fields.
     140             : template <size_t VolumeDim, typename LimiterLambda>
     141           1 : bool apply_limiter_to_characteristic_fields_in_all_directions(
     142             :     const gsl::not_null<Scalar<DataVector>*> mass_density_cons,
     143             :     const gsl::not_null<tnsr::I<DataVector, VolumeDim>*> momentum_density,
     144             :     const gsl::not_null<Scalar<DataVector>*> energy_density,
     145             :     const Mesh<VolumeDim>& mesh,
     146             :     const EquationsOfState::EquationOfState<false, 2>& equation_of_state,
     147             :     const LimiterLambda& prepare_and_apply_limiter,
     148             :     const bool compute_char_transformation_numerically = false) {
     149             :   // Temp variables for calculations
     150             :   // There are quite a few tensors in this allocation because in general we need
     151             :   // to preserve the input (cons field) tensors until they are overwritten at
     152             :   // the very end of the computation... then we need additional buffers for the
     153             :   // char fields, the limited cons fields, and the accumulated cons fields for
     154             :   // averaging.
     155             :   //
     156             :   // Possible optimization: specialize the 1D case which doesn't do average so
     157             :   // doesn't need as many buffers
     158             :   Variables<tmpl::list<
     159             :       NewtonianEuler::Tags::VMinus, NewtonianEuler::Tags::VMomentum<VolumeDim>,
     160             :       NewtonianEuler::Tags::VPlus, ::Tags::TempScalar<0>,
     161             :       ::Tags::TempI<0, VolumeDim>, ::Tags::TempScalar<1>, ::Tags::TempScalar<2>,
     162             :       ::Tags::TempI<1, VolumeDim>, ::Tags::TempScalar<3>>>
     163             :       temp_buffer(mesh.number_of_grid_points());
     164             :   auto& char_v_minus = get<NewtonianEuler::Tags::VMinus>(temp_buffer);
     165             :   auto& char_v_momentum =
     166             :       get<NewtonianEuler::Tags::VMomentum<VolumeDim>>(temp_buffer);
     167             :   auto& char_v_plus = get<NewtonianEuler::Tags::VPlus>(temp_buffer);
     168             :   auto& temp_mass_density_cons = get<::Tags::TempScalar<0>>(temp_buffer);
     169             :   auto& temp_momentum_density = get<::Tags::TempI<0, VolumeDim>>(temp_buffer);
     170             :   auto& temp_energy_density = get<::Tags::TempScalar<1>>(temp_buffer);
     171             :   auto& accumulate_mass_density_cons = get<::Tags::TempScalar<2>>(temp_buffer);
     172             :   auto& accumulate_momentum_density =
     173             :       get<::Tags::TempI<1, VolumeDim>>(temp_buffer);
     174             :   auto& accumulate_energy_density = get<::Tags::TempScalar<3>>(temp_buffer);
     175             : 
     176             :   // Initialize the accumulating tensors
     177             :   get(accumulate_mass_density_cons) = 0.;
     178             :   for (size_t i = 0; i < VolumeDim; ++i) {
     179             :     accumulate_momentum_density.get(i) = 0.;
     180             :   }
     181             :   get(accumulate_energy_density) = 0.;
     182             : 
     183             :   // Cellwise means, used in computing the cons/char transformations
     184             :   const auto mean_density =
     185             :       Scalar<double>{mean_value(get(*mass_density_cons), mesh)};
     186             :   const auto mean_momentum = [&momentum_density, &mesh]() {
     187             :     tnsr::I<double, VolumeDim> result{};
     188             :     for (size_t i = 0; i < VolumeDim; ++i) {
     189             :       result.get(i) = mean_value(momentum_density->get(i), mesh);
     190             :     }
     191             :     return result;
     192             :   }();
     193             :   const auto mean_energy =
     194             :       Scalar<double>{mean_value(get(*energy_density), mesh)};
     195             : 
     196             :   bool some_component_was_limited = false;
     197             : 
     198             :   // Loop over directions, then compute chars w.r.t. this direction and limit
     199             :   for (size_t d = 0; d < VolumeDim; ++d) {
     200             :     const auto unit_vector = [&d]() {
     201             :       auto components = make_array<VolumeDim>(0.);
     202             :       components[d] = 1.;
     203             :       return tnsr::i<double, VolumeDim>(components);
     204             :     }();
     205             :     const auto right_and_left =
     206             :         NewtonianEuler::Limiters::right_and_left_eigenvectors(
     207             :             mean_density, mean_momentum, mean_energy, equation_of_state,
     208             :             unit_vector, compute_char_transformation_numerically);
     209             :     const auto& right = right_and_left.first;
     210             :     const auto& left = right_and_left.second;
     211             : 
     212             :     // Transform tensors to characteristics
     213             :     NewtonianEuler::Limiters::characteristic_fields(
     214             :         make_not_null(&char_v_minus), make_not_null(&char_v_momentum),
     215             :         make_not_null(&char_v_plus), *mass_density_cons, *momentum_density,
     216             :         *energy_density, left);
     217             : 
     218             :     // Transform neighbor data and apply limiter
     219             :     const bool some_component_was_limited_with_this_unit_vector =
     220             :         prepare_and_apply_limiter(make_not_null(&char_v_minus),
     221             :                                   make_not_null(&char_v_momentum),
     222             :                                   make_not_null(&char_v_plus), left);
     223             : 
     224             :     some_component_was_limited =
     225             :         some_component_was_limited_with_this_unit_vector or
     226             :         some_component_was_limited;
     227             : 
     228             :     // Transform back to conserved variables. But skip the transformation if no
     229             :     // limiting occured with this unit vector.
     230             :     if (some_component_was_limited_with_this_unit_vector) {
     231             :       NewtonianEuler::Limiters::conserved_fields_from_characteristic_fields(
     232             :           make_not_null(&temp_mass_density_cons),
     233             :           make_not_null(&temp_momentum_density),
     234             :           make_not_null(&temp_energy_density), char_v_minus, char_v_momentum,
     235             :           char_v_plus, right);
     236             :     } else {
     237             :       temp_mass_density_cons = *mass_density_cons;
     238             :       temp_momentum_density = *momentum_density;
     239             :       temp_energy_density = *energy_density;
     240             :     }
     241             : 
     242             :     // Add to running sum for averaging
     243             :     const double one_over_dim = 1.0 / static_cast<double>(VolumeDim);
     244             :     get(accumulate_mass_density_cons) +=
     245             :         one_over_dim * get(temp_mass_density_cons);
     246             :     for (size_t i = 0; i < VolumeDim; ++i) {
     247             :       accumulate_momentum_density.get(i) +=
     248             :           one_over_dim * temp_momentum_density.get(i);
     249             :     }
     250             :     get(accumulate_energy_density) += one_over_dim * get(temp_energy_density);
     251             :   }  // for loop over dimensions
     252             : 
     253             :   *mass_density_cons = accumulate_mass_density_cons;
     254             :   *momentum_density = accumulate_momentum_density;
     255             :   *energy_density = accumulate_energy_density;
     256             :   return some_component_was_limited;
     257             : }
     258             : 
     259             : }  // namespace NewtonianEuler::Limiters

Generated by: LCOV version 1.14