model.hh

Go to the documentation of this file.

#pragma once
 
// external libraries
#include "DiFfRG/common/config_tree.hh"
#include <deal.II/base/point.h>
#include <deal.II/base/tensor.h>
 
// standard library
#include <cmath>
#include <limits>
#include <optional>
#include <vector>
 
// DiFfRG
#include <DiFfRG/discretization/common/affine_constraint_metadata.hh>
#include <DiFfRG/discretization/common/solution_sample.hh>
#include <DiFfRG/model/ad.hh>
#include <DiFfRG/model/component_descriptor.hh>
#include <DiFfRG/model/fv_boundaries.hh>
#include <DiFfRG/model/numflux.hh>
 
namespace DiFfRG
{
  using namespace dealii;
 
  namespace def
  {
    inline constexpr int no_wave_speed = -1;
 
    template <typename Model, typename Components_> class AbstractModel
    {
      Model &asImp()
      {
        static_assert(
            std::is_base_of_v<AbstractModel<Model, Components_>, Model>,
            "AbstractModel<Model, Components>: Model must inherit from AbstractModel<Model, Components> (CRTP). "
            "Check that your model class passes itself as the first template argument.");
        return static_cast<Model &>(*this);
      }
      const Model &asImp() const
      {
        static_assert(
            std::is_base_of_v<AbstractModel<Model, Components_>, Model>,
            "AbstractModel<Model, Components>: Model must inherit from AbstractModel<Model, Components> (CRTP). "
            "Check that your model class passes itself as the first template argument.");
        return static_cast<const Model &>(*this);
      }
 
    protected:
      Components_ m_components;
      auto &components() { return m_components; }
 
    public:
      const auto &get_components() const { return m_components; }
      using Components = Components_;
 
      template <int dim, typename Vector> void initial_condition(const Point<dim> &x, Vector &u_i) const = delete;
 
      template <int dim, typename NumberType, typename Vector, typename Vector_dot, size_t n_fe_functions>
      void mass([[maybe_unused]] std::array<NumberType, n_fe_functions> &m_i, [[maybe_unused]] const Point<dim> &x,
                [[maybe_unused]] const Vector &u_i, const Vector_dot &dt_u_i) const
      {
        for (uint i = 0; i < n_fe_functions; ++i)
          m_i[i] = dt_u_i[i];
      }
 
      template <int dim, typename NumberType, size_t n_fe_functions>
      void mass(std::array<std::array<NumberType, n_fe_functions>, n_fe_functions> &m_ij,
                [[maybe_unused]] const Point<dim> &x) const
      {
        for (uint i = 0; i < n_fe_functions; ++i)
          for (uint j = 0; j < n_fe_functions; ++j)
            m_ij[i][j] = 0.;
        for (uint i = 0; i < n_fe_functions; ++i)
          m_ij[i][i] = 1.;
      }
 
      template <int dim, typename NumberType, typename Solutions, size_t n_fe_functions>
      void flux([[maybe_unused]] std::array<Tensor<1, dim, NumberType>, n_fe_functions> &F_i,
                [[maybe_unused]] const Point<dim> &x, [[maybe_unused]] const Solutions &sol) const
      {
      }
 
      template <int dim, typename NumberType, typename Solutions, size_t n_fe_functions>
      void diffusion_flux([[maybe_unused]] std::array<Tensor<1, dim, NumberType>, n_fe_functions> &F_i,
                          [[maybe_unused]] const Point<dim> &x, [[maybe_unused]] const Solutions &sol) const
      {
      }
 
      template <int dim, typename NumberType, typename Solutions, size_t n_fe_functions>
      void source([[maybe_unused]] std::array<NumberType, n_fe_functions> &s_i, [[maybe_unused]] const Point<dim> &x,
                  [[maybe_unused]] const Solutions &sol) const
      {
      }
 
      template <uint dim> std::vector<bool> differential_components() const
      {
        std::vector<bool> differential_components(Model::Components::count_fe_functions(), false);
 
        // First we need two reference solutions u_i and dt_u_i, which we then both fill with 1.s
        std::array<double, Model::Components::count_fe_functions()> u_i{{}};
        std::array<double, Model::Components::count_fe_functions()> dt_u_i{{}};
        for (uint i = 0; i < Model::Components::count_fe_functions(); ++i) {
          u_i[i] = 1.;
          dt_u_i[i] = 1.;
        }
        // Set the point to be at 1. in all directions
        Point<dim> x;
        for (uint i = 0; i < dim; ++i)
          x[i] = 1.;
        // Get the mass function m_i
        std::array<double, Model::Components::count_fe_functions()> m_i{{}};
        asImp().mass(m_i, x, u_i, dt_u_i);
 
        // Now we check which components are differential by changing dt_u_i slightly and checking whether the mass
        // function changes.
        for (uint i = 0; i < Model::Components::count_fe_functions(); ++i) {
          dt_u_i[i] = 1. + 1e-1;
          std::array<double, Model::Components::count_fe_functions()> m_i_new{{}};
          asImp().mass(m_i_new, x, u_i, dt_u_i);
          dt_u_i[i] = 1.;
          for (uint j = 0; j < Model::Components::count_fe_functions(); ++j)
            if (!is_close(m_i[j], m_i_new[j])) differential_components[j] = true;
        }
 
        return differential_components;
      }
 
      template <size_t n_fe_functions> void wave_speed_blocks(std::array<int, n_fe_functions> &blocks) const
      {
        blocks.fill(0);
      }
 
 
 
      template <typename Vector> void initial_condition_variables([[maybe_unused]] Vector &v_a) const
      {
        // Just to avoid warnings
      }
 
      template <typename Vector, typename Solution>
      void dt_variables([[maybe_unused]] Vector &r_a, [[maybe_unused]] const Solution &sol) const
      {
        // Just to avoid warnings
      }
 
 
 
      template <int dim, typename Vector, typename Solutions>
      void extract([[maybe_unused]] Vector &result, [[maybe_unused]] const Point<dim> &x,
                   [[maybe_unused]] const Solutions &sol) const
      {
      }
 
 
 
      template <uint dependent, int dim, typename NumberType, typename Vector, size_t n_fe_functions_dep>
      void ldg_flux([[maybe_unused]] std::array<Tensor<1, dim, NumberType>, n_fe_functions_dep> &F,
                    [[maybe_unused]] const Point<dim> &x, [[maybe_unused]] const Vector &u) const
      {
      }
 
      template <uint dependent, int dim, typename NumberType, typename Vector, size_t n_fe_functions_dep>
      void ldg_source([[maybe_unused]] std::array<NumberType, n_fe_functions_dep> &s,
                      [[maybe_unused]] const Point<dim> &x, [[maybe_unused]] const Vector &u) const
      {
      }
 
      template <int dim, typename NumberType, typename Solutions_s, typename Solutions_n>
      void face_indicator([[maybe_unused]] std::array<NumberType, 2> &indicator,
                          [[maybe_unused]] const Tensor<1, dim> &normal, [[maybe_unused]] const Point<dim> &p,
                          [[maybe_unused]] const Solutions_s &sol_s, [[maybe_unused]] const Solutions_n &sol_n) const
      {
      }
 
      template <int dim, typename NumberType, typename Solution>
      void cell_indicator([[maybe_unused]] NumberType &indicator, [[maybe_unused]] const Point<dim> &p,
                          [[maybe_unused]] const Solution &sol) const
      {
      }
 
      template <int dim, typename Vector>
      std::array<double, dim> EoM([[maybe_unused]] const Point<dim> &x, const Vector &u) const
      {
        return std::array<double, dim>{{u[0]}};
      }
 
      static constexpr bool extract_uses_potential = true;
 
      template <int dim, typename Vector>
      std::array<double, dim> raw_potential_gradient([[maybe_unused]] const Point<dim> &x, const Vector &u) const
      {
        std::array<double, dim> gradient{};
        for (uint d = 0; d < dim && d < u.size(); ++d)
          gradient[d] = u[d];
        return gradient;
      }
 
      template <int dim, typename Vector> Point<dim> EoM_postprocess(const Point<dim> &EoM, const Vector &) const
      {
        return EoM;
      }
 
      template <typename FUN, typename DataOut> void readouts_multiple(FUN &helper, DataOut &) const
      {
        helper(
            "primary", [&](const auto &x, const auto &u_i) { return asImp().EoM(x, u_i); }, // chiral EoM
            [&](auto &output, const auto &x, const auto &sol) { asImp().readouts(output, x, sol); });
      }
 
      template <int dim, typename DataOut, typename Solutions>
      void readouts([[maybe_unused]] DataOut &output, [[maybe_unused]] const Point<dim> &x,
                    [[maybe_unused]] const Solutions &sol) const
      {
      }
 
      template <typename Constraints, typename Context>
      void affine_constraints(Constraints &constraints, const Context &context) const
      {
        if constexpr (requires(const Model &model, Constraints &constraint_matrix,
                               const Context &affine_constraint_context) {
                        model.apply_boundary_affine_constraints(constraint_matrix, affine_constraint_context);
                      })
          asImp().apply_boundary_affine_constraints(constraints, context);
        if constexpr (requires(const Model &model, Constraints &constraint_matrix,
                               const Context &affine_constraint_context) {
                        model.apply_affine_constraints(constraint_matrix, affine_constraint_context);
                      })
          asImp().apply_affine_constraints(constraints, context);
      }
    };
 
    namespace internal
    {
      template <typename> inline constexpr bool dependent_false_v = false;
 
      template <FixedString component_name, typename Model, int dim>
      double origin_constraint_coordinate(const Point<dim> &point)
      {
        if constexpr (dim == 1) {
          return point[0];
        } else {
          if constexpr (requires {
                          Model::template OriginConstraintCoordinate<component_name>::signed_coordinate(point);
                        }) {
            return Model::template OriginConstraintCoordinate<component_name>::signed_coordinate(point);
          } else {
            static_assert(dependent_false_v<Model>,
                          "Multidimensional origin affine-constraint helpers require Model::"
                          "OriginConstraintCoordinate<component_name>::signed_coordinate(point).");
            return 0.0;
          }
        }
      }
 
      template <FixedString component_name, typename Model, typename Context>
      std::vector<types::global_dof_index> select_origin_candidates([[maybe_unused]] const Context &context,
                                                                    const auto &view)
      {
        constexpr int dim = Context::dimension;
        static_assert(dim == 1 || dim == 2,
                      "Origin affine-constraint helpers currently support only one- and two-dimensional domains.");
        static_assert(Context::template component_size<component_name>() == 1,
                      "Origin affine-constraint helpers require a scalar FE-function component.");
 
        double best_abs_coordinate = std::numeric_limits<double>::infinity();
        bool has_non_negative_best = false;
 
        for (uint i = 0; i < view.dofs.n_elements(); ++i) {
          const double coordinate = origin_constraint_coordinate<component_name, Model>(view.points[i]);
          const double abs_coordinate = std::abs(coordinate);
          if (abs_coordinate < best_abs_coordinate) {
            best_abs_coordinate = abs_coordinate;
            has_non_negative_best = coordinate >= 0.0;
          } else if (abs_coordinate == best_abs_coordinate && coordinate >= 0.0) {
            has_non_negative_best = true;
          }
        }
 
        std::vector<types::global_dof_index> candidates;
        for (uint i = 0; i < view.dofs.n_elements(); ++i) {
          const double coordinate = origin_constraint_coordinate<component_name, Model>(view.points[i]);
          const double abs_coordinate = std::abs(coordinate);
          if (abs_coordinate != best_abs_coordinate) continue;
          if ((coordinate >= 0.0) != has_non_negative_best) continue;
          candidates.push_back(view.dofs.nth_index_in_set(i));
        }
 
        return candidates;
      }
 
      template <FixedString component_name, typename Context>
      std::optional<types::global_dof_index> select_origin_candidate(const Context &context, const auto &view)
      {
        constexpr int dim = Context::dimension;
        static_assert(dim == 1, "select_origin_candidate supports only one-dimensional domains; use "
                                "select_origin_candidates for multi-dimensional domains.");
 
        const auto candidates = select_origin_candidates<component_name, void>(context, view);
        if (candidates.empty()) return std::nullopt;
        return candidates.front();
      }
    } // namespace internal
 
    template <FixedString component_name, typename Model> class ConstrainOriginBoundaryPointToZero
    {
    public:
      template <typename Constraints, typename Context>
      void apply_boundary_affine_constraints(Constraints &constraints, const Context &context) const
      {
        const auto candidates = internal::select_origin_candidates<component_name, Model>(
            context, context.template boundary<component_name>());
        for (const auto dof : candidates) {
          constraints.add_line(dof);
          constraints.set_inhomogeneity(dof, 0.0);
        }
      }
    };
 
    template <FixedString component_name, typename Model> class ConstrainOriginSupportPointToZero
    {
    public:
      template <typename Constraints, typename Context>
      void apply_affine_constraints(Constraints &constraints, const Context &context) const
      {
        const auto candidates = internal::select_origin_candidates<component_name, Model>(
            context, context.template support<component_name>());
        for (const auto dof : candidates) {
          constraints.add_line(dof);
          constraints.set_inhomogeneity(dof, 0.0);
        }
      }
    };
 
    template <typename Model> class NoAffineConstraints
    {
    };
 
    class Time
    {
    public:
      void set_time(double t);
      const double &get_time() const;
 
    protected:
      // Initialized, like fRG::t below: a model is routinely built and read before the timestepper
      // first calls set_time(), and an indeterminate t is a bug that only shows up once the memory
      // happens to be dirty.
      double t = 0.;
    };
 
    class fRG
    {
    public:
      fRG(double Lambda);
 
      fRG(const ConfigTree &config);
 
      void set_time(double t);
 
      const double &get_time() const;
 
    protected:
      const double Lambda;
      double t = 0., k = 0., k2 = 0., k3 = 0., k4 = 0., k5 = 0., k6 = 0.;
      bool time_initialized = false;
    };
  } // namespace def
} // namespace DiFfRG