Add DirichletAnalytic BC to the LDG system SecondOrderScalarWave

src/Evolution/Systems/SecondOrderScalarWave/BoundaryConditions/CMakeLists.txt

@@ -5,11 +5,14 @@ spectre_target_sources(
   ${LIBRARY}
   PRIVATE
   BoundaryCondition.cpp
+  DirichletAnalytic.cpp
   )
 
 spectre_target_headers(
   ${LIBRARY}
   INCLUDE_DIRECTORY ${CMAKE_SOURCE_DIR}/src
   HEADERS
   BoundaryCondition.hpp
+  DirichletAnalytic.hpp
+  Factory.hpp
   )

src/Evolution/Systems/SecondOrderScalarWave/BoundaryConditions/DirichletAnalytic.cpp

@@ -0,0 +1,91 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "Evolution/Systems/SecondOrderScalarWave/BoundaryConditions/DirichletAnalytic.hpp"
+
+#include <cstddef>
+#include <memory>
+#include <pup.h>
+
+#include "DataStructures/TaggedTuple.hpp"
+#include "Evolution/Systems/SecondOrderScalarWave/Tags.hpp"
+#include "PointwiseFunctions/AnalyticSolutions/WaveEquation/Factory.hpp"
+#include "Utilities/CallWithDynamicType.hpp"
+#include "Utilities/GenerateInstantiations.hpp"
+
+namespace SecondOrderScalarWave::BoundaryConditions {
+template <size_t Dim>
+DirichletAnalytic<Dim>::DirichletAnalytic(const DirichletAnalytic& rhs)
+    : BoundaryCondition<Dim>{dynamic_cast<const BoundaryCondition<Dim>&>(rhs)},
+      analytic_prescription_(rhs.analytic_prescription_->get_clone()) {}
+
+template <size_t Dim>
+DirichletAnalytic<Dim>& DirichletAnalytic<Dim>::operator=(
+    const DirichletAnalytic& rhs) {
+  if (&rhs == this) {
+    return *this;
+  }
+  analytic_prescription_ = rhs.analytic_prescription_->get_clone();
+  return *this;
+}
+
+template <size_t Dim>
+DirichletAnalytic<Dim>::DirichletAnalytic(CkMigrateMessage* const msg)
+    : BoundaryCondition<Dim>(msg) {}
+
+template <size_t Dim>
+DirichletAnalytic<Dim>::DirichletAnalytic(
+    std::unique_ptr<evolution::initial_data::InitialData> analytic_prescription)
+    : analytic_prescription_(std::move(analytic_prescription)) {}
+
+template <size_t Dim>
+std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
+DirichletAnalytic<Dim>::get_clone() const {
+  return std::make_unique<DirichletAnalytic>(*this);
+}
+
+template <size_t Dim>
+void DirichletAnalytic<Dim>::pup(PUP::er& p) {
+  BoundaryCondition<Dim>::pup(p);
+  p | analytic_prescription_;
+}
+
+template <size_t Dim>
+std::optional<std::string> DirichletAnalytic<Dim>::dg_ghost(
+    const gsl::not_null<Scalar<DataVector>*> psi,
+    const gsl::not_null<Scalar<DataVector>*> pi,
+    const gsl::not_null<tnsr::i<DataVector, Dim, Frame::Inertial>*> phi,
+    const std::optional<
+        tnsr::I<DataVector, Dim, Frame::Inertial>>& /*face_mesh_velocity*/,
+    const tnsr::i<DataVector, Dim, Frame::Inertial>& /*normal_covector*/,
+    const tnsr::I<DataVector, Dim, Frame::Inertial>& coords,
+    const double time) const {
+  using boundary_value_tags = tmpl::list<Tags::Psi, Tags::Pi, Tags::Phi<Dim>>;
+  auto boundary_values = call_with_dynamic_type<
+      tuples::tagged_tuple_from_typelist<boundary_value_tags>,
+      Solutions::all_solutions<Dim>>(
+      analytic_prescription_.get(),
+      [&coords, &time](const auto* const analytic_solution) {
+        return analytic_solution->variables(coords, time,
+                                            boundary_value_tags{});
+      });
+  *psi = get<Tags::Psi>(boundary_values);
+  *pi = get<Tags::Pi>(boundary_values);
+  *phi = get<Tags::Phi<Dim>>(boundary_values);
+
+  return std::nullopt;
+}
+
+template <size_t Dim>
+// NOLINTNEXTLINE
+PUP::able::PUP_ID DirichletAnalytic<Dim>::my_PUP_ID = 0;
+
+#define DIM(data) BOOST_PP_TUPLE_ELEM(0, data)
+
+#define INSTANTIATION(r, data) template class DirichletAnalytic<DIM(data)>;
+
+GENERATE_INSTANTIATIONS(INSTANTIATION, (1, 2, 3))
+
+#undef INSTANTIATION
+#undef DIM
+}  // namespace SecondOrderScalarWave::BoundaryConditions

src/Evolution/Systems/SecondOrderScalarWave/BoundaryConditions/DirichletAnalytic.hpp

@@ -0,0 +1,101 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#pragma once
+
+#include <cstddef>
+#include <memory>
+#include <optional>
+#include <string>
+
+#include "DataStructures/DataVector.hpp"
+#include "DataStructures/Tensor/Tensor.hpp"
+#include "Evolution/BoundaryConditions/Type.hpp"
+#include "Evolution/Systems/SecondOrderScalarWave/BoundaryConditions/BoundaryCondition.hpp"
+#include "Options/String.hpp"
+#include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp"
+#include "Utilities/Gsl.hpp"
+#include "Utilities/Serialization/CharmPupable.hpp"
+#include "Utilities/TMPL.hpp"
+
+/// \cond
+namespace PUP {
+class er;
+}  // namespace PUP
+namespace Tags {
+struct Time;
+}  // namespace Tags
+namespace domain::Tags {
+template <size_t Dim, typename Frame>
+struct Coordinates;
+}  // namespace domain::Tags
+/// \endcond
+
+namespace SecondOrderScalarWave::BoundaryConditions {
+/*!
+ * \brief Sets Dirichlet boundary conditions from a prescribed analytic
+ * solution.
+ *
+ * Fills the exterior (ghost) values of the evolved variables \f$\Psi\f$ and
+ * \f$\Pi\f$ as well as the LDG auxiliary variable \f$\Phi_i\f$ from the
+ * analytic solution evaluated at the boundary.
+ */
+template <size_t Dim>
+class DirichletAnalytic final : public BoundaryCondition<Dim> {
+ public:
+  /// \brief What analytic solution to prescribe.
+  struct AnalyticPrescription {
+    static constexpr Options::String help =
+        "What analytic solution to prescribe.";
+    using type = std::unique_ptr<evolution::initial_data::InitialData>;
+  };
+
+  using options = tmpl::list<AnalyticPrescription>;
+
+  static constexpr Options::String help{
+      "DirichletAnalytic boundary conditions setting the value of Psi, Pi, and"
+      " Phi to the analytic solution."};
+
+  DirichletAnalytic() = default;
+  DirichletAnalytic(DirichletAnalytic&&) = default;
+  DirichletAnalytic& operator=(DirichletAnalytic&&) = default;
+  DirichletAnalytic(const DirichletAnalytic&);
+  DirichletAnalytic& operator=(const DirichletAnalytic&);
+  ~DirichletAnalytic() override = default;
+
+  explicit DirichletAnalytic(
+      std::unique_ptr<evolution::initial_data::InitialData>
+          analytic_prescription);
+
+  explicit DirichletAnalytic(CkMigrateMessage* msg);
+
+  WRAPPED_PUPable_decl_base_template(
+      domain::BoundaryConditions::BoundaryCondition, DirichletAnalytic);
+
+  auto get_clone() const -> std::unique_ptr<
+      domain::BoundaryConditions::BoundaryCondition> override;
+
+  static constexpr evolution::BoundaryConditions::Type bc_type =
+      evolution::BoundaryConditions::Type::Ghost;
+
+  void pup(PUP::er& p) override;
+
+  using dg_interior_evolved_variables_tags = tmpl::list<>;
+  using dg_interior_temporary_tags =
+      tmpl::list<domain::Tags::Coordinates<Dim, Frame::Inertial>>;
+  using dg_gridless_tags = tmpl::list<::Tags::Time>;
+
+  std::optional<std::string> dg_ghost(
+      gsl::not_null<Scalar<DataVector>*> psi,
+      gsl::not_null<Scalar<DataVector>*> pi,
+      gsl::not_null<tnsr::i<DataVector, Dim, Frame::Inertial>*> phi,
+      const std::optional<
+          tnsr::I<DataVector, Dim, Frame::Inertial>>& /*face_mesh_velocity*/,
+      const tnsr::i<DataVector, Dim, Frame::Inertial>& /*normal_covector*/,
+      const tnsr::I<DataVector, Dim, Frame::Inertial>& coords,
+      double time) const;
+
+ private:
+  std::unique_ptr<evolution::initial_data::InitialData> analytic_prescription_;
+};
+}  // namespace SecondOrderScalarWave::BoundaryConditions

src/Evolution/Systems/SecondOrderScalarWave/BoundaryConditions/Factory.hpp

@@ -0,0 +1,19 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#pragma once
+
+#include <cstddef>
+
+#include "Domain/BoundaryConditions/Periodic.hpp"
+#include "Evolution/Systems/SecondOrderScalarWave/BoundaryConditions/BoundaryCondition.hpp"
+#include "Evolution/Systems/SecondOrderScalarWave/BoundaryConditions/DirichletAnalytic.hpp"
+#include "Utilities/TMPL.hpp"
+
+namespace SecondOrderScalarWave::BoundaryConditions {
+/// Typelist of standard BoundaryConditions
+template <size_t Dim>
+using standard_boundary_conditions =
+    tmpl::list<DirichletAnalytic<Dim>,
+               domain::BoundaryConditions::Periodic<BoundaryCondition<Dim>>>;
+}  // namespace SecondOrderScalarWave::BoundaryConditions

src/Evolution/Systems/SecondOrderScalarWave/CMakeLists.txt

@@ -33,6 +33,8 @@ target_link_libraries(
   DomainBoundaryConditions
   LinearOperators
   Utilities
+  PRIVATE
+  WaveEquationSolutions
   )
 
 add_subdirectory(BoundaryConditions)

tests/Unit/Evolution/Systems/SecondOrderScalarWave/BoundaryConditions/DirichletAnalytic.py

@@ -0,0 +1,69 @@
+# Distributed under the MIT License.
+# See LICENSE.txt for details.
+
+import numpy as np
+
+_amplitude = 0.9
+_width = 0.6
+_profile_center = 0.0
+_center = np.asarray([1.1, 0.1, -0.9])
+_wave_vector = np.asarray([0.1, 1.1, 2.1])
+
+
+def _u(coords, time):
+    dim = len(coords)
+    k = _wave_vector[:dim]
+    omega = np.sqrt(k.dot(k))
+    return k.dot(coords - _center[:dim]) - omega * time
+
+
+def _profile(u):
+    return _amplitude * np.exp(-((u - _profile_center) ** 2) / _width**2)
+
+
+def _dprofile(u):
+    return (
+        (-2.0 * _amplitude / _width**2)
+        * (u - _profile_center)
+        * np.exp(-((u - _profile_center) ** 2) / _width**2)
+    )
+
+
+def error(
+    face_mesh_velocity,
+    outward_directed_normal_covector,
+    coords,
+    time,
+):
+    return None
+
+
+def psi(
+    face_mesh_velocity,
+    outward_directed_normal_covector,
+    coords,
+    time,
+):
+    return _profile(_u(coords, time))
+
+
+def pi(
+    face_mesh_velocity,
+    outward_directed_normal_covector,
+    coords,
+    time,
+):
+    dim = len(coords)
+    k = _wave_vector[:dim]
+    omega = np.sqrt(k.dot(k))
+    return omega * _dprofile(_u(coords, time))
+
+
+def phi(
+    face_mesh_velocity,
+    outward_directed_normal_covector,
+    coords,
+    time,
+):
+    dim = len(coords)
+    return _wave_vector[:dim] * _dprofile(_u(coords, time))