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 const
ant:
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
.
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.