Add python test function for boundary corrections #2700

nilsdeppe · 2020-12-21T22:54:31Z

Proposed changes

Add a generic function that does a random value test of the boundary corrections for DG. This means that the actual tests for boundary corrections are quite short. For example, the Rusanov solver for Burgers is:

// Distributed under the MIT License.
// See LICENSE.txt for details.

#include "Framework/TestingFramework.hpp"

#include <cstddef>
#include <array>
#include <string>

#include "Evolution/Systems/Burgers/BoundaryCorrections/Factory.hpp"
#include "Evolution/Systems/Burgers/BoundaryCorrections/Rusanov.hpp"
#include "Evolution/Systems/Burgers/System.hpp"
#include "Framework/SetupLocalPythonEnvironment.hpp"
#include "Framework/TestCreation.hpp"
#include "Helpers/Evolution/DiscontinuousGalerkin/BoundaryCorrections.hpp"
#include "NumericalAlgorithms/Spectral/Mesh.hpp"
#include "NumericalAlgorithms/Spectral/Spectral.hpp"

SPECTRE_TEST_CASE("Unit.Burgers.Rusanov", "[Unit][Burgers]") {
  pypp::SetupLocalPythonEnvironment local_python_env{
      "Evolution/Systems/Burgers/BoundaryCorrections"};

  TestHelpers::evolution::dg::test_boundary_correction_conservation<
      Burgers::System>(
      Burgers::BoundaryCorrections::Rusanov{},
      Mesh<0>{1, Spectral::Basis::Legendre, Spectral::Quadrature::Gauss}, {});

  TestHelpers::evolution::dg::test_boundary_correction_with_python<
      Burgers::System>(
      "Rusanov",
      {{"dg_package_data_u", "dg_package_data_normal_dot_flux",
        "dg_package_data_abs_char_speed"}},
      {{"dg_boundary_terms_u"}}, Burgers::BoundaryCorrections::Rusanov{},
      Mesh<0>{1, Spectral::Basis::Legendre, Spectral::Quadrature::Gauss}, {});

  const auto rusanov = TestHelpers::test_factory_creation<
      Burgers::BoundaryCorrections::BoundaryCorrection>("Rusanov:");

  TestHelpers::evolution::dg::test_boundary_correction_with_python<
      Burgers::System>(
      "Rusanov",
      {{"dg_package_data_u", "dg_package_data_normal_dot_flux",
        "dg_package_data_abs_char_speed"}},
      {{"dg_boundary_terms_u"}},
      dynamic_cast<const Burgers::BoundaryCorrections::Rusanov&>(*rusanov),
      Mesh<0>{1, Spectral::Basis::Legendre, Spectral::Quadrature::Gauss}, {});
}

Most work/code comes in from the python implementations. There's a bit more needed if an equation of state (or some other non-DataVector-wrapping class) needs to be converted. Here's Newtonian Euler as an example:

// Distributed under the MIT License.
// See LICENSE.txt for details.

#include "Framework/TestingFramework.hpp"

#include <cstddef>
#include <array>
#include <string>

#include "Evolution/Systems/NewtonianEuler/BoundaryCorrections/Factory.hpp"
#include "Evolution/Systems/NewtonianEuler/BoundaryCorrections/Rusanov.hpp"
#include "Evolution/Systems/NewtonianEuler/System.hpp"
#include "Framework/SetupLocalPythonEnvironment.hpp"
#include "Framework/TestCreation.hpp"
#include "Helpers/Evolution/DiscontinuousGalerkin/BoundaryCorrections.hpp"
#include "NumericalAlgorithms/Spectral/Mesh.hpp"
#include "NumericalAlgorithms/Spectral/Spectral.hpp"
#include "PointwiseFunctions/Hydro/EquationsOfState/PolytropicFluid.hpp"
#include "PointwiseFunctions/Hydro/EquationsOfState/IdealFluid.hpp"
#include "PointwiseFunctions/Hydro/Tags.hpp"
#include "Utilities/TMPL.hpp"
#include "Utilities/TaggedTuple.hpp"

namespace {
struct ConvertPolytropic {
  using unpacked_container = bool;
  using packed_container = EquationsOfState::PolytropicFluid<false>;
  using packed_type = bool;

  static inline unpacked_container unpack(
      const packed_container& /*packed*/,
      const size_t /*grid_point_index*/) noexcept {
    return true;
  }

  [[noreturn]] static inline void pack(
      const gsl::not_null<packed_container*> /*packed*/,
      const unpacked_container& /*unpacked*/,
      const size_t /*grid_point_index*/) {
    ERROR("Should not be converting an EOS from an unpacked to a packed type");
  }

  static inline size_t get_size(const packed_container& /*packed*/) noexcept {
    return 1;
  }
};

struct ConvertIdeal {
  using unpacked_container = bool;
  using packed_container = EquationsOfState::IdealFluid<false>;
  using packed_type = bool;

  static inline unpacked_container unpack(
      const packed_container& /*packed*/,
      const size_t /*grid_point_index*/) noexcept {
    return false;
  }

  [[noreturn]] static inline void pack(
      const gsl::not_null<packed_container*> /*packed*/,
      const unpacked_container& /*unpacked*/,
      const size_t /*grid_point_index*/) {
    ERROR("Should not be converting an EOS from an unpacked to a packed type");
  }

  static inline size_t get_size(const packed_container& /*packed*/) noexcept {
    return 1;
  }
};

struct DummyInitialData {
  using argument_tags = tmpl::list<>;
  struct source_term_type {
    using sourced_variables = tmpl::list<>;
    using argument_tags = tmpl::list<>;
  };
};

template <size_t Dim, typename EosTag>
void test(const size_t num_pts,
          const tuples::TaggedTuple<EosTag>& volume_data) {
  TestHelpers::evolution::dg::test_boundary_correction_conservation<
      NewtonianEuler::System<Dim, typename EosTag::type, DummyInitialData>>(
      NewtonianEuler::BoundaryCorrections::Rusanov<Dim>{},
      Mesh<Dim - 1>{num_pts, Spectral::Basis::Legendre,
                    Spectral::Quadrature::Gauss},
      volume_data);

  TestHelpers::evolution::dg::test_boundary_correction_with_python<
      NewtonianEuler::System<Dim, typename EosTag::type, DummyInitialData>,
      tmpl::list<ConvertPolytropic, ConvertIdeal>>(
      "Rusanov",
      {{"dg_package_data_mass_density", "dg_package_data_momentum_density",
        "dg_package_data_energy_density",
        "dg_package_data_normal_dot_flux_mass_density",
        "dg_package_data_normal_dot_flux_momentum_density",
        "dg_package_data_normal_dot_flux_energy_density",
        "dg_package_data_abs_char_speed"}},
      {{"dg_boundary_terms_mass_density", "dg_boundary_terms_momentum_density",
        "dg_boundary_terms_energy_density"}},
      NewtonianEuler::BoundaryCorrections::Rusanov<Dim>{},
      Mesh<Dim - 1>{num_pts, Spectral::Basis::Legendre,
                    Spectral::Quadrature::Gauss},
      volume_data);

  const auto rusanov = TestHelpers::test_factory_creation<
      NewtonianEuler::BoundaryCorrections::BoundaryCorrection<Dim>>("Rusanov:");

  TestHelpers::evolution::dg::test_boundary_correction_with_python<
      NewtonianEuler::System<Dim, typename EosTag::type, DummyInitialData>,
      tmpl::list<ConvertPolytropic, ConvertIdeal>>(
      "Rusanov",
      {{"dg_package_data_mass_density", "dg_package_data_momentum_density",
        "dg_package_data_energy_density",
        "dg_package_data_normal_dot_flux_mass_density",
        "dg_package_data_normal_dot_flux_momentum_density",
        "dg_package_data_normal_dot_flux_energy_density",
        "dg_package_data_abs_char_speed"}},
      {{"dg_boundary_terms_mass_density", "dg_boundary_terms_momentum_density",
        "dg_boundary_terms_energy_density"}},
      dynamic_cast<const NewtonianEuler::BoundaryCorrections::Rusanov<Dim>&>(
          *rusanov),
      Mesh<Dim - 1>{num_pts, Spectral::Basis::Legendre,
                    Spectral::Quadrature::Gauss},
      volume_data);
}
}  // namespace

SPECTRE_TEST_CASE("Unit.NewtonianEuler.Rusanov", "[Unit][Evolution]") {
  pypp::SetupLocalPythonEnvironment local_python_env{
      "Evolution/Systems/NewtonianEuler/BoundaryCorrections"};
  test<1>(1, tuples::TaggedTuple<hydro::Tags::EquationOfState<
                 EquationsOfState::PolytropicFluid<false>>>{
                 EquationsOfState::PolytropicFluid<false>{1.0e-3, 2.0}});
  test<2>(5, tuples::TaggedTuple<hydro::Tags::EquationOfState<
                 EquationsOfState::PolytropicFluid<false>>>{
                 EquationsOfState::PolytropicFluid<false>{1.0e-3, 2.0}});
  test<3>(5, tuples::TaggedTuple<hydro::Tags::EquationOfState<
                 EquationsOfState::PolytropicFluid<false>>>{
                 EquationsOfState::PolytropicFluid<false>{1.0e-3, 2.0}});

  test<1>(
      1, tuples::TaggedTuple<
             hydro::Tags::EquationOfState<EquationsOfState::IdealFluid<false>>>{
             EquationsOfState::IdealFluid<false>{1.3}});
  test<2>(
      5, tuples::TaggedTuple<
             hydro::Tags::EquationOfState<EquationsOfState::IdealFluid<false>>>{
             EquationsOfState::IdealFluid<false>{1.3}});
  test<3>(
      5, tuples::TaggedTuple<
             hydro::Tags::EquationOfState<EquationsOfState::IdealFluid<false>>>{
             EquationsOfState::IdealFluid<false>{1.3}});
}

Upgrade instructions

Code review checklist

The code is documented and the documentation renders correctly. Run
make doc to generate the documentation locally into BUILD_DIR/docs/html.
Then open index.html.
The code follows the stylistic and code quality guidelines listed in the
code review guide.
The PR lists upgrade instructions and is labeled bugfix or
major new feature if appropriate.

Further comments

moxcodes · 2020-12-22T20:12:00Z

tests/Unit/Evolution/DiscontinuousGalerkin/Test_BoundaryCorrectionsHelper.cpp


 namespace {
+struct VolumeDouble {
+  double t;


perhaps value as the member name?

moxcodes · 2020-12-22T20:16:51Z

tests/Unit/Evolution/DiscontinuousGalerkin/Test_BoundaryCorrectionsHelper.cpp

      /*mesh_velocity*/,
      const boost::optional<Scalar<DataVector>>& normal_dot_mesh_velocity,
-      const double volume_double) const noexcept {
+      const VolumeDoubleType volume_double_in) const noexcept {


I'm sure this will be clearer in the documentation for the specific boundary corrections when you start putting them in, but for my understanding, can you explain the role of this volume_double in these packaging functions?

It's just to test that quantities from the volume (i.e. DataBox) can be passed in. A practical example is needing the equation of state (which can't be projected onto the boundary) in order to calculate the speed of sound, which is needed to compute the characteristic speeds. The Newtonian Euler system uses this pattern, for example.

moxcodes · 2020-12-22T20:21:29Z

tests/Unit/Evolution/DiscontinuousGalerkin/Test_BoundaryCorrectionsHelper.cpp

-template <size_t Dim>
-struct Correction {
+template <size_t Dim, typename VolumeDoubleType>
+struct Correction final {


This Correction object is just for the test helper, but presumably its interface will be the same as the forthcoming individual boundary correction structs -- will the boundary corrections need to store member data? If not, can their functions be made static?

Some boundary corrections (not ones we are using) exist that have parameters that would be reasonable to tune in an input file. I think it's easier to have all the functions be non-static for uniformity and ease-of-understanding for new users. There are also classes of boundary corrections like HLL, HLLE, HLLEM, HLLC, etc. (there's something like 5-8 of these that I know of, but I'm sure there are even more) that could reasonably all be one class with an input file option for which one to use. This helps keep down the amount of boiler plate code (class, inheritence, function declarations, files, etc.) to support many different forms of nearly identical boundary corrects.

moxcodes · 2020-12-22T20:38:56Z

tests/Unit/Helpers/Evolution/DiscontinuousGalerkin/BoundaryCorrections.hpp

+          CAPTURE(python_module);
+          CAPTURE(
+              gsl::at(python_dg_package_data_functions, function_name_index));
+          const auto cxx_result = pypp::call<ResultType, ConversionClassList>(


seems slightly misleading to call this cxx_result as it's the thing that results from calling the python functions and then converting them to c++ objects.

Agreed, should be named python_result!

moxcodes · 2020-12-22T20:49:38Z

tests/Unit/Helpers/Evolution/DiscontinuousGalerkin/BoundaryCorrections.hpp

+          const auto cxx_result = pypp::call<ResultType, ConversionClassList>(
+              python_module,
+              gsl::at(python_dg_package_data_functions, function_name_index),
+              get<FaceTags>(fields_on_face)..., unit_normal_covector,


The interface to the python function seems potentially confusing -- the ordering of the arguments is not quite what I would expect (same order, but variables pack-expanded).

I'm not sure what you mean. In the C++ code the Variables classes are also expanded before being passed in.

Sorry, I think I was just mistaken with my previous comment -- what I was saying was that I would expect the same order but pack-expanded, but thought that you were passing in a different order. I think I just misunderstood what FaceTags would hold, so I'm now happy with this as-is

moxcodes · 2020-12-22T20:50:37Z

tests/Unit/Helpers/Evolution/DiscontinuousGalerkin/BoundaryCorrections.hpp

+ * You can convert custom types using the `ConversionClassList` template
+ * parameter, which is then passed to `pypp::call()`. This allows you to, e.g.,
+ * convert an equation of state into an array locally in a test file.
+ */


It would be helpful to add notes about the assumed relationship between python functions and the package functions tag list, and the System:variables_tag tag list - In particular, you require that there be a python function for each packaged field, and one for each tag in the system variables. The interface of the python function should also be specified -- the order of arguments and relationships to the tag collections where the Variables need to be pack-expanded.

nilsdeppe · 2021-01-01T22:55:47Z

I've pushed a fixup and rebased on develop. Please let me know if I've missed or misunderstood any comments (or anything else!).

fmahebert

Once @moxcodes gives okay to squash, you can also squash these small things in

fmahebert · 2021-01-04T22:24:22Z

tests/Unit/Helpers/Evolution/DiscontinuousGalerkin/BoundaryCorrections.hpp

+ *   in `dg_package_field_tags`
+ * - There must be one python function to compute the boundary correction for
+ *   each tag in `System::variables_tag`
+ * - The arguments to the python functions for computing the packaged data is


fmahebert · 2021-01-04T22:38:29Z

tests/Unit/Helpers/Evolution/DiscontinuousGalerkin/BoundaryCorrections.hpp

+          make_not_null(&gen), make_not_null(&dist), used_for_size);
+  const auto unit_normal_covector = make_with_random_values<
+      tnsr::i<DataVector, FaceDim + 1, Frame::Inertial>>(
+      make_not_null(&gen), make_not_null(&dist), used_for_size);


You call this a unit normal, but as far as I can tell it has a random magnitude (including, rarely, 0 — should this be checked against too?). Please either normalize the vector, rename the variable, or add a comment.

moxcodes · 2021-01-07T22:23:22Z

looks good -- go ahead and squash.

nilsdeppe · 2021-01-07T23:59:44Z

Rebased on develop and squashed, including @fmahebert's requested changes. (I decided to just normalize the vector)

nilsdeppe · 2021-01-08T00:22:00Z

Rebased again and squashed in boost::optional -> std::optional changes in the moving mesh code

nilsdeppe requested review from fmahebert and moxcodes December 21, 2020 22:54

moxcodes reviewed Dec 22, 2020

View reviewed changes

nilsdeppe force-pushed the boundary_correction_random branch from 4c194d3 to 60035cc Compare January 1, 2021 22:55

fmahebert reviewed Jan 4, 2021

View reviewed changes

nilsdeppe force-pushed the boundary_correction_random branch from 60035cc to 3db3111 Compare January 7, 2021 23:58

Add python test function for boundary corrections

bd58e69

nilsdeppe force-pushed the boundary_correction_random branch from 3db3111 to bd58e69 Compare January 8, 2021 00:21

fmahebert approved these changes Jan 8, 2021

View reviewed changes

moxcodes approved these changes Jan 8, 2021

View reviewed changes

kidder approved these changes Jan 10, 2021

View reviewed changes

kidder merged commit 169ff89 into sxs-collaboration:develop Jan 10, 2021

nilsdeppe deleted the boundary_correction_random branch January 12, 2021 04:36

Add python test function for boundary corrections #2700

Add python test function for boundary corrections #2700

Uh oh!

Conversation

nilsdeppe commented Dec 21, 2020

Proposed changes

Upgrade instructions

Code review checklist

Further comments

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

nilsdeppe commented Jan 1, 2021

Uh oh!

fmahebert left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

moxcodes commented Jan 7, 2021

Uh oh!

nilsdeppe commented Jan 7, 2021

Uh oh!

nilsdeppe commented Jan 8, 2021

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants