sandbox/Antoonvh/robin.c

Victor Gustave Robin could have looked like Peter Gustav Lejeune Dirichlet (or Claude Louis Marie Henri Navier)

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'

Looks OK (script)

#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));
}