/**
# Double projection

This option for the [centered Navier--Stokes solver](centered.h) is
inspired by [Almgren et al., 2000](#almgren200). These authors first
recall that, while for the exact projection method, the definition of
the pressure is unambiguous, in the case of the *approximate*
projection method (used in [centered.h]()), several pressures can be
defined.

Exploiting the different properties of these different definitions may
be useful in some cases.

## Standard approximate projection

We first summarise the standard approximate projection scheme, as
implemented in [centered.h]().

1. Advection/viscosity
$$
 \frac{u^{\star} - u^n}{\Delta t} = (u \cdot \nabla u)^{n + 1 / 2} +
	\alpha \nabla \cdot (\mu \nabla D)^{\star}
$$
2. Acceleration
$$
  u^{\star}_f = \overline{u^{\star}} + \Delta ta_f
$$
3. Projection
$$
  u^{n + 1}_f = P (u^{\star}_f, p^{n + 1})
$$
4. Centered pressure gradient correction
$$
  g^{n + 1} = \overline{a_f - \alpha_f \nabla p^{n + 1}}
$$
$$
  u^{n + 1} = u^{\star} + \Delta tg^{n + 1}
$$

with $P(u,p)$ the projection operator and the overline designating 
either cell-center to cell-face averaging or reciprocally.

## Double approximate projection

As its name indicates this scheme adds an extra projection step and
defines two pressures: the standard one used to project the
face-centered velocity field $u_f$ (renamed $\mathbf{\delta p}$
below), and $p^{n+1}$ used to compute the cell-centered pressure
gradient $g$. This new pressure is obtained by projecting only the
*update* (i.e. its evolution in time) to the centered velocity
field. Note that in the case of an exact projection these two
projections are identical since the divergence of the velocity field
at the start of the timestep is zero.

The scheme can be summarised as:

1. Advection/viscosity
$$
 \frac{u^{\star} - u^n}{\Delta t} = (u \cdot \nabla u)^{n + 1 / 2} +
	\alpha \nabla \cdot (\mu \nabla D)^{\star} 
	\mathbf{+ g^n = A^{n+1/2} + g^n}
$$
$$
 \mathbf{A_f = \overline{A^{n+1/2}} + \Delta ta_f}
$$
2. Acceleration
$$
  u^{\star}_f = \overline{u^{\star}}
$$
3. Projection
$$
  u^{n + 1}_f = P (u^{\star}_f, \mathbf{\delta p})
$$
4. Approximate projection
$$
  u^{n + 1} = u^{\star} - \Delta t \overline{\alpha_f \nabla \delta p}
$$
5. Second projection
$$
\mathbf{P(A_f,p^n+1)}
$$
$$
  g^{n + 1} = \overline{a_f - \alpha_f \nabla p^{n + 1}}
$$

where the additions to the previous scheme are highlighted in bold.

Why is this useful? The new pressure does not feel the *history* of
divergence of the centered velocity field. This is useful in
particular when this history includes the noise induced by adaptive
mesh refinement.

The cost to pay is however significant since an extra (potentially
expensive) projection is required.

## Implementation

We need a (temporary) field to store the update $A_f$ and a
(permanent) field to store the projection pressure $\delta p$. */

face vector Af;
scalar dp[];

/**
We make heavy use of the [event inheritance mechanism](/Basilisk
C#event-inheritance). All the events below are first defined in
[centered.h]().

At the beginning of the timestep (i.e. before advection), we store the
(interpolated) value of the initial velocity field in $A_f$. */

event advection_term (i++)
{
  Af = new face vector;
  foreach_face()
    Af.x[] = - fm.x[]*face_value (u.x, 0);
}

/**
After the advection and diffusion terms have been added to $u$, we
recover the update by adding the new face-interpolated value of the
velocity field to the initial face velocity, and add the acceleration
i.e. we perform step 1 above:
$$
 A_f = \overline{A^{n+1/2}} + \Delta ta_f
$$
*/

face vector ab;

event acceleration (i++)
{
  foreach_face()
    Af.x[] += fm.x[]*(face_value (u.x, 0) + dt*a.x[]);

  /**
  We also add the centered gradient $\mathbf{g^n}$ to the centered
  velocity field. */
  
  correction (dt);

  /**
  Step 2 above (i.e. $u_f^{\star} = \overline{u^{\star}}$) is
  performed by the acceleration event of the [centered
  solver](centered.h#acceleration), but we need to reset the
  acceleration to zero. */

  ab = a;
  a = zerof;
}

/**
The projection step 3 is also performed by the [centered
solver](centered.h#projection). We want to store the resulting
pressure in $dp$ rather than $p$, so we swap the two fields, before
performing the projection. */

event projection (i++)
{
  scalar_clone (dp, p);
  swap (scalar, p, dp);
}

/**
Step 4 is performed by the [centered
solver](centered.h#projection). Step 5 is done at the end of the
timestep. We first restore the fields modified above, then perform the
second projection and compute the corresponding centered gradient
$g^{n+1}$. */

event end_timestep (i++)
{
  swap (scalar, p, dp);
  a = ab;
  // this could be optimised since we do not use Af
  mgp = project (Af, p, alpha, dt, mgp.nrelax);
  delete ((scalar *){Af});
  centered_gradient (p, g);
}

/**
## References

~~~bib
@article{almgren2000,
  title={Approximate projection methods: Part I. Inviscid analysis},
  author={Almgren, Ann S and Bell, John B and Crutchfield, William Y},
  journal={SIAM Journal on Scientific Computing},
  volume={22},
  number={4},
  pages={1139--1159},
  year={2000},
  publisher={SIAM},
  url={https://epubs.siam.org/doi/pdf/10.1137/S1064827599357024}
}
~~~
*/