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