# A solver for the Korteweg-De Vries Equation

\displaystyle \frac{\partial c}{\partial t} = 6c\frac{\partial c}{\partial x} - \frac{\partial^3 c}{\partial x^3}

When we integrate over a finite volume (V) the evolution equation for the cell averaged value of c (C) reads \frac{\partial C}{\partial t} = \frac{1}{V}\Sigma F_i\cdot A_i, where the right-hand side concerns a summation of fluxes over the cell’s face areas (A_i). For F we employ the divergence theorem and write:

\displaystyle F_i = 3c_i^2-\left(\frac{\partial^2 c}{\partial x^2}\right)_i

These fluxes can be aproximated from the spatial structure of C(x) in the neighborhood of a cell’s face employing a discretized approach and the numerical solution of C can be advanced in time, using a discrete timestep (dt) according to the popular fourth-order Runge-Kutta time-integration method.

#define BGHOSTS 2
#include "runge-kutta.h"
#include "utils.h"

This function computes the tendencies with second order spatial accuracy. This requires two layers of ghost cells:

static void KdV2 (scalar * s, double six, scalar *kl){
boundary(s);
for (int j = 0; j < list_len(s); j++){
scalar m = s[j];
scalar k = kl[j];
face vector flx[];
foreach_face()
flx.x[] = 0.75*(six/6)*sq((m[] + m[-1])) - (m[-2] - m[-1] - m[] + m[1])/(2.*sq(Delta));
boundary_flux({flx});
foreach(){
k[] = 0;
foreach_dimension()
k[] += (flx.x[1] - flx.x[])/Delta;
}
}
}

The forward-in-time integrator is prone to instabilities and since the user may not be aware of any relevant details we take some care here: The discretized timestep Dt can be compared against the grid-cell size (\Delta) and the tendency from each term. For an equation of the form:

\displaystyle \partial_tc = ac\partial_xc

with a a constant with units [a] = LT^{-1}[c]^{-1}. thus we formulate a Courant number:

\displaystyle Co = \frac{a\mathrm{Dt}c}{\Delta}

Similarly for an equation of the form,

\displaystyle \partial_tc = b\partial_x^3c,

with b a constant with units [b] = L^3T^{-1} we define the KdV number (KdVn):

\displaystyle KdVn = \frac{\mathrm{Dt}b}{\Delta^3}.

\mathrm{Dt} should be chosen such that both Co and KdVn remain within limited bounds. By default we use:

double Co = 0.8;
double KdVn = 0.5;

If the user-requested timestep dt>\mathrm{Dt}, the time integration is split into it equal parts to ensure a limited value of Co and KdVn. The function below provides the described user interface and returns the number of iterations (it). The input is a list of scalar fields and the time step dt.

int KdV(scalar * s, double dt){
double six = 6.;
double delt = L0/(double)(1 << depth());
double maxs = 1E-30;
for (scalar m in s){
double max = fabs(statsf(m).max);
if (max > maxs)
max = maxs;
}
double Dtmin = min(Co*delt/(six*maxs), KdVn*cube(delt));
int it = 1;
if (Dtmin < dt)
it = (int)((dt/(Dtmin))+0.99);
double Dt = dt/((double)(it));
for (int iter = 0; iter < it; iter++)
runge_kutta(s, six, Dt, KdV2, 4);
return it;
}