sandbox/Antoonvh/robin.c

Mixed boundary conditions

From Yong Hui’s page we know that for a boundary condition of the form,

a \mathtt{s} + b \frac{\partial{\mathtt{s}}}{\partial \mathbf{n}} = c, The discretized version for 2b\neq -a\Delta reads, 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.

s[\mathtt{ghost}] = A\cdot \mathtt{dirichlet}(D) + B\cdot\mathtt{neumann}(N).

Expanding the macros gives:

s[\mathtt{ghost}] = 2AD - As[\quad]+ B\Delta N + Bs[ \quad] = \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:

2D = \frac{2 c \Delta}{2b + a\Delta},\rightarrow D = \frac{c \Delta}{2b + a\Delta}.

Next, for the second, angle-bracketed-term equation,

B - 1 = \frac{2b - a\Delta }{2b + a\Delta} \rightarrow 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;

the dimensional constraints below are not compatible

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