src/grid/tree-mpi.h
// this can be used to control the outputs of debug_mpi()
int debug_iteration = -1;
void debug_mpi (FILE * fp1);
typedef struct {
CacheLevel * halo; // ghost cell indices for each level
void * buf; // MPI buffer
MPI_Request r; // MPI request
int depth; // the maximum number of levels
int pid; // the rank of the PE
int maxdepth; // the maximum depth for this PE (= depth or depth + 1)
} Rcv;
typedef struct {
Rcv * rcv;
char * name;
int npid;
} RcvPid;
typedef struct {
RcvPid * rcv, * snd;
} SndRcv;
typedef struct {
Boundary parent;
SndRcv mpi_level, mpi_level_root, restriction;
Array * send, * receive; // which pids do we send to/receive from
} MpiBoundary;
static void cache_level_init (CacheLevel * c)
{
c->p = NULL;
c->n = c->nm = 0;
}
static void rcv_append (Point point, Rcv * rcv)
{
if (level > rcv->depth) {
qrealloc (rcv->halo, level + 1, CacheLevel);
for (int j = rcv->depth + 1; j <= level; j++)
cache_level_init (&rcv->halo[j]);
rcv->depth = level;
}
cache_level_append (&rcv->halo[level], point);
if (level > rcv->maxdepth)
rcv->maxdepth = level;
}
void rcv_print (Rcv * rcv, FILE * fp, const char * prefix)
{
for (int l = 0; l <= rcv->depth; l++)
if (rcv->halo[l].n > 0)
foreach_cache_level(rcv->halo[l], l)
fprintf (fp, "%s%g %g %g %d %d\n", prefix, x, y, z, rcv->pid, level);
}
static void rcv_free_buf (Rcv * rcv)
{
if (rcv->buf) {
prof_start ("rcv_pid_receive");
MPI_Wait (&rcv->r, MPI_STATUS_IGNORE);
free (rcv->buf);
rcv->buf = NULL;
prof_stop();
}
}
static void rcv_destroy (Rcv * rcv)
{
rcv_free_buf (rcv);
for (int i = 0; i <= rcv->depth; i++)
if (rcv->halo[i].n > 0)
free (rcv->halo[i].p);
free (rcv->halo);
}
static RcvPid * rcv_pid_new (const char * name)
{
RcvPid * r = qcalloc (1, RcvPid);
r->name = strdup (name);
return r;
}
static Rcv * rcv_pid_pointer (RcvPid * p, int pid)
{
assert (pid >= 0 && pid < npe());
int i;
for (i = 0; i < p->npid; i++)
if (pid == p->rcv[i].pid)
break;
if (i == p->npid) {
qrealloc (p->rcv, ++p->npid, Rcv);
Rcv * rcv = &p->rcv[p->npid-1];
rcv->pid = pid;
rcv->depth = rcv->maxdepth = 0;
rcv->halo = qmalloc (1, CacheLevel);
rcv->buf = NULL;
cache_level_init (&rcv->halo[0]);
}
return &p->rcv[i];
}
static void rcv_pid_append (RcvPid * p, int pid, Point point)
{
rcv_append (point, rcv_pid_pointer (p, pid));
}
static void rcv_pid_append_pids (RcvPid * p, Array * pids)
{
// appends the pids of @p to @pids without duplication
for (int i = 0; i < p->npid; i++) {
int pid = p->rcv[i].pid, j, * a;
for (j = 0, a = pids->p; j < pids->len/sizeof(int); j++,a++)
if (*a == pid)
break;
if (j == pids->len/sizeof(int))
array_append (pids, &pid, sizeof(int));
}
}
void rcv_pid_write (RcvPid * p, const char * name)
{
for (int i = 0; i < p->npid; i++) {
Rcv * rcv = &p->rcv[i];
char fname[80];
sprintf (fname, "%s-%d-%d", name, pid(), rcv->pid);
FILE * fp = fopen (fname, "w");
rcv_print (rcv, fp, "");
fclose (fp);
}
}
static void rcv_pid_print (RcvPid * p, FILE * fp, const char * prefix)
{
for (int i = 0; i < p->npid; i++)
rcv_print (&p->rcv[i], fp, prefix);
}
static void rcv_pid_destroy (RcvPid * p)
{
for (int i = 0; i < p->npid; i++)
rcv_destroy (&p->rcv[i]);
free (p->rcv);
free (p->name);
free (p);
}
static Boundary * mpi_boundary = NULL;
#define BOUNDARY_TAG(level) (level)
#define COARSEN_TAG(level) ((level) + 64)
#define REFINE_TAG() (128)
#define MOVED_TAG() (256)
void debug_mpi (FILE * fp1);
static void apply_bc (Rcv * rcv, scalar * list, scalar * listv,
vector * listf, int l, MPI_Status s)
{
double * b = rcv->buf;
foreach_cache_level(rcv->halo[l], l) {
for (scalar s in list) {
memcpy (&s[], b, sizeof(double)*s.block);
b += s.block;
}
for (vector v in listf)
foreach_dimension() {
memcpy (&v.x[], b, sizeof(double)*v.x.block);
b += v.x.block;
if (*b != nodata && allocated(1))
memcpy (&v.x[1], b, sizeof(double)*v.x.block);
b += v.x.block;
}
for (scalar s in listv) {
for (int i = 0; i <= 1; i++)
for (int j = 0; j <= 1; j++)
#if dimension == 3
for (int k = 0; k <= 1; k++) {
if (*b != nodata && allocated(i,j,k))
memcpy (&s[i,j,k], b, sizeof(double)*s.block);
b += s.block;
}
#else // dimension == 2
{
if (*b != nodata && allocated(i,j))
memcpy (&s[i,j], b, sizeof(double)*s.block);
b += s.block;
}
#endif // dimension == 2
}
}
size_t size = b - (double *) rcv->buf;
free (rcv->buf);
rcv->buf = NULL;
int rlen;
MPI_Get_count (&s, MPI_DOUBLE, &rlen);
if (rlen != size) {
fprintf (stderr,
"rlen (%d) != size (%ld), %d receiving from %d at level %d\n"
"Calling debug_mpi(NULL)...\n"
"Aborting...\n",
rlen, size, pid(), rcv->pid, l);
fflush (stderr);
debug_mpi (NULL);
MPI_Abort (MPI_COMM_WORLD, -2);
}
}
#ifdef TIMEOUT
static struct {
int count, source, tag;
const char * name;
} RecvArgs;
static void timeout (int sig, siginfo_t *si, void *uc)
{
fprintf (stderr,
"ERROR MPI_Recv \"%s\" (count = %d, source = %d, tag = %d):\n"
"Time out\n"
"Calling debug_mpi(NULL)...\n"
"Aborting...\n",
RecvArgs.name, RecvArgs.count, RecvArgs.source, RecvArgs.tag);
fflush (stderr);
debug_mpi (NULL);
MPI_Abort (MPI_COMM_WORLD, -1);
}
#endif // TIMEOUT
static void mpi_recv_check (void * buf, int count, MPI_Datatype datatype,
int source, int tag,
MPI_Comm comm, MPI_Status * status,
const char * name)
{
#ifdef TIMEOUT
RecvArgs.count = count;
RecvArgs.source = source;
RecvArgs.tag = tag;
RecvArgs.name = name;
timer_t timerid;
extern double t;
if (t > TIMEOUT) {
struct sigaction sa;
sa.sa_flags = SA_SIGINFO;
sa.sa_sigaction = timeout;
sigemptyset(&sa.sa_mask);
assert (sigaction(SIGRTMIN, &sa, NULL) != -1);
struct sigevent sev;
sev.sigev_notify = SIGEV_SIGNAL;
sev.sigev_signo = SIGRTMIN;
sev.sigev_value.sival_ptr = &timerid;
assert (timer_create(CLOCK_REALTIME, &sev, &timerid) != -1);
struct itimerspec its;
its.it_value.tv_sec = 10;
its.it_value.tv_nsec = 0;
its.it_interval.tv_sec = 0;
its.it_interval.tv_nsec = 0;
assert (timer_settime(timerid, 0, &its, NULL) != -1);
}
#endif // TIMEOUT
int errorcode = MPI_Recv (buf, count, datatype, source, tag, comm, status);
if (errorcode != MPI_SUCCESS) {
char string[MPI_MAX_ERROR_STRING];
int resultlen;
MPI_Error_string (errorcode, string, &resultlen);
fprintf (stderr,
"ERROR MPI_Recv \"%s\" (count = %d, source = %d, tag = %d):\n%s\n"
"Calling debug_mpi(NULL)...\n"
"Aborting...\n",
name, count, source, tag, string);
fflush (stderr);
debug_mpi (NULL);
MPI_Abort (MPI_COMM_WORLD, -1);
}
#ifdef TIMEOUT
if (t > TIMEOUT)
assert (timer_delete (timerid) != -1);
#endif
}
trace
static int mpi_waitany (int count, MPI_Request array_of_requests[], int *indx,
MPI_Status *status)
{
return MPI_Waitany (count, array_of_requests, indx, status);
}
static int list_lenb (scalar * list) {
int len = 0;
for (scalar s in list)
len += s.block;
return len;
}
static int vectors_lenb (vector * list) {
int len = 0;
for (vector v in list)
len += v.x.block;
return len;
}
static void rcv_pid_receive (RcvPid * m, scalar * list, scalar * listv,
vector * listf, int l)
{
if (m->npid == 0)
return;
prof_start ("rcv_pid_receive");
int len = list_lenb (list) + 2*dimension*vectors_lenb (listf) +
(1 << dimension)*list_lenb (listv);
MPI_Request r[m->npid];
Rcv * rrcv[m->npid]; // fixme: using NULL requests should be OK
int nr = 0;
for (int i = 0; i < m->npid; i++) {
Rcv * rcv = &m->rcv[i];
if (l <= rcv->depth && rcv->halo[l].n > 0) {
assert (!rcv->buf);
rcv->buf = malloc (sizeof (double)*rcv->halo[l].n*len);
#if 0
fprintf (stderr, "%s receiving %d doubles from %d level %d\n",
m->name, rcv->halo[l].n*len, rcv->pid, l);
fflush (stderr);
#endif
#if 1 /* initiate non-blocking receive */
MPI_Irecv (rcv->buf, rcv->halo[l].n*len, MPI_DOUBLE, rcv->pid,
BOUNDARY_TAG(l), MPI_COMM_WORLD, &r[nr]);
rrcv[nr++] = rcv;
#else /* blocking receive (useful for debugging) */
MPI_Status s;
mpi_recv_check (rcv->buf, rcv->halo[l].n*len, MPI_DOUBLE, rcv->pid,
BOUNDARY_TAG(l), MPI_COMM_WORLD, &s, "rcv_pid_receive");
apply_bc (rcv, list, listf, listv, l, s);
#endif
}
}
/* non-blocking receives (does nothing when using blocking receives) */
if (nr > 0) {
int i;
MPI_Status s;
mpi_waitany (nr, r, &i, &s);
while (i != MPI_UNDEFINED) {
Rcv * rcv = rrcv[i];
assert (l <= rcv->depth && rcv->halo[l].n > 0);
assert (rcv->buf);
apply_bc (rcv, list, listv, listf, l, s);
mpi_waitany (nr, r, &i, &s);
}
}
prof_stop();
}
trace
static void rcv_pid_wait (RcvPid * m)
{
/* wait for completion of send requests */
for (int i = 0; i < m->npid; i++)
rcv_free_buf (&m->rcv[i]);
}
static void rcv_pid_send (RcvPid * m, scalar * list, scalar * listv,
vector * listf, int l)
{
if (m->npid == 0)
return;
prof_start ("rcv_pid_send");
int len = list_lenb (list) + 2*dimension*vectors_lenb (listf) +
(1 << dimension)*list_lenb (listv);
/* send ghost values */
for (int i = 0; i < m->npid; i++) {
Rcv * rcv = &m->rcv[i];
if (l <= rcv->depth && rcv->halo[l].n > 0) {
assert (!rcv->buf);
rcv->buf = malloc (sizeof (double)*rcv->halo[l].n*len);
double * b = rcv->buf;
foreach_cache_level(rcv->halo[l], l) {
for (scalar s in list) {
memcpy (b, &s[], sizeof(double)*s.block);
b += s.block;
}
for (vector v in listf)
foreach_dimension() {
memcpy (b, &v.x[], sizeof(double)*v.x.block);
b += v.x.block;
if (allocated(1))
memcpy (b, &v.x[1], sizeof(double)*v.x.block);
else
*b = nodata;
b += v.x.block;
}
for (scalar s in listv) {
for (int i = 0; i <= 1; i++)
for (int j = 0; j <= 1; j++)
#if dimension == 3
for (int k = 0; k <= 1; k++) {
if (allocated(i,j,k))
memcpy (b, &s[i,j,k], sizeof(double)*s.block);
else
*b = nodata;
b += s.block;
}
#else // dimension == 2
{
if (allocated(i,j))
memcpy (b, &s[i,j], sizeof(double)*s.block);
else
*b = nodata;
b += s.block;
}
#endif // dimension == 2
}
}
#if 0
fprintf (stderr, "%s sending %d doubles to %d level %d\n",
m->name, rcv->halo[l].n*len, rcv->pid, l);
fflush (stderr);
#endif
MPI_Isend (rcv->buf, (b - (double *) rcv->buf),
MPI_DOUBLE, rcv->pid, BOUNDARY_TAG(l), MPI_COMM_WORLD,
&rcv->r);
}
}
prof_stop();
}
static void rcv_pid_sync (SndRcv * m, scalar * list, int l)
{
scalar * listr = NULL, * listv = NULL;
vector * listf = NULL;
for (scalar s in list)
if (!is_constant(s) && s.block > 0) {
if (s.face)
listf = vectors_add (listf, s.v);
else if (s.restriction == restriction_vertex)
listv = list_add (listv, s);
else
listr = list_add (listr, s);
}
rcv_pid_send (m->snd, listr, listv, listf, l);
rcv_pid_receive (m->rcv, listr, listv, listf, l);
rcv_pid_wait (m->snd);
free (listr);
free (listf);
free (listv);
}
static void snd_rcv_destroy (SndRcv * m)
{
rcv_pid_destroy (m->rcv);
rcv_pid_destroy (m->snd);
}
static void snd_rcv_init (SndRcv * m, const char * name)
{
char s[strlen(name) + 5];
strcpy (s, name);
strcat (s, ".rcv");
m->rcv = rcv_pid_new (s);
strcpy (s, name);
strcat (s, ".snd");
m->snd = rcv_pid_new (s);
}
static void mpi_boundary_destroy (Boundary * b)
{
MpiBoundary * m = (MpiBoundary *) b;
snd_rcv_destroy (&m->mpi_level);
snd_rcv_destroy (&m->mpi_level_root);
snd_rcv_destroy (&m->restriction);
array_free (m->send);
array_free (m->receive);
free (m);
}
trace
static void mpi_boundary_level (const Boundary * b, scalar * list, int l)
{
MpiBoundary * m = (MpiBoundary *) b;
rcv_pid_sync (&m->mpi_level, list, l);
rcv_pid_sync (&m->mpi_level_root, list, l);
}
trace
static void mpi_boundary_restriction (const Boundary * b, scalar * list, int l)
{
MpiBoundary * m = (MpiBoundary *) b;
rcv_pid_sync (&m->restriction, list, l);
}
void mpi_boundary_new()
{
mpi_boundary = (Boundary *) qcalloc (1, MpiBoundary);
mpi_boundary->destroy = mpi_boundary_destroy;
mpi_boundary->level = mpi_boundary_level;
mpi_boundary->restriction = mpi_boundary_restriction;
MpiBoundary * mpi = (MpiBoundary *) mpi_boundary;
snd_rcv_init (&mpi->mpi_level, "mpi_level");
snd_rcv_init (&mpi->mpi_level_root, "mpi_level_root");
snd_rcv_init (&mpi->restriction, "restriction");
mpi->send = array_new();
mpi->receive = array_new();
add_boundary (mpi_boundary);
}
static FILE * fopen_prefix (FILE * fp, const char * name, char * prefix)
{
if (fp) {
sprintf (prefix, "%s-%d ", name, pid());
return fp;
}
else {
strcpy (prefix, "");
char fname[80];
if (debug_iteration >= 0)
sprintf (fname, "%s-%d-%d", name, debug_iteration, pid());
else
sprintf (fname, "%s-%d", name, pid());
return fopen (fname, "w");
}
}
void debug_mpi (FILE * fp1)
{
void output_cells_internal (FILE * fp);
char prefix[80];
FILE * fp;
// cleanup
if (fp1 == NULL) {
char name[80];
sprintf (name, "halo-%d", pid()); remove (name);
sprintf (name, "cells-%d", pid()); remove (name);
sprintf (name, "faces-%d", pid()); remove (name);
sprintf (name, "vertices-%d", pid()); remove (name);
sprintf (name, "neighbors-%d", pid()); remove (name);
sprintf (name, "mpi-level-rcv-%d", pid()); remove (name);
sprintf (name, "mpi-level-snd-%d", pid()); remove (name);
sprintf (name, "mpi-level-root-rcv-%d", pid()); remove (name);
sprintf (name, "mpi-level-root-snd-%d", pid()); remove (name);
sprintf (name, "mpi-restriction-rcv-%d", pid()); remove (name);
sprintf (name, "mpi-restriction-snd-%d", pid()); remove (name);
sprintf (name, "mpi-border-%d", pid()); remove (name);
sprintf (name, "exterior-%d", pid()); remove (name);
sprintf (name, "depth-%d", pid()); remove (name);
sprintf (name, "refined-%d", pid()); remove (name);
}
// local halo
fp = fopen_prefix (fp1, "halo", prefix);
for (int l = 0; l < depth(); l++)
foreach_halo (prolongation, l)
foreach_child()
fprintf (fp, "%s%g %g %g %d\n", prefix, x, y, z, level);
if (!fp1)
fclose (fp);
if (!fp1) {
fp = fopen_prefix (fp1, "cells", prefix);
output_cells_internal (fp);
fclose (fp);
}
fp = fopen_prefix (fp1, "faces", prefix);
foreach_face()
fprintf (fp, "%s%g %g %g %d\n", prefix, x, y, z, level);
if (!fp1)
fclose (fp);
fp = fopen_prefix (fp1, "vertices", prefix);
foreach_vertex()
fprintf (fp, "%s%g %g %g %d\n", prefix, x, y, z, level);
if (!fp1)
fclose (fp);
fp = fopen_prefix (fp1, "neighbors", prefix);
foreach() {
int n = 0;
foreach_neighbor(1)
if (is_refined(cell))
n++;
fprintf (fp, "%s%g %g %g %d\n", prefix, x, y, z, cell.neighbors);
assert (cell.neighbors == n);
}
if (!fp1)
fclose (fp);
MpiBoundary * mpi = (MpiBoundary *) mpi_boundary;
fp = fopen_prefix (fp1, "mpi-level-rcv", prefix);
rcv_pid_print (mpi->mpi_level.rcv, fp, prefix);
if (!fp1)
fclose (fp);
fp = fopen_prefix (fp1, "mpi-level-root-rcv", prefix);
rcv_pid_print (mpi->mpi_level_root.rcv, fp, prefix);
if (!fp1)
fclose (fp);
fp = fopen_prefix (fp1, "mpi-restriction-rcv", prefix);
rcv_pid_print (mpi->restriction.rcv, fp, prefix);
if (!fp1)
fclose (fp);
fp = fopen_prefix (fp1, "mpi-level-snd", prefix);
rcv_pid_print (mpi->mpi_level.snd, fp, prefix);
if (!fp1)
fclose (fp);
fp = fopen_prefix (fp1, "mpi-level-root-snd", prefix);
rcv_pid_print (mpi->mpi_level_root.snd, fp, prefix);
if (!fp1)
fclose (fp);
fp = fopen_prefix (fp1, "mpi-restriction-snd", prefix);
rcv_pid_print (mpi->restriction.snd, fp, prefix);
if (!fp1)
fclose (fp);
fp = fopen_prefix (fp1, "mpi-border", prefix);
foreach_cell() {
if (is_border(cell))
fprintf (fp, "%s%g %g %g %d %d %d\n",
prefix, x, y, z, level, cell.neighbors, cell.pid);
else
continue;
if (is_leaf(cell))
continue;
}
if (!fp1)
fclose (fp);
fp = fopen_prefix (fp1, "exterior", prefix);
foreach_cell() {
if (!is_local(cell))
fprintf (fp, "%s%g %g %g %d %d %d %d\n",
prefix, x, y, z, level, cell.neighbors,
cell.pid, cell.flags & leaf);
#if 0
else if (is_active(cell) && !is_border(cell))
continue;
if (is_leaf(cell))
continue;
#endif
}
if (!fp1)
fclose (fp);
fp = fopen_prefix (fp1, "depth", prefix);
fprintf (fp, "depth: %d %d\n", pid(), depth());
fprintf (fp, "======= mpi_level.snd ======\n");
RcvPid * snd = mpi->mpi_level.snd;
for (int i = 0; i < snd->npid; i++)
fprintf (fp, "%d %d %d\n", pid(), snd->rcv[i].pid, snd->rcv[i].maxdepth);
fprintf (fp, "======= mpi_level.rcv ======\n");
snd = mpi->mpi_level.rcv;
for (int i = 0; i < snd->npid; i++)
fprintf (fp, "%d %d %d\n", pid(), snd->rcv[i].pid, snd->rcv[i].maxdepth);
if (!fp1)
fclose (fp);
fp = fopen_prefix (fp1, "refined", prefix);
foreach_cache (tree->refined)
fprintf (fp, "%s%g %g %g %d\n", prefix, x, y, z, level);
if (!fp1)
fclose (fp);
}
static void snd_rcv_free (SndRcv * p)
{
char name[strlen(p->rcv->name) + 1];
strcpy (name, p->rcv->name);
rcv_pid_destroy (p->rcv);
p->rcv = rcv_pid_new (name);
strcpy (name, p->snd->name);
rcv_pid_destroy (p->snd);
p->snd = rcv_pid_new (name);
}
static bool is_root (Point point)
{
if (is_refined(cell))
foreach_child()
if (is_local(cell))
return true;
return false;
}
// see src/figures/prolongation.svg
static bool is_local_prolongation (Point point, Point p)
{
#if dimension == 2
struct { int x, y; } dp = {p.i - point.i, p.j - point.j};
#elif dimension == 3
struct { int x, y, z; } dp = {p.i - point.i, p.j - point.j, p.k - point.k};
#endif
foreach_dimension() {
if (dp.x == 0 && (is_refined(neighbor(-1)) || is_refined(neighbor(1))))
return true;
if (is_refined(neighbor(dp.x)))
return true;
}
return false;
}
#define is_remote(cell) (cell.pid >= 0 && cell.pid != pid())
static void append_pid (Array * pids, int pid)
{
for (int i = 0, * p = (int *) pids->p; i < pids->len/sizeof(int); i++, p++)
if (*p == pid)
return;
array_append (pids, &pid, sizeof(int));
}
static int locals_pids (Point point, Array * pids)
{
if (is_leaf(cell)) { // prolongation
if (is_local(cell)) {
Point p = point;
foreach_neighbor(1) {
if (is_remote(cell) &&
(is_refined(cell) || is_local_prolongation (point, p)))
append_pid (pids, cell.pid);
if (is_refined(cell))
foreach_child()
if (is_remote(cell))
append_pid (pids, cell.pid);
}
}
}
else
foreach_neighbor(1) {
if (is_remote(cell))
append_pid (pids, cell.pid);
if (is_refined(cell))
foreach_child()
if (is_remote(cell))
append_pid (pids, cell.pid);
}
return pids->len/sizeof(int);
}
static int root_pids (Point point, Array * pids)
{
foreach_child()
if (is_remote(cell))
append_pid (pids, cell.pid);
return pids->len/sizeof(int);
}
// turns on tree_check() and co with outputs controlled by the condition
// #define DEBUGCOND (pid() >= 300 && pid() <= 400 && t > 0.0784876)
// turns on tree_check() and co without any outputs
// #define DEBUGCOND false
static void rcv_pid_row (RcvPid * m, int l, int * row)
{
for (int i = 0; i < npe(); i++)
row[i] = 0;
for (int i = 0; i < m->npid; i++) {
Rcv * rcv = &m->rcv[i];
if (l <= rcv->depth && rcv->halo[l].n > 0)
row[rcv->pid] = rcv->halo[l].n;
}
}
void check_snd_rcv_matrix (SndRcv * sndrcv, const char * name)
{
int maxlevel = depth();
mpi_all_reduce (maxlevel, MPI_INT, MPI_MAX);
int * row = qmalloc (npe(), int);
for (int l = 0; l <= maxlevel; l++) {
int status = 0;
if (pid() == 0) { // master
// send/receive matrix i.e.
// send[i][j] = number of points sent/received from i to j
int ** send = matrix_new (npe(), npe(), sizeof(int));
int ** receive = matrix_new (npe(), npe(), sizeof(int));
rcv_pid_row (sndrcv->snd, l, row);
MPI_Gather (row, npe(), MPI_INT, &send[0][0], npe(), MPI_INT, 0,
MPI_COMM_WORLD);
rcv_pid_row (sndrcv->rcv, l, row);
MPI_Gather (row, npe(), MPI_INT, &receive[0][0], npe(), MPI_INT, 0,
MPI_COMM_WORLD);
int * astatus = qmalloc (npe(), int);
for (int i = 0; i < npe(); i++)
astatus[i] = 0;
for (int i = 0; i < npe(); i++)
for (int j = 0; j < npe(); j++)
if (send[i][j] != receive[j][i]) {
fprintf (stderr, "%s: %d sends %d to %d at level %d\n",
name, i, send[i][j], j, l);
fprintf (stderr, "%s: %d receives %d from %d at level %d\n",
name, j, receive[j][i], i, l);
fflush (stderr);
for (int k = i - 2; k <= i + 2; k++)
if (k >= 0 && k < npe())
astatus[k] = 1;
for (int k = j - 2; k <= j + 2; k++)
if (k >= 0 && k < npe())
astatus[k] = 1;
}
MPI_Scatter (astatus, 1, MPI_INT, &status, 1, MPI_INT, 0, MPI_COMM_WORLD);
free (astatus);
matrix_free (send);
matrix_free (receive);
}
else { // slave
rcv_pid_row (sndrcv->snd, l, row);
MPI_Gather (row, npe(), MPI_INT, NULL, npe(), MPI_INT, 0, MPI_COMM_WORLD);
rcv_pid_row (sndrcv->rcv, l, row);
MPI_Gather (row, npe(), MPI_INT, NULL, npe(), MPI_INT, 0, MPI_COMM_WORLD);
MPI_Scatter (NULL, 1, MPI_INT, &status, 1, MPI_INT, 0, MPI_COMM_WORLD);
}
if (status) {
fprintf (stderr,
"check_snd_rcv_matrix \"%s\" failed\n"
"Calling debug_mpi(NULL)...\n"
"Aborting...\n",
name);
fflush (stderr);
debug_mpi (NULL);
MPI_Abort (MPI_COMM_WORLD, -3);
}
}
free (row);
}
static bool has_local_child (Point point)
{
foreach_child()
if (is_local(cell))
return true;
return false;
}
trace
void mpi_boundary_update_buffers()
{
if (npe() == 1)
return;
prof_start ("mpi_boundary_update_buffers");
MpiBoundary * m = (MpiBoundary *) mpi_boundary;
SndRcv * mpi_level = &m->mpi_level;
SndRcv * mpi_level_root = &m->mpi_level_root;
SndRcv * restriction = &m->restriction;
snd_rcv_free (mpi_level);
snd_rcv_free (mpi_level_root);
snd_rcv_free (restriction);
static const unsigned short used = 1 << user;
foreach_cell() {
if (is_active(cell) && !is_border(cell))
/* We skip the interior of the local domain.
Note that this can be commented out in case of suspicion that
something is wrong with border cell tagging. */
continue;
if (cell.neighbors) {
// sending
Array pids = {NULL, 0, 0};
int n = locals_pids (point, &pids);
if (n) {
foreach_child()
if (is_local(cell))
for (int i = 0, * p = (int *) pids.p; i < n; i++, p++)
rcv_pid_append (mpi_level->snd, *p, point);
free (pids.p);
}
// receiving
bool locals = false;
if (is_leaf(cell)) { // prolongation
if (is_remote(cell)) {
Point p = point;
foreach_neighbor(1)
if ((is_local(cell) &&
(is_refined(cell) || is_local_prolongation (point, p))) ||
is_root(point)) {
locals = true; break;
}
}
}
else
foreach_neighbor(1)
if (is_local(cell) || is_root(point)) {
locals = true; break;
}
if (locals)
foreach_child()
if (is_remote(cell))
rcv_pid_append (mpi_level->rcv, cell.pid, point),
cell.flags |= used;
// root cells
if (!is_leaf(cell)) {
// sending
if (is_local(cell)) {
Array pids = {NULL, 0, 0};
// root cell
int n = root_pids (point, &pids);
if (n) {
foreach_neighbor()
for (int i = 0, * p = (int *) pids.p; i < n; i++, p++)
if (cell.pid >= 0 && cell.pid != *p)
rcv_pid_append (mpi_level_root->snd, *p, point);
// restriction (remote root)
for (int i = 0, * p = (int *) pids.p; i < n; i++, p++)
rcv_pid_append (restriction->snd, *p, point);
free (pids.p);
}
}
// receiving
else if (is_remote(cell)) {
bool root = false;
foreach_child()
if (is_local(cell)) {
root = true; break;
}
if (root) {
int pid = cell.pid;
foreach_neighbor()
if (is_remote(cell))
rcv_pid_append (mpi_level_root->rcv, pid, point),
cell.flags |= used;
// restriction (remote root)
rcv_pid_append (restriction->rcv, pid, point);
}
}
}
}
// restriction (remote siblings/children)
if (level > 0) {
if (is_local(cell)) {
// sending
Array pids = {NULL, 0, 0};
if (is_remote(aparent(0)))
append_pid (&pids, aparent(0).pid);
int n = root_pids (parent, &pids);
if (n) {
for (int i = 0, * p = (int *) pids.p; i < n; i++, p++)
rcv_pid_append (restriction->snd, *p, point);
free (pids.p);
}
}
else if (is_remote(cell)) {
// receiving
if (is_local(aparent(0)) || has_local_child (parent))
rcv_pid_append (restriction->rcv, cell.pid, point);
}
}
}
/* we remove unused cells
we do a breadth-first traversal from fine to coarse, so that
coarsening of unused cells can proceed fully. */
static const unsigned short keep = 1 << (user + 1);
for (int l = depth(); l >= 0; l--)
foreach_cell()
if (level == l) {
if (level > 0 && (cell.pid < 0 || is_local(cell) || (cell.flags & used)))
aparent(0).flags |= keep;
if (is_refined(cell) && !(cell.flags & keep))
coarsen_cell (point, NULL);
cell.flags &= ~(used|keep);
continue; // level == l
}
/* we update the list of send/receive pids */
m->send->len = m->receive->len = 0;
rcv_pid_append_pids (mpi_level->snd, m->send);
rcv_pid_append_pids (mpi_level_root->snd, m->send);
rcv_pid_append_pids (mpi_level->rcv, m->receive);
rcv_pid_append_pids (mpi_level_root->rcv, m->receive);
prof_stop();
#if DEBUG_MPI
debug_mpi (NULL);
#endif
#ifdef DEBUGCOND
extern double t; NOT_UNUSED(t);
check_snd_rcv_matrix (mpi_level, "mpi_level");
check_snd_rcv_matrix (mpi_level_root, "mpi_level_root");
check_snd_rcv_matrix (restriction, "restriction");
if (DEBUGCOND)
debug_mpi (NULL);
tree_check();
#endif
}
trace
void mpi_boundary_refine (scalar * list)
{
prof_start ("mpi_boundary_refine");
MpiBoundary * mpi = (MpiBoundary *) mpi_boundary;
/* Send refinement cache to each neighboring process. */
Array * snd = mpi->send;
MPI_Request r[2*snd->len/sizeof(int)];
int nr = 0;
for (int i = 0, * dest = snd->p; i < snd->len/sizeof(int); i++,dest++) {
int len = tree->refined.n;
MPI_Isend (&tree->refined.n, 1, MPI_INT, *dest,
REFINE_TAG(), MPI_COMM_WORLD, &r[nr++]);
if (len > 0)
MPI_Isend (tree->refined.p, sizeof(Index)/sizeof(int)*len,
MPI_INT, *dest, REFINE_TAG(), MPI_COMM_WORLD, &r[nr++]);
}
/* Receive refinement cache from each neighboring process.
fixme: use non-blocking receives */
Array * rcv = mpi->receive;
Cache rerefined = {NULL, 0, 0};
for (int i = 0, * source = rcv->p; i < rcv->len/sizeof(int); i++,source++) {
int len;
mpi_recv_check (&len, 1, MPI_INT, *source, REFINE_TAG(),
MPI_COMM_WORLD, MPI_STATUS_IGNORE,
"mpi_boundary_refine (len)");
if (len > 0) {
Index p[len];
mpi_recv_check (p, sizeof(Index)/sizeof(int)*len,
MPI_INT, *source, REFINE_TAG(),
MPI_COMM_WORLD, MPI_STATUS_IGNORE,
"mpi_boundary_refine (p)");
Cache refined = {p, len, len};
foreach_cache (refined)
if (level <= depth() && allocated(0)) {
if (is_leaf(cell)) {
bool neighbors = false;
foreach_neighbor()
if (allocated(0) && (is_active(cell) || is_local(aparent(0)))) {
neighbors = true; break;
}
// refine the cell only if it has local neighbors
if (neighbors)
refine_cell (point, list, 0, &rerefined);
}
}
}
}
/* check that caches were received OK and free ressources */
if (nr)
MPI_Waitall (nr, r, MPI_STATUSES_IGNORE);
/* update the refinement cache with "re-refined" cells */
free (tree->refined.p);
tree->refined = rerefined;
prof_stop();
/* if any cell has been re-refined, we repeat the process to take care
of recursive refinements induced by the 2:1 constraint */
mpi_all_reduce (rerefined.n, MPI_INT, MPI_SUM);
if (rerefined.n)
mpi_boundary_refine (list);
for (scalar s in list)
s.dirty = true;
}
static void check_depth()
{
#if DEBUG_MPI
int max = 0;
foreach_cell_all()
if (level > max)
max = level;
if (depth() != max) {
FILE * fp = fopen ("layer", "w");
fprintf (fp, "depth() = %d, max = %d\n", depth(), max);
for (int l = 0; l <= depth(); l++) {
Layer * L = tree->L[l];
fprintf (fp, "Layer level = %d, nc = %d, len = %d\n", l, L->nc, L->len);
for (int i = 0; i < L->len; i++)
if (L->m[i]) {
fprintf (fp, " i = %d, refcount = %d\n", i,
*((int *)(((char *)L->m[i]) + L->len*sizeof(char *))));
for (int j = 0; j < L->len; j++)
if (L->m[i][j]) {
fprintf (fp, " j = %d\n", j);
fprintf (fp, "point %g %g %d\n",
X0 + (i - 1.5)*L0/(1 << l), Y0 + (j - 1.5)*L0/(1 << l),
l);
}
}
}
fclose (fp);
fp = fopen ("colls", "w");
output_cells_internal (fp);
fclose (fp);
assert (false);
}
#endif
}
typedef struct {
int refined, leaf;
} Remote;
#define REMOTE() ((Remote *)&val(remote,0))
trace
void mpi_boundary_coarsen (int l, int too_fine)
{
if (npe() == 1)
return;
check_depth();
assert (sizeof(Remote) == sizeof(double));
scalar remote[];
foreach_cell() {
if (level == l) {
if (is_local(cell)) {
REMOTE()->refined = is_refined(cell);
REMOTE()->leaf = is_leaf(cell);
}
else {
REMOTE()->refined = true;
REMOTE()->leaf = false;
}
continue;
}
if (is_leaf(cell))
continue;
}
mpi_boundary_level (mpi_boundary, {remote}, l);
foreach_cell() {
if (level == l) {
if (!is_local(cell)) {
if (is_refined(cell) && !REMOTE()->refined)
coarsen_cell_recursive (point, NULL);
else if (is_leaf(cell) && cell.neighbors && REMOTE()->leaf) {
int pid = cell.pid;
foreach_child()
cell.pid = pid;
}
}
continue;
}
if (is_leaf(cell))
continue;
}
check_depth();
if (l > 0) {
foreach_cell() {
if (level == l) {
remote[] = is_local(cell) ? cell.neighbors : 0;
continue;
}
if (is_leaf(cell))
continue;
}
mpi_boundary_level (mpi_boundary, {remote}, l);
foreach_cell() {
if (level == l)
if (!is_local(cell) && is_local(aparent(0)) && remote[]) {
aparent(0).flags &= ~too_fine;
continue;
}
if (is_leaf(cell))
continue;
}
}
}
static void flag_border_cells()
{
foreach_cell() {
if (is_active(cell)) {
short flags = cell.flags & ~border;
foreach_neighbor() {
if (!is_local(cell) || (level > 0 && !is_local(aparent(0)))) {
flags |= border; break;
}
// root cell
if (is_refined_check())
foreach_child()
if (!is_local(cell)) {
flags |= border; break;
}
if (flags & border)
break;
}
cell.flags = flags;
}
else {
cell.flags &= ~border;
// continue; // fixme
}
if (is_leaf(cell)) {
if (cell.neighbors) {
foreach_child()
cell.flags &= ~border;
if (is_border(cell)) {
bool remote = false;
foreach_neighbor (GHOSTS/2)
if (!is_local(cell)) {
remote = true; break;
}
if (remote)
foreach_child()
cell.flags |= border;
}
}
continue;
}
}
}
static int balanced_pid (long index, long nt, int nproc)
{
long ne = max(1, nt/nproc), nr = nt % nproc;
int pid = index < nr*(ne + 1) ?
index/(ne + 1) :
nr + (index - nr*(ne + 1))/ne;
return min(nproc - 1, pid);
}
// static partitioning: only used for tests
trace
void mpi_partitioning()
{
prof_start ("mpi_partitioning");
long nt = 0;
foreach (serial)
nt++;
/* set the pid of each cell */
long i = 0;
tree->dirty = true;
foreach_cell_post (is_active (cell))
if (is_active (cell)) {
if (is_leaf (cell)) {
cell.pid = balanced_pid (i++, nt, npe());
if (cell.neighbors > 0) {
int pid = cell.pid;
foreach_child()
cell.pid = pid;
}
if (!is_local(cell))
cell.flags &= ~active;
}
else {
cell.pid = child(0).pid;
bool inactive = true;
foreach_child()
if (is_active(cell)) {
inactive = false; break;
}
if (inactive)
cell.flags &= ~active;
}
}
flag_border_cells();
prof_stop();
mpi_boundary_update_buffers();
}
void restore_mpi (FILE * fp, scalar * list1)
{
long index = 0, nt = 0, start = ftell (fp);
scalar size[], * list = list_concat ({size}, list1);;
long offset = sizeof(double)*list_len(list);
// read local cells
static const unsigned short set = 1 << user;
scalar * listm = is_constant(cm) ? NULL : (scalar *){fm};
foreach_cell()
if (balanced_pid (index, nt, npe()) <= pid()) {
unsigned flags;
if (fread (&flags, sizeof(unsigned), 1, fp) != 1) {
fprintf (stderr, "restore(): error: expecting 'flags'\n");
exit (1);
}
for (scalar s in list) {
double val;
if (fread (&val, sizeof(double), 1, fp) != 1) {
fprintf (stderr, "restore(): error: expecting scalar\n");
exit (1);
}
if (s.i != INT_MAX)
s[] = val;
}
if (level == 0)
nt = size[];
cell.pid = balanced_pid (index, nt, npe());
cell.flags |= set;
if (!(flags & leaf) && is_leaf(cell)) {
if (balanced_pid (index + size[] - 1, nt, npe()) < pid()) {
fseek (fp, (sizeof(unsigned) + offset)*(size[] - 1), SEEK_CUR);
index += size[];
continue;
}
refine_cell (point, listm, 0, NULL);
}
index++;
if (is_leaf(cell))
continue;
}
// read non-local neighbors
fseek (fp, start, SEEK_SET);
index = 0;
foreach_cell() {
unsigned flags;
if (fread (&flags, sizeof(unsigned), 1, fp) != 1) {
fprintf (stderr, "restore(): error: expecting 'flags'\n");
exit (1);
}
if (cell.flags & set)
fseek (fp, offset, SEEK_CUR);
else {
for (scalar s in list) {
double val;
if (fread (&val, sizeof(double), 1, fp) != 1) {
fprintf (stderr, "restore(): error: expecting a scalar\n");
exit (1);
}
if (s.i != INT_MAX)
s[] = val;
}
cell.pid = balanced_pid (index, nt, npe());
if (is_leaf(cell) && cell.neighbors) {
int pid = cell.pid;
foreach_child()
cell.pid = pid;
}
}
if (!(flags & leaf) && is_leaf(cell)) {
bool locals = false;
foreach_neighbor(1)
if ((cell.flags & set) && (is_local(cell) || is_root(point))) {
locals = true; break;
}
if (locals)
refine_cell (point, listm, 0, NULL);
else {
fseek (fp, (sizeof(unsigned) + offset)*(size[] - 1), SEEK_CUR);
index += size[];
continue;
}
}
index++;
if (is_leaf(cell))
continue;
}
/* set active flags */
foreach_cell_post (is_active (cell)) {
cell.flags &= ~set;
if (is_active (cell)) {
if (is_leaf (cell)) {
if (cell.neighbors > 0) {
int pid = cell.pid;
foreach_child()
cell.pid = pid;
}
if (!is_local(cell))
cell.flags &= ~active;
}
else if (!is_local(cell)) {
bool inactive = true;
foreach_child()
if (is_active(cell)) {
inactive = false; break;
}
if (inactive)
cell.flags &= ~active;
}
}
}
flag_border_cells();
mpi_boundary_update (list);
free (list);
}z_indexing(): fills index with the Z-ordering index.
If leaves is true only leaves are indexed, otherwise all active cells are indexed.
On the master process (pid() == 0), the function returns the (global) maximum index (and -1 on all other processes).
On a single processor, we would just need something like (for leaves)
double i = 0;
foreach()
index[] = i++;In parallel, this is a bit more difficult.
trace
double z_indexing (scalar index, bool leaves)
{We first compute the size of each subtree.
scalar size[];
subtree_size (size, leaves);The maximum index value is the size of the entire tree (i.e. the value of size in the root cell on the master process) minus one.
double maxi = -1.;
if (pid() == 0)
foreach_level(0)
maxi = size[] - 1.;fixme: doc
foreach_level(0)
index[] = 0;
for (int l = 0; l < depth(); l++) {
boundary_iterate (restriction, {index}, l);
foreach_cell() {
if (level == l) {
if (is_leaf(cell)) {
if (is_local(cell) && cell.neighbors) {
int i = index[];
foreach_child()
index[] = i;
}
}
else { // not leaf
bool loc = is_local(cell);
if (!loc)
foreach_child()
if (is_local(cell)) {
loc = true; break;
}
if (loc) {
int i = index[] + !leaves;
foreach_child() {
index[] = i;
i += size[];
}
}
}
continue; // level == l
}
if (is_leaf(cell))
continue;
}
}
boundary_iterate (restriction, {index}, depth());
return maxi;
}