conservation.h

Go to the documentation of this file.

/** @file conservation.h
 */
/**
# A generic solver for systems of conservation laws
 
Using the ideas of [Kurganov and Tadmor,
2000](references.bib#kurganov2000) it is possible to write a generic
solver for systems of conservation laws of the form
\f[
\partial_t\left(\begin{array}{c}
    s_i\\
    \mathbf{v}_j\\
 \end{array}\right) + \nabla\cdot\left(\begin{array}{c}
    \mathbf{F}_i\\
    \mathbf{T}_j\\
 \end{array}\right) = 0
\f]
where \f$s_i\f$ is a list of scalar fields, \f$\mathbf{v}_j\f$ a list of
vector fields and \f$\mathbf{F}_i\f$, \f$\mathbf{T}_j\f$ are the corresponding
vector (resp. tensor) fluxes. 
 
Note that the [Saint-Venant solver](saint-venant.h) is a particular
case of this generic algorithm.
 
The user must provide the lists of conserved scalar and vector fields
*/
 
extern scalar * scalars;
extern vector * vectors;
 
/**
as well as a function which, given the state of each quantity,
returns the fluxes and the minimum/maximum eigenvalues
(i.e. `eigenvalue[0]`/`eigenvalue[1]`) of the corresponding Riemann
problem. */
 
void flux (const double * state, double * flux, double * eigenvalue);
 
/**
## Time-integration
 
### Setup
 
Time integration will be done with a generic
[predictor-corrector](predictor-corrector.h) scheme. */
 
#include "predictor-corrector.h"
 
/**
The generic time-integration scheme in `predictor-corrector.h` needs
to know which fields are updated i.e. all the scalars + the components
of all the vector fields. It also needs a function to compute the
time-derivatives of the evolving variables. */
 
scalar * evolving;
double update_conservation (scalar * conserved, scalar * updates, double dtmax);
 
/** @brief Event: defaults (i = 0) */
void event_defaults (void) 
{
  evolving = list_concat (scalars, (scalar *) vectors);
  update = update_conservation;
 
  /**
  We switch to a pure minmod limiter by default for increased
  robustness. */
  
  theta = 1.;
 
  /**
  With the MUSCL scheme we use the CFL depends on the dimension of the
  problem. */
 
  if (CFL > 1./dimension)
    CFL = 1./dimension;
  
  /**
  On trees we need to replace the default bilinear
  refinement/prolongation with linear so that reconstructed values
  also use slope limiting. */
 
  #if TREE
  for (int _s = 0; _s < 1; _s++) /* scalar in evolving */ {
    s.refine = refine_linear;
    set_prolongation (s, refine_linear);
    set_restriction (s, restriction_volume_average);
  }
  #endif
}
 
/**
At the end of the run we need to free the list (to avoid a memory
leak). */
 
event cleanup (i = end) free (evolving);
 
/**
User initialisation happens here. */
 
event init (i = 0);
 
/**
### Computing fluxes
 
The core of the central-upwind scheme (see e.g. section 3.1 of
[Kurganov & Levy, 2002](references.bib#kurganov2002)) is the
approximate solution of the Riemann problem given by the left and
right states to get the fluxes `f`. */
 
static double riemann (const double * right, const double * left,
           double Delta, double * f, int len, 
           double dtmax)
{
  double fr[len], fl[len], er[2], el[2];
  flux (right, fr, er);
  flux (left,  fl, el);
  double ap = max(er[1], el[1]); ap = max(ap, 0.);
  double am = min(er[0], el[0]); am = min(am, 0.);
  double a = max(ap, -am); 
 
  if (a > 0.) {
    for (int i = 0; i < len; i++)
      f[i] = (ap*fl[i] - am*fr[i] + ap*am*(right[i] - left[i]))/(ap - am);
    double dt = CFL*Delta/a;
    if (dt < dtmax)
      dtmax = dt;
  }
  else
    for (int i = 0; i < len; i++)
      f[i] = 0.;
  return dtmax;
}
 
double update_conservation (scalar * conserved, scalar * updates, double dtmax)
{
  /**
  The gradients of each quantity are stored in a list of dynamically-allocated
  fields. First-order reconstruction is used for the gradient fields. */
 
  vector * slopes = NULL;
  for (int _s = 0; _s < 1; _s++) /* scalar in conserved */ {
    vector slope = {0} /* new vector */;
    for (int _d = 0; _d < dimension; _d++) {
      slope.x.gradient = zero;
      #if TREE
      slope.x.prolongation = refine_linear;
      #endif
    }
    slopes = vectors_append (slopes, slope);
  }
  gradients (conserved, slopes);
 
  /**
  We allocated fields for storing fluxes for each scalar and vector
  quantity. */
 
  vector * lflux = NULL;
  int len = list_len (conserved);
  for (int _s = 0; _s < 1; _s++) /* scalar in conserved */ {
    vector f1 = {0} /* new vector */;
    lflux = vectors_append (lflux, f1);
  }
 
  /**
  The predictor-corrector scheme treats all fields as scalars (stored in
  the `conserved` list). We need to recover vector and tensor quantities
  from these lists. To do so, knowing the number of scalar fields, we
  split the scalar list into a list of scalars and a list of vectors. */
  
  int scalars_len = list_len (scalars);
 
  scalar * scalars = list_copy (conserved);
  if (scalars) scalars[scalars_len].i = -1;
  vector * vectors = vectors_from_scalars (&conserved[scalars_len]);
  
  /**
  We then do the same for the gradients i.e. split the list of vectors
  into a list of vectors and a list of tensors. */
 
  vector * scalar_slopes = vectors_copy (slopes);
  if (scalar_slopes) scalar_slopes[scalars_len] = (vector){{-1}};
  tensor * vector_slopes = tensors_from_vectors (&slopes[scalars_len]);
 
  /**
  And again for the fluxes. */
  
  vector * scalar_fluxes = vectors_copy (lflux);
  if (scalar_fluxes) scalar_fluxes[scalars_len] = (vector){{-1}};
  tensor * vector_fluxes = tensors_from_vectors (&lflux[scalars_len]);
 
  /**
  We are ready to compute the fluxes through each face of the domain. */
  
  for (int _i = 0; _i < _N; _i++) /* foreach_face */) {
 
    /**
    #### Left/right state reconstruction 
    
    We use the central values of each scalar/vector quantity and the
    pre-computed gradients to compute the left and right states. */
    
    double r[len], l[len]; // right/left Riemann states
    double f[len];         // fluxes for each conserved quantity
    double dx = Delta/2.;
    int i = 0;
    scalar s;
    vector g;
    for (s,g in scalars,scalar_slopes) {
      r[i] = s[] - dx*g.x[];
      l[i++] = s[-1] + dx*g.x[-1];
    }
    vector v;
    tensor t;
    for (v,t in vectors,vector_slopes) {      
      r[i] = v.x[] - dx*t.x.x[];
      l[i++] = v.x[-1] + dx*t.x.x[-1];
      #if dimension > 1
        r[i] = v.y[] - dx*t.y.x[];
  l[i++] = v.y[-1] + dx*t.y.x[-1];
      #endif
      #if dimension > 2
        r[i] = v.z[] - dx*t.z.x[];
  l[i++] = v.z[-1] + dx*t.z.x[-1];
      #endif
    }
 
    /**
    #### Riemann problem
    
    We then call the generic Riemann solver and store the resulting fluxes
    in the pre-allocated fields. */
    
    dtmax = riemann (r, l, Delta*cm[]/fm.x[], f, len, dtmax);
    i = 0;
    for (int _fs = 0; _fs < 1; _fs++) /* vector in scalar_fluxes */
      fs.x[] = fm.x[]*f[i++];
    for (int _fv = 0; _fv < 1; _fv++) /* tensor in vector_fluxes */ {
      fv.x.x[] = fm.x[]*f[i++];
      #if dimension > 1
        fv.y.x[] = fm.x[]*f[i++];
      #endif
      #if dimension > 2
        fv.z.x[] = fm.x[]*f[i++];
      #endif
    }
  }
 
  /**
  #### Update
  
  The update for each scalar quantity is the divergence of the fluxes. */
  
  for (int _i = 0; _i < _N; _i++) /* foreach */ {
    scalar ds;
    vector f;
    for (ds,f in updates,lflux) {
      ds[] = 0.;
      for (int _d = 0; _d < dimension; _d++)
  ds[] += (f.x[] - f.x[1])/(cm[]*Delta);
    }
  }
 
  /**
  #### Cleanup
  
  We finally deallocate the memory used to store lists and gradient
  fields. */
  
  free (scalars);
  free (vectors);
  free (scalar_slopes);
  free (vector_slopes);
  free (scalar_fluxes);
  free (vector_fluxes);
  delete ((scalar *) slopes);
  free (slopes);
  delete ((scalar *) lflux);
  free (lflux);
  
  return dtmax;
}