Logo Search packages:      
Sourcecode: ecs version File versions  Download package

cs_partition.c

/*============================================================================
 *  Domain partitioning for Code_Saturne 1.4
 *============================================================================*/

/*
  This file is part of the Code_Saturne Preprocessor, element of the
  Code_Saturne CFD tool.

  Copyright (C) 2007-2008 EDF S.A., France

  contact: saturne-support@edf.fr

  The Code_Saturne Preprocessor is free software; you can redistribute it
  and/or modify it under the terms of the GNU General Public License
  as published by the Free Software Foundation; either version 2 of
  the License, or (at your option) any later version.

  The Code_Saturne Preprocessor is distributed in the hope that it will be
  useful, but WITHOUT ANY WARRANTY; without even the implied warranty
  of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with the Code_Saturne Preprocessor; if not, write to the
  Free Software Foundation, Inc.,
  51 Franklin St, Fifth Floor,
  Boston, MA  02110-1301  USA
*/

/* Detect version of C used (C89 or C99) */

#if !defined(__STDC_VERSION__)
#  define __STDC_VERSION__ 1989
#endif

/* Include configuration file */

#include "ecs_config.h"

/*----------------------------------------------------------------------------
 * Standard C library headers
 *----------------------------------------------------------------------------*/

#include <assert.h>
#include <locale.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#if (__STDC_VERSION__ >= 199901L)
#include <stdint.h>
#endif

/*----------------------------------------------------------------------------
 * Internationalization macros
 *----------------------------------------------------------------------------*/

#if defined(ENABLE_NLS)

#include <libintl.h>
#define _(String) gettext(String)
#define gettext_noop(String) String
#define N_(String) gettext_noop(String)

#else

#define _(String) String
#define N_(String) String
#define textdomain(Domain)
#define bindtextdomain(Package, Directory)

#endif

/*----------------------------------------------------------------------------
 * METIS library headers
 *----------------------------------------------------------------------------*/

#if defined(HAVE_METIS)

#ifdef __cplusplus
extern "C" {
#endif

/*
  Caution:
  --------
  File included by <metis.h> for Metis 4.0 lead to compilation
  errors on Linux in C99 mode due to redeclaration of functions from
  <stdlib.h>, so we include our own prototypes here.
*/

typedef int idxtype;
void METIS_PartGraphRecursive(int *, idxtype *, idxtype *, idxtype *, idxtype *,
                              int *, int *, int *, int *, int *, idxtype *);
void METIS_PartGraphKway(int *, idxtype *, idxtype *, idxtype *, idxtype *,
                         int *, int *, int *, int *, int *, idxtype *);

#ifdef __cplusplus
}
#endif

#endif /* defined(HAVE_METIS) */

/*----------------------------------------------------------------------------
 * BFT library headers
 *----------------------------------------------------------------------------*/

#include <bft_error.h>
#include <bft_file.h>
#include <bft_mem.h>
#include <bft_mem_usage.h>
#include <bft_printf.h>
#include <bft_timer.h>

/*----------------------------------------------------------------------------
 *  Local headers
 *----------------------------------------------------------------------------*/

/*---------------------------------------------------------------------------*/

#ifdef __cplusplus
extern "C" {
#if 0
} /* Fake brace to force Emacs auto-indentation back to column 0 */
#endif
#endif /* __cplusplus */

/*============================================================================
 * Local Type Definitions
 *============================================================================*/

typedef int cs_int_t;

typedef struct {

  bft_file_t     *f;              /* Associated file structure pointer */

  size_t          header_size;    /* Header default size */
  size_t          header_align;   /* Header alignment */
  size_t          body_align;     /* Body alignment */

  size_t          buffer_size;    /* Current size of header buffer */
  unsigned char  *buffer;         /* Header buffer */

  size_t          n_values;       /* Number of values in section header */
  size_t          type_size;      /* Size of current type */
  char           *name;           /* Pointer to name field in section header */
  char           *type_name;      /* Pointer to type field in section header */
  void           *data;           /* Pointer to data in section header */

} _cs_io_t;

/*============================================================================
 * Static global variables
 *============================================================================*/

/*============================================================================
 * Private function definitions
 *============================================================================*/

/*----------------------------------------------------------------------------
 * Display the distribution of cells per partition
 *
 * parameters:
 *   n_cells   --> number of cells
 *   n_parts   --> number of partitions
 *   part      --> cell partition number
 *----------------------------------------------------------------------------*/

static void
_cell_part_histogram(size_t      n_cells,
                     int         n_parts,
                     const   int part[])
{
  int i, k;
  size_t j, n_min, n_max;
  double step;

  size_t *n_part_cells;
  size_t count[10];
  int n_steps = 10;

  if (n_cells <=1 || n_parts <= 1) /* Should never happen */
    return;

  bft_printf(_("  Number of cells per domain (histogramm):\n"));

  BFT_MALLOC(n_part_cells, n_parts, size_t);

  for (i = 0; i < n_parts; i++)
    n_part_cells[i] = 0;

  for (j = 0; j < n_cells; j++)
    n_part_cells[part[j] - 1] += 1;

  /* Compute min and max */

  n_min = n_part_cells[0];
  n_max = n_part_cells[0];

  for (i = 1; i < n_parts; i++) {
    if (n_part_cells[i] > n_max)
      n_max = n_part_cells[i];
    else if (n_part_cells[i] < n_min)
      n_min = n_part_cells[i];
  }

  /* Define axis subdivisions */

  for (i = 0; i < n_steps; i++)
    count[i] = 0;

  if (n_max - n_min > 0) {

    if (n_max-n_min < (size_t)n_steps)
      n_steps = n_max-n_min > 0 ? n_max-n_min : 1;

    step = (double)(n_max - n_min) / n_steps;

    /* Loop on partitions */

    for (i = 0; i < n_parts; i++) {

      /* Associated subdivision */

      for (j = 0, k = 1; k < n_steps; j++, k++) {
        if (n_part_cells[i] < n_min + k*step)
          break;
      }
      count[j] += 1;

    }

    for (i = 0, j = 1; i < n_steps - 1; i++, j++)
      bft_printf("    [ %10d ; %10d [ = %10d\n",
                 (int)(n_min + i*step),
                 (int)(n_min + j*step),
                 (int)(count[i]));

    bft_printf("    [ %10d ; %10d ] = %10d\n\n",
               (int)(n_min + (n_steps - 1)*step),
               (int)n_max,
               (int)(count[n_steps - 1]));

  }

  else { /* if (n_max == n_min) */
    bft_printf("    [ %10d ; %10d ] = %10d\n\n",
               (int)(n_min), (int)n_max, (int)n_parts);
  }

  BFT_FREE(n_part_cells);
}

/*----------------------------------------------------------------------------
 * Convert a buffer of type uint64_t to size_t
 *
 * parameters:
 *   buf <-- buffer
 *   val --> array to which values are converted
 *   n   <-- number of values to convert
 *----------------------------------------------------------------------------*/

static void
_convert_to_size(const unsigned char  buf[],
                 size_t               val[],
                 size_t               n)
{
  size_t i;

#if (__STDC_VERSION__ >= 199901L)

  for (i = 0; i < n; i++)
    val[i] = ((const uint64_t *)buf)[i];

#else

  if (sizeof(size_t) == 8) {
    for (i = 0; i < n; i++)
      val[i] = ((const size_t *)buf)[i];
  }
  else if (sizeof(unsigned long long) == 8) {
    for (i = 0; i < n; i++)
      val[i] = ((const unsigned long long *)buf)[i];
  }
  else
    bft_error(__FILE__, __LINE__, 0,
              _("Compilation configuration / porting error:\n"
                "Unable to determine a 64-bit unsigned int type.\n"
                "size_t is %d bits, unsigned long long %d bits"),
              sizeof(size_t)*8, sizeof(unsigned long long)*8);

#endif
}

/*----------------------------------------------------------------------------
 * Open input file, testing magic string for type.
 *
 * parameters:
 *   filename     <-- file name
 *
 * returns:
 *   File metadata structure
 *----------------------------------------------------------------------------*/

static _cs_io_t
_open_input(const char *filename)
{
  size_t alignments[3];
  char header_buf[65];
  char expected_content[] = "Face-based mesh definition, R0";

  _cs_io_t inp;

  inp.f = NULL;
  inp.header_size = 0;
  inp.header_align = 0;
  inp.body_align = 0;
  inp.buffer_size = 0;
  inp.buffer = NULL;

  inp.n_values = 0;
  inp.type_size = 0;
  inp.name = NULL;
  inp.type_name = NULL;
  inp.data = NULL;

  /* Open file */

  bft_printf(_("\nOpening input file: \"%s\"\n\n"), filename) ;

  bft_printf_flush();

  inp.f = bft_file_open(filename,
                        BFT_FILE_MODE_READ,
                        BFT_FILE_TYPE_BINARY);

  bft_file_set_big_endian(inp.f);

  /* Read "magic string" */

  bft_file_read(header_buf, 1, 64, inp.f);

  if (strncmp(header_buf, "Code_Saturne I/O, BE, R0", 64) != 0) {
    header_buf[64] = '\0';
    bft_error(__FILE__, __LINE__, 0,
              _("File format of \"%s\" is not recognized:\n"
                "First %d bytes: \"%s\"."),
              filename, 64, header_buf);
  }

  bft_file_read(header_buf, 1, 64, inp.f);

  if (strncmp(header_buf, expected_content, 64) != 0) {
    header_buf[64] = '\0';
    bft_error(__FILE__, __LINE__, 0,
              _("Content type \"%s\" not expected.\n"
                "Expected: \"%s\"."),
              header_buf, expected_content);
  }

  bft_file_read(header_buf, 8, 3, inp.f);

  _convert_to_size((unsigned char*)header_buf, alignments, 3);

  inp.header_size = alignments[0];
  inp.header_align = alignments[1];
  inp.body_align = alignments[2];

  inp.buffer_size = inp.header_size;
  BFT_MALLOC(inp.buffer, inp.buffer_size, unsigned char);

  /* Finish */

  return inp;
}

/*----------------------------------------------------------------------------
 * Close input file
 *
 * parameters:
 *   f <-> pointer to file object
 *----------------------------------------------------------------------------*/

static void
_close_input(_cs_io_t *inp)
{
  if (inp != NULL) {
    if (inp->f != NULL) {
      bft_printf(_("\nClosing input: \"%s\"\n\n"),
                 bft_file_get_name(inp->f)) ;
      inp->f = bft_file_free(inp->f);
    }
    inp->header_size = 0;
    inp->header_align = 0;
    inp->body_align = 0;
    inp->buffer_size = 0;
    BFT_FREE(inp->buffer);
  }
}

/*----------------------------------------------------------------------------
 * Read a record of type integers and convert it to cs_int_t if necessary.
 *
 * parameters:
 *   inp      <-- pointer to input object
 *
 * returns:
 *   pointer to data read (must be deallocated later).
 *----------------------------------------------------------------------------*/

static cs_int_t *
_read_int_t(_cs_io_t  *inp)
{
  size_t i;
  size_t n_elts = inp->n_values;
  cs_int_t  *retval = NULL;

  /* Position read pointer if necessary */

  {
    bft_file_off_t offset = bft_file_tell(inp->f);
    size_t ba = inp->body_align;
    offset += (ba - (offset % ba)) % ba;
    bft_file_seek(inp->f, offset, BFT_FILE_SEEK_SET);
  }

  /* Read data, convert if necessary */

  BFT_MALLOC(retval, n_elts, cs_int_t);

  if (sizeof(cs_int_t) == inp->type_size)
    bft_file_read(retval, inp->type_size, n_elts, inp->f);

  else {

    size_t int_size = inp->type_size;
    void *_buf = NULL;

    BFT_MALLOC(_buf, n_elts*int_size, unsigned char);

    bft_file_read(_buf, int_size, n_elts, inp->f);

    if (sizeof(int) == int_size) {
      int *buf = _buf;
      for (i = 0; i < n_elts; i++)
        retval[i] = buf[i];
    }
    else if (sizeof(long) == int_size) {
      long *buf = _buf;
      for (i = 0; i < n_elts; i++)
        retval[i] = buf[i];
    }
    else if (sizeof(long long) == int_size) {
      long long *buf = _buf;
      for (i = 0; i < n_elts; i++)
        retval[i] = buf[i];
    }
    else if (sizeof(short) == int_size) {
      short *buf = _buf;
      for (i = 0; i < n_elts; i++)
        retval[i] = buf[i];
    }
    else
      bft_error(__FILE__, __LINE__, 0,
                _("Compilation configuration / porting error:\n"
                  "Unable to find an integer type of size %d bits:\n"
                  "int is %d bits, long %d, long long %d, short %d."),
                int_size*8, sizeof(int)*8, sizeof(long)*8,
                sizeof(long long)*8, sizeof(short)*8);

    BFT_FREE(_buf);

  }

  return retval;
}

/*----------------------------------------------------------------------------
 * Read section header.
 *
 * parameters:
 *   inp <-- pointer to input object
 *
 * returns:
 *   number of bytes in section body
 *----------------------------------------------------------------------------*/

static size_t
_read_section_header(_cs_io_t  *inp)
{
  int type_name_error = 0;
  size_t body_size = 0;
  size_t header_vals[6];
  unsigned int_endian = 0;

  *((char *)(&int_endian)) = '\1'; /* Determine if we are little-endian */

  assert(inp != NULL);
  assert(inp->f != NULL);
  assert(inp->buffer != NULL);

  /* Position read pointer if necessary */
  /*------------------------------------*/

  {
    bft_file_off_t offset = bft_file_tell(inp->f);
    size_t ha = inp->header_align;
    offset += (ha - (offset % ha)) % ha;
    bft_file_seek(inp->f, offset, BFT_FILE_SEEK_SET);
  }

  /* Read header */
  /*-------------*/

  bft_file_read(inp->buffer, 1, inp->header_size, inp->f);

  if (int_endian == 1)
    bft_file_swap_endian(inp->buffer, inp->buffer, 8, 6);

  _convert_to_size(inp->buffer, header_vals, 6);

  if (header_vals[0] > inp->header_size) {

    if (header_vals[0] > inp->buffer_size) {
      while (header_vals[0] > inp->buffer_size)
        inp->buffer_size *=2;
      BFT_REALLOC(inp->buffer, inp->buffer_size, unsigned char);
    }

    bft_file_read(inp->buffer + inp->header_size,
                  1,
                  header_vals[0] - inp->header_size,
                  inp->f);

  }

  /* Set pointers to data fields */

  inp->n_values = header_vals[1];
  inp->type_size = 0;
  inp->data = NULL;
  inp->type_name = (char *)(inp->buffer + 48);
  inp->name = (char *)(inp->buffer + 56);

  if (header_vals[1] > 0 && inp->type_name[7] == 'e')
    inp->data = inp->buffer + 56 + header_vals[5];

  inp->type_size = 0;

  if (inp->n_values > 0) {

    /* Check type name and compute size of data */

    if (inp->type_name[0] == 'c') {
      if (inp->type_name[1] != ' ')
        type_name_error = 1;
      else
        inp->type_size = 1;
    }
    else if (   inp->type_name[0] == 'i'
             || inp->type_name[0] == 'u'
             || inp->type_name[0] == 'r') {

      if (inp->type_name[1] == '4')
        inp->type_size = 4;
      else if (inp->type_name[1] == '8')
        inp->type_size = 8;
      else
        type_name_error = 1;

    }
    else
      type_name_error = 1;

    if (type_name_error)
      bft_error(__FILE__, __LINE__, 0,
                _("Type \"%s\" is not known\n"
                  "Known types: \"c \", \"i4\", \"i8\", \"u4\", \"u8\", "
                  "\"r4\", \"r8\"."), inp->type_name);

    else if (inp->data == NULL)
      body_size = inp->type_size*inp->n_values;

    else if (int_endian == 1 && inp->type_size > 1)
      bft_file_swap_endian(inp->data,
                           inp->data,
                           inp->type_size,
                           inp->n_values);
  }

  return body_size;
}

/*----------------------------------------------------------------------------
 * Update face -> cells connectivity with periodic faces info
 *
 * parameters:
 *   n_couples    <-- number of periodic face couples
 *   face_couples <-- periodic face couples
 *   face_cells   <-> face->cells connectivity
 *----------------------------------------------------------------------------*/

static void
_add_periodic_faces(size_t           n_couples,
                    const cs_int_t   face_couples[],
                    cs_int_t         face_cells[])
{
  size_t i, face_id_1, face_id_2;

  for (i = 0; i < n_couples; i++) {

    face_id_1 = face_couples[i*2    ] - 1;
    face_id_2 = face_couples[i*2 + 1] - 1;

    if (face_cells[face_id_1 * 2] == 0) {

      assert(face_cells[face_id_1*2 + 1] != 0);
      assert(face_cells[face_id_2*2]     != 0);
      assert(face_cells[face_id_2*2 + 1] == 0);

      face_cells[face_id_1*2]     = face_cells[face_id_2*2];
      face_cells[face_id_2*2 + 1] = face_cells[face_id_1*2 + 1];

    }
    else {

      assert(face_cells[face_id_1*2 + 1] == 0);
      assert(face_cells[face_id_2*2]     == 0);
      assert(face_cells[face_id_2*2 + 1] != 0);

      face_cells[face_id_1*2 + 1] = face_cells[face_id_2*2 + 1];
      face_cells[face_id_2*2]     = face_cells[face_id_1*2];

    }

  }

}

/*----------------------------------------------------------------------------
 * Read input file
 *
 * parameters:
 *   no_perio   <-- ignore periodicity information if 1
 *   n_cells    --> number of cells in mesh
 *   n_faces    --> number of cells in mesh
 *   face_cells --> face->cells connectivity
 *----------------------------------------------------------------------------*/

static void
_read_input(int          no_perio,
            size_t      *n_cells,
            size_t      *n_faces,
            cs_int_t   **face_cells)
{
  size_t data_size;

  _cs_io_t  inp;

  int read_type;
  int eof_reached = 0;

  const char *read_type_name[3] = {N_("Read:   "),
                                   N_("Ignored:"),
                                   N_("Reached:")};

  /* Initialization */

  *n_cells = 0;
  *n_faces = 0;
  *face_cells = NULL;

  /* Open preprocessor input */

  inp = _open_input("preprocessor_output");

  assert(inp.f  != NULL);

  while (!eof_reached) {

    /* Read section header */

    read_type = 0;

    data_size = _read_section_header(&inp);

    if (!strcmp(inp.name, "n_cells")) {
      assert(!strcmp(inp.type_name, "u8"));
      assert(inp.n_values == 1 && inp.data != NULL);
      _convert_to_size(inp.data, n_cells, 1);
    }

    else if (!strcmp(inp.name, "n_faces")) {
      assert(!strcmp(inp.type_name, "u8"));
      assert(inp.n_values == 1 && inp.data != NULL);
      _convert_to_size(inp.data, n_faces, 1);
    }

    /* Main face -> cells connectivity */

    else if (!strcmp(inp.name, "face_cells")) {
      assert(   !strcmp(inp.type_name, "i8")
             || !strcmp(inp.type_name, "i4"));
      if (inp.n_values > 0)
        *face_cells = _read_int_t(&inp);
    }

    /* Additional periodicity connectivity */

    else if (!no_perio && !strncmp(inp.name,
                                   "periodicity_faces_",
                                   strlen("periodicity_faces_"))) {

      assert(   !strcmp(inp.type_name, "i8")
             || !strcmp(inp.type_name, "i4"));

      if (inp.n_values > 0) {

        cs_int_t  *perio_faces = NULL;

        if (face_cells == NULL)
          bft_error(__FILE__, __LINE__, 0,
                    _("Non-empty section \"%s\" encountered before\n"
                      "\"face_cells\" or with empty face->cells connectivity."),
                    inp.name);

        perio_faces = _read_int_t(&inp);

        _add_periodic_faces(inp.n_values/2,
                            perio_faces,
                            *face_cells);

        BFT_FREE(perio_faces);

      }
    }

    /* Ignored sections */

    else if (!strcmp(inp.name, "EOF")) {
      read_type = 2;
      eof_reached = 1;
    }

    else {
      read_type = 1;
      if (inp.n_values > 0 && inp.data == NULL) {
        bft_file_off_t offset = bft_file_tell(inp.f);
        size_t ha = inp.header_align;
        offset += (ha - (offset % ha)) % ha;
        offset += inp.n_values*inp.type_size;
        bft_file_seek(inp.f, offset, BFT_FILE_SEEK_SET);
      }
    }

    /* Print info */

    if (inp.n_values > 0) {
      bft_printf(_("  %s \"%-32s\"; Type: \"%s\"; Size: %lu\n"),
                 _(read_type_name[read_type]),
                 inp.name, inp.type_name, (unsigned long)(inp.n_values));
    }
    else {
      bft_printf(_("  %s \"%-32s\"\n"),
                 _(read_type_name[read_type]), inp.name);
    }

  }

  /* Clean-up */

  _close_input(&inp);
}

#if defined(HAVE_METIS)

/*----------------------------------------------------------------------------
 * Build cell -> cell connectivity
 *
 * parameters:
 *   n_cells        <-- number of cells in mesh
 *   n_faces        <-- number of cells in mesh
 *   face_cells     <-- face->cells connectivity
 *   cell_idx       --> cell->cells index
 *   cell_neighbors --> cell->cells connectivity
 *----------------------------------------------------------------------------*/

static void
_metis_cell_cells(size_t      n_cells,
                  size_t      n_faces,
                  cs_int_t   *face_cells,
                  idxtype   **cell_idx,
                  idxtype   **cell_neighbors)
{
  size_t i, id_1, id_2;

  cs_int_t  cell_id[2];

  idxtype  *n_neighbors;
  idxtype  *_cell_idx;
  idxtype  *_cell_neighbors ;

  /* Count and allocate arrays */

  BFT_MALLOC(n_neighbors, n_cells, idxtype);

  for (i = 0; i < n_cells; i++)
    n_neighbors[i] = 0;

  for (i = 0; i < n_faces; i++) {

    cell_id[0] = face_cells[i*2    ];
    cell_id[1] = face_cells[i*2 + 1];

    if (cell_id[0] == 0 || cell_id[1] == 0 || cell_id[0] == cell_id[1])
      continue;

    id_1 = cell_id[0] - 1;
    id_2 = cell_id[1] - 1;

    n_neighbors[id_1] += 1;
    n_neighbors[id_2] += 1;

  }

  BFT_MALLOC(_cell_idx, n_cells + 1, idxtype);

  _cell_idx[0] = 1 ;

  for (i = 0; i < n_cells; i++)
    _cell_idx[i + 1] = _cell_idx[i] + n_neighbors[i];

  BFT_MALLOC(_cell_neighbors, _cell_idx[n_cells], idxtype);

  for (i = 0; i < n_cells; i++)
    n_neighbors[i] = 0;

  for (i = 0; i < n_faces; i++) {

    cell_id[0] = face_cells[i*2    ];
    cell_id[1] = face_cells[i*2 + 1];

    if (cell_id[0] == 0 || cell_id[1] == 0 || cell_id[0] == cell_id[1])
      continue;

    id_1 = cell_id[0] - 1;
    id_2 = cell_id[1] - 1;

    _cell_neighbors[_cell_idx[id_1] + n_neighbors[id_1] - 1] = id_2 + 1;
    _cell_neighbors[_cell_idx[id_2] + n_neighbors[id_2] - 1] = id_1 + 1;

    n_neighbors[id_1] += 1;
    n_neighbors[id_2] += 1;

  }

  BFT_FREE(n_neighbors);

  *cell_idx = _cell_idx;
  *cell_neighbors = _cell_neighbors;

}

/*----------------------------------------------------------------------------
 * Compute partition using Metis
 *
 * parameters:
 *   n_cells       <-- number of cells in mesh
 *   n_parts       <-- number of partitions
 *   cell_cell_idx <-- cell->cells index
 *   cell_cell     <-- cell->cells connectivity
 *   cell_part     --> cell partition
 *----------------------------------------------------------------------------*/

static void
_part_metis(size_t    n_cells,
            int       n_parts,
            idxtype  *cell_idx,
            idxtype  *cell_neighbors,
            int      *cell_part)
{
  double  start_time[2], end_time[2];

  int     wgtflag    = 0; /* No weighting for faces or cells */
  int     numflag    = 1; /* 1 to n numbering (Fortran type) */
  int     options[5] = {0, 3, 1, 1, 0}; /* By default if options[0] = 0 */
  int     edgecut    = 0; /* <-- Number of faces on partition */

  int       _n_cells = n_cells;
  idxtype  *_cell_part = NULL;

  start_time[0] = bft_timer_wtime();
  start_time[1] = bft_timer_cpu_time();

  if (sizeof(idxtype) == sizeof(int))
    _cell_part = cell_part;

  else
    BFT_MALLOC(_cell_part, n_cells, idxtype);

  if (n_parts < 8) {

    bft_printf(_("Partitioning %d cells to %d domains"
                 " (METIS_PartGraphRecursive).\n"),
               (int)n_cells, n_parts);

    METIS_PartGraphRecursive(&_n_cells,
                             cell_idx,
                             cell_neighbors,
                             NULL,       /* vwgt:   cell weights */
                             NULL,       /* adjwgt: face weights */
                             &wgtflag,
                             &numflag,
                             &n_parts,
                             options,
                             &edgecut,
                             _cell_part);

  }

  else {

    bft_printf(_("Partitioning %d cells to %d domains"
                 " (METIS_PartGraphKway).\n"),
               (int)n_cells, n_parts);

    METIS_PartGraphKway(&_n_cells,
                        cell_idx,
                        cell_neighbors,
                        NULL,       /* vwgt:   cell weights */
                        NULL,       /* adjwgt: face weights */
                        &wgtflag,
                        &numflag,
                        &n_parts,
                        options,
                        &edgecut,
                        _cell_part);

  }

  end_time[0] = bft_timer_wtime();
  end_time[1] = bft_timer_cpu_time();

  bft_printf(_("\n"
               "  Total number of faces on parallel boundaries: %d\n"
               "  wall-clock time: %f s; CPU time: %f s\n\n"),
             (int)edgecut,
             (double)(end_time[0] - start_time[0]),
             (double)(end_time[1] - start_time[1]));

  if (sizeof(idxtype) != sizeof(int)) {
    size_t i;
    for (i = 0; i < n_cells; i++)
      cell_part[i] = _cell_part[i];
    BFT_FREE(_cell_part);
  }

  _cell_part_histogram(n_cells, n_parts, cell_part);

}

#endif /* defined(HAVE_METIS) */

/*----------------------------------------------------------------------------
 * Write padding zeroes to ensure following alignement if necessary
 *
 * parameters:
 *   outp      <-> pointer to output object
 *   alignment <-- required alignement
 *----------------------------------------------------------------------------*/

static void
_write_pad(_cs_io_t  *outp,
           size_t     alignment)
{
  bft_file_off_t offset = bft_file_tell(outp->f);
  size_t pad_size = (alignment - (offset % alignment)) % alignment;

  if (pad_size > 0) {

    char padding[128] = "";
    size_t rem_size = pad_size;

    memset(padding, 0, sizeof(padding));

    while (rem_size > 0) {
      size_t write_size = sizeof(padding);
      if (rem_size < write_size)
        write_size = rem_size;
      bft_file_write(padding, 1, write_size, outp->f);
      rem_size -= write_size;
    }

  }
}

/*----------------------------------------------------------------------------
 * Convert a buffer of type size_t to uint64_t
 *
 * parameters:
 *   buf --> buffer
 *   val <-- array from which values are converted
 *   n   <-- number of values to convert
 *----------------------------------------------------------------------------*/

static void
_convert_from_size(unsigned char  buf[],
                   const size_t   val[],
                   size_t         n)
{
  size_t i;

#if (__STDC_VERSION__ >= 199901L)

  for (i = 0; i < n; i++)
    ((uint64_t *)buf)[i] = val[i];

#else

  if (sizeof(size_t) == 8) {
    for (i = 0; i < n; i++)
      ((size_t *)buf)[i] = val[i];
  }
  else if (sizeof(unsigned long long) == 8) {
    for (i = 0; i < n; i++)
      ((unsigned long long *)buf)[i] = val[i];
  }
  else
    bft_error(__FILE__, __LINE__, 0,
              _("Compilation configuration / porting error:\n"
                "Unable to determine a 64-bit unsigned int type.\n"
                "size_t is %d bits, unsigned long long %d bits"),
              sizeof(size_t)*8, sizeof(unsigned long long)*8);

#endif
}

/*----------------------------------------------------------------------------
 * Open output file.
 *
 * parameters:
 *   n_ranks <-- number of ranks corresonding to output file
 *
 * returns:
 *   opened file structure wrapper
 *----------------------------------------------------------------------------*/

static _cs_io_t
_open_output(int n_ranks)
{
  size_t i;
  int n_ranks_size;

  size_t alignments[3];
  char header_buf[65];
  char header[64] = "";
  char *filename = NULL;

  _cs_io_t outp;

  outp.f = NULL;
  outp.header_size = 96;
  outp.header_align = 64;
  outp.body_align = 64;
  outp.buffer_size = 0;
  outp.buffer = NULL;

  outp.n_values = 0;
  outp.type_size = 0;
  outp.name = NULL;
  outp.type_name = NULL;
  outp.data = NULL;

  /* Build file name */

  for (i = n_ranks, n_ranks_size = 1;
       i >= 10;
       i /= 10, n_ranks_size += 1);

  BFT_MALLOC(filename,
             strlen("domain_number_") + n_ranks_size + 1,
             char);

  sprintf(filename, "domain_number_%d", n_ranks);

  /* Create file */

  bft_printf(_("  Opening output: %s\n"), filename);

  bft_printf_flush();

  outp.f = bft_file_open(filename,
                         BFT_FILE_MODE_WRITE,
                         BFT_FILE_TYPE_BINARY);

  BFT_FREE(filename);

  bft_file_set_big_endian(outp.f);

  /* Write header and comment information */

  memset(header, 0, 64);

  /* BE stands for big-endian, allowing for future
     native file mode generation, using BE/LE */

  strncpy(header, "Code_Saturne I/O, BE, R0", 63);

  bft_file_write(header, 1, 64, outp.f);

  memset(header, 0, 64);
  strncpy(header, "Domain partitioning, R0", 63);

  bft_file_write(header, 1, 64, outp.f);

  alignments[0] = outp.header_size;
  alignments[1] = outp.header_align;
  alignments[2] = outp.body_align;

  _convert_from_size((unsigned char*)header_buf, alignments, 3);

  bft_file_write(header_buf, 8, 3, outp.f);

  /* Allocate header buffer for future use */

  outp.buffer_size = outp.header_size;
  BFT_MALLOC(outp.buffer, outp.buffer_size, unsigned char);

  /* Finish */

  return outp;
}

/*----------------------------------------------------------------------------
 * Write section header.
 *
 * parameters:
 *   outp        <-> pointer to output object
 *   name        <-- section name
 *   n_vals      <-- number of associated values
 *   location_id <-- optional location id, or 0
 *   size_val    <-- optional embedded size value (1 value), or NULL
 *----------------------------------------------------------------------------*/

static void
_write_section_header(_cs_io_t      *outp,
                      const char    *name,
                      size_t         n_vals,
                      size_t         location_id,
                      const size_t  *size_val)
{
  size_t name_size = 0, name_pad_size = 0;
  size_t header_vals[6];
  size_t type_size = sizeof(cs_int_t);
  unsigned int_endian = 0;

  *((char *)(&int_endian)) = '\1'; /* Determine if we are little-endian */

  assert(outp != NULL);
  assert(outp->f != NULL);
  assert(outp->buffer != NULL);

  /* Position write pointer if necessary */
  /*-------------------------------------*/

  _write_pad(outp, outp->header_align);

  /* Prepare header */
  /*----------------*/

  header_vals[0] = 56;
  header_vals[1] = n_vals;
  header_vals[2] = location_id;
  header_vals[3] = 0;
  if (n_vals > 0)
    header_vals[4] = 1;
  else
    header_vals[4] = 0;

  name_size = strlen(name);
  name_pad_size = 8-(name_size%8); /* At least 1 NULL
                                      character with this rule */

  header_vals[5] = name_size + name_pad_size;
  header_vals[0] += (name_size + name_pad_size);

  if (n_vals == 1 && size_val != NULL) {
    header_vals[0] += 8;
    type_size = 8;
  }

  if (header_vals[0] > outp->buffer_size) {
    while (header_vals[0] > outp->buffer_size)
      outp->buffer_size *=2;
    BFT_REALLOC(outp->buffer, outp->buffer_size, unsigned char);
  }

  memset(outp->buffer, 0, outp->buffer_size);

  _convert_from_size(outp->buffer, header_vals, 6);

  if (int_endian == 1)
    bft_file_swap_endian(outp->buffer, outp->buffer, 8, 6);

  /* Element type name */

  outp->type_name = (char *)outp->buffer + 48;

  if (n_vals > 0) {
    outp->type_name[0] = 'i';
    if (size_val != NULL) {
      outp->type_name[7] = 'e';
      outp->type_name[1] = '8';
    }
    else {
      switch(sizeof(int)) {
      case 4:
        outp->type_name[1] = '4';
        break;
      case 8:
        outp->type_name[1] = '8';
        break;
      default:
        assert(sizeof(int) == 4 || sizeof(int) == 8);
      }
    }
  }

  /* Section name */

  outp->name = (char *)outp->buffer + 56;

  strcpy(outp->name, name);

  /* Optional embedded size */

  if (n_vals == 1 && size_val != NULL) {

    size_t bshift = header_vals[0] - 8;

    _convert_from_size(outp->buffer + bshift, size_val, 1);

    if (int_endian == 1)
      bft_file_swap_endian(outp->buffer + bshift,
                           outp->buffer + bshift, 8, 1);

  }

  /* Write prepared header to file */
  /*-------------------------------*/

  bft_file_write(outp->buffer, 1, header_vals[0], outp->f);
}

/*----------------------------------------------------------------------------
 * Close output file
 *
 * parameters:
 *   f <-> pointer to file object
 *----------------------------------------------------------------------------*/

static void
_close_output(_cs_io_t *outp)
{
  if (outp != NULL) {
    if (outp->f != NULL) {
      bft_printf(_("\n  Closing output: \"%s\"\n\n"),
                 bft_file_get_name(outp->f)) ;
      outp->f = bft_file_free(outp->f);
    }
    outp->header_size = 0;
    outp->header_align = 0;
    outp->body_align = 0;
    outp->buffer_size = 0;
    BFT_FREE(outp->buffer);
  }
}

/*----------------------------------------------------------------------------
 * Write a cell domain array to the output file.
 *
 * parameters:
 *   outp        <-> pointer to output object
 *   n           <-- number of me
 *   n_vals      <-- number of associated values
 *   location_id <-- optional location id, or 0
 *   size_val    <-- optional embedded size value (1 value), or NULL
 *----------------------------------------------------------------------------*/

static void
_write_domain(_cs_io_t    *outp,
              size_t       n,
              const int   *domain_num)
{
  assert(outp != NULL);
  assert(outp->f != NULL);

  /* Write header */

  _write_section_header(outp, "cell:domain number", n, 1, NULL);

  /* Position write pointer if necessary, then write */

  _write_pad(outp, outp->body_align);

  bft_file_write(domain_num, sizeof(int), n, outp->f);
}

/*----------------------------------------------------------------------------
 * Write output file.
 *
 * parameters:
 *   n_cells    <-- number of cells corresonding to mesh
 *   n_ranks    <-- number of ranks corresonding to output file
 *   domain_num <-- domain number array (size: n_cells)
 *----------------------------------------------------------------------------*/

static void
_write_output(size_t      n_cells,
              int         n_ranks,
              const int  *domain_num)
{
  size_t _n_cells = n_cells;
  size_t _n_ranks = n_ranks;

  _cs_io_t outp = _open_output(n_ranks);

  _write_section_header(&outp, "n_cells", 1, 1, &_n_cells);
  _write_section_header(&outp, "n_ranks", 1, 0, &_n_ranks);

  _write_domain(&outp, n_cells, domain_num);

  _close_output(&outp);
}

/*----------------------------------------------------------------------------
 * Print usage and exit.
 *
 * parameters:
 *   arg_0     <-- name of executable as given by argv[0]
 *   exit_code <-- EXIT_SUCCESS or EXIT_FAILURE
 *----------------------------------------------------------------------------*/

static void
_usage(const char  *arg_0,
       int          exit_code)
{
  bft_printf
    (_("\n"
       "  Usage: %s <n_ranks> [<n_ranks_2> ... <n_ranks_n>] [-s]\n\n"
       "  Compute domain partition(s) of Code_Saturne Preprocessor output.\n\n"
       "  <n_ranks> number of ranks of destination partition (> 1).\n\n"
       "  -no-perio ignore periodicity information if present.\n"
       "  -no-write do not write output.\n"
       "  -h        this message.\n\n"),
     arg_0);

  exit(exit_code);
}

/*----------------------------------------------------------------------------
 * Read command line arguments.
 *
 * parameters:
 *   argc     <-- number of command line arguments
 *   argv     <-- array of command line arguments
 *   no_write --> do not write output if 1
 *   no_perio --> ignore periodicity information if 1
 *   n_parts  --> number of partitionings required
 *   n_ranks  --> array of required ranks per partitioning
 *----------------------------------------------------------------------------*/

static void
_read_args(int               argc,
           char            **argv,
           int              *no_write,
           int              *no_perio,
           int              *n_parts,
           int             **n_ranks)
{
  int _n_ranks;

  int i = 1;

  /* Initialize return arguments */

  *no_write = 0;
  *no_perio  = 0;
  *n_parts   = 0;
  *n_ranks   = NULL;

  /* Parse and check command line */

  if (argc < 2)
    _usage(argv[0], EXIT_FAILURE);

  while (i < argc) {

    if (strcmp(argv[i], "-h") == 0)
      _usage(argv[0], EXIT_SUCCESS);

    else if (strcmp(argv[i], "-no-write") == 0)
      *no_write = 1;

    else if (strcmp(argv[i], "-no-perio") == 0)
      *no_perio = 1;

    else {
      _n_ranks = atoi(argv[i]);

      if (_n_ranks <= 1)
        _usage(argv[0], EXIT_FAILURE);

      else {
        BFT_REALLOC(*n_ranks, *n_parts + 1, int);
        (*n_ranks)[*n_parts] = _n_ranks;
        *n_parts += 1;
      }

    }

    i++;
  }

  if (*n_parts < 1)
    _usage(argv[0], EXIT_FAILURE);

  /* At this point, command line seems correct */

  bft_printf(_("\n"
               "  .------------------------------.\n"
               "  |   Code_Saturne Partitioner   |\n"
               "  `------------------------------'\n"));

  bft_printf(_("\n  External library:"));
#if defined(HAVE_METIS)
  bft_printf(" METIS 4.0");
#endif
  bft_printf("\n");
}

/*============================================================================
 * Public function definitions
 *============================================================================*/

int
main (int argc, char *argv[])
{
  double     start_time[2];
  size_t     n_cells = 0, n_faces = 0;
  cs_int_t  *face_cells = NULL;

  int no_write = 0;
  int no_perio = 0;
  int n_parts  = 0;
  int *n_ranks = NULL;

  if (getenv("LANG") != NULL)
     setlocale(LC_ALL,"");
  else
     setlocale(LC_ALL,"C");
  setlocale(LC_NUMERIC,"C");

#if defined(ENABLE_NLS)
  bindtextdomain (PACKAGE, LOCALEDIR);
  textdomain(PACKAGE);
#endif

  /* Initialize memory management */

#if defined(ECS_ARCH_Linux)
  bft_mem_usage_set_options(BFT_MEM_USAGE_TRACK_PR_SIZE |
                            BFT_MEM_USAGE_TRACK_ALLOC_SIZE);
#endif
  bft_mem_usage_init();

  bft_mem_init(getenv("CS_PARTITION_MEM_LOG"));

  /* Initialize timers */

  start_time[0] = bft_timer_wtime();
  start_time[1] = bft_timer_cpu_time();

  /* Parse command line arguments */

  _read_args(argc, argv, &no_write, &no_perio, &n_parts, &n_ranks);

  /* Build cell-cell connectivity information */

  _read_input(no_perio, &n_cells, &n_faces, &face_cells);

#if defined(HAVE_METIS)

  {
    int i;
    int *cell_part = NULL;
    idxtype *cell_idx = NULL, *cell_neighbors = NULL;

    _metis_cell_cells(n_cells,
                      n_faces,
                      face_cells,
                      &cell_idx,
                      &cell_neighbors);

    if (face_cells != NULL)
      BFT_FREE(face_cells);

    /* Compute partitionings */

    BFT_MALLOC(cell_part, n_cells, int);

    for (i = 0; i < n_parts; i++) {

      _part_metis(n_cells,
                  n_ranks[i],
                  cell_idx,
                  cell_neighbors,
                  cell_part);

      if (!no_write)
        _write_output(n_cells, n_ranks[i], cell_part);

    }

    BFT_FREE(cell_part);

    BFT_FREE(cell_idx);
    BFT_FREE(cell_neighbors);
  }

#endif

  /* Finalize */

  BFT_FREE(n_ranks);

  /* Timing and Memory use summary */

  {
    int    n_div;
    double f_size;
    double end_time[2];
    size_t pr_size;
    static char  unit_prefix[] = {'K', 'M', 'G'} ;

    bft_printf(_("Memory use and timing summary:\n\n")) ;

    pr_size = bft_mem_usage_max_pr_size() ;

    if (pr_size > 0) {
      f_size = pr_size ;
      for (n_div = 0 ; f_size > 1024. && n_div < 2 ; n_div++)
        f_size /= 1024. ;
      bft_printf(_("  Total memory:    %15.3f %cb\n\n"),
                 f_size, unit_prefix[n_div]);
    }

    /* Timers */

    end_time[0] = bft_timer_wtime();
    end_time[1] = bft_timer_cpu_time();

    bft_printf(_("  Wall-clock time: %15.3f s\n"
                 "  CPU time:        %15.3f s\n"),
               (double)(end_time[0] - start_time[0]),
               (double)(end_time[1] - start_time[1]));

  }

  bft_printf(_("\n"
               "  .------------------------.\n"
               "  |   Partitioner finish   |\n"
               "  `------------------------'\n\n"));

  bft_mem_end();

  exit(EXIT_SUCCESS);
}

/*----------------------------------------------------------------------------*/

#ifdef __cplusplus
}
#endif /* __cplusplus */

Generated by  Doxygen 1.6.0   Back to index