parteleser.c

Go to the documentation of this file.

/*
    This file is part of magpar.

    Copyright (C) 2002-2010 Werner Scholz

    www:   http://www.magpar.net/
    email: magpar(at)magpar.net

    magpar 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.

    magpar 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 magpar; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

static char const Id[] = "$Id: parteleser.c 2962 2010-02-04 19:50:44Z scholz $\n\n";
static char const Td[] = "$Today: " __FILE__ "  " __DATE__ "  "  __TIME__  " $\n\n";

#define PREALLOC_DG 25

#include "init.h"
#include "util/util.h"

#ifdef METIS
EXTERN_C_BEGIN
#include "metis.h"
EXTERN_C_END
#endif

int MetisPartition(GridData *gdata,int parts)
{
  MagparFunctionLogBegin;

  int rank,size;
  ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);

  if (rank) {
    MagparFunctionLogReturn(0);
  }

  PetscPrintf(PETSC_COMM_WORLD,
    "Meshpartitioning using METIS for %i processors\n",
    parts
  );

#ifdef METIS
  /* directly use METIS_PartGraphKWay, since we do already have ia and ja.
     Now we do it twice, but it does not make much difference since Metis is
     so fast.
  */
/* partition elements
    METIS_PartMeshDual(
*/
/* partition nodes */
  int t_edgecut;
  int zero=0,two=2;
  METIS_PartMeshNodal(
    &gdata->n_ele,
    &gdata->n_vert,
    gdata->elevert,
    &two,
    &zero,
    &parts,
    &t_edgecut,
    gdata->elenewproc,
    gdata->vertnewproc
  );
#else
  SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Metis library not linked in, but it is required for volume mesh partitioning.\n");
#endif

  /* create trivial local->global mappings */
  for (int i=0; i<gdata->n_vert; i++) {
    gdata->vertl2g[i]=i;
  }
  for (int i=0; i<gdata->n_ele; i++) {
    gdata->elel2g[i]=i;
  }

  /* Create a global reordering (mapping) of the vertices and elements
     by permuting the trivial mapping
  */
  ascat(&gdata->vertl2g,gdata->n_vert,1,gdata->vertnewproc,parts,0);
  ascat(&gdata->elel2g,gdata->n_ele,1,gdata->elenewproc,parts,0);

  MagparFunctionLogReturn(0);
}

/*
   Given the grid data, determines a new partitioning of
   the CELLS (elements) to reduce the number of cut edges between
   cells (elements).
*/

int DataPartitionElementsSer(GridData *gdata)
{
  MagparFunctionLogBegin;

  int rank,size;
  ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);

  if (rank) {
    MagparFunctionLogReturn(0);
  }

  PetscTruth flg;
  int dometispartition;
  ierr = PetscOptionsGetInt(PETSC_NULL,"-metispartition",&dometispartition,&flg);CHKERRQ(ierr);
  if (flg!=PETSC_TRUE) {
    dometispartition=0;
    PetscPrintf(PETSC_COMM_WORLD,
      "Option -metispartition not found, using default value %i\n",
      dometispartition
    );
  }

  if (dometispartition<0) {
    SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"-metispartition invalid! Exiting!\n");
  }
  else if (dometispartition==0) {
    PetscPrintf(PETSC_COMM_WORLD,"METIS not used for mesh partitioning.\n");
    MagparFunctionLogReturn(0);
  }

  if (dometispartition==1) {
    dometispartition=size;
  }
  if (dometispartition==1) {
    PetscPrintf(PETSC_COMM_WORLD,"Mesh partitioning not performed since not necessary on one processor.\n");
    MagparFunctionLogReturn(0);
  }
  if (dometispartition<size) {
    PetscPrintf(PETSC_COMM_WORLD,"%i partitions for %i processors: some processors unused!\n",dometispartition,size);
  }

  ierr = MetisPartition(gdata,dometispartition);CHKERRQ(ierr);
  ierr = CheckPartition(gdata,dometispartition);CHKERRQ(ierr);

  /* remap to real number of processors */

  if (dometispartition>size) {
    PetscPrintf(PETSC_COMM_WORLD,
      "Remapping %i partitions onto %i processors.\n",
      dometispartition,size
    );
    for (int i=0;i<gdata->n_ele;i++) {
      gdata->elenewproc[i]= distint(dometispartition,size,gdata->elenewproc[i]);
    }
    for (int i=0;i<gdata->n_vert;i++) {
      gdata->vertnewproc[i]=distint(dometispartition,size,gdata->vertnewproc[i]);
    }
  }

  /* Permute data according to new numbering */
  IS isvert,isele;
  ierr = ISCreateGeneral(PETSC_COMM_SELF,gdata->n_vert,gdata->vertl2g,&isvert);CHKERRQ(ierr);
  ierr = ISSetPermutation(isvert);CHKERRQ(ierr);
  ierr = ISCreateGeneral(PETSC_COMM_SELF,gdata->n_ele,gdata->elel2g,&isele);CHKERRQ(ierr);
  ierr = ISSetPermutation(isele);CHKERRQ(ierr);

  /* CheckPartition done above */
  ierr = PermuteData(gdata,isvert,isele);CHKERRQ(ierr);

  ierr = ISDestroy(isvert);CHKERRQ(ierr);
  ierr = ISDestroy(isele);CHKERRQ(ierr);

  MagparFunctionLogReturn(0);
}