rf4t.c

Higher-order accurate face-value prolongation

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:

   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:

   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)

5th order solenoidal refinement (script)

Divergence:

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