SparseGridsSimulation Class Reference

Class describing sparse grids simulations. More...

#include <sparse_grids_simulation.h>

Inheritance diagram for SparseGridsSimulation:

Public Member Functions
	SparseGridsSimulation ()
	~SparseGridsSimulation ()
int	define (int, int)
int	ReadInputs ()
int	Initialize ()
int	InitialSolution ()
int	InitializeBoundaries ()
int	InitializeImmersedBoundaries ()
int	InitializePhysics ()
int	InitializePhysicsData ()
int	InitializeSolvers ()
int	InitializationWrapup ()
int	Solve ()
void	WriteErrors (double, double)
bool	isDefined () const
int	mpiCommDup ()
void	usePetscTS (PetscBool a_flag)
int	SolvePETSc ()
Public Member Functions inherited from Simulation
	Simulation ()
virtual	~Simulation ()

Protected Member Functions
int	SanityChecks ()
int	ComputeSGDimsAndCoeffs ()
int	ComputeProcessorDistribution (ProcDistribution &, const GridDimensions &)
void	GetCTGridSizes (const int, std::vector< GridDimensions > &)
int	InitializeBarebones (SimulationObject *)
int	InitializeSolversBarebones (SimulationObject *)
int	CleanupBarebones (SimulationObject *)
int	SetSolverParameters (SimulationObject &, const GridDimensions &, const ProcDistribution &, const SimulationObject &, const int, const int)
void	OutputSolution (double)
void	CombinationTechnique (SimulationObject *const)
void	interpolate (SimulationObject *const, const SimulationObject *const)
void	interpolate (const GridDimensions &, double **const, const SimulationObject *const)
void	interpolateGrid (SimulationObject *const, const SimulationObject *const)
void	coarsenGrid1D (const GridDimensions &, const GridDimensions &, const double *const, double *const, int)
void	refineGrid1D (const GridDimensions &, const GridDimensions &, const double *const, double *const, int)
void	fillGhostCells (const GridDimensions &, const int, double *const, const int)
void	CalculateError ()
void	computeError (SimulationObject &, double *const)
bool	isPowerOfTwo (int x)
double	cfunc (int a, int b)
int	factorial (int a)
void	allocateGridArrays (const GridDimensions &a_dim, double **const a_x, const int a_ngpt=0)
void	allocateDataArrays (const GridDimensions &a_dim, const int a_nvars, double **const a_u, const int a_ngpt=0)

Protected Attributes
bool	m_is_defined
int	m_nsims_sg
int	m_ndims
int	m_rank
int	m_nproc
int	m_interp_order
std::vector< bool >	m_is_periodic
int	m_write_sg_solutions
int	m_print_sg_errors
SimulationObject *	m_sim_fg
std::vector< SimulationObject >	m_sims_sg
int	m_n_fg
int	m_imin
std::vector< SGCTElem >	m_combination
std::vector< ProcDistribution >	m_iprocs

Detailed Description

Class describing sparse grids simulations.

Class describing sparse grids simulations This class contains all data and functions needed to run a sparse grids simulation.

This class contains all data and functions needed to run sparse grids simulations.

Definition at line 38 of file sparse_grids_simulation.h.

Constructor & Destructor Documentation

SparseGridsSimulation ( )

inline

Constructor

Definition at line 43 of file sparse_grids_simulation.h.

    {
      m_is_defined = false;

      m_nsims_sg = -1;
      m_ndims = -1;
      m_rank = -1;
      m_nproc = -1;

      m_sim_fg = NULL;
      m_sims_sg.clear();

      m_n_fg = -1;
      m_imin = -1;

      m_interp_order = 6;

      m_is_periodic.clear();
      m_combination.clear();
    }

~SparseGridsSimulation ( )

inline

Destructor

Definition at line 65 of file sparse_grids_simulation.h.

    {
      int err;
      /* clean up full grid object */
      err = CleanupBarebones(m_sim_fg);
      delete m_sim_fg;
      /* clean up and delete sparse grids objects */
      err = Cleanup((void*) m_sims_sg.data(), m_nsims_sg);
      if (err) {
        printf( "Error: CleanUp() returned with status %d on process %d.\n",
                err, m_rank );
      }
      m_sims_sg.clear();
      /* done */
      return;
    }

Member Function Documentation

int define ( int  a_rank, int  a_nproc  )

virtual

Define this object

Define the sparse grids simulation object - here, only the full grid simulation object SparseGridsSimulation::m_sim_fg is created.

This function also reads sparse grids inputs from the file sparse_grids.inp. Rank 0 reads in the inputs and broadcasts them to all the processors.

The format of sparse_grids.inp is as follows:

begin
    <keyword>   <value>
    <keyword>   <value>
    ...
    <keyword>   <value>
end

where the list of keywords and their type are:

Keyword name	Type	Variable	Default value
log2_imin	int	SparseGridsSimulation::m_imin	2
interp_order	int	SparseGridsSimulation::m_interp_order	6
write_sg_solution	char[]	SparseGridsSimulation::m_write_sg_solutions	"no" (0)
write_sg_errors	char[]	SparseGridsSimulation::m_print_sg_errors	"no" (0)

Parameters

a_rank	MPI rank of this process
a_nproc	Total number of MPI ranks

Implements Simulation.

Definition at line 33 of file SparseGridsDefine.cpp.

{
  m_rank = a_rank;
  m_nproc = a_nproc;

  if (m_sim_fg != NULL) {
    fprintf(stderr, "Error in SparseGridsSimulation::define() -\n");
    fprintf(stderr, "  m_sim_fg not NULL!\n");
    exit(1);
  }

  /* default values */
  m_imin = 2;
  m_write_sg_solutions = 0;
  m_print_sg_errors = 0;
  m_interp_order = 6;

  if (!m_rank) {

    FILE *in;
    in = fopen(_SPARSEGRIDS_SIM_INP_FNAME_,"r");

    if (!in) {

      fprintf(stderr, "Error in SparseGridsSimulation::Define() -\n");
      fprintf(stderr, "  %s file not found.\n", _SPARSEGRIDS_SIM_INP_FNAME_);
      exit(1);

    } else {

      int ferr;
      char word[_MAX_STRING_SIZE_];
      ferr = fscanf(in,"%s", word); if (ferr != 1) return(1);

      if (std::string(word) == "begin") {

        while (std::string(word) != "end") {

          ferr = fscanf(in,"%s",word); if (ferr != 1) return(1);

          if (std::string(word) == "log2_imin") {

            ferr = fscanf(in,"%d",&m_imin); if (ferr != 1) return(1);

          } else if (std::string(word) == "interp_order") {

            ferr = fscanf(in,"%d",&m_interp_order); if (ferr != 1) return(1);

          } else if (std::string(word) == "write_sg_solutions") {

            char answer[_MAX_STRING_SIZE_];
            ferr = fscanf(in,"%s",answer); if (ferr != 1) return(1);
            m_write_sg_solutions = (strcmp(answer,"yes") ? 0 : 1);

          } else if (std::string(word) == "write_sg_errors") {

            char answer[_MAX_STRING_SIZE_];
            ferr = fscanf(in,"%s",answer); if (ferr != 1) return(1);
            m_print_sg_errors = (strcmp(answer,"yes") ? 0 : 1);

          } else if (std::string(word) != "end") {

            char useless[_MAX_STRING_SIZE_];
            ferr = fscanf(in,"%s",useless);
            printf("Warning: keyword %s in file \"%s\" with value %s not recognized or extraneous. Ignoring.\n",
                    _SPARSEGRIDS_SIM_INP_FNAME_, word, useless );

          }

          if (ferr != 1) return(1);
        }

      } else {
        fprintf( stderr, "Error: Illegal format in file \"%s\". Word read is: %s\n",
                 _SPARSEGRIDS_SIM_INP_FNAME_, word);
        return 1;
      }

      fclose(in);

    }

    /* print useful stuff to screen */
    printf("Sparse grids inputs:\n");
    printf("  log2 of minimum grid size:  %d\n", m_imin);
    printf("  interpolation order:  %d\n", m_interp_order);
    printf( "  write sparse grids solutions?  %s\n",
            ( m_write_sg_solutions == 1 ? "yes" : "no" ) );
  }

#ifndef serial
  MPI_Bcast(&m_imin,1,MPI_INT,0,MPI_COMM_WORLD);
  MPI_Bcast(&m_interp_order,1,MPI_INT,0,MPI_COMM_WORLD);
  MPI_Bcast(&m_write_sg_solutions,1,MPI_INT,0,MPI_COMM_WORLD);
  MPI_Bcast(&m_print_sg_errors,1,MPI_INT,0,MPI_COMM_WORLD);
#endif

  m_sim_fg = new SimulationObject;
  /* the following are deliberately set to junk values
   * to ensure they are never used */
  m_sim_fg->solver.my_idx = -1;
  m_sim_fg->solver.nsims = -1;
  /* these are okay */
  m_sim_fg->mpi.rank = m_rank;
  m_sim_fg->mpi.nproc = m_nproc;

  if (!m_rank) {
    printf("Allocated full grid simulation object(s).\n");
  }

  m_is_defined = true;
  return 0;
}

int ReadInputs ( )

inlinevirtual

Read solver inputs for the full grid simulation object

Implements Simulation.

Definition at line 86 of file sparse_grids_simulation.h.

    {
      int retval = ::ReadInputs(  (void*) m_sim_fg,
                                  1, 
                                  m_rank );
      return retval;
    }

int Initialize ( )

virtual

Initialize the simulation

Initialize all the stuff required for a sparse grids simulation. Before this function is called, the solver inputs for the full grid object has already been read, and therefore, the global dimensions and processor decompositions are available. This function does the following:

Implements Simulation.

Definition at line 15 of file SparseGridsInitialize.cpp.

{
  int ierr;

  /* find out the number of spatial dimensions */
  m_ndims = m_sim_fg->solver.ndims;

  /* make sure full grid simulation object is allocated */
  if (m_sim_fg == NULL) {
    fprintf(stderr, "Error in SparseGridsSimulation::Initialize()\n");
    fprintf(stderr, "  m_sim_fg is NULL on rank %d!\n", m_rank);
    return 1;
  }

  /* set full grid processor distribution based on grid size and 
   * number of MPI ranks */
  GridDimensions dim_fg(m_ndims);
  for (int d=0; d<m_ndims; d++) dim_fg[d] = m_sim_fg->solver.dim_global[d];
  ProcDistribution iproc_fg;
  ComputeProcessorDistribution(iproc_fg, dim_fg);
  for (int d=0; d<m_ndims; d++) m_sim_fg->mpi.iproc[d] = iproc_fg[d];
  MPIBroadcast_integer( m_sim_fg->mpi.iproc,
                        m_ndims,
                        0,
                        &(m_sim_fg->mpi.world)); CHECKERR(ierr);

  ::WriteInputs( (void*) m_sim_fg, 1, m_rank);
  if (!m_rank) {
    printf("Processor distribution for full grid object: ");
    for (int d=0; d<m_ndims; d++) printf(" %3d", iproc_fg[d]);
    printf("\n");
  }
 
  if (!m_rank) printf("Initializing sparse grids...\n");
  if (!m_rank) printf("  Number of spatial dimensions: %d\n", m_ndims);

  /* some sanity checks */
  ierr = SanityChecks();
  if (ierr) return ierr;

  /* compute the sizes of grids that go into the combination technique */
  if (!m_rank) {
    printf("  Computing sparse grids dimensions...\n");
  }
  ierr = ComputeSGDimsAndCoeffs();
  if (ierr) return ierr;
  m_nsims_sg = m_combination.size();
  if (!m_rank) {
    printf( "  Number of sparse grid domains in combination technique: %d\n",
            m_nsims_sg );
  }
  if (m_nsims_sg == 0) {
    fprintf(stderr, "Error in SparseGridsSimulation::Initialize()\n");
    fprintf(stderr, "  ComputeSGDimsAndCoeffs() returned empty vector!\n");
    return 1;
  }

  /* compute processor distributions for the sparse grids */
  if (!m_rank) {
    printf("  Computing processor decompositions...\n");
  }
  m_iprocs.resize(m_nsims_sg);
  for (int i = 0; i < m_nsims_sg; i++) {
    ierr = ComputeProcessorDistribution(  m_iprocs[i],
                                          m_combination[i]._dim_ );

    if (ierr) return ierr;
  }

  /* Print some stuff to screen */
  if (!m_rank) {
    printf("Sparse Grids: Combination technique grids sizes and coefficients are:-\n");
    for (int i = 0; i < m_nsims_sg; i++) {
      printf("  %3d: dim = (", i);
      for (int d=0; d<m_ndims; d++) printf(" %4d ", m_combination[i]._dim_[d]);
      printf("),  coeff = %+1.2e, ", m_combination[i]._coeff_);
      printf("iproc = (");
      for (int d=0; d<m_ndims; d++) printf(" %4d ", m_iprocs[i][d]);
      printf(")\n");
    }
  }

  /* allocate full grid data structures (only what is needed) */
  ierr = InitializeBarebones(m_sim_fg);
  if (ierr) return ierr;

  /* create sparse grids simulation objects */
  m_sims_sg.resize(m_nsims_sg);
  /* set the solver parameters for the sparse grids sim objects */
  for (int i = 0; i < m_nsims_sg; i++) {
#ifndef serial
    MPI_Comm_dup(MPI_COMM_WORLD, &(m_sims_sg[i].mpi.world));
#endif
    ierr = SetSolverParameters( m_sims_sg[i],
                                m_combination[i]._dim_,
                                m_iprocs[i],
                                *m_sim_fg,
                                i,
                                m_nsims_sg );
    if (ierr) return ierr;
  }

  /* allocate their data structures (everything) */
  ierr = ::Initialize( (void*) m_sims_sg.data(), m_nsims_sg);
  if (ierr) return ierr;

  /* compute and report total NDOFs for sparse grids */
  long ndof_sg = 0;
  for (int i = 0; i < m_nsims_sg; i++) {
    ndof_sg += m_sims_sg[i].solver.npoints_global;
  }
  long ndof_fg = m_sim_fg->solver.npoints_global;
  if (!m_rank) {
    printf("Total number of DOFs:-\n");
    printf("  using sparse grids: %d\n", (int)ndof_sg);
    printf("  using conventional grid: %d\n", (int)ndof_fg);
  }

  /* done */
  return 0;
}

int InitialSolution ( )

virtual

Read initial solution

Reads in the initial solution for the full grid, saves it in the full grid object. Then interpolates (coarsens) the grid coordinates to all the sparge grids. The solution is not yet interpolated because boundary information is not yet available. It is done in SparseGridsSimulation::InitializationWrapup().

Implements Simulation.

Definition at line 17 of file SparseGridsInitialSolution.cpp.

{
  int ierr;

  /* first, read in the full grid initial solution */
  int retval = ::InitialSolution( (void*) m_sim_fg, 1 );
  if (retval) return retval;

  /* now, interpolate to all the sparse grids */
  for (int n = 0; n < m_nsims_sg; n++) {

    if (!m_rank) {
      printf("Interpolating grid coordinates to sparse grids domain %d.\n", n);
    }
    interpolateGrid( &m_sims_sg[n], m_sim_fg );

    HyPar* solver = &(m_sims_sg[n].solver);
    MPIVariables* mpi = &(m_sims_sg[n].mpi);

    int nvars = solver->nvars;
    int ghosts = solver->ghosts;
    int *dim_local = solver->dim_local;

    /* calculate dxinv */
    {
      int offset = 0;
      for (int d = 0; d < m_ndims; d++) {
        for (int i = 0; i < dim_local[d]; i++) {
          solver->dxinv[i+offset+ghosts] 
            = 2.0 / (   solver->x[i+1+offset+ghosts]
                      - solver->x[i-1+offset+ghosts]  );
        }
        offset += (dim_local[d] + 2*ghosts);
      }
    }

    /* exchange MPI-boundary values of dxinv between processors */
    {
      int offset = 0;
      for (int d = 0; d < m_ndims; d++) {
        ierr = MPIExchangeBoundaries1D( mpi,
                                        &(solver->dxinv[offset]),
                                        dim_local[d],
                                        ghosts,
                                        d,
                                        m_ndims );
        if (ierr) return ierr;
        offset += (dim_local[d] + 2*ghosts);
      }
    }

    /* fill in ghost values of dxinv at physical boundaries by extrapolation */
    {
      int offset = 0;
      for (int d = 0; d < m_ndims; d++) {
        double *dxinv = &(solver->dxinv[offset]);
        int    *dim = dim_local;
        if (mpi->ip[d] == 0) {
          /* fill left boundary along this dimension */
          for (int i = 0; i < ghosts; i++) dxinv[i] = dxinv[ghosts];
        }
        if (mpi->ip[d] == mpi->iproc[d]-1) {
          /* fill right boundary along this dimension */
          for (int i = dim[d]+ghosts; i < dim[d]+2*ghosts; i++) {
            dxinv[i] = dxinv[dim[d]+ghosts-1];
          }
        }
        offset  += (dim[d] + 2*ghosts);
      }
    }

  }

  return 0; 
}

int InitializeBoundaries ( )

inlinevirtual

Initialize the boundary conditions of the simulation

Implements Simulation.

Definition at line 101 of file sparse_grids_simulation.h.

    {
      int retval = ::InitializeBoundaries(  (void*) m_sims_sg.data(),
                                            m_nsims_sg );
      if (m_nsims_sg > 0) {
        m_is_periodic.resize(m_ndims);
        for (int d=0; d<m_ndims; d++) {
          m_is_periodic[d] = ( m_sims_sg[0].solver.isPeriodic[d] == 1 ? true : false);
          m_sim_fg->solver.isPeriodic[d] = m_sims_sg[0].solver.isPeriodic[d];
        }
      }
      return retval;
    }

int InitializeImmersedBoundaries ( )

inlinevirtual

Initialize the immersed boundary conditions of the simulation

Implements Simulation.

Definition at line 116 of file sparse_grids_simulation.h.

    {
      int retval = ::InitializeImmersedBoundaries(  (void*) m_sims_sg.data(),
                                                    m_nsims_sg );
      return retval;
    }

int InitializePhysics ( )

inlinevirtual

Initialize the physics of the simulations

Implements Simulation.

Definition at line 124 of file sparse_grids_simulation.h.

    {
      int retval = ::InitializePhysics( (void*) m_sims_sg.data(),
                                        m_nsims_sg );
      return retval;
    }

int InitializePhysicsData ( )

inlinevirtual

Initialize the physics of the simulations

Implements Simulation.

Definition at line 132 of file sparse_grids_simulation.h.

    {
      for (int ns = 0; ns < m_sims_sg.size(); ns++) {
        int retval = ::InitializePhysicsData( (void*) &(m_sims_sg[ns]),
                                              -1, 
                                              m_nsims_sg,
                                              m_sim_fg->solver.dim_global);
        if (retval) {
          fprintf(stderr, "Error in SparseGridsSimulation::InitializePhysicsData()\n");
          fprintf(stderr, "  InitializePhysicsData returned with error code %d on rank %d.\n",
                  retval, m_sims_sg[ns].mpi.rank);
          return retval;
        }
      }
      return 0;
    }

int InitializeSolvers ( )

inlinevirtual

Initialize the numerical solvers of the simulations

Implements Simulation.

Definition at line 150 of file sparse_grids_simulation.h.

    {
      int retval = ::InitializeSolvers( (void*) m_sims_sg.data(),
                                        m_nsims_sg );
      InitializeSolversBarebones(m_sim_fg);

      /* some modifications to output filename roots */
      for (int i=0; i<m_nsims_sg; i++) {
        strcpy(m_sims_sg[i].solver.op_fname_root, "op_sg");
        strcpy(m_sims_sg[i].solver.aux_op_fname_root, "ts0_sg");
      }
      strcpy(m_sim_fg->solver.op_fname_root, "op_fg");
      strcpy(m_sim_fg->solver.aux_op_fname_root, "ts0_fg");
      return retval;
    }

int InitializationWrapup ( )

virtual

Wrap up initializations

Finish all the initializations: Interpolate the initial solution on the full grid to all the sparse grids, now that the boundary conditions are available.

Reimplemented from Simulation.

Definition at line 15 of file SparseGridsInitializationWrapup.cpp.

{
  int ierr;

  for (int n = 0; n < m_nsims_sg; n++) {

    if (!m_rank) {
      printf("Interpolating initial solution to sparse grids domain %d.\n", n);
    }
    interpolate( &m_sims_sg[n], m_sim_fg );

    HyPar* solver = &(m_sims_sg[n].solver);
    MPIVariables* mpi = &(m_sims_sg[n].mpi);

    int nvars = solver->nvars;
    int ghosts = solver->ghosts;
    int *dim_local = solver->dim_local;

    /* exchange MPI-boundary values of u between processors */
    MPIExchangeBoundariesnD(  m_ndims,
                              nvars,
                              dim_local,
                              ghosts,
                              mpi,
                              solver->u  );

    /* calculate volume integral of the initial solution */
    ierr = ::VolumeIntegral(solver->VolumeIntegralInitial,
                            solver->u,
                            solver,
                            mpi );
    if (ierr) return ierr;

    if (!mpi->rank) {
      printf("  Volume integral of the initial solution on sparse grids domain %d:\n", n);
      for (int d=0; d<nvars; d++) {
        printf("    %2d:  %1.16E\n",d,solver->VolumeIntegralInitial[d]);
      }
    }

    /* Set initial total boundary flux integral to zero */
    _ArraySetValue_( solver->TotalBoundaryIntegral, nvars, 0 );

  }

  return 0; 
}

int Solve ( )

virtual

Run the simulation using native time integrators

This function integrates the semi-discrete ODE (obtained from discretizing the PDE in space) using natively implemented time integration methods. It initializes the time integration object, iterates the simulation for the required number of time steps, and calculates the errors. After the specified number of iterations, it writes out some information to the screen and the solution to a file.

Implements Simulation.

Definition at line 17 of file SparseGridsSolve.cpp.

{
  int tic     = 0;

  /* write out iblank to file for visualization */
  if (m_write_sg_solutions == 1) {
    for (int ns = 0; ns < m_nsims_sg; ns++) {
      if (m_sims_sg[ns].solver.flag_ib) {
  
        char fname_root[_MAX_STRING_SIZE_] = "iblank_sg";
        if (m_nsims_sg > 1) {
          char index[_MAX_STRING_SIZE_];
          GetStringFromInteger(ns, index, (int)log10((m_nsims_sg)+1));
          strcat(fname_root, "_");
          strcat(fname_root, index);
        }
  
        WriteArray( m_sims_sg[ns].solver.ndims,
                    1,
                    m_sims_sg[ns].solver.dim_global,
                    m_sims_sg[ns].solver.dim_local,
                    m_sims_sg[ns].solver.ghosts,
                    m_sims_sg[ns].solver.x,
                    m_sims_sg[ns].solver.iblank,
                    &(m_sims_sg[ns].solver),
                    &(m_sims_sg[ns].mpi),
                    fname_root );
      }
    }
  }
  if (m_sim_fg->solver.flag_ib) {

    char fname_root[_MAX_STRING_SIZE_] = "iblank";
    WriteArray( m_sim_fg->solver.ndims,
                1,
                m_sim_fg->solver.dim_global,
                m_sim_fg->solver.dim_local,
                m_sim_fg->solver.ghosts,
                m_sim_fg->solver.x,
                m_sim_fg->solver.iblank,
                &(m_sim_fg->solver),
                &(m_sim_fg->mpi),
                fname_root );
  }

  /* Define and initialize the time-integration object */
  TimeIntegration TS;
  if (!m_rank) printf("Setting up time integration.\n");
  TimeInitialize((void*)m_sims_sg.data(), m_nsims_sg, m_rank, m_nproc, &TS);

  if (!m_rank) {
    printf( "Solving in time (from %d to %d iterations)\n",
            TS.restart_iter,TS.n_iter);
  }

  for (TS.iter = TS.restart_iter; TS.iter < TS.n_iter; TS.iter++) {

    /* Write initial solution to file if this is the first iteration */
    if (!TS.iter) OutputSolution(TS.waqt);

    /* Call pre-step function */
    TimePreStep  (&TS);

    /* Step in time */
    TimeStep     (&TS);

    /* Call post-step function */
    TimePostStep (&TS);

    /* Print information to screen */
    TimePrintStep(&TS);
    tic++;

    /* Write intermediate solution to file */
    if ((TS.iter+1)%m_sims_sg[0].solver.file_op_iter == 0) { 
      OutputSolution(TS.waqt);
      tic = 0; 
    }

  }

  /* write a final solution file, if last iteration did not write one */
  if (tic || (!TS.n_iter)) OutputSolution(TS.waqt);

  if (!m_rank) {
    printf("Completed time integration (Final time: %f).\n",TS.waqt);
    if (m_nsims_sg > 1) printf("\n");
  }

  /* calculate error if exact solution has been provided */
  CalculateError();

  TimeCleanup(&TS);

  return(0);
}

void WriteErrors ( double  solver_runtime, double  main_runtime  )

virtual

Write simulation errors and wall times to file

Writes out the errors and other data for the sparse grids simulation.

Parameters

solver_runtime	Measured runtime of solver
main_runtime	Measured total runtime

Implements Simulation.

Definition at line 14 of file SparseGridsWriteErrors.cpp.

{
  if (!m_rank) {

    /* Write sparse grids stuff, if asked for */
    if (m_print_sg_errors == 1) {
      for (int n = 0; n < m_nsims_sg; n++) {
  
        char  err_fname[_MAX_STRING_SIZE_],
              cons_fname[_MAX_STRING_SIZE_],
              fc_fname[_MAX_STRING_SIZE_];
  
        strcpy(err_fname, "errors");
        strcpy(cons_fname,"conservation");
        strcpy(fc_fname,  "function_counts");
  
        if (m_nsims_sg > 1) {
  
          strcat(err_fname,"_");
          strcat(cons_fname,"_");
          strcat(fc_fname,"_");
  
          char index[_MAX_STRING_SIZE_];
          GetStringFromInteger(n, index, (int)log10(m_nsims_sg)+1);
  
          strcat(err_fname,index);
          strcat(cons_fname,index);
          strcat(fc_fname,index);
        }
  
        strcat(err_fname,".dat");
        strcat(cons_fname,".dat");
        strcat(fc_fname,".dat");
  
        FILE *out; 
        /* write out solution errors and wall times to file */
        int d;
        out = fopen(err_fname,"w");
        for (d=0; d<m_sims_sg[n].solver.ndims; d++) fprintf(out,"%4d ",m_sims_sg[n].solver.dim_global[d]);
        for (d=0; d<m_sims_sg[n].solver.ndims; d++) fprintf(out,"%4d ",m_sims_sg[n].mpi.iproc[d]);
        fprintf(out,"%1.16E  ",m_sims_sg[n].solver.dt);
        fprintf(out,"%1.16E %1.16E %1.16E   ",m_sims_sg[n].solver.error[0],m_sims_sg[n].solver.error[1],m_sims_sg[n].solver.error[2]);
        fprintf(out,"%1.16E %1.16E\n",solver_runtime,main_runtime);
        fclose(out);
        /* write out conservation errors to file */
        out = fopen(cons_fname,"w");
        for (d=0; d<m_sims_sg[n].solver.ndims; d++) fprintf(out,"%4d ",m_sims_sg[n].solver.dim_global[d]);
        for (d=0; d<m_sims_sg[n].solver.ndims; d++) fprintf(out,"%4d ",m_sims_sg[n].mpi.iproc[d]);
        fprintf(out,"%1.16E  ",m_sims_sg[n].solver.dt);
        for (d=0; d<m_sims_sg[n].solver.nvars; d++) fprintf(out,"%1.16E ",m_sims_sg[n].solver.ConservationError[d]);
        fprintf(out,"\n");
        fclose(out);
        /* write out function call counts to file */
        out = fopen(fc_fname,"w");
        fprintf(out,"%d\n",m_sims_sg[n].solver.n_iter);
        fprintf(out,"%d\n",m_sims_sg[n].solver.count_hyp);
        fprintf(out,"%d\n",m_sims_sg[n].solver.count_par);
        fprintf(out,"%d\n",m_sims_sg[n].solver.count_sou);
  #ifdef with_petsc
        fprintf(out,"%d\n",m_sims_sg[n].solver.count_RHSFunction);
        fprintf(out,"%d\n",m_sims_sg[n].solver.count_IFunction);
        fprintf(out,"%d\n",m_sims_sg[n].solver.count_IJacobian);
        fprintf(out,"%d\n",m_sims_sg[n].solver.count_IJacFunction);
  #endif
        fclose(out);
  
        /* print solution errors, conservation errors, and wall times to screen */
        if (m_sims_sg[n].solver.error[0] >= 0) {
          printf("Computed errors for sparse grids domain %d:\n", n);
          printf("  L1   Error: %1.16E\n",m_sims_sg[n].solver.error[0]);
          printf("  L2   Error: %1.16E\n",m_sims_sg[n].solver.error[1]);
          printf("  Linf Error: %1.16E\n",m_sims_sg[n].solver.error[2]);
        }
        if (!strcmp(m_sims_sg[n].solver.ConservationCheck,"yes")) {
          printf("Conservation Errors:\n");
          for (d=0; d<m_sims_sg[n].solver.nvars; d++) printf("\t%1.16E\n",m_sims_sg[n].solver.ConservationError[d]);
          printf("\n");
        }
  
      }
    }

    /* First write stuff for the full grid solution */
    {
      char  err_fname[_MAX_STRING_SIZE_];
      strcpy(err_fname,"errors_fg");
      strcat(err_fname,".dat");

      FILE *out; 
      /* write out solution errors and wall times to file */
      int d;
      out = fopen(err_fname,"w");
      for (d=0; d<m_sim_fg->solver.ndims; d++) fprintf(out,"%4d ",m_sim_fg->solver.dim_global[d]);
      for (d=0; d<m_sim_fg->solver.ndims; d++) fprintf(out,"%4d ",m_sim_fg->mpi.iproc[d]);
      fprintf(out,"%1.16E  ",m_sim_fg->solver.dt);
      fprintf(out,"%1.16E %1.16E %1.16E   ",m_sim_fg->solver.error[0],m_sim_fg->solver.error[1],m_sim_fg->solver.error[2]);
      fprintf(out,"%1.16E %1.16E\n",solver_runtime,main_runtime);
      fclose(out);

      /* print solution errors, conservation errors, and wall times to screen */
      if (m_sim_fg->solver.error[0] >= 0) {
        printf("Computed errors for full grid solution:\n");
        printf("  L1   Error: %1.16E\n",m_sim_fg->solver.error[0]);
        printf("  L2   Error: %1.16E\n",m_sim_fg->solver.error[1]);
        printf("  Linf Error: %1.16E\n",m_sim_fg->solver.error[2]);
      }
    }

    /* report wall times */
    printf("Solver runtime (in seconds): %1.16E\n",solver_runtime);
    printf("Total  runtime (in seconds): %1.16E\n",main_runtime);

  }

  return;
}

bool isDefined ( ) const

inlinevirtual

Return whether object is defined or not

Implements Simulation.

Definition at line 176 of file sparse_grids_simulation.h.

{ return m_is_defined; }

int mpiCommDup ( )

inlinevirtual

Create duplicate MPI communicators

Implements Simulation.

Definition at line 180 of file sparse_grids_simulation.h.

    {
      if (!m_sim_fg) {
        fprintf(stderr, "Error in SparseGridsSimulation::mpiCommDup()\n");
        fprintf(stderr, "  m_sim_fg is not allocated on rank %d!\n", m_rank);
        return 1;
      }
      MPI_Comm_dup(MPI_COMM_WORLD, &(m_sim_fg->mpi.world));
      return 0;
    }

void usePetscTS ( PetscBool  a_flag )

inlinevirtual

Set flag whether to use PETSc time integration

Implements Simulation.

Definition at line 194 of file sparse_grids_simulation.h.

    {
      m_sim_fg->solver.use_petscTS = a_flag;
    }

int SolvePETSc ( )

inlinevirtual

Run the simulation using PETSc time integrators

Implements Simulation.

Definition at line 200 of file sparse_grids_simulation.h.

    {
      int retval = ::SolvePETSc(  (void*) m_sims_sg.data(),
                                  m_nsims_sg,
                                  m_rank,
                                  m_nproc );
      return retval;
    }

int SanityChecks ( )

protected

Some basic sanity checks

Do some basic sanity checks to make sure that a sparse grids simulation can be carried out for this setup.

For the full grid:

• Check if grid sizes are same along all dimensions.
• Check if grid size is a power of 2.
• Check if number of MPI ranks are same along all dimensions.
• Check if number of MPI ranks is a power of 2.

Definition at line 18 of file SparseGridsSanityChecks.cpp.

{
  int *dim_global_fg = m_sim_fg->solver.dim_global;
  int *iproc_fg = m_sim_fg->mpi.iproc;

  /* check if grid sizes are same along all dimensions */
  {
    bool flag = true;
    for (int d=0; d<m_ndims; d++) {
      flag = flag && (dim_global_fg[d] == dim_global_fg[0]);
    }
    if (!flag) {
      fprintf(stderr, "Error in SparseGridsSimulation::SanityChecks()\n");
      fprintf(stderr, "  full grid dimensions are not equal in all dimensions.\n");
      return 1;
    }
  }

  /* check if grid size is a power of 2 */
  {
    bool flag = isPowerOfTwo(dim_global_fg[0]);
    if (!flag) {
      fprintf(stderr, "Error in SparseGridsSimulation::SanityChecks()\n");
      fprintf(stderr, "  full grid dimensions are not a power of 2.\n");
      return 1;
    }
  }
  
  return 0;
}

int ComputeSGDimsAndCoeffs ( )

protected

Compute the number, coefficients, and dimensions of all the sparse grids that go into the combination technique.

Definition at line 24 of file SparseGridsMiscFuncs.cpp.

{
  m_n_fg = log2(m_sim_fg->solver.dim_global[0]);
  int d = m_ndims;

  double c1 = 1;
  m_combination.clear();
  for (int s = 1; s <= d; s++) {

    double coeff = c1 * cfunc(d-1,s-1);

    std::vector<GridDimensions> dims(0);
    GetCTGridSizes((m_n_fg+(m_ndims-1)*(m_imin-1))+d-s, dims);
    if (dims.size() == 0) {
      fprintf(stderr, "Error in SparseGridsSimulation::ComputeSGDimsAndCoeffs()\n");
      fprintf(stderr, "  GetCTGridSize() returned empty vector!\n");
      m_combination.clear();
      return 1;
    }
    for (int i=0; i<dims.size(); i++) {
      m_combination.push_back(std::pair<double,GridDimensions>(coeff,dims[i]));
    }

    c1 *= -1;

  }

  return 0;
}

int ComputeProcessorDistribution ( ProcDistribution &  a_iprocs, const GridDimensions &  a_dim  )

protected

Compute the load-balanced MPI ranks distributions

Compute the load-balanced processor distribution for a given grid size and the total number of MPI ranks available

Parameters

a_iprocs	Processor distribution to compute
a_dim	Grid dimensions

Definition at line 85 of file SparseGridsMiscFuncs.cpp.

{
  a_iprocs.resize(m_ndims);
  /* get the normal vector for the grid dimensions */
  std::vector<double> dvec;
  StdVecOps::createNormalVector(dvec, a_dim);

  /* calculate the maximum number of MPI ranks in each dimension */
  int max_procs[m_ndims], min_procs[m_ndims];
  for (int i=0; i<m_ndims; i++) {
    max_procs[i] = a_dim[i]/_MIN_GRID_PTS_PER_PROC_;
    min_procs[i] = 1;
  }
  int max_nproc; _ArrayProduct1D_(max_procs, m_ndims, max_nproc);
  if (max_nproc < m_nproc) {
    fprintf(stderr, "Error in SparseGridsSimulation::ComputeProcessorDistribution() - rank %d\n", m_rank);
    fprintf(stderr, "  Number of MPI ranks greater than the maximum number of MPI ranks that can be used.\n");
    fprintf(stderr, "  Please re-run with %d MPI ranks.\n", max_nproc);
    return 1;
  }

  /* find all the processor distributions that are okay, i.e., their product
   * is the total number of MPI ranks #SparseGridsSimulation::m_nproc */
  std::vector<ProcDistribution> iproc_candidates(0);
  int iproc[m_ndims], ubound[m_ndims], done = 0;
  for (int d=0; d<m_ndims; d++) ubound[d] = max_procs[d]+1;
  _ArraySetValue_(iproc, m_ndims, 1);
  while (!done) {
    int prod; _ArrayProduct1D_(iproc, m_ndims, prod);
    if (prod == m_nproc) {
      ProcDistribution iproc_vec(m_ndims);
      for (int d = 0; d < m_ndims; d++) iproc_vec[d] = iproc[d];
      iproc_candidates.push_back(iproc_vec);
    }
    _ArrayIncrementIndexWithLBound_(m_ndims,ubound,min_procs,iproc,done);
  }
  if (iproc_candidates.size() == 0) {
    fprintf(stderr, "Error in SparseGridsSimulation::ComputeProcessorDistribution() - rank %d\n", m_rank);
    fprintf(stderr, "  Unable to fine candidate iprocs!\n");
    return 1;
  }

  /* find the candidate that is closest to the normalized dim vector */
  double min_norm = DBL_MAX;
  for (int i = 0; i<iproc_candidates.size(); i++) {
    std::vector<double> pvec;
    StdVecOps::createNormalVector(pvec, iproc_candidates[i]);
    double norm = StdVecOps::compute2Norm(pvec, dvec);
    if (norm < min_norm) {
      min_norm = norm;
      a_iprocs = iproc_candidates[i];
    }
  }

  return 0;
}

void GetCTGridSizes ( const int  a_N, std::vector< GridDimensions > &  a_dims  )

protected

Compute all grids such that the base2 logs of the size in each dimension adds up to the input argument

Parameters

a_N	Desired sum of the base2 logs of grid sizes
a_dims	Vector of grid sizes

Definition at line 56 of file SparseGridsMiscFuncs.cpp.

{
  a_dims.clear();

  int index[m_ndims], ubounds[m_ndims], lbounds[m_ndims];
  _ArraySetValue_(index,  m_ndims, m_imin);
  _ArraySetValue_(ubounds, m_ndims, a_N);
  _ArraySetValue_(lbounds, m_ndims, m_imin);

  int done = 0;
  while (!done) {
    int sum; _ArraySum1D_(index, sum, m_ndims);
    if (sum == a_N) {
      std::vector<int> tmp(m_ndims);
      for (int d=0; d<m_ndims; d++) {
        raiseto_int( tmp[d], 2,index[d] );
      }
      a_dims.push_back(tmp);
    }
    _ArrayIncrementIndexWithLBound_(m_ndims,ubounds,lbounds,index,done);
  }

  return;
}

int InitializeBarebones ( SimulationObject *  simobj )

protected

Initialize a "barebones" simulation object

A barebones initialization function for a simulation object. It will allocate data holders for only the stuff needed for storing a solution; this object cannot be used for an actual simulation. This function does the following:

• initializes the values for MPI variables
• allocates memory for arrays to store full grid solution

Parameters

simobj

simulation objects of type SimulationObject

Definition at line 194 of file SparseGridsMiscFuncs.cpp.

{
  /* Set this to "true" since this is a barebones initialization */
  simobj->is_barebones = 1;

  HyPar* solver = &(simobj->solver);
  MPIVariables* mpi = &(simobj->mpi);

  /* allocations */
  mpi->ip             = (int*) calloc (m_ndims,sizeof(int));
  mpi->is             = (int*) calloc (m_ndims,sizeof(int));
  mpi->ie             = (int*) calloc (m_ndims,sizeof(int));
  mpi->bcperiodic     = (int*) calloc (m_ndims,sizeof(int));
  solver->dim_local   = (int*) calloc (m_ndims,sizeof(int));
  solver->isPeriodic  = (int*) calloc (m_ndims,sizeof(int));

#ifndef serial
  /* Domain partitioning */
  int total_proc = 1;
  for (int i=0; i<m_ndims; i++) total_proc *= mpi->iproc[i];

  /* calculate ndims-D rank of each process (ip[]) from rank in MPI_COMM_WORLD */
  IERR MPIRanknD( m_ndims,
                  m_rank,
                  mpi->iproc,
                  mpi->ip ); CHECKERR(ierr);

  /* calculate local domain sizes along each dimension */
  for (int i=0; i<m_ndims; i++) {
    solver->dim_local[i] = MPIPartition1D( solver->dim_global[i],
                                           mpi->iproc[i],
                                           mpi->ip[i] );
  }

  /* calculate local domain limits in terms of global domain */
  IERR MPILocalDomainLimits(  m_ndims,
                              m_rank,
                              &(simobj->mpi),
                              solver->dim_global,
                              mpi->is,
                              mpi->ie  );
  CHECKERR(ierr);

  /* create sub-communicators for parallel computations 
   * along grid lines in each dimension */
   IERR MPICreateCommunicators(m_ndims,&(simobj->mpi)); CHECKERR(ierr);

  /* initialize periodic BC flags to zero */
  for (int i=0; i<solver->ndims; i++) mpi->bcperiodic[i] = 0;
  
  /* create communication groups */
  IERR MPICreateIOGroups(&(simobj->mpi)); CHECKERR(ierr);

#else

  for (int i=0; i<m_ndims; i++) {
    mpi->ip[i]            = 0;
    solver->dim_local[i]  = solver->dim_global[i];
    mpi->iproc[i]         = 1;
    mpi->is[i]            = 0;
    mpi->ie[i]            = solver->dim_local[i];
    mpi->bcperiodic[i]    = 0;
  }

#endif

  solver->npoints_global 
    = solver->npoints_local 
    = solver->npoints_local_wghosts
    = 1;
  for (int i=0; i<m_ndims; i++) {
    solver->npoints_global *= solver->dim_global[i];
    solver->npoints_local *= solver->dim_local [i];
    solver->npoints_local_wghosts *= (solver->dim_local[i]+2*solver->ghosts);
  }

  /* Allocations */
  if (!m_rank) printf("Allocating data arrays for full grid.\n");
  solver->index = (int*) calloc (m_ndims,sizeof(int));
  solver->stride_with_ghosts = (int*) calloc (solver->ndims,sizeof(int));
  solver->stride_without_ghosts = (int*) calloc (solver->ndims,sizeof(int));
  int accu1 = 1, accu2 = 1;
  for (int i=0; i<solver->ndims; i++) {
    solver->stride_with_ghosts[i]    = accu1;
    solver->stride_without_ghosts[i] = accu2;
    accu1 *= (solver->dim_local[i]+2*solver->ghosts);
    accu2 *=  solver->dim_local[i];
  }


  int size;
  /* state variables */
  size = 1;
  for (int i=0; i<m_ndims; i++) size *= (solver->dim_local[i]+2*solver->ghosts);
  solver->u = (double*) calloc (solver->nvars*size,sizeof(double));

  /* grid */
  size = 0;
  for (int i=0; i<m_ndims; i++) size += (solver->dim_local[i]+2*solver->ghosts);
  solver->x     = (double*) calloc (size,sizeof(double));
  solver->dxinv = (double*) calloc (size,sizeof(double));
    
  /* allocate MPI send/receive buffer arrays */
  int bufdim[solver->ndims], maxbuf = 0;
  for (int d = 0; d < solver->ndims; d++) {
    bufdim[d] = 1;
    for (int i = 0; i < solver->ndims; i++) {
      if (i == d) bufdim[d] *= solver->ghosts;
      else        bufdim[d] *= solver->dim_local[i];
    }
    if (bufdim[d] > maxbuf) maxbuf = bufdim[d];
  }
  maxbuf *= solver->nvars;
  mpi->maxbuf  = maxbuf;
  mpi->sendbuf = (double*) calloc (2*solver->ndims*maxbuf,sizeof(double));
  mpi->recvbuf = (double*) calloc (2*solver->ndims*maxbuf,sizeof(double));

  solver->VolumeIntegral        = (double*) calloc (solver->nvars,sizeof(double));
  solver->VolumeIntegralInitial = (double*) calloc (solver->nvars,sizeof(double));
  solver->TotalBoundaryIntegral = (double*) calloc (solver->nvars,sizeof(double));
  solver->ConservationError     = (double*) calloc (solver->nvars,sizeof(double));

  for (int i=0; i<solver->nvars; i++) solver->ConservationError[i] = -1;

  return(0);
}

int InitializeSolversBarebones ( SimulationObject *  sim )

protected

Initialize some function pointers for a "barebones" simulation object

This function initializes all some function pointers needed by a barebones simulation object.

Parameters

sim	simulation objects of type SimulationObject

Definition at line 324 of file SparseGridsMiscFuncs.cpp.

{
  if (sim->is_barebones != 1) {
    fprintf(stderr, "Error in SparseGridsSimulation::InitializeSolversBarebones() - \n");
    fprintf(stderr, "  simulation object is not barebones type.\n");
  }

  int ns;

  HyPar        *solver   = &(sim->solver);
  MPIVariables *mpi      = &(sim->mpi);

  solver->ParabolicFunction         = NULL;
  solver->SecondDerivativePar       = NULL;
  solver->FirstDerivativePar        = NULL;
  solver->InterpolateInterfacesPar  = NULL;

  solver->interp                = NULL;
  solver->compact               = NULL;
  solver->lusolver              = NULL;
  solver->SetInterpLimiterVar   = NULL;
  solver->flag_nonlinearinterp  = 0;
  solver->time_integrator = NULL;
  solver->msti = NULL;
  
  /* Solution output function */
  solver->WriteOutput    = NULL; /* default - no output */
  solver->filename_index = NULL;
  if (!strcmp(solver->output_mode,"serial")) {
    solver->index_length = 5;
    solver->filename_index = (char*) calloc (solver->index_length+1,sizeof(char));
    int i; for (i=0; i<solver->index_length; i++) solver->filename_index[i] = '0';
    solver->filename_index[solver->index_length] = (char) 0;
    if (!strcmp(solver->op_file_format,"text")) {
      solver->WriteOutput = WriteText;
      strcpy(solver->solnfilename_extn,".dat");
    } else if (!strcmp(solver->op_file_format,"tecplot2d")) {
      solver->WriteOutput = WriteTecplot2D;
      strcpy(solver->solnfilename_extn,".dat");
    } else if (!strcmp(solver->op_file_format,"tecplot3d")) {
      solver->WriteOutput = WriteTecplot3D;
      strcpy(solver->solnfilename_extn,".dat");
    } else if ((!strcmp(solver->op_file_format,"binary")) || (!strcmp(solver->op_file_format,"bin"))) {
      solver->WriteOutput = WriteBinary;
      strcpy(solver->solnfilename_extn,".bin");
    } else if (!strcmp(solver->op_file_format,"none")) {
      solver->WriteOutput = NULL;
    } else {
      fprintf(stderr,"Error (domain %d): %s is not a supported file format.\n",
              ns, solver->op_file_format);
      return(1);
    }
    if ((!strcmp(solver->op_overwrite,"no")) && solver->restart_iter) {
      /* if it's a restart run, fast-forward the filename */
      int t;
      for (t=0; t<solver->restart_iter; t++) 
        if ((t+1)%solver->file_op_iter == 0) IncrementFilenameIndex(solver->filename_index,solver->index_length);
    }
  } else if (!strcmp(solver->output_mode,"parallel")) {
    if (!strcmp(solver->op_file_format,"none")) solver->WriteOutput = NULL;
    else {
      /* only binary file writing supported in parallel mode */
      /* use post-processing scripts to convert              */
      solver->WriteOutput = WriteBinary;
      strcpy(solver->solnfilename_extn,".bin");
    }
  } else {

    fprintf(stderr,"Error (domain %d): %s is not a supported output mode.\n",
            ns, solver->output_mode);
    fprintf(stderr,"Should be \"serial\" or \"parallel\".    \n");
    return(1);

  }

  /* Solution plotting function */
  strcpy(solver->plotfilename_extn,".png");
#ifdef with_python
  solver->py_plt_func = NULL;
  solver->py_plt_func_args = NULL;
  {
    char python_plotting_fname[_MAX_STRING_SIZE_] = "plotSolution";
    PyObject* py_plot_name = PyUnicode_DecodeFSDefault(python_plotting_fname);
    PyObject* py_plot_module = PyImport_Import(py_plot_name);
    Py_DECREF(py_plot_name);
    if (py_plot_module) {
      solver->py_plt_func = PyObject_GetAttrString(py_plot_module, "plotSolution");
      if (!solver->py_plt_func) {
        if (!mpi->rank) {
          printf("Unable to load plotSolution function from Python module.\n");
        }
      } else {
        if (!mpi->rank) {
          printf("Loaded Python module for plotting.\n");
          printf("Loaded plotSolution function from Python module.\n");
        }
      }
    } else {
      if (!mpi->rank) {
        printf("Unable to load Python module for plotting.\n");
      }
    }
  }
#endif

  return 0;
}

int CleanupBarebones ( SimulationObject *  sim )

protected

Cleanup (deallocate) a "barebones" simulation object

Cleanup (deallocation) function for a barebones simulation object that has been allocated using SparseGridsSimulation::InitializeBarebones().

Parameters

sim	simulation objects of type SimulationObject

Definition at line 146 of file SparseGridsMiscFuncs.cpp.

{
  if (sim->is_barebones == 0) {
    fprintf(stderr, "Error in SparseGridsSimulation::CleanupBarebones()\n");
    fprintf(stderr, "  Simulation object is not a barebones one.\n");
    return 1;
  }

  HyPar* solver = &(sim->solver);
  MPIVariables* mpi = &(sim->mpi);

  /* Free the communicators created */
  IERR MPIFreeCommunicators(solver->ndims,mpi); CHECKERR(ierr);

  /* These variables are allocated in Initialize.c */
  free(solver->dim_global);
  free(solver->dim_global_ex);
  free(solver->dim_local);
  free(solver->index);
  free(solver->isPeriodic);
  free(solver->u);
  free(solver->x);
  free(solver->dxinv);
  free(mpi->iproc);
  free(mpi->ip);
  free(mpi->is);
  free(mpi->ie);
  free(mpi->bcperiodic);
  free(mpi->sendbuf);
  free(mpi->recvbuf);
  free(solver->VolumeIntegral);
  free(solver->VolumeIntegralInitial);
  free(solver->TotalBoundaryIntegral);
  free(solver->ConservationError);
  free(solver->stride_with_ghosts);
  free(solver->stride_without_ghosts);
  if (solver->filename_index) free(solver->filename_index);

  return(0);
}

int SetSolverParameters ( SimulationObject &  a_dst_sim, const GridDimensions &  a_dim_global, const ProcDistribution &  a_iproc, const SimulationObject &  a_src_sim, const int  a_idx, const int  a_nsims  )

protected

Set solver parameters for a simulation object

Set the solver parameters (stuff that is usually read in from solver.inp) for a simulation object, where the global grid sizes and processor distribution are specified, and all other parameters are same as a source simulation object

Parameters

a_dst_sim	Simulation object
a_dim_global	Specified global grid sizes
a_iproc	Specified processor distibution
a_src_sim	Source simulation object
a_idx	Index or serial number
a_nsims	Total number of simulations

Definition at line 11 of file SparseGridsSetSolverParameters.cpp.

{
  a_dst_sim.solver.my_idx = a_idx;
  a_dst_sim.solver.nsims = a_nsims;

  a_dst_sim.mpi.rank = m_rank;
  a_dst_sim.mpi.nproc = m_nproc;

  a_dst_sim.solver.ndims = a_src_sim.solver.ndims;
  a_dst_sim.solver.nvars = a_src_sim.solver.nvars;
  a_dst_sim.solver.ghosts = a_src_sim.solver.ghosts;

  a_dst_sim.solver.dim_global = (int*) calloc (a_dst_sim.solver.ndims,sizeof(int));
  a_dst_sim.mpi.iproc         = (int*) calloc (a_dst_sim.solver.ndims,sizeof(int));
  for (int d = 0; d < a_dst_sim.solver.ndims; d++) {
    a_dst_sim.solver.dim_global[d] = a_dim_global[d];
    a_dst_sim.mpi.iproc[d] = a_iproc[d];
  }

  a_dst_sim.solver.n_iter       = a_src_sim.solver.n_iter;
  a_dst_sim.solver.restart_iter = a_src_sim.solver.restart_iter;

  strcpy(a_dst_sim.solver.time_scheme, a_src_sim.solver.time_scheme);
  strcpy(a_dst_sim.solver.time_scheme_type, a_src_sim.solver.time_scheme_type);
  strcpy(a_dst_sim.solver.spatial_scheme_hyp, a_src_sim.solver.spatial_scheme_hyp);
  strcpy(a_dst_sim.solver.SplitHyperbolicFlux, a_src_sim.solver.SplitHyperbolicFlux);
  strcpy(a_dst_sim.solver.interp_type, a_src_sim.solver.interp_type);
  strcpy(a_dst_sim.solver.spatial_type_par, a_src_sim.solver.spatial_type_par);
  strcpy(a_dst_sim.solver.spatial_scheme_par, a_src_sim.solver.spatial_scheme_par);

  a_dst_sim.solver.dt = a_src_sim.solver.dt;

  strcpy(a_dst_sim.solver.ConservationCheck, a_src_sim.solver.ConservationCheck);

  a_dst_sim.solver.screen_op_iter = a_src_sim.solver.screen_op_iter;
  a_dst_sim.solver.file_op_iter = a_src_sim.solver.file_op_iter;

  strcpy(a_dst_sim.solver.op_file_format, a_src_sim.solver.op_file_format);
  strcpy(a_dst_sim.solver.ip_file_type, a_src_sim.solver.ip_file_type);

  strcpy(a_dst_sim.solver.input_mode, a_src_sim.solver.input_mode);
  strcpy(a_dst_sim.solver.output_mode, a_src_sim.solver.output_mode);
  a_dst_sim.mpi.N_IORanks = a_src_sim.mpi.N_IORanks;

  strcpy(a_dst_sim.solver.op_overwrite, a_src_sim.solver.op_overwrite);
  strcpy(a_dst_sim.solver.plot_solution, a_src_sim.solver.plot_solution);
  strcpy(a_dst_sim.solver.model, a_src_sim.solver.model);
  strcpy(a_dst_sim.solver.ib_filename, a_src_sim.solver.ib_filename);

  a_dst_sim.solver.flag_ib = a_src_sim.solver.flag_ib;

#ifdef with_petsc
  a_dst_sim.solver.use_petscTS = a_src_sim.solver.use_petscTS;
#endif

  return(0);
}

void OutputSolution ( double  a_time )

protected

Output solutions to file

Write solutions to file

Parameters

a_time

simulation time

Definition at line 11 of file SparseGridsOutputSolution.cpp.

{
  /* if asked for, write individual sparse grids solutions */
  if (m_write_sg_solutions == 1) {
    for (int ns = 0; ns < m_nsims_sg; ns++) {
      if (m_sims_sg[ns].solver.PhysicsOutput) {
        m_sims_sg[ns].solver.PhysicsOutput( &(m_sims_sg[ns].solver),
                                            &(m_sims_sg[ns].mpi),
                                            a_time );
      }
    }
    ::OutputSolution((void*)m_sims_sg.data(), m_nsims_sg, a_time);
  }

  /* Combine the sparse grids solutions to full grid */
  CombinationTechnique(m_sim_fg);

  /* Write the full grid solution */
  ::OutputSolution((void*)m_sim_fg, 1, a_time);

  return;
}

void CombinationTechnique ( SimulationObject * const  a_sim )

protected

Combination technique

Implements the combination technique where the solutions from all the sparse grids are combined to give a higher-resolution solution.

The sparse grids domains may have different processor layouts, so this combination is carried out on rank 0, and the solution is distributed.

Parameters

a_sim

target simulation object on which to combine

Definition at line 23 of file SparseGridsCombinationTechnique.cpp.

{
  double** u_sg = (double**) calloc (m_nsims_sg, sizeof(double*));
  std::vector<double> coeffs(m_nsims_sg, 0.0);
  for (int n=0; n<m_nsims_sg; n++) {
    u_sg[n] = m_sims_sg[n].solver.u;
    coeffs[n] = m_combination[n]._coeff_;
  }

  ::CombineSolutions( m_sims_sg.data(),
                      u_sg,
                      m_nsims_sg,
                      a_sim,
                      a_sim->solver.u,
                      coeffs.data() );

  /* done */
  return;
}

void interpolate ( SimulationObject * const  a_dst, const SimulationObject * const  a_src  )

protected

Interpolate data from one simulation object to another

Interpolate data from one grid to another. Note that along each dimension, the ratio of the number of grid points in the source grid and that in the destination grid must be an integer power of 2 (negative or positive).

Parameters

a_dst	Destination object
a_src	Source object

Definition at line 15 of file SparseGridsInterpolate.cpp.

{
  /* 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 interpolate ( const GridDimensions &  a_dim_dst, double ** const  a_u_dst, const SimulationObject * const  a_src  )

protected

Interpolate data from a simulation to a global C-array

Interpolate data from one grid to another of a desired resolution. Note that along each dimension, the ratio of the number of grid points in the source grid and that in the destination grid must be an integer power of 2 (negative or positive).

The incoming pointer must be NULL. After this function is executed, it will point to a chunk of memory with the interpolated solution. This is the global solution.

Since the source and destination grids may have different processor layouts, the interpolation is carried out on rank 0 by allocating global arrays on it. Therefore, this function is not scalable.

Parameters

a_dim_dst	grid dimensions to interpolate to
a_u_dst	pointer to array containing interpolated data
a_src	Source object

Definition at line 63 of file SparseGridsInterpolate.cpp.

{
  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;
}

void interpolateGrid ( SimulationObject * const  a_dst, const SimulationObject * const  a_src  )

protected

Interpolate data from one simulation object to another

Interpolate grid coordinates from one grid to another. Note that along each dimension, the ratio of the number of grid points in the source grid and that in the destination grid must be an integer power of 2 (negative or positive).

Since the source and destination grids may have different processor layouts, the interpolation is carried out on rank 0 by allocating global arrays on it. Therefore, this function is not scalable.

Parameters

a_dst	Destination object
a_src	Source object

Definition at line 20 of file SparseGridsInterpolateGrid.cpp.

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

  /* gather the source on rank 0 */
  double* xg_src = NULL;
  if (!m_rank) allocateGridArrays(dim_src, &xg_src);
  {
    int offset_global, offset_local;
    offset_global = offset_local = 0;
    for (int d=0; d<m_ndims; d++) {
      MPIGatherArray1D(  (void*) &(a_src->mpi),
                         (a_src->mpi.rank ? NULL : &xg_src[offset_global]),
                         &(a_src->solver.x[offset_local+a_src->solver.ghosts]),
                         a_src->mpi.is[d],
                         a_src->mpi.ie[d],
                         a_src->solver.dim_local[d],
                         0 );
      offset_global += a_src->solver.dim_global[d];
      offset_local  += a_src->solver.dim_local [d] + 2*a_src->solver.ghosts;
    }
  }

  /* now do the interpolation, dimension-by-dimension */
  double *xg_dst = NULL;
  if (!m_rank) {

    GridDimensions dim_to(m_ndims,0);
    GridDimensions dim_from(m_ndims,0);

    double* x_from;
    double* x_to;

    dim_to = dim_src;
    x_to = xg_src;
    x_from = NULL;

    for (int dir = 0; dir < m_ndims; dir++) {

      dim_from = dim_to;
      dim_to[dir] = dim_dst[dir];

      if (dim_from[dir] == dim_to[dir]) continue;

      double fac = (dim_to[dir] > dim_from[dir] ? 
                        (double)dim_to[dir]/(double)dim_from[dir] 
                      : (double)dim_from[dir]/(double)dim_to[dir] );
      if (!isPowerOfTwo((int)fac)) {
        fprintf(stderr,"Error in SparseGridsSimulation::interpolate() - \n");
        fprintf(stderr,"  refinement/coarsening factor not a power of 2!\n");
        exit(1);
      }

      if (x_from != NULL) free(x_from);
      x_from = x_to;

      allocateGridArrays(dim_to, &x_to);
      if (dim_to[dir] < dim_from[dir]) {
        coarsenGrid1D(dim_from, dim_to, x_from, x_to, dir);
      } else {
        refineGrid1D(dim_from, dim_to, x_from, x_to, dir);
      }
  
    }

    /* dim_to should be equal to dim_dst now */
    for (int d = 0; d < m_ndims; d++) {
      if (dim_to[d] != dim_dst[d]) {
        fprintf(stderr,"Error in SparseGridsSimulation::interpolate() - \n");
        fprintf(stderr,"  dim_to[%d] != dim_dst[%d]!\n", d, d);
        exit(1);
      }
    }

    if (x_from != NULL) free(x_from);
    xg_dst = x_to;

  }

  /* partition the destination */
  int offset_global = 0, offset_local = 0;
  for (int d=0; d<m_ndims; d++) {
    MPIPartitionArray1D(  (void*) &(a_dst->mpi),
                          (a_dst->mpi.rank ? NULL : &xg_dst[offset_global]),
                          &(a_dst->solver.x[offset_local+a_dst->solver.ghosts]),
                          a_dst->mpi.is[d],
                          a_dst->mpi.ie[d],
                          a_dst->solver.dim_local[d],
                          0);
    offset_global += a_dst->solver.dim_global[d];
    offset_local  += a_dst->solver.dim_local [d] + 2*a_dst->solver.ghosts;
  }

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

  /* exchange MPI-boundary values of x between processors */
  {
    int offset = 0;
    for (int d = 0; d < m_ndims; d++) {
      MPIExchangeBoundaries1D(  (void*) &(a_dst->mpi),
                                &(a_dst->solver.x[offset]),
                                a_dst->solver.dim_local[d],
                                a_dst->solver.ghosts,
                                d,
                                m_ndims);
      offset += (a_dst->solver.dim_local[d] + 2*a_dst->solver.ghosts);
    }
  }
  /* fill in ghost values of x at physical boundaries by extrapolation */
  {
    int offset = 0;
    for (int d = 0; d < m_ndims; d++) {
      double* X       = &(a_dst->solver.x[offset]);
      int*    dim     = a_dst->solver.dim_local;
      int     ghosts  = a_dst->solver.ghosts;
      if (a_dst->mpi.ip[d] == 0) {
        /* fill left boundary along this dimension */
        for (int i = 0; i < ghosts; i++) {
          int delta = ghosts - i;
          X[i] = X[ghosts] + ((double) delta) * (X[ghosts]-X[ghosts+1]);
        }
      }
      if (a_dst->mpi.ip[d] == a_dst->mpi.iproc[d]-1) {
        /* fill right boundary along this dimension */
        for (int i = dim[d]+ghosts; i < dim[d]+2*ghosts; i++) {
          int delta = i - (dim[d]+ghosts-1);
          X[i] =  X[dim[d]+ghosts-1] 
                  + ((double) delta) * (X[dim[d]+ghosts-1]-X[dim[d]+ghosts-2]);
        }
      }
      offset  += (dim[d] + 2*ghosts);
    }
  }
  return;
}

void coarsenGrid1D ( const GridDimensions &  a_dim_src, const GridDimensions &  a_dim_dst, const double * const  a_x_src, double * const  a_x_dst, int  a_dir  )

protected

Coarsen along a given dimension

Coarsen along a given dimension - the source and destination must have the same sizes along all other dimensions.

Note that the grid must not have any ghost points!

Parameters

a_dim_src	Grid size of source data
a_dim_dst	Grid size of destination data
a_x_src	Source grid
a_x_dst	Destination grid
a_dir	Dimension along which to coarsen

Definition at line 168 of file SparseGridsInterpolateGrid.cpp.

{
  for (int d = 0; d < m_ndims; d++) {
    if ((d != a_dir) && (a_dim_src[d] != a_dim_dst[d])) {
      fprintf(stderr, "Error in SparseGridsSimulation::coarsenGrid1D() -\n");
      fprintf(stderr, " a_dim_src[%d] != a_dim_dst[%d]\n", d, d);
      exit(1);
    }
  }

  int n_src = a_dim_src[a_dir];
  int n_dst = a_dim_dst[a_dir];
  if (n_dst > n_src) {
    fprintf(stderr, "Error in SparseGridsSimulation::coarsenGrid1D() -\n");
    fprintf(stderr, " destination grid is finer than source grid along a_dir!\n");
    exit(1);
  }

  double fac = ((double) n_src) / ((double) n_dst);
  int stride = (int) fac;
  if (std::abs(((double)stride)-fac) > _MACHINE_ZERO_) {
    fprintf(stderr, "Error in SparseGridsSimulation::coarsenGrid1D() -\n");
    fprintf(stderr, "  non-integer coarsening factor!\n");
    exit(1);
  }

  const double* x_src = a_x_src;
  double* x_dst = a_x_dst;

  for (int d = 0; d < m_ndims; d++) {

    if (d == a_dir) {

      for (int i = 0; i < a_dim_dst[d]; i++) {
        double avg = 0;
        for (int j=stride*i; j<stride*(i+1); j++) {
          if (j >= n_src) {
            fprintf(stderr, "Error in SparseGridsSimulation::coarsenGrid1D() -\n");
            fprintf(stderr, "  j >= n_src!\n");
            exit(1);
          }
          avg += x_src[j];
        }
        avg /= (double) stride;
        x_dst[i] = avg;
      }

    } else {

      _ArrayCopy1D_(x_src, x_dst, a_dim_dst[d]);
      
    }

    x_src += a_dim_src[d];
    x_dst += a_dim_dst[d];
  }
  
  return;
}

void refineGrid1D ( const GridDimensions &  a_dim_src, const GridDimensions &  a_dim_dst, const double * const  a_x_src, double * const  a_x_dst, int  a_dir  )

protected

Refine along a given dimension

Refine along a given dimension - the source and destination must have the same sizes along all other dimensions.

Note that the grid must not have any ghost points!

Parameters

a_dim_src	Grid size of source data
a_dim_dst	Grid size of destination data
a_x_src	Source grid
a_x_dst	Destination grid
a_dir	Dimension along which to refine

Definition at line 238 of file SparseGridsInterpolateGrid.cpp.

{
  for (int d = 0; d < m_ndims; d++) {
    if ((d != a_dir) && (a_dim_src[d] != a_dim_dst[d])) {
      fprintf(stderr, "Error in SparseGridsSimulation::coarsenGrid1D() -\n");
      fprintf(stderr, " a_dim_src[%d] != a_dim_dst[%d]\n", d, d);
      exit(1);
    }
  }

  int n_src = a_dim_src[a_dir];
  int n_dst = a_dim_dst[a_dir];
  if (n_dst < n_src) {
    fprintf(stderr, "Error in SparseGridsSimulation::refineGrid1D() -\n");
    fprintf(stderr, "  destination grid is coarser than source grid along a_dir!\n");
    exit(1);
  }

  double fac = ((double) n_dst) / ((double) n_src);
  int stride = (int) fac;
  if (std::abs(((double)stride)-fac) > _MACHINE_ZERO_) {
    fprintf(stderr, "Error in SparseGridsSimulation::refineGrid1D() -\n");
    fprintf(stderr, "  non-integer refinement factor!\n");
    exit(1);
  }

  int offset = 0;
  for (int d = 0; d < a_dir; d++) offset += a_dim_src[d];
  const double* x_src = a_x_src + offset;
  double* x_dst = a_x_dst + offset;

  fprintf(stderr, "Error in SparseGridsSimulation::refineGrid1D() -\n");
  fprintf(stderr, "  NOT YET IMPLEMENTED! Why do you need this?\n");
  exit(1);
  
  return;
}

void fillGhostCells ( const GridDimensions &  a_dim, const int  a_ngpt, double * const  a_u, const int  a_nvars  )

protected

Fill ghost cells for interpolation

Fill the ghost cells of a solution. Note that the solution array must be a global one (not one that is partitioned among MPI ranks). This is not a parallel operation (it will execute independently on multiple MPI ranks, if called by multiple processes).

For periodicity along any dimension, the ghost cells are filled appropriately from the other side of the domain. Otherwise, the interior data is extrapolated by a 4th order polynomial (assuming uniform grid spacing).

Parameters

a_dim	grid dimensions of solution
a_ngpt	number of ghost cells
a_u	solution array
a_nvars	number of vector components of the solution

Definition at line 18 of file SparseGridsFillGhostCells.cpp.

{
  if (m_is_periodic.size() != m_ndims) {
    fprintf(stderr, "Error in SparseGridsSimulation::fillGhostCells() -\n");
    fprintf(stderr, "m_is_periodic.size() != m_ndims \n");
    exit(1);
  }

  for (int d = 0; d < m_ndims; d++) {

    int bounds[m_ndims];
    _ArrayCopy1D_(a_dim, bounds, m_ndims);
    bounds[d] = a_ngpt;

    int index[m_ndims];
    _ArraySetValue_(index, m_ndims, 0);

    if (m_is_periodic[d]) {
      
      /* periodic case */

      int done = 0;
      while (!done) {

        {
          /* low end - face = 1 */

          int p_gpt = 0, 
              p_int = 0;

          int index_gpt[m_ndims];
          _ArrayCopy1D_(index, index_gpt, m_ndims);
          index_gpt[d] -= a_ngpt;
          _ArrayIndex1D_(m_ndims, a_dim, index_gpt, a_ngpt, p_gpt);

          int index_int[m_ndims];
          _ArrayCopy1D_(index, index_int, m_ndims);
          index_int[d] += (a_dim[d]-a_ngpt);
          _ArrayIndex1D_(m_ndims, a_dim, index_int, a_ngpt, p_int);

          _ArrayCopy1D_((a_u+a_nvars*p_int), (a_u+a_nvars*p_gpt), a_nvars);
        }

        {
          /* high end - face = -1 */

          int p_gpt = 0, 
              p_int = 0;

          int index_gpt[m_ndims];
          _ArrayCopy1D_(index, index_gpt, m_ndims);
          index_gpt[d] += a_dim[d];
          _ArrayIndex1D_(m_ndims, a_dim, index_gpt, a_ngpt, p_gpt);

          int index_int[m_ndims];
          _ArrayCopy1D_(index, index_int, m_ndims);
          _ArrayIndex1D_(m_ndims, a_dim, index_int, a_ngpt, p_int);

          _ArrayCopy1D_((a_u+a_nvars*p_int), (a_u+a_nvars*p_gpt), a_nvars);
        }

        _ArrayIncrementIndex_(m_ndims, bounds, index, done);

      }

    } else {

      /* not periodic - extrapolate */

      int done = 0;
      while (!done) {

        {
          /* low end - face = 1 */

          int p_gpt = 0, 
              p_int_0 = 0,
              p_int_1 = 0,
              p_int_2 = 0,
              p_int_3 = 0;

          int index_gpt[m_ndims];
          _ArrayCopy1D_(index, index_gpt, m_ndims);
          index_gpt[d] -= a_ngpt;
          _ArrayIndex1D_(m_ndims, a_dim, index_gpt, a_ngpt, p_gpt);

          int index_int[m_ndims];
          _ArrayCopy1D_(index, index_int, m_ndims);

          index_int[d] = 0;
          _ArrayIndex1D_(m_ndims, a_dim, index_int, a_ngpt, p_int_0);
          index_int[d]++;
          _ArrayIndex1D_(m_ndims, a_dim, index_int, a_ngpt, p_int_1);
          index_int[d]++;
          _ArrayIndex1D_(m_ndims, a_dim, index_int, a_ngpt, p_int_2);
          index_int[d]++;
          _ArrayIndex1D_(m_ndims, a_dim, index_int, a_ngpt, p_int_3);

          double alpha = - (double) (a_ngpt - index[d]);
          double c0 = -((-2.0 + alpha)*(-1.0 + alpha)*alpha)/6.0;
          double c1 = ((-2.0 + alpha)*(-1.0 + alpha)*(1.0 + alpha))/2.0;
          double c2 = (alpha*(2.0 + alpha - alpha*alpha))/2.0;
          double c3 = (alpha*(-1.0 + alpha*alpha))/6.0;

          for (int v = 0; v < a_nvars; v++) {
    
            a_u[p_gpt*a_nvars+v] =    c0 * a_u[p_int_0*a_nvars+v]
                                    + c1 * a_u[p_int_1*a_nvars+v]
                                    + c2 * a_u[p_int_2*a_nvars+v]
                                    + c3 * a_u[p_int_3*a_nvars+v];
    
          }

        }

        {
          /* high end - face = -1 */

          int p_gpt = 0, 
              p_int_0 = 0,
              p_int_1 = 0,
              p_int_2 = 0,
              p_int_3 = 0;

          int index_gpt[m_ndims];
          _ArrayCopy1D_(index, index_gpt, m_ndims);
          index_gpt[d] += a_dim[d];
          _ArrayIndex1D_(m_ndims, a_dim, index_gpt, a_ngpt, p_gpt);

          int index_int[m_ndims];
          _ArrayCopy1D_(index, index_int, m_ndims);

          index_int[d] = a_dim[d]-1;
          _ArrayIndex1D_(m_ndims, a_dim, index, a_ngpt, p_int_0);
          index_int[d]--;
          _ArrayIndex1D_(m_ndims, a_dim, index, a_ngpt, p_int_1);
          index_int[d]--;
          _ArrayIndex1D_(m_ndims, a_dim, index, a_ngpt, p_int_2);
          index_int[d]--;
          _ArrayIndex1D_(m_ndims, a_dim, index, a_ngpt, p_int_3);

          double alpha = - (double) (index[d]+1);
          double c0 = -((-2.0 + alpha)*(-1.0 + alpha)*alpha)/6.0;
          double c1 = ((-2.0 + alpha)*(-1.0 + alpha)*(1.0 + alpha))/2.0;
          double c2 = (alpha*(2.0 + alpha - alpha*alpha))/2.0;
          double c3 = (alpha*(-1.0 + alpha*alpha))/6.0;

          for (int v = 0; v < a_nvars; v++) {
    
            a_u[p_gpt*a_nvars+v] =    c0 * a_u[p_int_0*a_nvars+v]
                                    + c1 * a_u[p_int_1*a_nvars+v]
                                    + c2 * a_u[p_int_2*a_nvars+v]
                                    + c3 * a_u[p_int_3*a_nvars+v];
    
          }

        }

        _ArrayIncrementIndex_(m_ndims, bounds, index, done);

      }

    }

  }

  return;
}

void CalculateError ( )

protected

Calculate errors

Calculates the error in the solution if the exact solution is available. The exact solution should be provided in the file "exact.inp" in the same format as the initial solution.

Definition at line 17 of file SparseGridsCalculateError.cpp.

{
  int exact_flag;
  double *uex = NULL;

  /* read in full grid exact solution, if available */
  {
    HyPar* solver = &(m_sim_fg->solver);
    MPIVariables* mpi = &(m_sim_fg->mpi);
    long size = solver->nvars * solver->npoints_local_wghosts;
    uex = (double*) calloc (size, sizeof(double));
  
    char fname_root[_MAX_STRING_SIZE_] = "exact";
    IERR ExactSolution( solver,
                        mpi,
                        uex,
                        fname_root,
                        &exact_flag ); CHECKERR(ierr);
  }

  for (int n=0; n<m_nsims_sg; n++) {
//    TimeError(  &(m_sims_sg[n].solver),
//                &(m_sims_sg[n].mpi),
//                uex );
  }

  if (!exact_flag) {

    /* No exact solution available */
    for (int n=0; n<m_nsims_sg; n++) {
      m_sims_sg[n].solver.error[0] 
        = m_sims_sg[n].solver.error[1] 
        = m_sims_sg[n].solver.error[2] 
        = -1;
    }
    m_sim_fg->solver.error[0] 
      = m_sim_fg->solver.error[1] 
      = m_sim_fg->solver.error[2] 
      = -1;

  } else {

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

    double *uex2 = NULL;

    if (m_print_sg_errors == 1) {
      long size =   m_sim_fg->solver.nvars 
                  * m_sim_fg->solver.npoints_local_wghosts;
      uex2 = (double*) calloc(size, sizeof(double));
      _ArrayCopy1D_(uex, uex2, size);
    }

    /* calculate error for full grid */
    computeError( *m_sim_fg, uex );

    /* calculate error for sparse grids */
    if (m_print_sg_errors == 1) {
      for (int n = 0; n < m_nsims_sg; n++) {
  
        GridDimensions dim_fg(m_ndims,0);
        StdVecOps::copyFrom(dim_fg, m_sim_fg->solver.dim_global, m_ndims);
  
        GridDimensions dim_sg(m_ndims,0);
        StdVecOps::copyFrom(dim_sg, m_sims_sg[n].solver.dim_global, m_ndims);
  
        /* assemble the global exact solution on full grid */
        double *uex_global_fg = NULL;
        if (!m_rank) {
          allocateDataArrays( dim_fg, 
                              m_sim_fg->solver.nvars, 
                              &uex_global_fg, 
                              m_sim_fg->solver.ghosts);
        }
        MPIGatherArraynDwGhosts( m_ndims,
                                 (void*) &(m_sim_fg->mpi),
                                 uex_global_fg,
                                 uex2,
                                 m_sim_fg->solver.dim_global,
                                 m_sim_fg->solver.dim_local,
                                 m_sim_fg->solver.ghosts,
                                 m_sim_fg->solver.nvars );
  
        /* interpolate to sparse grid - 
         * this will delete the full grid array*/
        double *uex_global_sg = NULL;
        if (!m_rank) {
          int ierr = ::InterpolateGlobalnDVar(  dim_sg.data(),
                                                &uex_global_sg,
                                                dim_fg.data(),
                                                uex_global_fg,
                                                m_sims_sg[n].solver.nvars,
                                                m_sims_sg[n].solver.ghosts,
                                                m_ndims,
                                                periodic_arr.data() );
          if (ierr) {
            fprintf(stderr,"InterpolateGlobalnDVar() returned with error!\n");
            exit(1);
          }
        }
  
        /* allocate local exact solution on this sparse grid */
        long size = m_sims_sg[n].solver.nvars 
                    * m_sims_sg[n].solver.npoints_local_wghosts;
        double* uex_sg = (double*) calloc(size, sizeof(double));
  
        /* partition the global exact solution to local on this sparse grid */
        MPIPartitionArraynDwGhosts( m_ndims,
                                    (void*) &(m_sims_sg[n].mpi),
                                    (m_rank ? NULL : uex_global_sg),
                                    uex_sg,
                                    m_sims_sg[n].solver.dim_global,
                                    m_sims_sg[n].solver.dim_local,
                                    m_sims_sg[n].solver.ghosts,
                                    m_sims_sg[n].solver.nvars );
  
        /* delete the global exact solution array */
        if (!m_rank) free(uex_global_sg);
  
        /* compute the error */
        computeError( m_sims_sg[n], uex_sg);
  
        /* delete the exact solution array */
        free(uex_sg);
      }
      free(uex2);
    }

  }

  free(uex);
  return;
}

void computeError ( SimulationObject &  a_sim, double * const  a_uex  )

protected

Compute error for a simulation object

compute errors for a particular simulation object

Parameters

a_sim	Simulation object
a_uex	Exact solution

Definition at line 155 of file SparseGridsCalculateError.cpp.

{
  static const double tolerance = 1e-15;

  HyPar* solver = &(a_sim.solver);
  MPIVariables* mpi = &(a_sim.mpi);

  if (a_uex == NULL) {
    fprintf(stderr, "Error in SparseGrids::computeError() -\n");
    fprintf(stderr, "  exact solution pointer in NULL.\n");
    exit(1);
  }

  /* calculate solution norms (for relative error) */
  double sum, global_sum;
  double solution_norm[3] = {0.0,0.0,0.0};
  /* L1 */
  sum = ArraySumAbsnD   (solver->nvars,solver->ndims,solver->dim_local,
                         solver->ghosts,solver->index,a_uex);
  global_sum = 0; MPISum_double(&global_sum,&sum,1,&mpi->world);
  solution_norm[0] = global_sum/((double)solver->npoints_global);
  /* L2 */
  sum = ArraySumSquarenD(solver->nvars,solver->ndims,solver->dim_local,
                         solver->ghosts,solver->index,a_uex);
  global_sum = 0; MPISum_double(&global_sum,&sum,1,&mpi->world);
  solution_norm[1] = sqrt(global_sum/((double)solver->npoints_global));
  /* Linf */
  sum = ArrayMaxnD      (solver->nvars,solver->ndims,solver->dim_local,
                         solver->ghosts,solver->index,a_uex);
  global_sum = 0; MPIMax_double(&global_sum,&sum,1,&mpi->world);
  solution_norm[2] = global_sum;

  /* compute error = difference between exact and numerical solution */
  long size = solver->nvars*solver->npoints_local_wghosts;
  _ArrayAXPY_(solver->u,-1.0,a_uex,size);

  /* calculate L1 norm of error */
  sum = ArraySumAbsnD   (solver->nvars,solver->ndims,solver->dim_local,
                         solver->ghosts,solver->index,a_uex);
  global_sum = 0; MPISum_double(&global_sum,&sum,1,&mpi->world);
  solver->error[0] = global_sum/((double)solver->npoints_global);

  /* calculate L2 norm of error */
  sum = ArraySumSquarenD(solver->nvars,solver->ndims,solver->dim_local,
                         solver->ghosts,solver->index,a_uex);
  global_sum = 0; MPISum_double(&global_sum,&sum,1,&mpi->world);
  solver->error[1] = sqrt(global_sum/((double)solver->npoints_global));

  /* calculate Linf norm of error */
  sum = ArrayMaxnD      (solver->nvars,solver->ndims,solver->dim_local,
                         solver->ghosts,solver->index,a_uex);
  global_sum = 0; MPIMax_double(&global_sum,&sum,1,&mpi->world);
  solver->error[2] = global_sum;

  /* 
    decide whether to normalize and report relative errors, 
    or report absolute errors.
  */
  if (    (solution_norm[0] > tolerance) 
      &&  (solution_norm[1] > tolerance) 
      &&  (solution_norm[2] > tolerance) ) {
    solver->error[0] /= solution_norm[0];
    solver->error[1] /= solution_norm[1];
    solver->error[2] /= solution_norm[2];
  }

  return;
}

bool isPowerOfTwo ( int  x )

inlineprotected

Checks if an integer is a power of 2

Definition at line 317 of file sparse_grids_simulation.h.

    {
      if (x == 0)  return false;

      while (x > 1) {
        if (x%2 != 0) return false;
        x /= 2;
      }
      return true;
    }

double cfunc ( int  a, int  b  )

inlineprotected

The combination/choose function from probability stuff

Definition at line 329 of file sparse_grids_simulation.h.

    {
      return ( ( (double)factorial(a) ) / ( (double) (factorial(b)*factorial(a-b)) ) );
    }

int factorial ( int  a )

inlineprotected

Factorial function

Definition at line 335 of file sparse_grids_simulation.h.

    {
      int retval = 1.0;
      for (int i=1; i<=a; i++) retval *= i;
      return retval;
    }

void allocateGridArrays ( const GridDimensions &  a_dim, double **const  a_x, const int  a_ngpt = 0  )

inlineprotected

Allocate grid array given grid dimension

Definition at line 343 of file sparse_grids_simulation.h.

    {
      GridDimensions dim_wghosts = a_dim;
      for (int i=0; i<dim_wghosts.size(); i++) {
        dim_wghosts[i] += (2*a_ngpt);
      }
      long size_x = StdVecOps::sum(dim_wghosts);
      (*a_x) = (double*) calloc (size_x, sizeof(double));
      for (int i=0; i<size_x; i++) (*a_x)[i] = 0.0;
      return;
    }

void allocateDataArrays ( const GridDimensions &  a_dim, const int  a_nvars, double **const  a_u, const int  a_ngpt = 0  )

inlineprotected

Allocate data arrays given grid dimension

Definition at line 358 of file sparse_grids_simulation.h.

    {
      GridDimensions dim_wghosts = a_dim;
      for (int i=0; i<dim_wghosts.size(); i++) {
        dim_wghosts[i] += (2*a_ngpt);
      }
      long size_u = a_nvars*StdVecOps::product(dim_wghosts);
      (*a_u) = (double*) calloc (size_u, sizeof(double));
      for (int i=0; i<size_u; i++) (*a_u)[i] = 0.0;
      return;
    }

Field Documentation

bool m_is_defined

protected

Boolean to show if this object is defined

Definition at line 212 of file sparse_grids_simulation.h.

int m_nsims_sg

protected

Number of sparse grids simulations

Definition at line 214 of file sparse_grids_simulation.h.

int m_ndims

protected

Number of spatial dimensions

Definition at line 215 of file sparse_grids_simulation.h.

int m_rank

protected

MPI rank of this process

Definition at line 216 of file sparse_grids_simulation.h.

int m_nproc

protected

Total number of MPI ranks

Definition at line 216 of file sparse_grids_simulation.h.

int m_interp_order

protected

Order of interpolation between grids (input - sparse_grids.inp )

Definition at line 220 of file sparse_grids_simulation.h.

std::vector<bool> m_is_periodic

protected

Periodicity along each dimension

Definition at line 222 of file sparse_grids_simulation.h.

int m_write_sg_solutions

protected

Write out the sparse grid solutions to file? (input - sparse_grids.inp )

Definition at line 225 of file sparse_grids_simulation.h.

int m_print_sg_errors

protected

Print and write out the sparse grid errors? (input - sparse_grids.inp )

Definition at line 228 of file sparse_grids_simulation.h.

SimulationObject* m_sim_fg

protected

full grid simulation object

Definition at line 230 of file sparse_grids_simulation.h.

std::vector<SimulationObject> m_sims_sg

protected

vector of sparse grids simulation objects

Definition at line 231 of file sparse_grids_simulation.h.

int m_n_fg

protected

Base2 log of the number of grid points along a dimension of the full grid

Definition at line 233 of file sparse_grids_simulation.h.

int m_imin

protected

Base2 log of the minimum number of grid points along any dimension (input - sparse_grids.inp )

Definition at line 234 of file sparse_grids_simulation.h.

std::vector<SGCTElem> m_combination

protected

Coefficients and grid dimensions for the combination technique

Definition at line 236 of file sparse_grids_simulation.h.

std::vector<ProcDistribution> m_iprocs

protected

MPI ranks along each dimension for the grids in the combination technique

Definition at line 237 of file sparse_grids_simulation.h.

---

The documentation for this class was generated from the following files:

• include/sparse_grids_simulation.h
• src/SparseGrids/SparseGridsCalculateError.cpp
• src/SparseGrids/SparseGridsCombinationTechnique.cpp
• src/SparseGrids/SparseGridsDefine.cpp
• src/SparseGrids/SparseGridsFillGhostCells.cpp
• src/SparseGrids/SparseGridsInitializationWrapup.cpp
• src/SparseGrids/SparseGridsInitialize.cpp
• src/SparseGrids/SparseGridsInitialSolution.cpp
• src/SparseGrids/SparseGridsInterpolate.cpp
• src/SparseGrids/SparseGridsInterpolateGrid.cpp
• src/SparseGrids/SparseGridsMiscFuncs.cpp
• src/SparseGrids/SparseGridsOutputSolution.cpp
• src/SparseGrids/SparseGridsSanityChecks.cpp
• src/SparseGrids/SparseGridsSetSolverParameters.cpp
• src/SparseGrids/SparseGridsSolve.cpp
• src/SparseGrids/SparseGridsWriteErrors.cpp