movedata.c

Go to the documentation of this file.

/*
    This file is part of magpar.

    Copyright (C) 2002-2009 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: movedata.c 2700 2009-08-07 21:07:13Z scholz $\n\n";
static char const Td[] = "$Today: " __FILE__ "  " __DATE__ "  "  __TIME__  " $\n\n";

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

int DataMoveData(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 (size==1) {
#ifdef HAVE_ELEVERTALL
    gdata->elevertall=gdata->elevert;
#endif

    int *ia,*ja;

    ierr = Mesh2Nodal(gdata->n_ele,gdata->n_vert,gdata->elevert,&ia,&ja);CHKERRQ(ierr);
    gdata->mesh2nodal_ia=ia;
    gdata->mesh2nodal_ja=ja;

    ierr = Mesh2Dual(gdata->n_ele,gdata->n_vert,gdata->elevert,&ia,&ja);CHKERRQ(ierr);
    gdata->mesh2dual_ia=ia;
    gdata->mesh2dual_ja=ja;

    MagparFunctionLogReturn(0);
  }

  ierr = MPI_Bcast(&gdata->n_prop,1,MPI_INT,0,PETSC_COMM_WORLD);CHKERRQ(ierr);
  if (rank) {
    ierr = PetscMalloc(gdata->n_prop*NP*sizeof(PetscReal),&gdata->propdat);CHKERRQ(ierr);
  }
  ierr = MPI_Bcast(gdata->propdat,gdata->n_prop*NP,MPIU_SCALAR,0,PETSC_COMM_WORLD);CHKERRQ(ierr);

  ierr = MPI_Bcast(&gdata->n_ele,1,MPI_INT,0,PETSC_COMM_WORLD);CHKERRQ(ierr);

#ifdef HAVE_ELEVERTALL
  ierr = PetscMalloc(NV*gdata->n_ele*sizeof(int),&gdata->elevertall);CHKERRQ(ierr);
  if (!rank) {
    ierr = PetscMemcpy(gdata->elevertall,gdata->elevert,NV*gdata->n_ele*sizeof(int));CHKERRQ(ierr)
  }
  ierr = MPI_Bcast(gdata->elevertall,NV*gdata->n_ele,MPI_INT,0,PETSC_COMM_WORLD);CHKERRQ(ierr);
#ifdef HAVE_M2IJ
  int *ia,*ja;
  if (!rank) {
    ierr = Mesh2Nodal(gdata->n_ele,gdata->n_vert,gdata->elevert,&ia,&ja);CHKERRQ(ierr);
    gdata->mesh2nodal_ia=ia;
    gdata->mesh2nodal_ja=ja;

    ierr = Mesh2Dual(gdata->n_ele,gdata->n_vert,gdata->elevert,&ia,&ja);CHKERRQ(ierr);
    gdata->mesh2dual_ia=ia;
    gdata->mesh2dual_ja=ja;
  }
  else {
    ierr = PetscMalloc((gdata->n_vert+1)*sizeof(int),&gdata->mesh2nodal_ia);CHKERRQ(ierr);
    ierr = PetscMalloc((gdata->n_ele+1)*sizeof(int),&gdata->mesh2dual_ia);CHKERRQ(ierr);
  }
  ierr = MPI_Bcast(gdata->mesh2nodal_ia,gdata->n_vert+1,MPI_INT,0,PETSC_COMM_WORLD);CHKERRQ(ierr);
  ierr = MPI_Bcast(gdata->mesh2dual_ia,gdata->n_ele+1,MPI_INT,0,PETSC_COMM_WORLD);CHKERRQ(ierr);

  if (rank) {
    ierr = PetscMalloc(gdata->mesh2nodal_ia[gdata->n_vert]*sizeof(int),&gdata->mesh2nodal_ja);CHKERRQ(ierr);
    ierr = PetscMalloc(gdata->mesh2dual_ia[gdata->n_ele]*sizeof(int),&gdata->mesh2dual_ja);CHKERRQ(ierr);
  }
  ierr = MPI_Bcast(gdata->mesh2nodal_ja,gdata->mesh2nodal_ia[gdata->n_vert],MPI_INT,0,PETSC_COMM_WORLD);CHKERRQ(ierr);
  ierr = MPI_Bcast(gdata->mesh2dual_ja,gdata->mesh2dual_ia[gdata->n_ele],MPI_INT,0,PETSC_COMM_WORLD);CHKERRQ(ierr);
#endif
#endif

  int t_gotele;
  t_gotele=ascat(&gdata->elevert,gdata->ln_ele,NV,gdata->elenewproc,size,1);
  int t_gotn;
  t_gotn=  ascat(&gdata->eleprop,gdata->ln_ele,1, gdata->elenewproc,size,1);assert(t_gotele==t_gotn);

  /* distribute boundary faces evenly over processors */
  ierr = MPI_Bcast(&gdata->n_vert_bnd,1,MPI_INT,0,PETSC_COMM_WORLD);CHKERRQ(ierr);
  ierr = MPI_Bcast(&gdata->n_bnd_fac, 1,MPI_INT,0,PETSC_COMM_WORLD);CHKERRQ(ierr);
  int *bndfacnewproc;
  ierr = PetscMalloc(gdata->n_bnd_fac*sizeof(int),&bndfacnewproc);CHKERRQ(ierr);
  if (!rank) {
    for (int i=0;i<gdata->n_bnd_fac;i++) {
      bndfacnewproc[i]=distint(gdata->n_bnd_fac,size,i);
    }
  }
  /* cannot move to facnb.c: partsurfser.c needs bndfacvert on proc 0 */
  int t_gotbndfac;
  t_gotbndfac=ascat(&gdata->bndfacvert,gdata->ln_bnd_fac,NN,bndfacnewproc,size,1);
  gdata->ln_bnd_fac=t_gotbndfac;
  ierr = PetscFree(bndfacnewproc);CHKERRQ(ierr);
  PetscSynchronizedPrintf(PETSC_COMM_WORLD,
    "<%i>Received local number of boundary faces: %i\n",
    rank,gdata->ln_bnd_fac
  );
  PetscSynchronizedFlush(PETSC_COMM_WORLD);
  int tsum;
  ierr=MPI_Allreduce((int*)&gdata->ln_bnd_fac,(int*)&tsum,1,MPI_INT,MPI_SUM,PETSC_COMM_WORLD);CHKERRQ(ierr);
  if (!rank) assert(gdata->n_bnd_fac==tsum);
  PetscPrintf(PETSC_COMM_WORLD,
    "Total number of boundary faces: %i\n",
    gdata->n_bnd_fac
  );

  /* Create a global reordering of the element numbers */
  int t_gotvert;
  t_gotvert=ascat(&gdata->vertl2g,gdata->ln_vert,1,gdata->vertnewproc,size,1);
  t_gotn=ascat(&gdata->elel2g,gdata->ln_ele,1,gdata->elenewproc,size,1);assert(t_gotele==t_gotn);

  /* broadcast vertxyz to all processors */
  ierr = MPI_Bcast(&gdata->n_vert,    1,MPI_INT,0,PETSC_COMM_WORLD);CHKERRQ(ierr);
  if (rank) {
    ierr = PetscMalloc(ND*gdata->n_vert*sizeof(PetscReal),&gdata->vertxyz);CHKERRQ(ierr);
  }
  ierr = MPI_Bcast(gdata->vertxyz,gdata->n_vert*ND,MPIU_SCALAR,0,PETSC_COMM_WORLD);CHKERRQ(ierr);

  /* broadcast all boundary vertex indicators to all processors
     (required for bmatrix)
  */
  if (rank) {
    ierr = PetscMalloc(gdata->n_vert*sizeof(int),&gdata->vertbndg2bnd);CHKERRQ(ierr);
  }
  ierr = MPI_Bcast(gdata->vertbndg2bnd,gdata->n_vert,MPI_INT,0,PETSC_COMM_WORLD);CHKERRQ(ierr);

  t_gotn=ascat(&gdata->vertprop,gdata->ln_vert,1,gdata->vertnewproc,size,1);assert(t_gotvert==t_gotn);

  /* scatter magnetization from the first processor
     into the distributed vector
  */

  Vec tv_M;
  ierr = VecCreate(PETSC_COMM_WORLD,&tv_M);CHKERRQ(ierr);
  ierr = VecSetSizes(tv_M,ND*t_gotvert,ND*gdata->n_vert);CHKERRQ(ierr);
  ierr = VecSetFromOptions(tv_M);CHKERRQ(ierr);
  ierr = VecSetBlockSize(tv_M,ND);CHKERRQ(ierr);

  IS tis_from,tis_to;
  if (!rank) {
    /* preserve the ordering since the data have already been permuted above */
    ierr = ISCreateStride(PETSC_COMM_SELF,ND*gdata->n_vert,0,1,&tis_from);CHKERRQ(ierr);
    ierr = ISCreateStride(PETSC_COMM_SELF,ND*gdata->n_vert,0,1,&tis_to);CHKERRQ(ierr);

  }
  else {
    PetscInt t_idx_from[1],t_idx_to[1];
    tis_from=NULL;
    tis_to=NULL;
    ISCreateGeneral(PETSC_COMM_SELF,0,t_idx_from,&tis_from);
    ISCreateGeneral(PETSC_COMM_SELF,0,t_idx_to,  &tis_to);
  }

  VecScatter t_scatter;
  VecScatterCreate(gdata->M,tis_from,tv_M,tis_to,&t_scatter);
  ISDestroy(tis_from);
  ISDestroy(tis_to);

  VecScatterBegin(t_scatter,gdata->M,tv_M,INSERT_VALUES,SCATTER_FORWARD);
  VecScatterEnd  (t_scatter,gdata->M,tv_M,INSERT_VALUES,SCATTER_FORWARD);
  VecDestroy(gdata->M);
  gdata->M=tv_M;

#ifdef ADDONS
  /* broadcast vertex coupling information to all processors */
  int j;
  if (!rank) {
    if (gdata->vertcoupl!=PETSC_NULL) {
      j=1;
    }
    else {
      j=0;
    }
  }
  ierr = MPI_Bcast(&j,1,MPI_INT,0,PETSC_COMM_WORLD);CHKERRQ(ierr);

  if (j) {
    if (rank) {
      ierr = PetscMalloc(gdata->n_vert*sizeof(int),&gdata->vertcoupl);CHKERRQ(ierr);
    }

    ierr = MPI_Bcast(gdata->vertcoupl,gdata->n_vert,MPI_INT,0,PETSC_COMM_WORLD);CHKERRQ(ierr);
  }
  else {
    gdata->vertcoupl=PETSC_NULL;
  }

  if (gdata->VH1!=PETSC_NULL) {
    Vec tv_H1;
    ierr = VecDuplicate(tv_M,&tv_H1);CHKERRQ(ierr);
    VecScatterBegin(t_scatter,gdata->VH1,tv_H1,INSERT_VALUES,SCATTER_FORWARD);
    VecScatterEnd  (t_scatter,gdata->VH1,tv_H1,INSERT_VALUES,SCATTER_FORWARD);
    VecDestroy(gdata->VH1);
    gdata->VH1=tv_H1;
  }
  if (gdata->VH2!=PETSC_NULL) {
    Vec tv_H2;
    ierr = VecDuplicate(tv_M,&tv_H2);CHKERRQ(ierr);
    VecScatterBegin(t_scatter,gdata->VH2,tv_H2,INSERT_VALUES,SCATTER_FORWARD);
    VecScatterEnd  (t_scatter,gdata->VH2,tv_H2,INSERT_VALUES,SCATTER_FORWARD);
    VecDestroy(gdata->VH2);
    gdata->VH2=tv_H2;
  }
#endif
  VecScatterDestroy(t_scatter);

  t_gotn=ascat(&gdata->elenewproc,gdata->ln_ele,1,gdata->elenewproc,size,1);assert(t_gotele==t_gotn);
  t_gotn=ascat(&gdata->vertnewproc,gdata->ln_vert,1,gdata->vertnewproc,size,1);assert(t_gotvert==t_gotn);

  gdata->ln_ele=t_gotele;
  gdata->ln_vert=t_gotvert;

  PetscSynchronizedPrintf(PETSC_COMM_WORLD,"<%i>Received local number of elements: %i\n",rank,gdata->ln_ele);
  PetscSynchronizedPrintf(PETSC_COMM_WORLD,"<%i>Received local number of vertices: %i\n",rank,gdata->ln_vert);
  PetscSynchronizedFlush(PETSC_COMM_WORLD);

  /* check correct data transfer */
  for (int i=0; i<gdata->ln_ele; i++)
  {
    assert(gdata->elenewproc[i]==rank);
  }
  for (int i=0; i<gdata->ln_vert; i++)
  {
    assert(gdata->vertnewproc[i]==rank);
  }
  PetscSynchronizedFlush(PETSC_COMM_WORLD);

  int vidlow, vidhigh, vlocsize;
  ierr = VecGetOwnershipRange(gdata->M,(PetscInt*)&vidlow,(PetscInt*)&vidhigh);CHKERRQ(ierr);
  ierr = VecGetLocalSize(gdata->M,(PetscInt*)&vlocsize);CHKERRQ(ierr);
  vidlow   /= ND;
  vidhigh  /= ND;
  vlocsize /= ND;
  assert(vidhigh-vidlow==vlocsize);
  assert(vlocsize==gdata->ln_vert);
  PetscSynchronizedPrintf(PETSC_COMM_WORLD,"<%i>Ownership(M): from vertex id %i to vertex id %i\n",rank,vidlow,vidhigh-1);
  PetscSynchronizedFlush(PETSC_COMM_WORLD);

  ierr = PetscFree(gdata->vertnewproc);CHKERRQ(ierr);
  gdata->vertnewproc=PETSC_NULL;
  ierr = PetscFree(gdata->elenewproc);CHKERRQ(ierr);
  gdata->elenewproc=PETSC_NULL;

  MagparFunctionLogReturn(0);
}