async_io.cc

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) NumeriCor GmbH
// SPDX-License-Identifier: LicenseRef-APL-1.1
 
#include "basics.h"
#include "sim_utils.h"
#include "fem.h"
#include "physics.h"
#include "async_io.h"
 
namespace opencarp {
namespace async {
 
void IO_poll_for_output(async_IO_queue & io_queue)
{
  int rank = get_rank();
 
#if 0
  int size = get_size();
 
  for(int pid = 0; pid < size; pid++) {
    if(rank == pid)
      printf("IO rank %d / %d: polling ..\n", rank+1, size);
 
    MPI_Barrier(PETSC_COMM_WORLD);
  }
#endif
 
  bool do_continue = true;
  int cmd_buff = 0;
  MPI_Status status;
 
  while(do_continue) {
    if(!rank)
      MPI_Recv(&cmd_buff, 1, MPI_INT, 0, ASYNC_TAG, user_globals::IO_Intercomm, &status);
 
    MPI_Bcast(&cmd_buff, 1, MPI_INT, 0, PETSC_COMM_WORLD);
 
    switch(cmd_buff) {
      case ASYNC_CMD_EXIT:
        do_continue = false;
        break;
 
      case ASYNC_CMD_REGISTER_OUTPUT:
        IO_register_output(io_queue);
        break;
 
      case ASYNC_CMD_OUTPUT:
        IO_do_output(io_queue);
        break;
 
      default: break;
    }
  }
 
  for(IGBheader* igb : io_queue.IGBs) {
    if(igb->fileptr()) fclose((FILE*)igb->fileptr());
    delete igb;
  }
 
#if 0
  for(int pid = 0; pid < size; pid++) {
    if(rank == pid)
      printf("IO rank %d / %d: exiting ..\n", rank+1, size);
 
    MPI_Barrier(PETSC_COMM_WORLD);
  }
#endif
}
 
void COMPUTE_send_exit_flag()
{
  int rank = get_rank();
 
  if(!rank) {
    int msg = ASYNC_CMD_EXIT;
    MPI_Send(&msg, 1, MPI_INT, 0, ASYNC_TAG, user_globals::IO_Intercomm);
  }
}
 
int  COMPUTE_register_output(const SF::vector<mesh_int_t> & idx,
                             const int dpn,
                             const char* name,
                             const char* units)
{
  intercomm_layout il;
  il.setup(user_globals::IO_Intercomm);
 
  long int loc_size = idx.size();
  SF::vector<long int> layout;
  layout_from_count(loc_size, layout, PETSC_COMM_WORLD);
 
  if(il.loc_rank == 0) {
    char header[2048];
 
    const timer_manager & tm = *user_globals::tm_manager;
    const int num_io = tm.timers[iotm_spacedt]->numIOs;
    const double dimt = tm.end - tm.start;
 
    snprintf(header, sizeof header, "%d %d %lf %s %s", dpn, num_io, dimt, name, units);
 
    int msg = ASYNC_CMD_REGISTER_OUTPUT;
    MPI_Send(&msg, 1, MPI_INT, 0, ASYNC_TAG, user_globals::IO_Intercomm);
 
    MPI_Send(header, 2048, MPI_CHAR, 0, ASYNC_TAG, user_globals::IO_Intercomm);
    MPI_Send(layout.data(), layout.size(), MPI_LONG, 0, ASYNC_TAG, user_globals::IO_Intercomm);
  }
 
  MPI_Barrier(PETSC_COMM_WORLD);
 
  int receive_rank = COMPUTE_get_receive_rank(il);
  MPI_Send(idx.data(), idx.size()*sizeof(mesh_int_t), MPI_BYTE, receive_rank,
           ASYNC_TAG, user_globals::IO_Intercomm);
 
  int id = -1;
 
  if(il.loc_rank == 0) {
    MPI_Status status;
    MPI_Recv(&id, 1, MPI_INT, 0, ASYNC_TAG, user_globals::IO_Intercomm, &status);
  }
 
  MPI_Bcast(&id, 1, MPI_INT, 0, PETSC_COMM_WORLD);
  return id;
}
 
void IO_prepare_sort(const SF::vector<mesh_int_t> & inp_idx,
                     SF::commgraph<size_t> & grph,
                     SF::vector<mesh_int_t> & perm_before_comm,
                     SF::vector<mesh_int_t> & perm_after_comm)
{
  int size = get_size(), rank = get_rank();
  SF::vector<mesh_int_t> idx(inp_idx);
 
  // determine global min and max indices
  mesh_int_t gmax = SF::global_max(idx, PETSC_COMM_WORLD);
  mesh_int_t gmin = SF::global_min(idx, PETSC_COMM_WORLD);
 
  // block size
  mesh_int_t bsize = (gmax - gmin) / size + 1;
 
  // distribute tuples uniquely and linearly ascending across the ranks ----------------
  SF::vector<mesh_int_t> dest(idx.size());
  SF::interval(perm_before_comm, 0, idx.size());
 
  // find a destination for every tuple
  for(size_t i=0; i<dest.size(); i++)
    dest[i] = (idx[i] - gmin) / bsize;
 
  // find permutation to sort tuples in the send buffer
  binary_sort_copy(dest, perm_before_comm);
 
  // fill send buffer
  SF::vector<mesh_int_t> snd_idx(idx.size());
  for(size_t i=0; i<perm_before_comm.size(); i++)
    snd_idx[i] = idx[perm_before_comm[i]];
 
  // communicate
  grph.configure(dest, PETSC_COMM_WORLD);
 
  size_t rsize = sum(grph.rcnt);
  SF::vector<mesh_int_t> recv_idx(rsize);
 
  MPI_Exchange(grph, snd_idx, recv_idx, PETSC_COMM_WORLD);
 
  // sort the received values locally
  SF::interval(perm_after_comm, 0, recv_idx.size());
  binary_sort_copy(recv_idx, perm_after_comm);
}
 
 
IGBheader* IO_open_igb(const int numIOs, const double dimt,
                       const size_t gsize, const int dpn,
                       const char* name, const char* units)
{
  IGBheader* pigb = new IGBheader();
  IGBheader & igb = *pigb;
  int       err    = 0;
 
  igb.x(gsize / dpn);
  igb.dim_x(igb.x()-1);
  igb.inc_x(1);
 
  igb.y(1); igb.z(1);
  igb.t(numIOs);
  igb.dim_t(dimt);
 
  switch(dpn) {
    default:
    case 1: igb.type(IGB_FLOAT); break;
    case 3: igb.type(IGB_VEC3_f); break;
    case 4: igb.type(IGB_VEC4_f); break;
    case 9: igb.type(IGB_VEC9_f); break;
  }
 
  igb.unites_x("um"); igb.unites_y("um"); igb.unites_z("um");
  igb.unites_t("ms");
  igb.unites(units);
 
  igb.inc_t(param_globals::spacedt);
 
  if(get_rank() == 0) {
    FILE_SPEC file = f_open(name, "wb");
    if(file != NULL) {
      igb.fileptr(file->fd);
      igb.write();
      delete file;
    }
    else err++;
  }
 
  err = get_global(err, MPI_SUM);
  if(err) {
    log_msg(0,5,0, "%s error: Could not set up data output! Aborting!", __func__);
    EXIT(1);
  }
 
  return pigb;
}
 
void IO_register_output(async_IO_queue & io_queue)
{
  set_dir(OUTPUT);
 
  intercomm_layout il;
  il.setup(user_globals::IO_Intercomm);
 
  char header[2048];
  SF::vector<long int> data_layout(il.rem_size+1);
  MPI_Status status;
 
  if(il.loc_rank == 0) {
    MPI_Recv(header, 2048, MPI_CHAR, 0, ASYNC_TAG, user_globals::IO_Intercomm, &status);
    MPI_Recv(data_layout.data(), data_layout.size(), MPI_LONG, 0, ASYNC_TAG, user_globals::IO_Intercomm, &status);
  }
 
  MPI_Bcast(header, 2048, MPI_CHAR, 0, PETSC_COMM_WORLD);
  MPI_Bcast(data_layout.data(), data_layout.size(), MPI_LONG, 0, PETSC_COMM_WORLD);
 
  int dpn = 0, numIOs = 0;
  double dimt = 0.0;
 
  char name_str[2048], units_str[2048];
  sscanf(header, "%d %d %lf %s %s", &dpn, &numIOs, &dimt, name_str, units_str);
 
  size_t gsize = data_layout[data_layout.size()-1];
  IGBheader* igb = IO_open_igb(numIOs, dimt, gsize, dpn, name_str, units_str);
 
  SF::vector<int> senders;
  IO_get_sender_ranks(il, senders);
 
  size_t num_recv = 0;
  for(int s : senders)
    num_recv += data_layout[s+1] - data_layout[s];
 
  SF::vector<mesh_int_t> idx_buff(num_recv);
 
  SF::vector<MPI_Request> req (senders.size());
  SF::vector<MPI_Status>  stat(senders.size());
 
  for(size_t i=0, dsp=0; i<senders.size(); i++) {
    int    send_rank = senders[i];
    size_t send_size = data_layout[send_rank+1] - data_layout[send_rank];
    MPI_Irecv(idx_buff.data() + dsp, send_size*sizeof(mesh_int_t), MPI_BYTE, send_rank,
              ASYNC_TAG, user_globals::IO_Intercomm, req.data()+i);
    dsp += send_size;
  }
 
  MPI_Waitall(senders.size(), req.data(), stat.data());
 
  SF::commgraph<size_t> cg;
  SF::vector<mesh_int_t> pafter, pbefore;
  IO_prepare_sort(idx_buff, cg, pbefore, pafter);
 
  int id = io_queue.add(igb, data_layout, cg, pbefore, pafter);
 
  if(il.loc_rank == 0)
    MPI_Send(&id, 1, MPI_INT, 0, ASYNC_TAG, user_globals::IO_Intercomm);
}
 
int  COMPUTE_get_receive_rank(const intercomm_layout & il)
{
  SF::vector<int> remote_dist;
  SF::divide(il.loc_size, il.rem_size, remote_dist);
 
  int remote_idx = 0;
  int rank_cnt = remote_dist[remote_idx];
 
  while(remote_idx < int(remote_dist.size()) && il.loc_rank >= rank_cnt) {
    remote_idx++;
    rank_cnt += remote_dist[remote_idx];
  }
 
  return remote_idx;
}
void IO_get_sender_ranks(const intercomm_layout & il, SF::vector<int> & sender)
{
  SF::vector<int> recv_cnt, recv_dsp, recv_ranks;
  SF::divide(il.rem_size, il.loc_size, recv_cnt);
  SF::dsp_from_cnt(recv_cnt, recv_dsp);
 
  sender.resize(0);
  sender.reserve(recv_cnt[il.loc_rank]);
 
  int start = recv_dsp[il.loc_rank], stop = recv_dsp[il.loc_rank+1];
  for(int i=start; i<stop; i++)
    sender.push_back(i);
}
 
void IO_sort_data(SF::vector<float> & data,
                  const SF::vector<mesh_int_t> & perm_b,
                  const SF::vector<mesh_int_t> & perm_a,
                  SF::commgraph<size_t> & cg)
{
  SF::vector<float> snd_data(data.size());
  for(size_t i=0; i<perm_b.size(); i++)
    snd_data[i] = data[perm_b[i]];
 
  SF::vector<float> recv_data(perm_a.size());
  SF::MPI_Exchange(cg, snd_data, recv_data, PETSC_COMM_WORLD);
 
  data.resize(perm_a.size());
  for(size_t i=0; i<perm_a.size(); i++)
    data[i] = recv_data[perm_a[i]];
}
 
void IO_do_output(async_IO_queue & io_queue)
{
  intercomm_layout il;
  il.setup(user_globals::IO_Intercomm);
 
  MPI_Status status;
  int id = -1;
 
  if(il.loc_rank == 0)
    MPI_Recv(&id, 1, MPI_INT, 0, ASYNC_TAG, user_globals::IO_Intercomm, &status);
  MPI_Bcast(&id, 1, MPI_INT, 0, PETSC_COMM_WORLD);
 
  assert(id > -1 && id < int(io_queue.IGBs.size()));
 
  IGBheader*               igb    = io_queue.IGBs[id];
  SF::vector<long int> &   layout = io_queue.layouts[id];
  SF::commgraph<size_t> &  cg     = io_queue.cg[id];
  SF::vector<mesh_int_t> & perm_b = io_queue.perm_b[id];
  SF::vector<mesh_int_t> & perm_a = io_queue.perm_a[id];
 
  SF::vector<int> senders;
  IO_get_sender_ranks(il, senders);
 
  size_t num_recv = 0;
  for(int s : senders)
    num_recv += layout[s+1] - layout[s];
 
  SF::vector<float> buff(num_recv);
  SF::vector<MPI_Request> req (senders.size());
  SF::vector<MPI_Status>  stat(senders.size());
 
  for(size_t i=0, dsp=0; i<senders.size(); i++) {
    int    send_rank = senders[i];
    size_t send_size = layout[send_rank+1] - layout[send_rank];
 
    MPI_Irecv(buff.data() + dsp, send_size, MPI_FLOAT, send_rank,
              ASYNC_TAG, user_globals::IO_Intercomm, req.data()+i);
    dsp += send_size;
  }
 
  MPI_Waitall(senders.size(), req.data(), stat.data());
  IO_sort_data(buff, perm_b, perm_a, cg);
  SF::root_write((FILE*)igb->fileptr(), buff, PETSC_COMM_WORLD);
}
 
void COMPUTE_do_output(SF_real* dat, const int lsize, const int IO_id)
{
  intercomm_layout il;
  il.setup(user_globals::IO_Intercomm);
 
  if(il.loc_rank == 0) {
    int msg = ASYNC_CMD_OUTPUT;
    MPI_Send(&msg, 1, MPI_INT, 0, ASYNC_TAG, user_globals::IO_Intercomm);
 
    msg = IO_id;
    MPI_Send(&msg, 1, MPI_INT, 0, ASYNC_TAG, user_globals::IO_Intercomm);
  }
 
  int receive_rank = COMPUTE_get_receive_rank(il);
  SF::vector<float> data;
  data.assign(dat, dat+lsize);
 
  MPI_Send(data.data(), data.size(), MPI_FLOAT, receive_rank, ASYNC_TAG, user_globals::IO_Intercomm);
}
 
void COMPUTE_do_output(SF_real* dat, const SF::vector<mesh_int_t> & idx, const int IO_id)
{
  intercomm_layout il;
  il.setup(user_globals::IO_Intercomm);
 
  if(il.loc_rank == 0) {
    int msg = ASYNC_CMD_OUTPUT;
    MPI_Send(&msg, 1, MPI_INT, 0, ASYNC_TAG, user_globals::IO_Intercomm);
 
    msg = IO_id;
    MPI_Send(&msg, 1, MPI_INT, 0, ASYNC_TAG, user_globals::IO_Intercomm);
  }
 
  int receive_rank = COMPUTE_get_receive_rank(il);
 
  SF::vector<float> data;
  data.resize(idx.size()); data.resize(0);
 
  for(mesh_int_t ii : idx)
    data.push_back(dat[ii]);
 
  MPI_Send(data.data(), data.size(), MPI_FLOAT, receive_rank, ASYNC_TAG, user_globals::IO_Intercomm);
}
 
 
}}