model.hh

Go to the documentation of this file.

#pragma once
 
// external libraries
#include "DiFfRG/common/json.hh"
#include <deal.II/base/point.h>
#include <deal.II/base/tensor.h>
 
// DiFfRG
#include <DiFfRG/model/ad.hh>
#include <DiFfRG/model/component_descriptor.hh>
#include <DiFfRG/model/numflux.hh>
 
namespace DiFfRG
{
  using namespace dealii;
 
  namespace def
  {
    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(std::array<NumberType, n_fe_functions> &m_i, const Point<dim> &x, const Vector &u_i,
                const Vector_dot &dt_u_i) const
      {
        // Just to avoid warnings
        (void)x;
        (void)u_i;
 
        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, const Point<dim> &x) const
      {
        // Just to avoid warnings
        (void)x;
 
        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(std::array<Tensor<1, dim, NumberType>, n_fe_functions> &F_i, const Point<dim> &x,
                const Solutions &sol) const
      {
        // Just to avoid warnings
        (void)F_i;
        (void)x;
        (void)sol;
      }
 
      template <int dim, typename NumberType, typename Solutions, size_t n_fe_functions>
      void source(std::array<NumberType, n_fe_functions> &s_i, const Point<dim> &x, const Solutions &sol) const
      {
        // Just to avoid warnings
        (void)s_i;
        (void)x;
        (void)sol;
      }
 
      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 <typename Vector> void initial_condition_variables(Vector &v_a) const
      {
        // Just to avoid warnings
        (void)v_a;
      }
 
      template <typename Vector, typename Solution> void dt_variables(Vector &r_a, const Solution &sol) const
      {
        // Just to avoid warnings
        (void)r_a;
        (void)sol;
      }
 
 
 
      template <int dim, typename Vector, typename Solutions>
      void extract(Vector &, const Point<dim> &, const Solutions &) const
      {
      }
 
 
 
      template <uint dependent, int dim, typename NumberType, typename Vector, size_t n_fe_functions_dep>
      void ldg_flux(std::array<Tensor<1, dim, NumberType>, n_fe_functions_dep> &F, const Point<dim> &x,
                    const Vector &u) const
      {
        (void)F;
        (void)x;
        (void)u;
      }
 
      template <uint dependent, int dim, typename NumberType, typename Vector, size_t n_fe_functions_dep>
      void ldg_source(std::array<NumberType, n_fe_functions_dep> &s, const Point<dim> &x, const Vector &u) const
      {
        (void)s;
        (void)x;
        (void)u;
      }
 
 
 
      template <int dim, typename NumberType, typename Solutions_s, typename Solutions_n>
      void face_indicator(std::array<NumberType, 2> & /*indicator*/, const Tensor<1, dim> & /*normal*/,
                          const Point<dim> & /*p*/, const Solutions_s & /*sol_s*/, const Solutions_n & /*sol_n*/) const
      {
      }
 
      template <int dim, typename NumberType, typename Solution>
      void cell_indicator(NumberType & /*indicator*/, const Point<dim> & /*p*/, const Solution & /*sol*/) const
      {
      }
 
 
      template <int dim, typename Vector> std::array<double, dim> EoM(const Point<dim> &x, const Vector &u) const
      {
        // Just to avoid warnings
        (void)x;
        return std::array<double, dim>{{u[0]}};
      }
 
      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([&](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(DataOut &output, const Point<dim> &x, const Solutions &sol) const
      {
        // Just to avoid warnings
        (void)output;
        (void)x;
        (void)sol;
      }
 
      template <int dim, typename Constraints>
      void affine_constraints(Constraints &constraints, const std::vector<IndexSet> &component_boundary_dofs,
                              const std::vector<std::vector<Point<dim>>> &component_boundary_points)
      {
        // Just to avoid warnings
        (void)constraints;
        (void)component_boundary_dofs;
        (void)component_boundary_points;
      }
 
    };
 
    class Time
    {
    public:
      void set_time(double t);
      const double &get_time() const;
 
    protected:
      double t;
    };
 
    class fRG
    {
    public:
      fRG(double Lambda);
 
      fRG(const JSONValue &json);
 
      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