Tutorial 5: Extractors, and reading the solution away from the EoM#

Every tutorial so far flowed a field-dependent potential and read the answer off at its minimum. That is usually what you want – the minimum is the physical state – but not always. This tutorial is about the case where it is not, and about the two pieces of DiFfRG that handle it: extractors, which carry data from the field-space solution into the ODE sector, and extractor_point, which lets the model choose where that data is read.

The full code is in Tutorials/tut5. There is no Mathematica here and no generated kernel: in \(d=3\) with a Litim regulator the loop integral is elementary, so the whole flow is four lines of C++ and the tutorial stays about the mechanism.

The physics: a kink in the potential#

We flow the \(O(N)\) effective potential \(V_k(\rho)\), \(\rho = \phi_i\phi_i/2\), through its derivative \(u = \partial_\rho V\). The curvature masses are

\[m_\pi^2 = u\,,\qquad m_\sigma^2 = u + 2\rho\,\partial_\rho u\,,\]

and with a Litim regulator \(R_k(q) = (k^2-q^2)\theta(k^2-q^2)\) the Wetterich equation collapses to something you can write down by hand. Inside the shell \(\partial_t R = 2k^2\) and the propagator is the constant \(1/(k^2+m^2)\); the shell has volume \(k^3/(6\pi^2)\) in three dimensions, so

\[\partial_t V_k = -\frac{k^5}{6\pi^2}\left[\frac{N-1}{k^2+m_\pi^2} + \frac{1}{k^2+m_\sigma^2}\right]\,.\]

Differentiating once in \(\rho\) turns this into a conservation law for \(u\), which is what the Kurganov-Tadmor finite-volume assembler solves. The Goldstone term depends on \(u\) alone and is the advective flux; the radial term carries \(\partial_\rho u\) through \(m_\sigma^2\) and is the diffusive one:

template <typename NT> NT loop(const NT m2) const { return powr<5>(k) / (6. * powr<2>(M_PI) * (k2 + m2)); }

Now tune the couplings close to a first-order transition – in parameter.toml, lambda2 = 0.04857, lambda4 = -0.5, lambda6 = 1.0. Away from such a point the potential is either plainly symmetric or plainly broken and \(u\) stays smooth. Close to it, the flow develops a kink: \(u\) stays continuous while \(\partial_\rho u\) jumps, at a small but finite \(\rho\). It separates the symmetric region around the origin from the broken branch beyond it.

This is where the EoM stops being the whole story. In this parameter range the global minimum is still at \(\rho = 0\), so everything a readout reports is a statement about the symmetric phase. The interesting physics – the branch that becomes the broken phase – lives on the far side of the kink.

Step 1: extractors, at the default point#

An extractor is a scalar the model reads off the field-space solution at one point per evaluation. It is the only bridge from the FE sector into the Variables sector: whatever extract writes is exactly what dt_variables sees.

using FEFunctionDesc = FEFunctionDescriptor<Scalar<"u">>;
using VariableDesc   = VariableDescriptor<Scalar<"A">>;
using ExtractorDesc  = ExtractorDescriptor<Scalar<"m2Sigma_x">, Scalar<"rho_x">>;
using Components     = ComponentDescriptor<FEFunctionDesc, VariableDesc, ExtractorDesc>;

template <typename NT, typename Solution>
void extract(std::array<NT, Components::count_extractors()> &extractors, const Point<dim> &x,
             const Solution &sol) const
{
  const auto rho = x[0];
  const auto u = get<"fe_functions">(sol)[idxf("u")];
  const auto du = get<"fe_derivatives">(sol)[idxf("u")][0];

  extractors[idxe("m2Sigma_x")] = u + 2. * rho * du;
  extractors[idxe("rho_x")] = rho;               // remember *where* we were
}

Storing rho_x looks redundant, and with the default behaviour it is: extract is called at the EoM, so rho_x is the EoM. Run the tutorial as it stands but with the extractor_point method commented out, and the two columns of the output agree exactly:

        k    rho_EoM   rho_extract     m2s(EoM)    m2s(extr)
  0.01832   0.000000      0.000000      0.01188      0.01188

That is the baseline. It is also the entire contract you get for free.

Step 2: why the EoM is sometimes the wrong place#

Look at the profile once the kink has formed (output.h5, group FE). \(u\) dips to a positive minimum around \(\rho \approx 0.08\) and never crosses zero – hence \(\rho_\mathrm{EoM} = 0\) – while a sharp step in \(\partial_\rho u\) sits near \(\rho \approx 5\times10^{-3}\). The readout at the origin cannot see it.

This is not an artefact of the tutorial’s toy. It is the situation that motivates the hook in practice: an LPA\('\) truncation flows a wave function \(Z_\phi\) from an anomalous dimension \(\eta_\phi\), and \(\eta_\phi\) is generated by the three-point vertex

\[\Gamma_{\pi\pi\sigma} = \delta_{ab}\,\sigma\,\partial_\rho^2 V\,,\qquad \sigma = \sqrt{2\rho}\,,\]

which vanishes identically at \(\rho = 0\). Evaluated at the EoM in this phase, \(\eta_\phi\) is exactly zero and the truncation collapses to LPA. Evaluated just above the kink, it is not.

Step 3: `extractor_point`#

A shock is not a pointwise property. The EoM machinery in discretization/common/eom.hh builds a scalar potential from the pointwise field Model::EoM(x, u) and minimises it, which can express “where does \(u\) vanish” but not “where does \(\partial_\rho u\) jump” – that is a statement about a cell and its neighbours. So the model is handed the whole discrete solution, as a DiFfRG::SolutionSample: one entry per active cell, at the cell centre, sorted by coordinate, carrying values, gradients and the local grid spacing. Building it is cheap by construction – for a finite-volume solution the values are the cell averages, read straight off the dofs, and the gradients are central differences. It is emphatically not the scheme’s limited reconstruction: that would mean a stencil fill over the whole mesh on every residual evaluation, and it stays where it belongs, in the flux path.

template <int d, typename NT> Point<d> extractor_point(const Point<d> &EoM, const SolutionSample<d, NT> &sample) const
{
  const auto shock = find_lowest_shock(sample);
  if (!shock) return EoM;
  return sample[std::min(*shock + prm.offset_cells, sample.size() - 1)].point;
}

Defining this method is the whole opt-in. There is deliberately no default implementation in DiFfRG::def::AbstractModel: its absence is what the DiFfRG::HasExtractorPoint concept detects, so models that do not need it never pay for building the sample.

The detector itself is ordinary numerics, and the tutorial’s version is cut down, but three of its choices are worth stating because they are easy to get wrong:

Measure the step in the slope, never a jump in \(u\). The kink is a weak discontinuity: \(u\) is continuous. In a finite-element discretization the jump in \(u\) across a node is \(O(\Delta\rho)\), so any threshold on it is really a threshold on the grid spacing, and the answer moves when you refine.
Require the step to be locally realised. The solver occasionally puts a small spike into \(\partial_\rho u\) with a smooth shoulder leading up to it. The shoulder alone tilts the medians on either side against each other and fakes a step, but only a fraction of it is realised across the transition zone itself. Comparing the variation actually realised there against the step is the test that separates real kinks (ratio \(\approx 1\)) from those artefacts (\(0.3\)–\(0.4\)). Note that the window this is measured over has to span the whole transition zone – measured too narrowly it undercounts a genuine kink and rejects it.
A change of grid spacing is a step in \(\partial_\rho u\) too. On the piecewise-uniform grid used here, a candidate whose window straddles a block boundary says nothing about the solution. SolutionSample carries cell_width, so this is a direct test rather than a hard cutoff in \(\rho\).

Run it now and the two columns come apart exactly when the kink appears, and stay apart:

        k    rho_EoM   rho_extract     m2s(EoM)    m2s(extr)
02237   0.000000      0.000000      0.01188      0.01188
02024   0.000000      0.005470      0.01188      0.00086   <-- kink
01832   0.000000      0.005270      0.01188      0.00047
01500   0.000000      0.004990      0.01188      0.00006
01228   0.000000      0.004770      0.01188      0.00007
00674   0.000000      0.004490      0.01188      0.00002

Detection switches on at \(k \approx 0.02\), when the kink becomes sharp enough to stand out, and then holds for every remaining frame. The kink moves inward as \(k\) falls, and \(m_\sigma^2\) read there is two to three orders of magnitude below its value at the origin.

The readouts did not move. rho_EoM is still \(0\) and m2s(EoM) is still the value at the origin, because readouts is evaluated at the EoM supplied by readouts_multiple, and extractor_point does not touch it. That separation is the point of the hook: the model reports where it should and feeds its Variables from where it must.

Note

extract runs on every residual and Jacobian evaluation, not once per output frame, so a detector that flips between cells inside a Newton step will cost you. If yours is marginal, the fix is to make it decide more stably – widen the windows, or cache the last accepted position and require a new candidate to beat it by a margin, as Examples/QuarkMesonLPAprime does for its EoM.

Warning

The extractor in this tutorial only feeds a diagnostic. The moment an extractor feeds the flux instead, its cost changes character, and not in the way you would guess.

An implicit timestepper builds its prediction from the assumption that the right-hand side varies smoothly in time. An extractor that is re-solved at every evaluation changes the right-hand side every time, and the predictor resets on each change – regardless of how small it is, because the finite-volume assembler does not put \(\partial(\text{extractor})/\partial u\) into its Jacobian. So the cost scales with how often the extractor moves, not by how much.

Measured on the LPA\('\) model this tutorial points to, over one unit of RG-time, in residual evaluations against LPA’s 202 for the same interval:

extractor refreshed	residual evaluations
every evaluation	2006
every 0.1 in RG-time	666
held fixed	~250

Held fixed it costs what LPA costs, so the extra term in the flux is not stiff at all – the whole expense is that it moves. If you hit this, refresh the extractor on a coarse RG-time grid and accept the lag, rather than reaching for the timestepper tolerances. Tightening those is the natural first instinct and it is the wrong lever.

Step 4: extractors driving a Variable#

Extractors exist to feed the ODE sector. This tutorial closes that loop with a plain diagnostic, accumulated along the flow at whatever point step 3 chose:

\[\partial_t A = m_\sigma^2(\rho_x)\,,\]

template <typename Vector, typename Solution> void dt_variables(Vector &r_a, const Solution &sol) const
{
  r_a[idxv("A")] = -get<"extractors">(sol)[idxe("m2Sigma_x")];   // residual convention: A_dot = -r_a
}

A has no field-space discretization at all – it is one number carried alongside the PDE. In tut5.cc that is reflected in the timestepper:

using TimeStepper = TimeStepperSUNDIALS_IDA_BoostABM<Assembler>;

IDA drives the stiff field-space sector adaptively while the Adams-Bashforth-Moulton stepper carries A explicitly at a fixed step, so its Jacobian is never needed. That stepper reads both [timestepping.implicit] and [timestepping.explicit] from the parameter file, and the constructor throws unless minimal_dt <= dt <= maximal_dt in the explicit block – so fill it in rather than relying on the defaults.

Warning

A here is a diagnostic, not a truncation. A real LPA\('\) would flow \(Z_\phi\) with \(\partial_t Z_\phi = -\eta_\phi Z_\phi\), computing \(\eta_\phi\) inside extract from the two-point function – and, because \(\eta_\phi\) sits inside its own regulator insertion, solving for it by fixed-point iteration there. See Examples/QuarkMesonLPAprime for that pattern.

What to take away#

extract is evaluated at extractor_point; readouts is evaluated at its own EoM. They are independent.
extractor_point is optional and detected by its presence. Omit it and nothing changes.
SolutionSample gives the model the whole discrete solution, which is what non-local features require. It samples cell centres, so its resolution is the mesh’s – which is the scale such features are defined at anyway.
The tolerance, not the grid, is what pins a kink position down: at rel_tol = 1e-7 it still drifts by several percent between grids, at 1e-8 it is stable.

Next steps#

Tutorial 3 generates the flow kernel symbolically instead of writing it by hand.
Examples/QuarkMesonLPAprime is a full LPA\('\) model: extractors feeding three Variables, with a self-consistent iteration for the anomalous dimensions.
The EoM search itself, and the /discretization/EoM_* settings that control it, are documented in Tutorial 1.