sandbox/yonghui/smalltest/robin.c

Robin condition test

See how to use the patch from Thomas.

My defintion of Robin condition is a little bit different \displaystyle a \boldsymbol{u} + b \frac{\partial{\boldsymbol{u}}}{\partial \boldsymbol{n}} = c Write in discretized form at boundary \displaystyle v[ ghost ] = \frac{2 c \Delta + (2b - a\Delta) v [\quad]}{2b + a\Delta} The code I used:

double aax = bbx = 1.;

@define robin(aax,bbx,ccx) ((2.*(ccx)*Delta/(2.*(bbx)+(aax)*Delta)) +((2.*(bbx)-(aax)*Delta)/(2.*(bbx)+(aax)*Delta))*val(_s,0,0,0))

@define robin_homogeneous() (((2.*(bbx)-(aax)*Delta)/(2.*(bbx)+(aax)*Delta))*val(_s,0,0,0))

Until now you have to define the values of aax and bbx in the code, in order to make sure robin_homogeneous() works.

#include "grid/multigrid.h"
#include "axi.h"
#include "navier-stokes/centered.h"
#include "navier-stokes/swirl.h"
#include "view.h"
#define ROBIN 1

Boundary conditions

u.n[top] = dirichlet(0);
u.t[top] = dirichlet(0);
#if ROBIN
w[top] = robin(aax,bbx,fc3);
#else
w[top]   = dirichlet(aax);
#endif

Setup

FILE *fp;
double fa3,fb,fc3,AAA;

int main()
{
  fp=fopen("file1","w");
  N = 64;
  periodic(right);
  const face vector muc[] = {.1,.1};
  stokes = true;
  mu = muc;
  DT = 0.1;
  run();
}

event init (i = 0) {
  // a b c
  aax = 2.;
  bbx = 2.;
  fc3 =	2.;
}

event logfile (i += 10;i<400)
{
  fprintf(stderr,"%d %g\n",i,t);
  if (i > 380) {
    double wwt,ppt;
    for (double xxx=0.01; xxx <= L0 ;xxx += 0.01){
      wwt = interpolate(w,L0/2.,xxx);
      ppt = interpolate(p,L0/2.,xxx);
      fprintf(fp,"%g %g %g\n",xxx,wwt,ppt);
    }
  }
}