# sandbox/bugs/flowrate.c

This is not a bug: periodic boundaries are indeed not “real” boundaries.

Here we try to compute the flowrate at the right boundary’s section of the tri-periodic (cubic) domain of a flow driven by a pressure drop.

Problem: when using the foreach_boundary(right) syntax (i.e with the compute_flowrate_at_boundary_right function defined bellow), there is no cell transversal and thus the flowrate remains 0.

When using a manual detection of the cells neighboring the right boundary, (i.e with compute_flowrate_basic_right function defined bellow) one can see that there is indeed a flowrate but this solution is not satisfactory as the location of the plane is not the boundary’s.

#define LEVEL 6
#include "grid/octree.h"
#include "navier-stokes/centered.h"

#define rhoval 1000.
#define muval 60.

#define mydt 0.00375
#define maxiter (20)

double deltau;
scalar un[];

int main() {

L0 = 1;

DT = mydt;

init_grid(1 << (LEVEL));

/* Tri-periodic flow driven by body force */

foreach_dimension()
periodic(top);

const face vector dp[] = {10./L0/rhoval, 0, 0};
a = dp;

run();
}

event init (i = 0) {
origin (0, 0, 0);

/* set dynamic viscosity */
const face vector muc[] = {muval, muval, muval};
mu = muc;

/* set density of the flow */
const scalar rhoc[] = rhoval;
rho = rhoc;

/* The flow is at rest initially. */
foreach(){
u.x[] = 0.;
un[] = u.x[];
}
}

void compute_flowrate_at_boundary_right(const double Vol, coord * vel) {

foreach_dimension()
vel->x = 0.;

boundary(all);
int k =0;
foreach_boundary(right){
k ++;
if (k == 1) printf("thread %d, x = %g",pid(),x);
foreach_dimension()
vel->x += sq(Delta)*u.x[];
}

foreach_dimension()
vel->x /= Vol;

#if _MPI
foreach_dimension(){
mpi_all_reduce(vel->x, MPI_DOUBLE, MPI_SUM);
}
#endif

}

void compute_flowrate_basic_right(const double Vol,  coord * supvelo){

Cache intDomain = {0};
Point lpoint;
double xboundary = L0;

// allocate cache for the integration domain (one slice of x, 2D domain)

foreach(){
if (x > (xboundary - Delta))
{
lpoint = locate(x, y, z);
cache_append(&intDomain, lpoint, 0);

}
}

cache_shrink(&intDomain);
printf("thread %d, Number of point in integration domain %d\n",pid(),intDomain.n);

foreach_dimension()
supvelo->x = 0.;

int k = 0;

foreach_cache(intDomain){
k ++;
if (k == 1) printf("thread %d, x = %g",pid(),x);
foreach_dimension()
supvelo->x += sq(Delta)*u.x[];
}

foreach_dimension()
supvelo->x /= Vol;

free(intDomain.p);

#if _MPI
foreach_dimension(){
mpi_all_reduce(supvelo->x, MPI_DOUBLE, MPI_SUM);
}
#endif

}

event flowrate(i++;i<maxiter){

coord flowrate;
double Volume = pow(L0,3);

compute_flowrate_at_boundary_right(Volume, &flowrate);
/* compute_flowrate_basic_right(Volume, &flowrate); */

fprintf(ferr,"%g %g %g %g %g\n",t, flowrate.x, flowrate.y, flowrate.z, Volume);

}
plot 'log' u 1:2 w l