Tutorial: convergence test
==========================

An essential property of numerical methods is their accuracy. With
standard approaches to differential equations, one may use
generic convergence tests to characterize algorithms and/or
implementations.

As an example usage of TurTLE, we show here how to characterize
convergence for the field solver, or the AdamsBashforth
particle solver.
The purpose of this tutorial is to introduce TurTLE users to nontrivial
workflows that involve keeping track of a number of simulations and
prescribing their initial conditions.

The field solver test is a fairly simple pursuit, showing basic
operations with the `turtle` launcher.

The particle solver test is also minimal within the methodological
constraints imposed by the use of the Adams Bashforth integration
method, leading to a relatively increased complexity.
The example shows how to account for a number of technical details: full
control over multi-step method initialization and full control of fluid
solver forcing parameters.


Navier-Stokes solution
----------------------

TurTLE implements two distinct solvers of the Navier-Stokes equations:
one for their vorticity formulation, and one for the more common velocity
formulation.
While there are inescapable differences in the two implementations, they
both use the same time-stepping algorithm (a 3rd order RungeKutta
method).

We test convergence properties for a "statistically stationary initial
condition", since the smoothness properties of the field will influence
convergence properties in general.
As an example of the variety of questions that we can address with
TurTLE, we also perform a direct comparison of the NSVE (vorticity
solver) and NSE (velocity solver) results.

A basic Python script that performs this convergence test is provided
under `examples/convergence/fields_temporal.py`, and it is discussed
below.
In brief, simulations are performed using the same initial conditions,
but different time steps.
Error estimates are computed as differences between the different
solutions.
The figure shows the L2 norm of these differences, normalized by
the root-mean-squared velocity magnitude.
Note that the actual value of this "relative error" is irrelevant
(it depends on total integration time, spatial resolution, etc), but the
plot shows that (1) solutions converge with the 3rd order of the
timestep and (2) the systematic difference between NSE and NSVE behaves
quite differently from the two individual error estimates.

.. image:: err_vs_dt_field.svg

Particle tracking errors
------------------------

There is a more complex mechanism that generates numerical errors for
particle tracking in this setting.
In particular, the accuracy of the field itself is relevant, as well as
the accuracy of the interpolation method coupled to the ODE solver.
For instance the default particle tracking solver is a 4th order
Adams Bashforth method in TurTLE, which seems at first unnecessary since the
fluid only provides a 3rd order solution.
The default interpolation method is that of cubic splines, which should
in principle limit ODE error convergence to 2nd order (because convergence
relies on existence of relevant Taylor expansion terms, and interpolated
fields only have a finite number of continuous derivatives).
This discussion is briefly continued below, but we will not do it justice
in this tutorial.

We present below a Python script that performs a particle tracking
convergence test, provided as
`examples/convergence/particles_temporal.py`.
The figure shows trajectory integration errors, in a fairly simple
setting, for two interpolation methods: cubic splines (I4O3) and 5th
order splines on a kernel of 8 points (I8O5).
Furthermore, we show both the maximum and the mean errors for the same
set of trajectories.
While both mean errors seem to converge with the 3rd order of the
timestep (consistent with the error of the velocity field), the maximum
error obtained for cubic splines converges with the second order ---
which is the expected behavior for particle tracking in a field that
only has one continuous derivative.

.. image:: err_vs_dt_particle.svg

The results would merit a longer and more careful analysis outside the
scope of this particular tutorial, so we limit ourselves to noting that
(1) different ODEs will have different numerical properties (i.e.
rotating particles, or inertial particles) and (2) interpolation errors
may in general grow with the Reynolds number, since interpolation
errors depend on values of gradients (however the distribution of
extreme gradient values also changes with Reynolds number, so the
overall behavior is complex).


Python code for fields
----------------------

.. include:: convergence_fields.rst


Python code for particles
-------------------------

.. include:: convergence_particles.rst