/**
# Higher-order accurate face-value prolongation  

 Here we test the methods for 4th order face refinement and Solenoidal
 refinement.
 */
#include "higher-order.h"
#include "poisson.h" //For projection()
#include "utils.h"   

face vector fv[];
scalar s[];

double max_divergence (face vector fv);//Prototype

int main() {
  L0 = 10;
  X0 = Y0 = -L0/2;
  foreach_dimension() {
    periodic (left);
    fv.x.prolongation = refine_face_4_x;
  }
  fv.x.refine = refine_face_solenoidal;
  
  /**
     First, the prolongation method is tested:
     
     ~~~gnuplot 4th order prolongation
     set logscale xy 2
     set xr [9:4000]
     set yr [1e-12:1e-2]
     set grid
     set xlabel 'N'
     set ylabel 'L_1 Error'
     set size square
     plot 'out' u 1:2, 1000*x**(-4)
     ~~~
  */
  for (int l = 4; l < 12; l++) {
    init_grid (1 << l);
    foreach_face() {
      coord f = {y, -x};
      fv.x[] =  f.x*exp(-sq(x) - sq(y));
    }
    boundary ((scalar*){fv});
    wavelet (fv.x, s);
    fprintf (stdout, "%d %g\n", N, normf(s).avg); 
  }
  /**
     Next, the Solenoidal refinement test is performed
     
     Accuracy:
     
     ~~~gnuplot 5th order solenoidal refinement
     set logscale xy 2
     set xr [9:4000]
     set yr [1e-13:1e-2]
     set grid
     set size square
     set xlabel 'N'
     set ylabel 'L_1 Error'
     plot 'log' u 1:2, 500000*x**(-5)
     ~~~
     
     Divergence:

     ~~~gnuplot The divergence is within the machine precision
     reset
     set logscale x
     set xr [9:4000]
     set yr [-1e-11:1e-11]
     set grid
     set xlabel 'N'
     set ylabel 'Max. absolute divergence'
     set size square
     plot 'log' u 1:4
     ~~~
     
     
  */
  NITERMIN = 25;
  for (int l = 4; l < 12; l++) {
    init_grid (1 << l);
    foreach_face() {
      coord f = {y, -x};
      fv.x[] =  f.x*exp(-sq(x) - sq(y));
    }
    project (fv, s);
    double md = max_divergence (fv); //approx zero
    double * fva = malloc (2*(sq(N))*sizeof(double));
    long int j = 0;
    foreach_face()
      fva[j++] = fv.x[]; //Array of reference values.
    
    unrefine (level == l - 1);
    refine (level < l);

    double err = 0;
    j = 0;
    foreach_face()
      err += fabs(fv.x[] - fva[j++]);
    fprintf (stderr, "%d %g %g %g\n",
	     N, err/sq(N), md, max_divergence (fv));
    free (fva);
  }
}

double max_divergence (face vector fv) {
  double max = 0; 
  foreach() {
    double d = 0;
    foreach_dimension()
      d += fv.x[1] - fv.x[];
    d /= Delta;
    if (fabs(d) > max)
      max = fabs(d);
  }
  return max;
}