sandbox/Antoonvh/robin.c
Mixed boundary conditions
From Yong Hui’s page we know that for a boundary condition of the form,
\displaystyle a \mathtt{s} + b \frac{\partial{\mathtt{s}}}{\partial \mathbf{n}} = c, The discretized version for 2b\neq -a\Delta reads, \displaystyle s [ \mathtt{ghost} ] = \left[ \frac{2 c \Delta}{2b + a\Delta} \right]+ \Bigg \langle \frac{2b - a\Delta }{2b + a\Delta} \Bigg \rangle s[\quad]
The ansatz is that this may be expressed as a linear mix of a dirichlet
and a Neumann
condition.
\displaystyle s[\mathtt{ghost}] = A\cdot \mathtt{dirichlet}(D) + B\cdot\mathtt{neumann}(N).
Expanding the macros gives:
\displaystyle s[\mathtt{ghost}] = 2AD - As[\quad]+ B\Delta N + Bs[ \quad] \displaystyle = \left[2AD + B\Delta N\right]+ \langle B-A\rangle s[\quad].
Each term at the rhs forms a seperate equation. Such that we have four unkowns (A,B,D,N) and only two equations. Such underdetermined system may have many solutions. Lets check if a solution exist for A = 1 and N = 0. The equation for the square-bracketed term then becomes:
\displaystyle 2D = \frac{2 c \Delta}{2b + a\Delta},\rightarrow \displaystyle D = \frac{c \Delta}{2b + a\Delta}.
Next, for the second, angle-bracketed-term equation,
\displaystyle B - 1 = \frac{2b - a\Delta }{2b + a\Delta} \rightarrow \displaystyle B = \frac{2b - a\Delta }{2b + a\Delta} + 1.
Wrapping it up (in many braces):
#define robin(a,b,c) ((dirichlet ((c)*Delta/(2*(b) + (a)*Delta))) + ((neumann (0))*((2*(b) - (a)*Delta)/(2*(b) + (a)*Delta) + 1.)))
Mind the devide-by-zero
risk.
A test
Vyaas Gururajan coded a neat test case. Special care must be taken, because it is not always critical.
set xlabel 'x'
set ylabel 'u'
set size square
set grid
plot 'out' u 1:3 w l lw 6 t 'Analytical Solution', '' w l lw 3 t 'Approximate Solution'
#include "grid/multigrid1D.h"
#include "poisson.h"
Here is the analaytical solution taken from Professor Fitzpatrick’s notes:
double alphaL = 2., betaL = -1., gammaL = -3.;
double alphaR = 1., betaR = 3., gammaR = -2.;
double analyticalU (double x) {
double d = alphaL*alphaR + alphaL*betaR - betaL*alphaR;
double g = (gammaL*(alphaR + betaR) - betaL*(gammaR - (alphaR+betaR)/3.))/d;
double h = (alphaL*(gammaR - (alphaR + betaR)/3.) - gammaL*alphaR)/d;
return (g + h*x + sq(x)/2. - sq(sq(x))/6.);
}
scalar u[], rhs[];
u[left] = robin (alphaL, -betaL, gammaL); //Mind the sign!
u[right] = robin (alphaR, betaR, gammaR);
int main() {
init_grid (1 << 7);
foreach()
rhs[] = 1. - 2.*sq(x);
poisson (u, rhs, tolerance = 1e-9);
foreach()
printf("%g %g %g\n",x, u[], analyticalU(x));
}