Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <cstddef>
       7             : #include <optional>
       8             : #include <utility>
       9             : 
      10             : #include "DataStructures/DataBox/PrefixHelpers.hpp"
      11             : #include "DataStructures/DataBox/Prefixes.hpp"
      12             : #include "DataStructures/Variables.hpp"
      13             : #include "Domain/TagsTimeDependent.hpp"
      14             : #include "Evolution/DgSubcell/Projection.hpp"
      15             : #include "Evolution/DgSubcell/SubcellOptions.hpp"
      16             : #include "Evolution/DgSubcell/Tags/CellCenteredFlux.hpp"
      17             : #include "Evolution/DgSubcell/Tags/DidRollback.hpp"
      18             : #include "Evolution/DgSubcell/Tags/Mesh.hpp"
      19             : #include "Evolution/DgSubcell/Tags/SubcellOptions.hpp"
      20             : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
      21             : #include "NumericalAlgorithms/Spectral/Parity.hpp"
      22             : #include "Utilities/Gsl.hpp"
      23             : #include "Utilities/TMPL.hpp"
      24             : 
      25             : namespace evolution::dg::subcell::fd {
      26             : /*!
      27             :  * \brief Mutator that wraps the system's `FluxMutator` to correctly set the
      28             :  * cell-centered fluxes on the subcell grid.
      29             :  *
      30             :  * Currently we only use high-order FD if the FD order was specified in the
      31             :  * input file. We will need to extend this to support adaptive-order in the
      32             :  * future. In that case we need to check if the FD reconstruction reports back
      33             :  * the order to use.
      34             :  */
      35             : template <typename System, typename FluxMutator, size_t Dim,
      36             :           bool ComputeOnlyOnRollback, typename Fr = Frame::Inertial>
      37           1 : struct CellCenteredFlux {
      38           0 :   using flux_variables = typename System::flux_variables;
      39           0 :   using variables = typename System::variables_tag::tags_list;
      40             : 
      41           0 :   using return_tags =
      42             :       tmpl::list<subcell::Tags::CellCenteredFlux<flux_variables, Dim>>;
      43           0 :   using argument_tags = tmpl::push_front<
      44             :       typename FluxMutator::argument_tags, subcell::Tags::SubcellOptions<Dim>,
      45             :       subcell::Tags::Mesh<Dim>, domain::Tags::Mesh<Dim>,
      46             :       domain::Tags::MeshVelocity<Dim, Frame::Inertial>,
      47             :       ::Tags::Variables<variables>, subcell::Tags::DidRollback>;
      48             : 
      49             :   template <typename... FluxTags, typename... Args>
      50           0 :   static void apply(
      51             :       const gsl::not_null<std::optional<Variables<tmpl::list<FluxTags...>>>*>
      52             :           cell_centered_fluxes,
      53             :       const subcell::SubcellOptions& subcell_options,
      54             :       const Mesh<Dim>& subcell_mesh, const Mesh<Dim>& dg_mesh,
      55             :       const std::optional<tnsr::I<DataVector, Dim>>& dg_mesh_velocity,
      56             :       const ::Variables<variables>& cell_centered_flux_vars,
      57             :       const bool did_rollback, Args&&... args) {
      58             :     if (did_rollback or not ComputeOnlyOnRollback) {
      59             :       if (subcell_options.finite_difference_derivative_order() !=
      60             :           ::fd::DerivativeOrder::Two) {
      61             :         if (not cell_centered_fluxes->has_value()) {
      62             :           (*cell_centered_fluxes) =
      63             :               Variables<db::wrap_tags_in<::Tags::Flux, flux_variables,
      64             :                                          tmpl::size_t<Dim>, Fr>>{
      65             :                   subcell_mesh.number_of_grid_points()};
      66             :         }
      67             :         FluxMutator::apply(
      68             :             make_not_null(&get<FluxTags>((*cell_centered_fluxes).value()))...,
      69             :             std::forward<Args>(args)...);
      70             :         if (dg_mesh_velocity.has_value()) {
      71             :           for (size_t i = 0; i < Dim; i++) {
      72             :             //
      73             :             // Project mesh velocity on face mesh. We only need the component
      74             :             // orthogonal to the face.
      75             :             const DataVector& cell_centered_mesh_velocity =
      76             :                 evolution::dg::subcell::fd::project(
      77             :                     dg_mesh_velocity.value().get(i), dg_mesh,
      78             :                     subcell_mesh.extents(),
      79             :                     i == 0 ? Spectral::Parity::Odd : Spectral::Parity::Even);
      80             : 
      81             :             tmpl::for_each<flux_variables>(
      82             :                 [&cell_centered_flux_vars, &cell_centered_mesh_velocity,
      83             :                  &cell_centered_fluxes, &i](auto tag_v) {
      84             :                   using tag = tmpl::type_from<decltype(tag_v)>;
      85             :                   using flux_tag =
      86             :                       ::Tags::Flux<tag, tmpl::size_t<Dim>, Frame::Inertial>;
      87             :                   using FluxTensor = typename flux_tag::type;
      88             :                   const auto& var = get<tag>(cell_centered_flux_vars);
      89             :                   auto& flux = get<flux_tag>(cell_centered_fluxes->value());
      90             :                   for (size_t storage_index = 0; storage_index < var.size();
      91             :                        ++storage_index) {
      92             :                     const auto tensor_index =
      93             :                         var.get_tensor_index(storage_index);
      94             :                     const auto flux_storage_index =
      95             :                         FluxTensor::get_storage_index(prepend(tensor_index, i));
      96             :                     flux[flux_storage_index] -=
      97             :                         cell_centered_mesh_velocity * var[storage_index];
      98             :                   }
      99             :                 });
     100             :           }
     101             :         }
     102             :       }
     103             :     }
     104             :   }
     105             : };
     106             : }  // namespace evolution::dg::subcell::fd