async_io.cc

Go to the documentation of this file.

// ----------------------------------------------------------------------------
// openCARP is an open cardiac electrophysiology simulator.
//
// Copyright (C) 2022 openCARP project
//
// This program is licensed under the openCARP Academic Public License (APL)
// v1.0: You can use and redistribute it and/or modify it in non-commercial
// academic environments under the terms of APL as published by the openCARP
// project v1.0, or (at your option) any later version. Commercial use requires
// a commercial license (info@opencarp.org).
//
// This program is distributed without any warranty; see the openCARP APL for
// more details.
//
// You should have received a copy of the openCARP APL along with this program
// and can find it online: http://www.opencarp.org/license
// ----------------------------------------------------------------------------

#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);
}


}}