PetscGlobalDOF.cpp

Go to the documentation of this file.

#ifdef with_petsc

#include <stdlib.h>
#include <vector>
#include <basic.h>
#include <arrayfunctions.h>
#include <mpivars_cpp.h>
#include <simulation_object.h>
#include <petscinterface.h>

static int ApplyPeriodicity(  int     dir,    
                              int     ndims,  
                              int     *size,  
                              int     ghosts, 
                              double  *phi     )
{
  int bounds[ndims], index1[ndims], index2[ndims], offset[ndims], 
      done, p1 = 0, p2 = 0;
  _ArrayCopy1D_(size,bounds,ndims); bounds[dir] = ghosts;

  done = 0; _ArraySetValue_(index1,ndims,0);
  while (!done) {
    _ArraySetValue_(offset,ndims,0); offset[dir] = -ghosts;
    _ArrayIndex1DWO_(ndims,size,index1,offset,ghosts,p1);
    _ArrayCopy1D_(index1,index2,ndims);
    index2[dir] = index1[dir] + size[dir]-ghosts;
    _ArrayIndex1D_(ndims,size,index2,ghosts,p2);

    phi[p1] = phi[p2];
    _ArrayIncrementIndex_(ndims,bounds,index1,done);
  }

  done = 0; _ArraySetValue_(index1,ndims,0);
  while (!done) {
    _ArraySetValue_(offset,ndims,0); offset[dir] = size[dir];
    _ArrayIndex1DWO_(ndims,size,index1,offset,ghosts,p1);
    _ArrayIndex1D_(ndims,size,index1,ghosts,p2);

    phi[p1] = phi[p2];
    _ArrayIncrementIndex_(ndims,bounds,index1,done);
  }
  return(0);
}

int PetscGlobalDOF(void* c )
{
  PETScContext* ctxt = (PETScContext*) c;
  SimulationObject* sim = (SimulationObject*) ctxt->simobj;
  int nsims = ctxt->nsims;

  ctxt->globalDOF.resize(nsims, nullptr);

  /* compute MPI offset */
  std::vector<int> local_sizes(ctxt->nproc ,0);
  local_sizes[ctxt->rank] = ctxt->npoints;
  MPIMax_integer(local_sizes.data(),local_sizes.data(),ctxt->nproc,&sim[0].mpi.world);

  int MPIOffset = 0;
  for (int i=0; i<ctxt->rank; i++) MPIOffset += local_sizes[i];

  int simOffset = 0;
  for (int ns = 0; ns < nsims; ns++) {

    HyPar* solver = &(sim[ns].solver);
    MPIVariables* mpi = &(sim[ns].mpi);

    int   *dim        = solver->dim_local,
          ndims       = solver->ndims,
          ghosts      = solver->ghosts,
          npoints     = solver->npoints_local,
          npoints_wg  = solver->npoints_local_wghosts,
          nv          = ndims + 1, i;

    ctxt->globalDOF[ns] = (double*) calloc( npoints_wg, sizeof(double) );
    _ArraySetValue_(ctxt->globalDOF[ns], npoints_wg, -1.0);

    for (int i = 0; i < npoints; i++) {
      int p = (ctxt->points[ns]+i*nv)[ndims];
      ctxt->globalDOF[ns][p] = (double) (i + simOffset + MPIOffset);
    }

    for (int i=0; i<ndims; i++) {
      if (solver->isPeriodic[i]) {
        ApplyPeriodicity(i,ndims,dim,ghosts,ctxt->globalDOF[ns]);
      }
    }
    MPIExchangeBoundariesnD(ndims,1,dim,ghosts,mpi,ctxt->globalDOF[ns]);

    simOffset += npoints;
  }

  return 0;
}

#endif