sandbox/Antoonvh/robin.c

    Victor Gustave Robin could have looked like Peter Gustav Lejeune Dirichlet

    Victor Gustave Robin could have looked like Peter Gustav Lejeune Dirichlet

    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)

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