SpECTRE Documentation Coverage Report
Current view: top level - Evolution/Systems/ScalarWave - UpwindPenaltyCorrection.cpp Hit Total Coverage
Commit: b1342d46f40e2d46bbd11d0cef68fd973031a24b Lines: 0 1 0.0 %
Date: 2020-09-24 20:24:42
Legend: Lines: hit not hit

          Line data    Source code
       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #include "Evolution/Systems/ScalarWave/UpwindPenaltyCorrection.hpp"
       5             : 
       6             : #include <algorithm>
       7             : #include <array>
       8             : 
       9             : #include "DataStructures/DataBox/PrefixHelpers.hpp"
      10             : #include "DataStructures/DataVector.hpp"
      11             : #include "DataStructures/Tags/TempTensor.hpp"
      12             : #include "DataStructures/Tensor/Tensor.hpp"
      13             : #include "DataStructures/Variables.hpp"  // IWYU pragma: keep
      14             : #include "Evolution/Systems/ScalarWave/System.hpp"
      15             : #include "Utilities/GenerateInstantiations.hpp"
      16             : #include "Utilities/Gsl.hpp"   // IWYU pragma: keep
      17             : #include "Utilities/TMPL.hpp"  // IWYU pragma: keep
      18             : 
      19             : /// \cond
      20             : namespace ScalarWave {
      21             : template <size_t Dim>
      22             : void UpwindPenaltyCorrection<Dim>::package_data(
      23             :     const gsl::not_null<Scalar<DataVector>*> packaged_char_speed_v_psi,
      24             :     const gsl::not_null<tnsr::i<DataVector, Dim, Frame::Inertial>*>
      25             :         packaged_char_speed_v_zero,
      26             :     const gsl::not_null<Scalar<DataVector>*> packaged_char_speed_v_plus,
      27             :     const gsl::not_null<Scalar<DataVector>*> packaged_char_speed_v_minus,
      28             :     const gsl::not_null<tnsr::i<DataVector, Dim, Frame::Inertial>*>
      29             :         packaged_char_speed_n_times_v_plus,
      30             :     const gsl::not_null<tnsr::i<DataVector, Dim, Frame::Inertial>*>
      31             :         packaged_char_speed_n_times_v_minus,
      32             :     const gsl::not_null<Scalar<DataVector>*> packaged_char_speed_gamma2_v_psi,
      33             :     const gsl::not_null<tnsr::a<DataVector, 3, Frame::Inertial>*>
      34             :         packaged_char_speeds,
      35             : 
      36             :     const Scalar<DataVector>& v_psi,
      37             :     const tnsr::i<DataVector, Dim, Frame::Inertial>& v_zero,
      38             :     const Scalar<DataVector>& v_plus, const Scalar<DataVector>& v_minus,
      39             :     const std::array<DataVector, 4>& char_speeds,
      40             :     const Scalar<DataVector>& constraint_gamma2,
      41             :     const tnsr::i<DataVector, Dim, Frame::Inertial>& interface_unit_normal)
      42             :     const noexcept {
      43             :   // Computes the contribution to the numerical flux from one side of the
      44             :   // interface.
      45             :   //
      46             :   // Note: when PenaltyFlux::operator() is called, an Element passes in its own
      47             :   // packaged data to fill the interior fields, and its neighbor's packaged data
      48             :   // to fill the exterior fields. This introduces a sign flip for each normal
      49             :   // used in computing the exterior fields.
      50             :   get(*packaged_char_speed_v_psi) = char_speeds[0] * get(v_psi);
      51             :   *packaged_char_speed_v_zero = v_zero;
      52             :   for (size_t i = 0; i < Dim; ++i) {
      53             :     packaged_char_speed_v_zero->get(i) *= char_speeds[1];
      54             :   }
      55             :   get(*packaged_char_speed_v_plus) = char_speeds[2] * get(v_plus);
      56             :   get(*packaged_char_speed_v_minus) = char_speeds[3] * get(v_minus);
      57             :   for (size_t d = 0; d < Dim; ++d) {
      58             :     packaged_char_speed_n_times_v_plus->get(d) =
      59             :         get(*packaged_char_speed_v_plus) * interface_unit_normal.get(d);
      60             :     packaged_char_speed_n_times_v_minus->get(d) =
      61             :         get(*packaged_char_speed_v_minus) * interface_unit_normal.get(d);
      62             :   }
      63             :   for (size_t i = 0; i < 4; ++i) {
      64             :     (*packaged_char_speeds)[i] = gsl::at(char_speeds, i);
      65             :   }
      66             :   get(*packaged_char_speed_gamma2_v_psi) =
      67             :       get(constraint_gamma2) * get(*packaged_char_speed_v_psi);
      68             : }
      69             : 
      70             : template <size_t Dim>
      71             : void UpwindPenaltyCorrection<Dim>::operator()(
      72             :     const gsl::not_null<Scalar<DataVector>*> pi_boundary_correction,
      73             :     const gsl::not_null<tnsr::i<DataVector, Dim, Frame::Inertial>*>
      74             :         phi_boundary_correction,
      75             :     const gsl::not_null<Scalar<DataVector>*> psi_boundary_correction,
      76             : 
      77             :     const Scalar<DataVector>& char_speed_v_psi_int,
      78             :     const tnsr::i<DataVector, Dim, Frame::Inertial>& char_speed_v_zero_int,
      79             :     const Scalar<DataVector>& char_speed_v_plus_int,
      80             :     const Scalar<DataVector>& char_speed_v_minus_int,
      81             :     const tnsr::i<DataVector, Dim, Frame::Inertial>&
      82             :         char_speed_normal_times_v_plus_int,
      83             :     const tnsr::i<DataVector, Dim, Frame::Inertial>&
      84             :         char_speed_normal_times_v_minus_int,
      85             :     const Scalar<DataVector>& char_speed_constraint_gamma2_v_psi_int,
      86             :     const tnsr::a<DataVector, 3, Frame::Inertial>& char_speeds_int,
      87             : 
      88             :     const Scalar<DataVector>& char_speed_v_psi_ext,
      89             :     const tnsr::i<DataVector, Dim, Frame::Inertial>& char_speed_v_zero_ext,
      90             :     const Scalar<DataVector>& char_speed_v_plus_ext,
      91             :     const Scalar<DataVector>& char_speed_v_minus_ext,
      92             :     const tnsr::i<DataVector, Dim, Frame::Inertial>&
      93             :         char_speed_minus_normal_times_v_plus_ext,
      94             :     const tnsr::i<DataVector, Dim, Frame::Inertial>&
      95             :         char_speed_minus_normal_times_v_minus_ext,
      96             :     const Scalar<DataVector>& char_speed_constraint_gamma2_v_psi_ext,
      97             :     const tnsr::a<DataVector, 3, Frame::Inertial>& char_speeds_ext)
      98             :     const noexcept {
      99             :   const size_t num_pts = char_speeds_int[0].size();
     100             :   Variables<tmpl::list<::Tags::TempScalar<0>, ::Tags::TempScalar<1>,
     101             :                        ::Tags::TempScalar<2>, ::Tags::TempScalar<3>,
     102             :                        ::Tags::TempScalar<4>, ::Tags::TempScalar<5>,
     103             :                        ::Tags::TempScalar<6>, ::Tags::TempScalar<7>>>
     104             :       buffer(num_pts);
     105             :   DataVector& weighted_lambda_psi_int = get(get<::Tags::TempScalar<0>>(buffer));
     106             :   weighted_lambda_psi_int = step_function(-char_speeds_int[0]);
     107             :   DataVector& weighted_lambda_psi_ext = get(get<::Tags::TempScalar<1>>(buffer));
     108             :   weighted_lambda_psi_ext = -step_function(char_speeds_ext[0]);
     109             : 
     110             :   DataVector& weighted_lambda_zero_int =
     111             :       get(get<::Tags::TempScalar<2>>(buffer));
     112             :   weighted_lambda_zero_int = step_function(-char_speeds_int[1]);
     113             :   DataVector& weighted_lambda_zero_ext =
     114             :       get(get<::Tags::TempScalar<3>>(buffer));
     115             :   weighted_lambda_zero_ext = -step_function(char_speeds_ext[1]);
     116             : 
     117             :   DataVector& weighted_lambda_plus_int =
     118             :       get(get<::Tags::TempScalar<4>>(buffer));
     119             :   weighted_lambda_plus_int = step_function(-char_speeds_int[2]);
     120             :   DataVector& weighted_lambda_plus_ext =
     121             :       get(get<::Tags::TempScalar<5>>(buffer));
     122             :   weighted_lambda_plus_ext = -step_function(char_speeds_ext[2]);
     123             : 
     124             :   DataVector& weighted_lambda_minus_int =
     125             :       get(get<::Tags::TempScalar<6>>(buffer));
     126             :   weighted_lambda_minus_int = step_function(-char_speeds_int[3]);
     127             :   DataVector& weighted_lambda_minus_ext =
     128             :       get(get<::Tags::TempScalar<7>>(buffer));
     129             :   weighted_lambda_minus_ext = -step_function(char_speeds_ext[3]);
     130             : 
     131             :   // D_psi = Theta(-lambda_psi^{ext}) lambda_psi^{ext} v_psi^{ext}
     132             :   //       - Theta(-lambda_psi^{int}) lambda_psi^{int} v_psi^{int}
     133             :   // where the unit normals on both sides point in the same direction, out
     134             :   // of the current element. Since lambda_psi from the neighbor is computing
     135             :   // with the normal vector pointing into the current element in the code,
     136             :   // we need to swap the sign of lambda_psi^{ext}. Theta is the heaviside step
     137             :   // function with Theta(0) = 0.
     138             :   psi_boundary_correction->get() =
     139             :       weighted_lambda_psi_ext * get(char_speed_v_psi_ext) -
     140             :       weighted_lambda_psi_int * get(char_speed_v_psi_int);
     141             : 
     142             :   get(*pi_boundary_correction) =
     143             :       0.5 * (weighted_lambda_plus_ext * get(char_speed_v_plus_ext) +
     144             :              weighted_lambda_minus_ext * get(char_speed_v_minus_ext)) +
     145             :       weighted_lambda_psi_ext * get(char_speed_constraint_gamma2_v_psi_ext)
     146             : 
     147             :       - 0.5 * (weighted_lambda_plus_int * get(char_speed_v_plus_int) +
     148             :                weighted_lambda_minus_int * get(char_speed_v_minus_int)) -
     149             :       weighted_lambda_psi_int * get(char_speed_constraint_gamma2_v_psi_int);
     150             : 
     151             :   for (size_t d = 0; d < Dim; ++d) {
     152             :     // Overall minus sign on ext because of normal vector is opposite direction.
     153             :     phi_boundary_correction->get(d) =
     154             :         -0.5 * (weighted_lambda_minus_ext *
     155             :                     char_speed_minus_normal_times_v_minus_ext.get(d) -
     156             :                 weighted_lambda_plus_ext *
     157             :                     char_speed_minus_normal_times_v_plus_ext.get(d)) +
     158             :         weighted_lambda_zero_ext * char_speed_v_zero_ext.get(d)
     159             : 
     160             :         - 0.5 * (weighted_lambda_plus_int *
     161             :                      char_speed_normal_times_v_plus_int.get(d) -
     162             :                  weighted_lambda_minus_int *
     163             :                      char_speed_normal_times_v_minus_int.get(d)) -
     164             :         weighted_lambda_zero_int * char_speed_v_zero_int.get(d);
     165             :   }
     166             : }
     167             : 
     168             : #define DIM(data) BOOST_PP_TUPLE_ELEM(0, data)
     169             : 
     170             : #define INSTANTIATION(_, data) \
     171             :   template class UpwindPenaltyCorrection<DIM(data)>;
     172             : 
     173             : GENERATE_INSTANTIATIONS(INSTANTIATION, (1, 2, 3))
     174             : 
     175             : #undef INSTANTIATION
     176             : #undef DIM
     177             : }  // namespace ScalarWave
     178             : /// \endcond

Generated by: LCOV version 1.14