profiles_foreach_region.h

Profile functions
- average_foreach_region(): Averages over a specified region.
- profile_foreach_region(): Generates profiles for a list of scalar fields and writes them to a file.
  - Usage

Profile functions

average_foreach_region(): Averages over a specified region.

This function calculates the average of a list of scalar fields within a specified region using foreach_region(). Fields are weighted by the field w. The function uses parallel reduction to accumulate values and weights, then normalizes the accumulated values.

\tilde{Q} \approx \begin{cases} \dfrac{\iint w(x,y,t) Q(x,y,t)~dx}{\iint w(x,y,t) ~dx } & \text{ if 2D} \\ \\ \dfrac{\iint w(x,y,t) Q(x,y,t)~dx~dy}{\iint w(x,y,t) ~dx~dy } & \text{ if 3D} \end{cases}

The arguments and their default values are:

list: Pointer to the list of scalar fields.
w: Scalar field containing the weights.
v: Array to store the calculated averages
hprof: Coordinate of the plane (default: 0.)
n: number of points along each dimension. Default is N.
linear: use first-order (default) or bilinear interpolation.
box: the lower-left and upper-right coordinates of the domain to consider. Default is the entire domain.

double average_foreach_region(scalar *list, scalar w, double *v, double hprof = 0., coord box[2] = {{X0, hprof}, {X0 + L0, hprof}}, coord n = {N, 1})
{

  int len = list_len(list); // Get length of the scalar list
  int m = 0;                // Counter for points within the y-range
  double a = 0 ;            // Accumulator for weights
  
  coord p;
  NOT_UNUSED(len);

  // Main loop: Iterate over the specified region
  foreach_region(p, box, n, reduction(+ : a) reduction(+ : m) reduction(+ : v[:len])){
    m++;

    // Calculate weight factor (cell volume)
    double b = 1. / Delta ;
    foreach_dimension(){
      b *= Delta;
    }
    
    if (w.i == unity.i){ // If no weights are provided    
      a += b;
    }
    else { // If weights are provided  
      a += interpolate_linear(point, w, p.x, p.y, p.z) * b;
    }

    
    // Accumulate weighted values for each scalar field      
    int k = 0;
    if (w.i == unity.i){ // If no weights are provided    
      for (scalar s in list){
        v[k++] += interpolate_linear(point, s, p.x, p.y, p.z) * b;
      }
    }
    else {// If weights are provided  
      for (scalar s in list){
        v[k++] += interpolate_linear(point, s, p.x, p.y, p.z) * interpolate_linear(point, w, p.x, p.y, p.z) * b;
      }
    }
  }

  // Normalize accumulated values
  for (int g = 0; g < len; g++){
    v[g] /= a;
  }

  // Return average cell size or its square root, depending on dimension
#if (dimension == 2)
  return a / (double)m;
#else
  return sqrt(a / (double)m);
#endif
}

profile_foreach_region(): Generates profiles for a list of scalar fields and writes them to a file.

This profiling function employs the same profiling strategy as the void profile() function presented here. However, this function utilizes foreach_region() and optional arguments.

A profile with n points between hmin and hmax is generated by averaging every field in list over m (or m*m if in 3D) equally spaced points. The resulting profile is then stored in a file named filename.

The arguments and their default values are:

list: List of scalar fields (default: all)
w: Weights field (default: unity)
filename: Output file name (default: profiles.asc)
hmin: Lower y (or z if 3D) coordinate (default: Y0 + _mindelta)
hmax: Upper y (or z if 3D) coordinate (default: Y0 + L0 - _mindelta)
rf: Reduction factor of query heights (default: 1)
n: number of points along the profile. Default is N.
m: number of points along each direction of sampled region. Default is N.

Usage

  scalar * list = {Y,u};
  profile_foreach_region(list, filename="profiles.asc");

which should write a file

#y@0	Y	u.x	u.y	
-0.499023438	1	0.203056828	-0.00047600522
-0.497070313	1	0.203056828	-0.00047600522
-0.495117188	1	0.203056828	-0.00047600522
-0.493164063	1	0.203056828	-0.00047600522
...

see, also example 1.

#define _mindelta (L0 / (double)(1 << (grid->maxdepth + 1)))
void profile_foreach_region(scalar *list = all,                                  
                            scalar w = unity,
                            char *filename = "profiles.asc",                                  
                            double hmin = Z0 + _mindelta,
                            double hmax = Z0 + L0 - _mindelta,
                            double xmin = X0,
                            double xmax = X0 + L0,
                            double ymin = Y0,
                            double ymax = Y0 + L0,
                            double rf = 1,
                            int n = N,
                            int m1 = N,
                            int m2 = N,
                            const char *mode = "a") {

  double deltahn = (hmax - hmin) / ((double)n - 0.99999999);

  FILE *fp = NULL;
  int len = list_len(list);
  
  // The primary worker (pid 0) handles file writing
  if (pid() == 0) {
    fp = fopen(filename, mode);
    if (fp == NULL) {
      perror(filename);
      exit(1);
    }

    // Write header
    fprintf(fp, "#y@%g\t", t);
    for (scalar s in list)
      fprintf(fp, "%s\t", s.name);
    fprintf(fp, "\n");
  }

  double hprof = hmin;
  
  // Iterate over different y-coordinates (in 2D) or z-coordinates (in 3D)
  while (hprof <= hmax) {
    // Define the region of interest
    #if dimension == 2
      coord box[2] = {{xmin, hprof}, {xmax, hprof}};
      coord nsamples = {m1, 1};
    #else
      coord box[2] = {{xmin, ymin, hprof}, {xmax, ymax, hprof}};
      coord nsamples = {m1, m2, 1};
    #endif

    // Calculate averages for the current region
    double aver[len];
    memset(&aver, 0, sizeof(aver));
    double deltah = average_foreach_region(list, w, aver, hprof, box, nsamples);

    // Write results to file (primary worker only)
    if (pid() == 0) {
      for (int k = 0; k < len; k++) {
        if (k == 0) {
          fprintf(fp, "%.9g\t%.9g", hprof, aver[k]);
        }
        else {
          fprintf(fp, "\t%.9g", aver[k]);
        }
      }
      fprintf(fp, "\n");
    }

    // Calculate next y- (or z-) coordinate
    deltah = deltahn / rf;
    hprof += rf * deltah;
  }

  // Close file (primary worker only)
  if (pid() == 0) {
    fprintf(fp, "\n");
    fflush(fp);
    if (fp != stdout)
      fclose(fp);
  }
}
#undef _mindelta