sandbox/tutorial_bgum2019/diffusion2.h

A diffusion-equation solver

For a species-field concentration c, the diffusion equations reads (or a version thereof),

\displaystyle \frac{\partial c}{\partial t} = \nabla \cdot \left( \kappa \nabla c \right).

With \kappa the diffusivity of the medium. Because of reasons, we wish to use a finite-volume discretization to solve for the temporal evolution of c(\mathbf{x},t). We recognize the flux vector (\mathbf{F}):

\displaystyle \frac{\partial c}{\partial t} + \nabla \cdot \mathbf{F} = 0,

as,

\displaystyle \mathbf{F} = -\kappa \nabla c.

Given the alligment of the faces, we can make a function that writes the fluxes (F) of c in a face vector field. We also request its user to provide the diffusivity on faces, which could be constant:

void flux_diffusion (scalar c, (const) face vector kappa, face vector F) {

Neighboring cells are possibly behind a boundary of some sort. As such, we need to compute the relevant solution values.

  boundary ({c});

By convention, a face separates a cell from its left, bottom, front neighbor in the x,y and z direction, respectively. The foreach_face() function will automatically rotate over all dimensions.

  foreach_face()
    F.x[] = -kappa.x[]*(c[] - c[-1])/Delta; //2nd order accurate gradient *estimation*
}

Using this function we can easily compute the tendency for c in a cell (j) with N faces:

\displaystyle \left(\frac{\partial c}{\partial t}\right)_j = \frac{1}{V_j} \sum_{i=1}^N \mathbf{F_i} \cdot \mathbf{A_i}.

Assuming that the ratio \frac{V_j}{A_i}=\Delta. Which is true for one-dimensional, square or cubic cells.

void tendency_from_flux (face vector F, scalar dc) {
  boundary_flux ({F});  //flux on levels
  foreach() {
    dc[] = 0;
    foreach_dimension() //Rotates over the dimensions
      dc[] +=  (F.x[] - F.x[1])/Delta; 
  }
}

With the tendency, the solution can be advanced in time with timestep dt.

void advance (scalar c, scalar dc, double dt) {
  foreach()
    c[] += dt*dc[];
}

Now we can construct a user-interface function, diffusion() that employs a fordward-Euler scheme.

void diffusion (scalar c, double dt, face vector kappa) {
  face vector F[];
  scalar dc[];
  flux_diffusion (c, kappa, F);
  tendency_from_flux (F, dc);
  advance (c, dc, dt);
}

At the cost of more computational effort and memory, we could also choose to use the mid-point rule to update the solution. It is second-order accurate in time!

void diffusion_midpoint (scalar c, double dt, face vector kappa) {
  face vector F[];
  scalar dc[], c_temp[];
  foreach()
    c_temp[] = c[];                 //create a scratch
  flux_diffusion(c_temp, kappa, F); 
  tendency_from_flux (F, dc);
  advance (c_temp, dc, dt/2.);      //advance to the mid point.
  flux_diffusion(c_temp, kappa, F); //Re-estimate the fluxes at the mid point,
  tendency_from_flux (F, dc);       //and the tendency there.
  advance (c, dc, dt);              //update the solution
}

Further,

We could use the Runke-Kutta schemes (upto 4th order) or a predictor corrector scheme.