Numerical Models#

In DiFfRG the computational setup is fully described by a numerical model. A model essentially is a description of a large set of coupled differential equations and some additional methods to handle data output.

In general, we have three components to any flow, which are FE functions \( u_i(x),\,i\in\{0,\dots,N_f-1\} \), variables \( v_a,\,a\in\{0,\dots,N_v-1\} \) and extractors \( e_b,\,b\in\{0,\dots,N_e-1\} \).

The latter two are just independent variables, whereas the FE functions depend additionaly on a field variable, \( u_i(x),\,x\in\mathbb{R}^d \). In other words, the FE functions explicitly live on a spatial discretization of the field space.

Defining a model#

Any model you define should at least be inherited from the abstract model class DiFfRG::def::AbstractModel to ensure that all necessary methods are at least defined to do nothing, i.e.

using namespace DiFfRG;
class MyModel : public def::AbstractModel<MyModel>, ...
{
  ...
};

Inside the class we can now overwrite all methods from DiFfRG::def::AbstractModel in order to implement the right system of flow equations.

Spatial discretization#

The FE functions usually correspond to expansion coefficients in a derivative expansion. As an example consider a bosonic theory as in this paper: Treating a purely bosonic theory in first-order derivative expansion, the effective action is given by

\[\large \Gamma_k[\phi] = \int_x \bigg(\frac{1}{2}Z(\rho)(\partial_\mu\phi)^2 + V(\rho) \bigg)\,,\]

where \( \rho = \phi^2 / 2 \). A flowing reparametrization of the field \( \phi \) is being performed and is given by

\[\large \dot\phi(x) = \frac{1}{2} \eta(\rho) \phi\,.\]

where we introduced the anomalous dimension \( \eta = \frac{\partial_{t_+} Z}{Z} \).

The flow is then fully parametrized in terms of FE functions

\[\begin{split}\begin{aligned}\large u_1(x) &= m^2(\rho) = \partial_\rho V(\rho)\,, \\\large u_2(x) &= \eta(\rho)\,, \end{aligned}\end{split}\]

where we also chose the field \( x = \rho \). We see here that the FE functions live on a spatial discretization of the d-dimensional field space \( \mathbb{R}^d \).

With the above ansatz one can quickly compute flow equations from the Wetterich equation,

\[\large k\partial_k \Gamma_k[\Phi] = \frac{1}{2}\text{Tr}\, G_{\alpha\beta}\,k\partial_k R^{\alpha\beta}\,.\]

We remark that the time

\[\large t = t_+ = \ln\left(\frac{\Lambda}{k}\right)\,,\]

as used in DiFfRG, is opposite in sign to the RG-time as defined in most literature, \(t_- = \ln\left(\frac{k}{\Lambda}\right)\). This is simply due to many time solvers not accepting negative time arguments.

In order to discretize the flow equations on a finite element space, the flow equations are expressed in the standard differential-algebraic form

\[\large m_i(\partial_t u_j, u_j, x) + \partial_x F_i(u_j, \partial_x u_j, \partial_x^2 u_j, e_b, v_a, x) + s_i(u_j, \partial_x u_j, \partial_x^2 u_j, e_b, v_a, x) = 0\,,\]

where \( m_i \) are called the mass functions, \( F_i \) the fluxes and \( s_i \) the sources. The latter two are functions of the FE functions, their derivatives, the field variable, the variables and the extractors.

In principle, the above system of equations can contain both equations containing the time derivatives, i.e. differential components, and equations without time derivatives, i.e. algebraic components. In order to solve the resulting DAEs one is currently restricted to the SUNDIALS IDA solver, which is however highly efficient and actually recommended for most cases.

Alternatively, the restricted formulation, allowing only for differential components,

\[\large m_{ij}(x) \partial_t u_j + \partial_x F_i(u_j, \partial_x u_j, \partial_x^2 u_j, e_b, v_a, x) + s_i(u_j, \partial_x u_j, \partial_x^2 u_j, e_b, v_a, x) = 0\,,\]

is used for all other provided ODE solvers, i.e. Runge-Kutta methods.

Note, that in the above definitions a change from \( t = t_+ \to t_- \) simply moves all terms onto the other side, i.e. when calculating the flow equations in the standard \(t_-\), one can still copy and paste everything without changing signs if the mass functions are simply \( m_i = \partial_{t_+} u_i = - \partial_{t_-} u_i \) (as is default).

The above components of standard form have direct analogues in the abstract numerical model which must be reimplemented for any system of flow equations. For actual implementation examples, especially regarding the template structure, please take a look at the models contained in the DiFfRG/models/ folder.

The relevant methods are also documented in DiFfRG::def::AbstractModel and read as follows:

The mass function \(m_i(\partial_t u_j, u_j, x)\) is implemented in the method

  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, const Vector_dot &u_dot) const;

Note, that the precise template structure is not important, the only important thing is that the types are consistent with the rest of the model. It is however necessary to leave at least the NumberType, Vector, and Vector_dot template parameters, as these can differ between calls.
The m_i argument is the resulting mass function \(m_i\), with \(N_f\) components. This method should fill the m_i argument with the desired structure of the flow equation. x is a d-dimensional array of field coordinates, and both u (~\(u_i(x)\)) and u_dot (~\(\partial_t u_i(x)\)) have \(N_f\) components.
The standard implementation of this method simply sets \(m_i = \partial_t u_i\).

.

If not using a DAE, the mass matrix \(m_{ij}(x)\) is implemented in the method

  template <int dim, typename NumberType, size_t n_fe_functions>
  void mass(std::array<std::array<NumberType, n_fe_functions>, M::Components::count_fe_functions()> &m_ij, const Point<dim> &x) const;

The m_ij argument is the resulting mass matrix \(m_{ij}\), with \(N_f\) components in each dimension. This method should fill the mass argument with the desired structure of the flow equation. x is a d-dimensional array of field coordinates. The standard implementation of this method simply sets \(m_{ij} = \delta_{ij}\).

.

The flux function \(F_i(u_j, \partial_x u_j, \partial_x^2 u_j, e_b, v_a, x)\) is given by

  template <int dim, typename NumberType, typename Solution, size_t n_fe_functions>
  void flux(std::array<Tensor<1, dim, NumberType>, n_fe_functions> &F_i, const Point<dim> &x, const Solution &sol) const;

Onve again, it is necessary to leave the NumberType and Solution templates, whereas the rest can be dropped. F_i has \(N_f\) components, x gives the coordinate in field space and sol contains all other arguments of the flux function. In practice, sol is a std::tuple<...> which contains 0. the array u_j

the array of arrays \(\partial_x u_j\)
the array of arrays of arrays \(\partial_x^2 u_j\)
the array of extractors \(e_b\) Lastly, the variables are communicated separately to the model, see the Other variables-section below The standard implementation of this method simply sets \(F_i = 0\).

.

The source function \(s_i(u_j, \partial_x u_j, \partial_x^2 u_j, e_b, v_a, x)\) is given by

  template <int dim, typename NumberType, typename Solution, size_t n_fe_functions>
  void source(std::array<NumberType, n_fe_functions> &s_i, const Point<dim> &x, const Solution &sol) const;

Again, it is necessary to leave the NumberType and Solution templates, whereas the rest can be dropped. s_i has \(N_f\) components, x gives the coordinate in field space and sol contains all other arguments of the flux function, with the layout as explained above in the flux case. The standard implementation of this method simply sets \(s_i = 0\).

Picking up the example from above, we can now sketch the implementation of the numerical model as follows:

using namespace DiFfRG;

using FEFunctionDesc = FEFunctionDescriptor<Scalar<"u">, Scalar<"v">>;
using Components = ComponentDescriptor<FEFunctionDesc>;
constexpr auto idxf = FEFunctionDesc{};

class MyModel : public def::AbstractModel<MyModel, Components>,
                public def::fRG,                    // this handles the fRG time
                ...
{
public:
  MyModel(const JSONValue& json) : def::fRG(json), prm(json) {}

  template <int dim, typename NumberType, typename Vector, typename Vector_dot>
  void mass(std::array<NumberType, Components::count_fe_functions()> &m_i, const Point<dim> &x, const Vector &u, const Vector_dot &u_dot) const
  {
    m_i[idxf("u")] = u_dot[idxf("u")];
    m_i[idxf("v")] = -u[idxf("v")];
  }

  template <int dim, typename NumberType, typename Solution>
  void flux(std::array<Tensor<1, dim, NumberType>, Components::count_fe_functions(0)> &F_i, const Point<dim> &x, const Solution &sol) const;
  {
    F_i[idxf("u")][0] = ...; // Flux of m^2
    F_i[idxf("v")][0] = ...; // Flux of eta
  }

  template <int dim, typename NumberType, typename Solution, typename M = Model>
  void source(std::array<NumberType, M::Components::count_fe_functions(0)> &s_i, const Point<dim> &x, const Solution &sol) const
  {
    s_i[idxf("u")] = ...; // Source of m^2
    s_i[idxf("v")] = ...; // Source of eta
  }
};

Assemblers and discretizations#

The actual numerical calculation of the flow equations (rather, their weak form) is done by the so-called assemblers. These are responsible for the actual discretization of the flow equations on the finite element space. In DiFfRG, we provide a set of assemblers for different discretizations, which are all derived from the abstract assembler class DiFfRG::AbstractAssembler.

Although all require at least the above interface methods, certain additional methods are required for certain discretizations. For example, the discontinuous Galerkin (DG) assemblers require the implementation of the numerical fluxes, and both discontinuous and continuous Galerkin (CG, also called simply FEM here) assemblers require the implementation of the boundary condition fluxes.

To understand the underlying numerics, see e.g. this review and also the excellent deal.ii tutorials.

For further reference, please refer to the documentation of the respective assemblers.

DiFfRG::DG::Assembler
DiFfRG::dDG::Assembler
DiFfRG::CG::Assembler

Underlying the assemblers are the actual discretizations, which are implemented in the DiFfRG::discretization namespace. These are responsible for the actual discretization of the field space, i.e. the construction of the finite element space. In DiFfRG, we provide a set of discretizations for different finite element spaces.

DiFfRG::DG::Discretization
DiFfRG::CG::Discretization

Running#

Putting everything together, we can write a straightforward main function to run the flow equations:

#include <DiFfRG/DiFfRG.hh>

using namespace DiFfRG;

// Make choices for types: the model, its discretization, the assembler and the timestepper
using Model = MyModel;
constexpr uint dim = Model::dim;
using Discretization = CG::Discretization<Model::Components, double, RectangularMesh<dim>>;
using VectorType = typename Discretization::VectorType;
using SparseMatrixType = typename Discretization::SparseMatrixType;
using Assembler = CG::Assembler<Discretization, Model>;
using TimeStepper = TimeStepperSUNDIALS_IDA<VectorType, SparseMatrixType, dim, UMFPack>;

int main(int argc, char *argv[])
{
  // Initialize DiFfRG (MPI + Kokkos) and read the parameter file / CLI overrides
  const auto config_helper = DiFfRG::Init(argc, argv).get_configuration_helper();
  const auto json = config_helper.get_json();

  // Define the objects needed to run the simulation
  Model model(json);
  RectangularMesh<dim> mesh(json);
  Discretization discretization(mesh, json);
  Assembler assembler(discretization, model, json);
  DataOutput<dim, VectorType> data_out(json);
  HAdaptivity mesh_adaptor(assembler, json);
  TimeStepper time_stepper(json, &assembler, &data_out, &mesh_adaptor);

  // Set up the initial condition
  FE::FlowingVariables initial_condition(discretization);
  initial_condition.interpolate(model);

  // Start the timestepping
  Timer timer;
  time_stepper.run(&initial_condition, 0., json.get_double("/timestepping/final_time"));

  // Print a bit of exit information to the logger.
  assembler.log("log");
  spdlog::get("log")->info("Simulation finished after " + time_format(timer.wall_time()));
  return 0;
}

Every object is constructed from the parsed json configuration, so all parameters (grid, FE order, tolerances, output, physics) are read from parameter.json. The timestepper is the SUNDIALS IDA solver, which is the recommended solver for most cases. If you use a discontinuous Galerkin discretization (DG/dDG/LDG) it is also necessary to supply a numerical flux, which can be done by modifying the numerical model as follows:

class MyModel : public def::AbstractModel<MyModel, Components>,
                public def::fRG,                    // this handles the fRG time
                public def::LLFFlux<MyModel>,        // use a LLF numflux
                public def::FlowBoundaries<MyModel>, // use Inflow/Outflow boundaries
                public def::AD<MyModel>              // define all jacobians per AD
{
  ...
};

Here, the local Lax-Friedrichs flux has been used for the numerical fluxes and the boundaries have been defined to be inflow-/outflow. We have also chosen to use the autodiff functionality for the calculation of the jacobians.

Other variables#

Besides FE functions, a model can carry variables \(v_a\) and extractors \(e_b\). These are degrees of freedom that do not live on the field-space discretization — they are plain values (or arrays of values) evolved alongside (or instead of) the FE functions. They are the natural representation for momentum-dependent truncations, where the flowing objects are dressing functions and couplings tabulated on a momentum grid rather than on the FEM field grid.

Variables and extractors are declared in the ComponentDescriptor next to the FE functions, using either a Scalar<"name"> (a single value) or a grid-valued FunctionND<"name", sizes...>:

using namespace DiFfRG;

static constexpr uint p_grid_size = 96;
using VariableDesc  = VariableDescriptor<Scalar<"m2">, FunctionND<"ZA", p_grid_size>>;
using ExtractorDesc = ExtractorDescriptor<Scalar<"observable">>;
// A pure variable system has an empty FEFunctionDescriptor<>:
using Components = ComponentDescriptor<FEFunctionDescriptor<>, VariableDesc, ExtractorDesc>;

constexpr auto idxv = VariableDesc{};   // idxv("ZA") + i indexes grid point i of ZA

The difference between the two: variables are evolved in RG time by their own flow equation, while extractors are auxiliary quantities computed from the current state and made available to the FE-function flux/source methods (useful for coupling a FEM sector to global quantities).

A model implements the following methods (from DiFfRG::def::AbstractModel, default to no-ops):

The initial condition of the variables,

template <typename Vector> void initial_condition_variables(Vector &v_a) const;

The flow of the variables \(\partial_t v_a\),

template <typename Vector, typename Solution> void dt_variables(Vector &r_a, const Solution &sol) const;

where r_a is the residual to fill and sol is a named tuple from which the current variables are obtained via get<"variables">(sol).

The extractors, filled from the current solution,

template <int dim, typename Vector, typename Solutions> void extract(Vector &e_b, const Point<dim> &x, const Solutions &sol) const;

A system that consists of variables only (no FE functions) is assembled by DiFfRG::Variables::Assembler (spatial dimension 0) and carried by DiFfRG::FlowingVariables; when FE functions are also present the two sectors are coupled and the FEM assembler handles both. Momentum grids are represented by coordinate systems (e.g. DiFfRG::LogarithmicCoordinates1D) and evaluated through interpolators (e.g. DiFfRG::SplineInterpolator1D); the flow kernels are generated as grid map integrators.

For a complete, worked momentum-dependent example, see Tutorial 4.