line_integrate.h

Computing a line integral
Option 1
Option 2
Option 3

Computing a line integral

We want to compute S, defined as

\displaystyle S(x, y_1,z) = \int_{y_1}^{\mathrm{top}} s \mathrm{d}y.

Altough this is easy on a cartesian grid (because of the particular allignment), a tree grid may pose a challenge. Here we present a few approaches.

Note that the equation implies that S(x, \mathrm{top}, z) = 0, giving a boundary condition at the top boundary:

S[top] = dirichlet (0.); // or just 0?

On this page we present a few iterative strategies that make use of multigrid acceleration.

#include "poisson.h"

Option 1

Solve:

\displaystyle \frac{\partial a}{\partial y} = -b,

via,

\displaystyle \frac{\partial^2 a}{\partial y^2} = -\frac{\partial b}{\partial y},

Using the Poisson solver with \alpha = \{0, 1, 0\}). Mind the proper boundary conditions. Note that this approach is not very practical because the finite TOLERANCE has an indirect relation to the error in S, and the BC on the derivative becomes very important.

struct Integrate_dy {
  scalar a, b;
  (const) scalar lambda;
  double tolerance;
  int nrelax, minlevel;
  scalar * res;
};

mgstats integrate_dy (struct Integrate_dy p) {
  scalar a = p.a, b = p.b;
  boundary ({b});
  scalar dbdy[];
  foreach()
    dbdy[] = -(b[0,1] - b[0,-1])/(2.*Delta);
  const face vector alpha[] = {0., 1., 0.};
  return poisson (a, dbdy, alpha, p.lambda, p.tolerance,
		  p.nrelax, p.minlevel, p.res);
}

Option 2

Alternatively, we discretize

\displaystyle \frac{\partial a}{\partial y} = -b,

with 2nd-order accuracy as;

(a[0,1] - a[0,-1])/2./Delta = -b[];

The residual function is:

static double residual_int (scalar * al, scalar * bl, scalar * resl, void * data) {
  scalar a = al[0], b = bl[0], res = resl[0];
  double maxres = 0.;
#if TREE
  /* conservative coarse/fine discretisation (2nd order)?? */
  face vector g[];
  foreach_face(y) {
    g.y[] = face_gradient_y (a, 0);
  }
  boundary_flux ({g});
  foreach (reduction(max:maxres)) {
    res[] = b[] + (g.y[0,1] + g.y[])/2.;
#else // !TREE
    foreach (reduction(max:maxres)) { 
      res[] = -b[] - (a[0,1] - a[0,-1])/(2.*Delta);
#endif // !TREE    
      if (fabs (res[]) > maxres) 
	maxres = fabs (res[]);
      //printf ("%g %g %g %g\n", x, y, res[], -b[]);
    }
    boundary (resl);
    printf ("%g \n", maxres);
    return maxres;
  }

This centered derivative formulation does not provide a suggestion for the relaxation of a[]. But perhaps it does for its neighbors. If they are not ghosts, we modify both so that the discrete system is statisfied.

a[0,1]  <- (a[0,1]  - b[]*2*Delta + a[0,-1])/2.;
a[0,-1] <- (a[0,-1] + b[]*2*Delta + a[0,1]) /2.;

We use a marker field to indicate ghosts.

scalar marker[];
marker[bottom] = nodata;
marker[top] = nodata;
 
static inline void refine_nodata (Point point, scalar s) {
  foreach_child()
    s[] = nodata;
}
 
 static void relax_int (scalar * al, scalar * bl, int l, void * data) {
   scalar a = al[0], b = bl[0];
#if JACOBI
  scalar c[];
#else
  scalar c = a;
#endif
  foreach_level_or_leaf (l) {
    if (marker[0,1] != nodata && marker[0,-1] != nodata) {
      double a01 = a[0,1]; 
      c[0,1]  = (a[0,1]  - (b[]*2*Delta + a[0,-1]))/2.;
      c[0,-1] = (a[0,-1] + (b[]*2*Delta + a01))    /2.;
      //c[0,1]  = (- b[]*2*Delta + a[0,-1]);
      //c[0,-1] = (  b[]*2*Delta + a01);
    
    }

If a neighbor is a ghosts, we only modify the other.

a[0, +/- 1] = -/+ b[]*2*Delta + a[0, -/+ 1];

    else if (marker[0,1] != nodata) {
      c[0,1]  = -b[]*2*Delta + a[0,-1];
    } else {
      c[0,-1] = b[]*2*Delta + a[0,1];
    }
  }
#if JACOBI
  foreach_level_or_leaf (l)
    a[] = (a[] + 2.*c[])/3.;
#endif
}

We also provinde a user-interface to this iterative Multigrid-accelerated method. Again, it does not work.

struct Integrate_dn {
  scalar a, b;
  (const) face vector alpha;
  (const) scalar lambda;
  double tolerance;
  int nrelax, minlevel;
  scalar * res;
};

mgstats integrate_dn (struct Integrate_dn p) {
  scalar a = p.a, b = p.b;
  marker.prolongation = refine_nodata;
  boundary ({marker});
  return mg_solve({a}, {b}, residual_int, relax_int,
		  &p, p.nrelax, p.res,
		  tolerance = p.tolerance, minlevel = p.minlevel);
}

Option 3

We may also discretize with first-order accuracy:

(a[0,1] - a[])/Delta = -b[];

The residual is;

  static double residual_1st (scalar * al, scalar * bl, scalar * resl, void * data) {
    scalar a = al[0], b = bl[0], res = resl[0];
    double maxres = 0.;
#if TREE
    /* conservative coarse/fine discretisation (2nd order)?? */
    face vector g[];
    foreach_face(y) {
      g.y[] = face_gradient_y (a, 0);
    }
    boundary_flux ({g});
  foreach (reduction(max:maxres)) {
    res[] = b[] + (g.y[0,1]);
#else // !TREE
    foreach (reduction(max:maxres)) { 
      res[] = b[] + (a[0,1] - a[0])/(Delta);
#endif // !TREE    
      if (fabs (res[]) > maxres) 
	maxres = fabs (res[]);
      //printf ("%g %g %g %g\n", x, y, res[], -b[]);
    }
    boundary (resl);
    return maxres;
  }

And the relaxation procedure:

  static void relax_1st (scalar * al, scalar * bl, int l, void * data) {
    scalar a = al[0], b = bl[0];
#if JACOBI
    scalar c[];
#else
    scalar c = a;
#endif
    foreach_level_or_leaf(l) 
      c[]  = b[]*Delta + a[0,1];
    
    
#if JACOBI
  foreach_level_or_leaf (l)
    a[] = (a[] + 2.*c[])/3.;
#endif
  }

Due to the biased nature of the stencil, the convergence properties become sensitive to the cell-iteration sequence. For cases with a top boundary condition for a, it appears we benefit from a top-down iterator for the relaxation.

Using this solution, this method works quite OK.

#include "topdown-iterator.h"
  
  trace
    mgstats integrate_1st (struct Integrate_dn p) {
    update_cache_f2();
    scalar a = p.a, b = p.b;
    return mg_solve({a}, {b}, residual_1st, relax_1st,
		    &p, p.nrelax, p.res,
		    tolerance = p.tolerance, minlevel = p.minlevel);
    update_cache_f();
  }