sandbox/Antoonvh/profile6.h
Flexible solution profiling functions
It may be used to generate profiles along arbitraty interfaces.
e.g.
profile (all, sqrt(sq(x - X0) + sq(y - Y0)), "r_prof"); writes azimutially averaged profiles of the solution along the radial direction (origin centered) to a file named `r_prof’.
Another e.g.
profile_equi (all, x, 10, "x_prof"); writes 10 y-z averaged solution values along the x direction a to a file named `x_prof’.
It works on trees, in serial, with MPI or with openmp if compiled correctly…
The functionalities rely on some helpful functions via fractions.h.
#include "fractions.h"The first function computes the averages of the scalar fields in a list, along the interface that is defined by the corresponding volume fraction field. The function returns the average grid size of the cells.
struct ob_av {
scalar * list; // List of scalar fields (mandatory)
double * v; // array for the values (mandatory)
scalar c; // Volume fractions for the interface (mandatory)
face vector cs; // Corresponding face fraction field (optional)
};
#if EMBED
double embed_interpolate2 (Point point, scalar s, coord p); //prototype
#endif
double interface_average (struct ob_av oa){
double area = 0, gs = 0;
int j = 0, g = 0, m = list_len (oa.list);
NOT_UNUSED(m);
scalar c = oa.c;
double * data = oa.v;
for (scalar s in oa.list)
oa.v[g++] = 0.;
foreach(reduction(+:area) reduction(+:j)
reduction(+:gs) reduction (+:data[:m])) {
if (c[] > 1e-5 && c[] < 1. - 1e-5) {
gs += Delta;
j++;
coord n, p[1];
if (oa.cs.x.i)
n = facet_normal (point, c, oa.cs);
else
n = mycs (point, c);
double alpha = plane_alpha (c[], n);
double ar = plane_area_center(n, alpha, p)/Delta;
foreach_dimension()
ar *= Delta;
area += ar;
int g = 0;
for (scalar s in oa.list) {
//s[];
//#if !EMBED
data[g++] += ar*interpolate(s, x + Delta*p->x,
y + Delta*p->y, z + Delta*p->z);
//#else
// data[g++] += ar*embed_interpolate2 (point, s, p[0]);
//#endif
}
}
}
g = 0;
if (area > 0){
for (scalar s in oa.list)
oa.v[g++] /= area;
}
if (j > 0)
return gs/(double)j;
else
return L0/100.;
}The function profiles uses the interface_average() function to loop over the domain. It therefore requires a phi vertex field.
struct prof {
scalar * list; // list of scalar field. The default is `all`
vertex scalar phi; // The phi coordinates field (default is: phi = y)
char * fname; // Optional file name
double max; // lower phi coordinate. The default is close to its minimum
double min; // upper phi coordinate. The default is close to its maximum
double rf; // reduction factor of query heights. The default is 1
FILE * fp; // File pointer, if `fname` is not provided. The default is `stdout`
int n; // Number of phi-values to query for equidistant profiles
char * func_str; // An identifier string (default is `phi`)
};
void profiles (struct prof p) {We start with setting a bunch of default settings. Such as the list includes all scalars, except phi.
if(!p.list)
p.list = all;
scalar * list = NULL ;
for (scalar s in p.list){
if (strcmp ("phi", s.name))
list = list_append (list, s);
}
vertex scalar phi;
if (p.phi.i)
phi = p.phi;For historic reasons, we make a vertical profile is phi is not provided.
else {
phi = new_vertex_scalar ("phi");
foreach_vertex ()
phi[] = y;
p.func_str = "y";
}The min and max coordinates are the extreme values, taken from the phi field with the addional maxdepth trick of Vincent Heusinkveld.
if (!p.min){
double mn = HUGE;
foreach_vertex(reduction(min:mn)){
if (phi[] < mn)
mn = phi[];
}
p.min = mn + L0/(double)(1 << (grid->maxdepth + 1));
}
if (!p.max){
double mx = -HUGE;
foreach_vertex(reduction(max:mx)){
if (phi[] > mx)
mx = phi[];
}
p.max = mx - L0/(double)(1 << (grid->maxdepth + 1));
}And by default there is no reduction factor and we print to the stdout file pointer.
if (!p.rf)
p.rf = 1;
if (!p.fname && !p.fp)
p.fp = stdout;
double dzn;
if (p.n)
dzn = (p.max - p.min) / ((double)p.n - 0.99999999);
if (!p.func_str)
p.func_str = "phi";
FILE * fp = p.fp;
char * file = p.fname;
if (pid()==0){
if (file && (fp = fopen (file, "w")) == NULL) {
perror (file);
exit (1);
}
assert (fp);
fprintf(fp,"#%s\t", p.func_str);
for(scalar s in list)
fprintf(fp,"%s\t",s.name);
fprintf(fp,"\n");
}
int len = list_len(list);
boundary(list);
double phi_p = p.min;
scalar s[];
face vector sf[];
while (phi_p < p.max){
double aver[len];
fractions (phi, s, sf, val = phi_p);
double dphi = interface_average (list, aver, s, sf);
if (pid() == 0){
int k = 0;
for (scalar s in list){
if (k == 0)
fprintf (fp, "%g\t%g", phi_p, aver[k]);
if (k > 0)
fprintf(fp,"\t%g",aver[k]);
k++;
if (k == len)
fprintf(fp,"\n");
}
}
if (p.n)
dphi = dzn/p.rf;Note that the units of phi and dphi may be incompatible.
Two addional user-interface functions
Because the phi field may be inconvienient, we make some user-friendly defenitions. Note that I do not know how to do optional arguments, so the user must provide all arguments. It relies on the default settings of the profiles() function above.
For grid-adaptive profiles use profile():
#define profile(list, fun, fname) do { \
vertex scalar phi[]; \
foreach_vertex() \
phi[] = fun; \
char * str = #fun; \
profiles (list, phi, fname, func_str = str); \
} while(0);And for equidistant profiles use profile_equi():
#define profile_equi(list, fun, nr, fname) do { \
vertex scalar phi[]; \
foreach_vertex() \
phi[] = fun; \
char * str = #fun; \
profiles (list, phi, fname, n = nr, func_str = str); \
} while(0);Embedded boundary
A 3D-compatible interpolate near an embedded boundary
#if EMBED
double embed_interpolate2 (Point point, scalar s, coord p) {
int i = sign(p.x), j = sign(p.y);
#if dimension == 2
if (cs[i] && cs[0,j] && cs[i,j])
// bilinear interpolation when all neighbors are defined
return ((s[]*(1. - fabs(p.x)) + s[i]*fabs(p.x))*(1. - fabs(p.y)) +
(s[0,j]*(1. - fabs(p.x)) + s[i,j]*fabs(p.x))*fabs(p.y));
#else //dimension == 3, see cartesian-common.h
int k = sign(p.z);
x = fabs(p.x); y = fabs(p.y); z = fabs(p.z);
/* trilinear interpolation */
if (cs[i] && cs[0,j] && cs[i,j] && cs[0,0,k] &&
cs[i,0,k] && cs[0,j,k] && cs[i,j,k]) {
return (((s[]*(1. - x) + s[i]*x)*(1. - y) +
(s[0,j]*(1. - x) + s[i,j]*x)*y)*(1. - z) +
((s[0,0,k]*(1. - x) + s[i,0,k]*x)*(1. - y) +
(s[0,j,k]*(1. - x) + s[i,j,k]*x)*y)*z);
}
#endif
else {
// linear interpolation with gradients biased toward the
// cells which are defined
double val = s[];
foreach_dimension() {
int i = sign(p.x);
if (cs[i])
val += fabs(p.x)*(s[i] - s[]);
else if (cs[-i])
val += fabs(p.x)*(s[] - s[-i]);
}
return val;
}
}
#endif