#include "grid/bitree.h"
#include "poisson.h"
#include "utils.h"

#define sol(x,c1,c2) ((c1*x)+c2+(pow(M_PI,0.5)*(x)*erf(x)/2)+(exp(-(x*x))/2))

double c1,c2,err;
scalar b[],a[];
mgstats mg;
FILE * fp1;
FILE * fp2;
char name[100];
double errzero;

double numintsol(double tol,double xi, double D,double c1,double c2){
  double into=5;
  double intn=10;
  double integral;
  double j=10;
  while ((fabs(into-intn)/fabs(intn))>tol){
    into=intn;
    intn=0.;
    integral=0;
    for (int m=0;m<j;m++){
      double xp=(xi-(D/2)+(D/(j*2))+(((double)m)*D/j));
      integral+=sol(xp,c1,c2);
    }
    intn=integral/j;
    j=j*2;
  }
  return intn;
}

double f(double xi,double Delta){
  return (pow(M_PI,0.5)*(erf(((xi)+(Delta/2)))-erf(((xi)-(Delta/2)))))/(2*Delta);
}

int main(){
  L0=10;
  X0=-L0/2;
  a[right]=dirichlet(0.);
  a[left]=neumann(0.);
  c1 = pow(M_PI,0.5)/2;
  c2 =-sol(5,c1,0);
  fp2=fopen("Convergencewithzeta","w");
  TOLERANCE=10; //Set large Tolerance
  int n=0;

  for (int j=0;j<20;j++){ 
    double zeta=1./pow(2.,(double)j);
    if (j==0) //Benchmark result
      zeta=0.;
    init_grid(1<<12);
    foreach(){
      a[]=0.;
      b[]=f(x,Delta);
    }
    while (adapt_wavelet({b},(double[]){zeta},12).nc){
      foreach(){
	a[]=0.;
	b[]=f(x,Delta);
      }
      boundary({a,b});
    }
    boundary(all);
    sprintf(name,"MGcycles%d.dat",j);
    fp1=fopen(name,"w");
    fprintf(ferr,"Grid j=%d\n",j);
    boundary({a,b});

    for (int i=0;i<10;i++){
      mg = poisson(a,b); //Solve the system
      err=0;
      double errz=0;
      n=0;
      foreach(){
	err+=fabs(a[]-numintsol(10e-5*pow(Delta,3.),x,Delta,c1,c2))*Delta;
	errz+=fabs(a[]-sol(x,c1,c2))*Delta;
	n++;
      }
      if(j==0)
	errzero=err;
      fprintf(fp1,"%d\t%d\t%g\t%g\t%g\t%d\t%g\t%g\n",i,mg.i,mg.resb,mg.resa,mg.sum,mg.nrelax,err,errz);
    }
    fclose(fp1);
    if(j>0)
      fprintf(fp2,"%d\t%g\t%g\t%d\n",j,zeta,fabs(err-errzero),n);
    fflush(fp2);
  }
  fclose(fp2);
}