SparseGridsInterpolate.cpp

Go to the documentation of this file.

#include <arrayfunctions.h>
#include <mathfunctions_cpp.h>
#include <std_vec_ops.h>
#include <sparse_grids_simulation.h>

void SparseGridsSimulation::interpolate(  SimulationObject* const       a_dst,  
                                          const SimulationObject* const a_src   
                                      )
{
  /* get the destination grid dimensions */
  GridDimensions dim_dst;
  StdVecOps::copyFrom(dim_dst, a_dst->solver.dim_global, m_ndims);

  /* get number of vector components */
  int nvars = a_src->solver.nvars;
  if (nvars != a_dst->solver.nvars) {
    fprintf(stderr, "Error in SparseGridsSimulation::interpolate(): unequal nvars\n");
    exit(1);
  }

  double *ug_dst = NULL;
  interpolate(dim_dst, &ug_dst, a_src); 

  /* partition the destination */
  MPIPartitionArraynDwGhosts( m_ndims,
                              (void*) &(a_dst->mpi),
                              (a_dst->mpi.rank ? NULL : ug_dst),
                              a_dst->solver.u,
                              a_dst->solver.dim_global,
                              a_dst->solver.dim_local,
                              a_dst->solver.ghosts,
                              nvars); 

  if (!m_rank) {
    free(ug_dst);
  }

  return;
}

void SparseGridsSimulation::interpolate(const GridDimensions& a_dim_dst, 
                                        double** const a_u_dst, 
                                        const SimulationObject* const a_src 
                                       )
{
  if ((*a_u_dst) != NULL) {
    fprintf(stderr, "Error in SparseGridsSimulation::interpolate() - \n");
    fprintf(stderr, "  a_u_dst is not NULL!\n");
    exit(1);
  }

  /* get the source grid dimensions */
  GridDimensions dim_src;
  StdVecOps::copyFrom(dim_src, a_src->solver.dim_global, m_ndims);

  /* get number of vector components and number of ghost points*/
  int nvars = a_src->solver.nvars;
  int ghosts = a_src->solver.ghosts;

  /* gather the source on rank 0 */
  double *ug_src = NULL;
  if (!m_rank) allocateDataArrays(dim_src, nvars, &ug_src, ghosts);
  MPIGatherArraynDwGhosts( m_ndims,
                           (void*) &(a_src->mpi),
                           ug_src,
                           a_src->solver.u,
                           a_src->solver.dim_global,
                           a_src->solver.dim_local,
                           ghosts,
                           nvars );

  std::vector<int> periodic_arr(m_ndims);
  for (int i=0; i<m_ndims; i++) {
    periodic_arr[i] = (m_is_periodic[i] ? 1 : 0);
  }

  if (!m_rank) {
    int ierr = ::InterpolateGlobalnDVar(  a_dim_dst.data(),
                                          a_u_dst,
                                          dim_src.data(),
                                          ug_src,
                                          nvars,
                                          ghosts,
                                          m_ndims,
                                          periodic_arr.data());
    if (ierr) {
      fprintf(stderr,"InterpolateGlobalnDVar() returned with error!\n");
      exit(1);
    }
  }

  return;
}