# The SAG equation

The SAG equation is a Fick diffusion equation of a specie $c$ $\frac{\partial c}{\partial t}={\nabla }^{2}c$ which can be solved with the reaction–diffusion solver.

``````#include "grid/multigrid.h"
#include "run.h"
#include "diffusion.h"``````

Concentration at time $t+dt$, $t$ and maximum difference between the two.

``````scalar c[], cold[];
double errmax;``````

The generic time loop needs a timestep. We will store the statistics on the diffusion solvers in `mgd`.

``````double dt;
mgstats mgd;``````

A “crystal growth” boundary condition of the form $\frac{\partial c}{\partial y}=bi\phantom{\rule{0.278em}{0ex}}c\phantom{\rule{0.278em}{0ex}}\phantom{\rule{0.278em}{0ex}}\mathrm{\text{on}}\phantom{\rule{0.278em}{0ex}}\phantom{\rule{0.278em}{0ex}}y=0$

is imposed for $-1 and $y=0$, where $bi$ is a kind of Biot number. This mixed condition can be discretised to second-order accuracy as

``(c[] - c[bottom])/Δ = bi*(c[] + c[bottom])/2.;``

The rest of the wall is covered with a mask so that no growth occurs i.e. $\frac{\partial c}{\partial y}=0$ (a Neumann boundary condition). Combining everything then gives

``````double bi = 1;
c[bottom] = fabs(x) < 1 ? neumann(0) : c[]*(2. - bi*Δ)/(2. + bi*Δ);``````

On the top boundary the species concentration is imposed.

``c[top] = dirichlet(1);``

And there is no flux on the right and left boundaries.

``````c[right] = neumann(0);
c[left]  = neumann(0);``````

## Parameters

The domain is the square box $\left[-5:5\right]×\left[0:10\right]$.

``````int main() {
L0 = 10.;
X0 = -L0/2;
N = 256;
errmax = 5.e-3;
run();
}``````

## Initial conditions

In the absence of a mask (i.e. with the Biot condition along the bottom boundary) an exact stationary solution is

``#define cexact(y) ((bi*y + 1)/(bi*L0 + 1))``

We use this as initial condition.

``````event init (i = 0) {
foreach() {
c[] = cexact(y);
cold[] = c[];
}
boundary ({c, cold});
}``````

## Time integration

``event integration (i++) {``

We first set the timestep according to the timing of upcoming events. We choose a maximum timestep of 0.2 which ensures the stability of the reactive terms for this example.

``  dt = dtnext (0.2);``

We use the diffusion solver to advance the system from $t$ to $t+dt$.

``  mgd = diffusion (c, dt);``

If the difference between c and cold is small enough, we stop.

``````  double dc = change (c, cold);
if (i > 0 && dc < errmax)
return 1; // stop
printf ("%g %g\n", t, dc);
}``````

## Results

``````event profiles (i += 10) {
FILE * fpx = fopen("cutx.txt", "w");
FILE * fpy = stderr;``````

The flux along $y$ is $\partial c/\partial y$

``````  scalar flux[];
foreach()
flux[] = (c[0,0] - c[0,-1])/Δ;
boundary ({flux});

for (double x = -L0/2 ; x < L0/2; x += L0/N) {
double y = x + L0/2;
fprintf (fpx, "%g %g %g\n",
x, interpolate (c, x, 0) , interpolate (flux, x, 0));
fprintf (fpy, "%g %g %g\n",
y, interpolate (c, 0, y) , interpolate (flux, 0, y));
}
fclose(fpx);
}``````

At the end of the simulation, we create an image of the error field, in PNG format.

``````event pictures (t = end) {
scalar  dce[];
foreach()
dce[] =  c[] - cexact(y);
boundary ({dce});
output_ppm (dce, file = "c.png", spread = 2, linear = true);
}``````

We compare the exact and computed solutions for a cross-section at $x=0$.

We also display the concentration flux divided by the reference flux $bi/\left(biL0+1\right)$.

And the error field.

## Bibliography

• N. Dupuis, J. Décobert, P.-Y. Lagrée, N. Lagay, D. Carpentier, F. Alexandre (2008): “Demonstration of planar thick InP layers by selective MOVPE”. Journal of Crystal Growth issue 23, 15 November 2008, Pages 4795-4798.

• N. Dupuis, J. Décobert, P.-Y. Lagrée , N. Lagay, C. Cuisin, F. Poingt, C. Kazmierski, A. Ramdane, A. Ougazzaden (2008): “Mask pattern interference in AlGaInAs MOVPE Selective Area Growth : experimental and modeling analysis”. Journal of Applied Physics 103, 113113 (2008).