# hub.darcs.net :: basilisk -> basilisk -> files

Basilisk source code

## root / src / diffusion.h.page

 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104  /** # Time-implicit discretisation of reaction--diffusion equations We want to discretise implicitly the reaction--diffusion equation $$\theta\partial_tf = \nabla\cdot(D\nabla f) + \beta f + r$$ where $\beta f + r$ is a reactive term, $D$ is the diffusion coefficient and $\theta$ can be a density term. Using a time-implicit backward Euler discretisation, this can be written $$\theta\frac{f^{n+1} - f^{n}}{dt} = \nabla\cdot(D\nabla f^{n+1}) + \beta f^{n+1} + r$$ Rearranging the terms we get $$\nabla\cdot(D\nabla f^{n+1}) + (\beta - \frac{\theta}{dt}) f^{n+1} = - \frac{\theta}{dt}f^{n} - r$$ This is a Poisson--Helmholtz problem which can be solved with a multigrid solver. */ #include "poisson.h" /** The parameters of the diffusion() function are a scalar field f, scalar fields r and $\beta$ defining the reactive term, the timestep dt and a face vector field containing the diffusion coefficient D. If D or $\theta$ are omitted they are set to one. If $\beta$ is omitted it is set to zero. Both D and $\beta$ may be constant fields. Note that the r, $\beta$ and $\theta$ fields will be modified by the solver. The function returns the statistics of the Poisson solver. */ struct Diffusion { // mandatory scalar f; double dt; // optional face vector D; // default 1 scalar r, beta; // default 0 scalar theta; // default 1 }; trace mgstats diffusion (struct Diffusion p) { /** If *dt* is zero we don't do anything. */ if (p.dt == 0.) { mgstats s = {0}; return s; } /** We define $f$ and $r$ for convenience. */ scalar f = p.f, r = automatic (p.r); /** We define a (possibly constant) field equal to $\theta/dt$. */ const scalar idt[] = - 1./p.dt; (const) scalar theta_idt = p.theta.i ? p.theta : idt; if (p.theta.i) { scalar theta_idt = p.theta; foreach() theta_idt[] *= idt[]; } /** We use r to store the r.h.s. of the Poisson--Helmholtz solver. */ if (p.r.i) foreach() r[] = theta_idt[]*f[] - r[]; else // r was not passed by the user foreach() r[] = theta_idt[]*f[]; /** If $\beta$ is provided, we use it to store the diagonal term $\lambda$. */ scalar lambda = theta_idt; if (p.beta.i) { scalar beta = p.beta; foreach() beta[] += theta_idt[]; lambda = beta; boundary ({lambda}); } /** Finally we solve the system. */ return poisson (f, r, p.D, lambda); }