sandbox/Antoonvh/KdV.h

    A solver for the Korteweg-De Vries Equation

    The Korteweg-De Vries(KdV) equation reads:

    \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;
    }