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