HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
SparseGridsSimulation Class Reference

Class describing sparse grids simulations. More...

#include <sparse_grids_simulation.h>

Inheritance diagram for SparseGridsSimulation:
Simulation

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
 
SimulationObjectm_sim_fg
 
std::vector< SimulationObjectm_sims_sg
 
int m_n_fg
 
int m_imin
 
std::vector< SGCTElemm_combination
 
std::vector< ProcDistributionm_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()

Constructor

Definition at line 43 of file sparse_grids_simulation.h.

44  {
45  m_is_defined = false;
46 
47  m_nsims_sg = -1;
48  m_ndims = -1;
49  m_rank = -1;
50  m_nproc = -1;
51 
52  m_sim_fg = NULL;
53  m_sims_sg.clear();
54 
55  m_n_fg = -1;
56  m_imin = -1;
57 
58  m_interp_order = 6;
59 
60  m_is_periodic.clear();
61  m_combination.clear();
62  }
std::vector< SGCTElem > m_combination
std::vector< SimulationObject > m_sims_sg
std::vector< bool > m_is_periodic

◆ ~SparseGridsSimulation()

~SparseGridsSimulation ( )
inline

Destructor

Definition at line 65 of file sparse_grids_simulation.h.

66  {
67  int err;
68  /* clean up full grid object */
70  delete m_sim_fg;
71  /* clean up and delete sparse grids objects */
72  err = Cleanup((void*) m_sims_sg.data(), m_nsims_sg);
73  if (err) {
74  printf( "Error: CleanUp() returned with status %d on process %d.\n",
75  err, m_rank );
76  }
77  m_sims_sg.clear();
78  /* done */
79  return;
80  }
std::vector< SimulationObject > m_sims_sg
int Cleanup(void *, int)
Definition: Cleanup.c:39
int CleanupBarebones(SimulationObject *)

Member Function Documentation

◆ define()

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_rankMPI rank of this process
a_nprocTotal number of MPI ranks

Implements Simulation.

Definition at line 33 of file SparseGridsDefine.cpp.

36 {
37  m_rank = a_rank;
38  m_nproc = a_nproc;
39 
40  if (m_sim_fg != NULL) {
41  fprintf(stderr, "Error in SparseGridsSimulation::define() -\n");
42  fprintf(stderr, " m_sim_fg not NULL!\n");
43  exit(1);
44  }
45 
46  /* default values */
47  m_imin = 2;
50  m_interp_order = 6;
51 
52  if (!m_rank) {
53 
54  FILE *in;
55  in = fopen(_SPARSEGRIDS_SIM_INP_FNAME_,"r");
56 
57  if (!in) {
58 
59  fprintf(stderr, "Error in SparseGridsSimulation::Define() -\n");
60  fprintf(stderr, " %s file not found.\n", _SPARSEGRIDS_SIM_INP_FNAME_);
61  exit(1);
62 
63  } else {
64 
65  int ferr;
66  char word[_MAX_STRING_SIZE_];
67  ferr = fscanf(in,"%s", word); if (ferr != 1) return(1);
68 
69  if (std::string(word) == "begin") {
70 
71  while (std::string(word) != "end") {
72 
73  ferr = fscanf(in,"%s",word); if (ferr != 1) return(1);
74 
75  if (std::string(word) == "log2_imin") {
76 
77  ferr = fscanf(in,"%d",&m_imin); if (ferr != 1) return(1);
78 
79  } else if (std::string(word) == "interp_order") {
80 
81  ferr = fscanf(in,"%d",&m_interp_order); if (ferr != 1) return(1);
82 
83  } else if (std::string(word) == "write_sg_solutions") {
84 
85  char answer[_MAX_STRING_SIZE_];
86  ferr = fscanf(in,"%s",answer); if (ferr != 1) return(1);
87  m_write_sg_solutions = (strcmp(answer,"yes") ? 0 : 1);
88 
89  } else if (std::string(word) == "write_sg_errors") {
90 
91  char answer[_MAX_STRING_SIZE_];
92  ferr = fscanf(in,"%s",answer); if (ferr != 1) return(1);
93  m_print_sg_errors = (strcmp(answer,"yes") ? 0 : 1);
94 
95  } else if (std::string(word) != "end") {
96 
97  char useless[_MAX_STRING_SIZE_];
98  ferr = fscanf(in,"%s",useless);
99  printf("Warning: keyword %s in file \"%s\" with value %s not recognized or extraneous. Ignoring.\n",
100  _SPARSEGRIDS_SIM_INP_FNAME_, word, useless );
101 
102  }
103 
104  if (ferr != 1) return(1);
105  }
106 
107  } else {
108  fprintf( stderr, "Error: Illegal format in file \"%s\". Word read is: %s\n",
110  return 1;
111  }
112 
113  fclose(in);
114 
115  }
116 
117  /* print useful stuff to screen */
118  printf("Sparse grids inputs:\n");
119  printf(" log2 of minimum grid size: %d\n", m_imin);
120  printf(" interpolation order: %d\n", m_interp_order);
121  printf( " write sparse grids solutions? %s\n",
122  ( m_write_sg_solutions == 1 ? "yes" : "no" ) );
123  }
124 
125 #ifndef serial
126  MPI_Bcast(&m_imin,1,MPI_INT,0,MPI_COMM_WORLD);
127  MPI_Bcast(&m_interp_order,1,MPI_INT,0,MPI_COMM_WORLD);
128  MPI_Bcast(&m_write_sg_solutions,1,MPI_INT,0,MPI_COMM_WORLD);
129  MPI_Bcast(&m_print_sg_errors,1,MPI_INT,0,MPI_COMM_WORLD);
130 #endif
131 
133  /* the following are deliberately set to junk values
134  * to ensure they are never used */
135  m_sim_fg->solver.my_idx = -1;
136  m_sim_fg->solver.nsims = -1;
137  /* these are okay */
138  m_sim_fg->mpi.rank = m_rank;
140 
141  if (!m_rank) {
142  printf("Allocated full grid simulation object(s).\n");
143  }
144 
145  m_is_defined = true;
146  return 0;
147 }
Structure defining a simulation.
int nsims
Definition: hypar.h:64
#define _MAX_STRING_SIZE_
Definition: basic.h:14
#define _SPARSEGRIDS_SIM_INP_FNAME_
int my_idx
Definition: hypar.h:61

◆ ReadInputs()

int ReadInputs ( )
inlinevirtual

Read solver inputs for the full grid simulation object

Implements Simulation.

Definition at line 86 of file sparse_grids_simulation.h.

87  {
88  int retval = ::ReadInputs( (void*) m_sim_fg,
89  1,
90  m_rank );
91  return retval;
92  }

◆ Initialize()

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.

16 {
17  int ierr;
18 
19  /* find out the number of spatial dimensions */
21 
22  /* make sure full grid simulation object is allocated */
23  if (m_sim_fg == NULL) {
24  fprintf(stderr, "Error in SparseGridsSimulation::Initialize()\n");
25  fprintf(stderr, " m_sim_fg is NULL on rank %d!\n", m_rank);
26  return 1;
27  }
28 
29  /* set full grid processor distribution based on grid size and
30  * number of MPI ranks */
31  GridDimensions dim_fg(m_ndims);
32  for (int d=0; d<m_ndims; d++) dim_fg[d] = m_sim_fg->solver.dim_global[d];
33  ProcDistribution iproc_fg;
34  ComputeProcessorDistribution(iproc_fg, dim_fg);
35  for (int d=0; d<m_ndims; d++) m_sim_fg->mpi.iproc[d] = iproc_fg[d];
37  m_ndims,
38  0,
39  &(m_sim_fg->mpi.world)); CHECKERR(ierr);
40 
41  ::WriteInputs( (void*) m_sim_fg, 1, m_rank);
42  if (!m_rank) {
43  printf("Processor distribution for full grid object: ");
44  for (int d=0; d<m_ndims; d++) printf(" %3d", iproc_fg[d]);
45  printf("\n");
46  }
47 
48  if (!m_rank) printf("Initializing sparse grids...\n");
49  if (!m_rank) printf(" Number of spatial dimensions: %d\n", m_ndims);
50 
51  /* some sanity checks */
52  ierr = SanityChecks();
53  if (ierr) return ierr;
54 
55  /* compute the sizes of grids that go into the combination technique */
56  if (!m_rank) {
57  printf(" Computing sparse grids dimensions...\n");
58  }
59  ierr = ComputeSGDimsAndCoeffs();
60  if (ierr) return ierr;
61  m_nsims_sg = m_combination.size();
62  if (!m_rank) {
63  printf( " Number of sparse grid domains in combination technique: %d\n",
64  m_nsims_sg );
65  }
66  if (m_nsims_sg == 0) {
67  fprintf(stderr, "Error in SparseGridsSimulation::Initialize()\n");
68  fprintf(stderr, " ComputeSGDimsAndCoeffs() returned empty vector!\n");
69  return 1;
70  }
71 
72  /* compute processor distributions for the sparse grids */
73  if (!m_rank) {
74  printf(" Computing processor decompositions...\n");
75  }
76  m_iprocs.resize(m_nsims_sg);
77  for (int i = 0; i < m_nsims_sg; i++) {
79  m_combination[i]._dim_ );
80 
81  if (ierr) return ierr;
82  }
83 
84  /* Print some stuff to screen */
85  if (!m_rank) {
86  printf("Sparse Grids: Combination technique grids sizes and coefficients are:-\n");
87  for (int i = 0; i < m_nsims_sg; i++) {
88  printf(" %3d: dim = (", i);
89  for (int d=0; d<m_ndims; d++) printf(" %4d ", m_combination[i]._dim_[d]);
90  printf("), coeff = %+1.2e, ", m_combination[i]._coeff_);
91  printf("iproc = (");
92  for (int d=0; d<m_ndims; d++) printf(" %4d ", m_iprocs[i][d]);
93  printf(")\n");
94  }
95  }
96 
97  /* allocate full grid data structures (only what is needed) */
99  if (ierr) return ierr;
100 
101  /* create sparse grids simulation objects */
102  m_sims_sg.resize(m_nsims_sg);
103  /* set the solver parameters for the sparse grids sim objects */
104  for (int i = 0; i < m_nsims_sg; i++) {
105 #ifndef serial
106  MPI_Comm_dup(MPI_COMM_WORLD, &(m_sims_sg[i].mpi.world));
107 #endif
108  ierr = SetSolverParameters( m_sims_sg[i],
109  m_combination[i]._dim_,
110  m_iprocs[i],
111  *m_sim_fg,
112  i,
113  m_nsims_sg );
114  if (ierr) return ierr;
115  }
116 
117  /* allocate their data structures (everything) */
118  ierr = ::Initialize( (void*) m_sims_sg.data(), m_nsims_sg);
119  if (ierr) return ierr;
120 
121  /* compute and report total NDOFs for sparse grids */
122  long ndof_sg = 0;
123  for (int i = 0; i < m_nsims_sg; i++) {
124  ndof_sg += m_sims_sg[i].solver.npoints_global;
125  }
126  long ndof_fg = m_sim_fg->solver.npoints_global;
127  if (!m_rank) {
128  printf("Total number of DOFs:-\n");
129  printf(" using sparse grids: %d\n", (int)ndof_sg);
130  printf(" using conventional grid: %d\n", (int)ndof_fg);
131  }
132 
133  /* done */
134  return 0;
135 }
int InitializeBarebones(SimulationObject *)
#define CHECKERR(ierr)
Definition: basic.h:18
#define _dim_
std::vector< int > GridDimensions
int MPIBroadcast_integer(int *, int, int, void *)
Definition: MPIBroadcast.c:23
int ndims
Definition: hypar.h:26
#define _coeff_
int npoints_global
Definition: hypar.h:40
std::vector< SGCTElem > m_combination
int ComputeProcessorDistribution(ProcDistribution &, const GridDimensions &)
MPI_Comm world
std::vector< SimulationObject > m_sims_sg
std::vector< int > ProcDistribution
int SetSolverParameters(SimulationObject &, const GridDimensions &, const ProcDistribution &, const SimulationObject &, const int, const int)
std::vector< ProcDistribution > m_iprocs
int WriteInputs(void *, int, int)
Definition: WriteInputs.c:15
int * dim_global
Definition: hypar.h:33

◆ InitialSolution()

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.

18 {
19  int ierr;
20 
21  /* first, read in the full grid initial solution */
22  int retval = ::InitialSolution( (void*) m_sim_fg, 1 );
23  if (retval) return retval;
24 
25  /* now, interpolate to all the sparse grids */
26  for (int n = 0; n < m_nsims_sg; n++) {
27 
28  if (!m_rank) {
29  printf("Interpolating grid coordinates to sparse grids domain %d.\n", n);
30  }
32 
33  HyPar* solver = &(m_sims_sg[n].solver);
34  MPIVariables* mpi = &(m_sims_sg[n].mpi);
35 
36  int nvars = solver->nvars;
37  int ghosts = solver->ghosts;
38  int *dim_local = solver->dim_local;
39 
40  /* calculate dxinv */
41  {
42  int offset = 0;
43  for (int d = 0; d < m_ndims; d++) {
44  for (int i = 0; i < dim_local[d]; i++) {
45  solver->dxinv[i+offset+ghosts]
46  = 2.0 / ( solver->x[i+1+offset+ghosts]
47  - solver->x[i-1+offset+ghosts] );
48  }
49  offset += (dim_local[d] + 2*ghosts);
50  }
51  }
52 
53  /* exchange MPI-boundary values of dxinv between processors */
54  {
55  int offset = 0;
56  for (int d = 0; d < m_ndims; d++) {
57  ierr = MPIExchangeBoundaries1D( mpi,
58  &(solver->dxinv[offset]),
59  dim_local[d],
60  ghosts,
61  d,
62  m_ndims );
63  if (ierr) return ierr;
64  offset += (dim_local[d] + 2*ghosts);
65  }
66  }
67 
68  /* fill in ghost values of dxinv at physical boundaries by extrapolation */
69  {
70  int offset = 0;
71  for (int d = 0; d < m_ndims; d++) {
72  double *dxinv = &(solver->dxinv[offset]);
73  int *dim = dim_local;
74  if (mpi->ip[d] == 0) {
75  /* fill left boundary along this dimension */
76  for (int i = 0; i < ghosts; i++) dxinv[i] = dxinv[ghosts];
77  }
78  if (mpi->ip[d] == mpi->iproc[d]-1) {
79  /* fill right boundary along this dimension */
80  for (int i = dim[d]+ghosts; i < dim[d]+2*ghosts; i++) {
81  dxinv[i] = dxinv[dim[d]+ghosts-1];
82  }
83  }
84  offset += (dim[d] + 2*ghosts);
85  }
86  }
87 
88  }
89 
90  return 0;
91 }
int nvars
Definition: hypar.h:29
void interpolateGrid(SimulationObject *const, const SimulationObject *const)
double * x
Definition: hypar.h:107
int MPIExchangeBoundaries1D(void *, double *, int, int, int, int)
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
std::vector< SimulationObject > m_sims_sg
int * dim_local
Definition: hypar.h:37
int ghosts
Definition: hypar.h:52
Structure of MPI-related variables.
double * dxinv
Definition: hypar.h:110

◆ InitializeBoundaries()

int InitializeBoundaries ( )
inlinevirtual

Initialize the boundary conditions of the simulation

Implements Simulation.

Definition at line 101 of file sparse_grids_simulation.h.

102  {
103  int retval = ::InitializeBoundaries( (void*) m_sims_sg.data(),
104  m_nsims_sg );
105  if (m_nsims_sg > 0) {
106  m_is_periodic.resize(m_ndims);
107  for (int d=0; d<m_ndims; d++) {
108  m_is_periodic[d] = ( m_sims_sg[0].solver.isPeriodic[d] == 1 ? true : false);
109  m_sim_fg->solver.isPeriodic[d] = m_sims_sg[0].solver.isPeriodic[d];
110  }
111  }
112  return retval;
113  }
std::vector< SimulationObject > m_sims_sg
std::vector< bool > m_is_periodic
int * isPeriodic
Definition: hypar.h:162

◆ InitializeImmersedBoundaries()

int InitializeImmersedBoundaries ( )
inlinevirtual

Initialize the immersed boundary conditions of the simulation

Implements Simulation.

Definition at line 116 of file sparse_grids_simulation.h.

117  {
118  int retval = ::InitializeImmersedBoundaries( (void*) m_sims_sg.data(),
119  m_nsims_sg );
120  return retval;
121  }
std::vector< SimulationObject > m_sims_sg

◆ InitializePhysics()

int InitializePhysics ( )
inlinevirtual

Initialize the physics of the simulations

Implements Simulation.

Definition at line 124 of file sparse_grids_simulation.h.

125  {
126  int retval = ::InitializePhysics( (void*) m_sims_sg.data(),
127  m_nsims_sg );
128  return retval;
129  }
std::vector< SimulationObject > m_sims_sg

◆ InitializePhysicsData()

int InitializePhysicsData ( )
inlinevirtual

Initialize the physics of the simulations

Implements Simulation.

Definition at line 132 of file sparse_grids_simulation.h.

133  {
134  for (int ns = 0; ns < m_sims_sg.size(); ns++) {
135  int retval = ::InitializePhysicsData( (void*) &(m_sims_sg[ns]),
136  -1,
137  m_nsims_sg,
139  if (retval) {
140  fprintf(stderr, "Error in SparseGridsSimulation::InitializePhysicsData()\n");
141  fprintf(stderr, " InitializePhysicsData returned with error code %d on rank %d.\n",
142  retval, m_sims_sg[ns].mpi.rank);
143  return retval;
144  }
145  }
146  return 0;
147  }
std::vector< SimulationObject > m_sims_sg
int * dim_global
Definition: hypar.h:33

◆ InitializeSolvers()

int InitializeSolvers ( )
inlinevirtual

Initialize the numerical solvers of the simulations

Implements Simulation.

Definition at line 150 of file sparse_grids_simulation.h.

151  {
152  int retval = ::InitializeSolvers( (void*) m_sims_sg.data(),
153  m_nsims_sg );
155 
156  /* some modifications to output filename roots */
157  for (int i=0; i<m_nsims_sg; i++) {
158  strcpy(m_sims_sg[i].solver.op_fname_root, "op_sg");
159  strcpy(m_sims_sg[i].solver.aux_op_fname_root, "ts0_sg");
160  }
161  strcpy(m_sim_fg->solver.op_fname_root, "op_fg");
162  strcpy(m_sim_fg->solver.aux_op_fname_root, "ts0_fg");
163  return retval;
164  }
int InitializeSolversBarebones(SimulationObject *)
char aux_op_fname_root[_MAX_STRING_SIZE_]
Definition: hypar.h:208
std::vector< SimulationObject > m_sims_sg
char op_fname_root[_MAX_STRING_SIZE_]
Definition: hypar.h:206

◆ InitializationWrapup()

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.

16 {
17  int ierr;
18 
19  for (int n = 0; n < m_nsims_sg; n++) {
20 
21  if (!m_rank) {
22  printf("Interpolating initial solution to sparse grids domain %d.\n", n);
23  }
25 
26  HyPar* solver = &(m_sims_sg[n].solver);
27  MPIVariables* mpi = &(m_sims_sg[n].mpi);
28 
29  int nvars = solver->nvars;
30  int ghosts = solver->ghosts;
31  int *dim_local = solver->dim_local;
32 
33  /* exchange MPI-boundary values of u between processors */
35  nvars,
36  dim_local,
37  ghosts,
38  mpi,
39  solver->u );
40 
41  /* calculate volume integral of the initial solution */
43  solver->u,
44  solver,
45  mpi );
46  if (ierr) return ierr;
47 
48  if (!mpi->rank) {
49  printf(" Volume integral of the initial solution on sparse grids domain %d:\n", n);
50  for (int d=0; d<nvars; d++) {
51  printf(" %2d: %1.16E\n",d,solver->VolumeIntegralInitial[d]);
52  }
53  }
54 
55  /* Set initial total boundary flux integral to zero */
56  _ArraySetValue_( solver->TotalBoundaryIntegral, nvars, 0 );
57 
58  }
59 
60  return 0;
61 }
int nvars
Definition: hypar.h:29
int MPIExchangeBoundariesnD(int, int, int *, int, void *, double *)
double * u
Definition: hypar.h:116
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
std::vector< SimulationObject > m_sims_sg
#define _ArraySetValue_(x, size, value)
int * dim_local
Definition: hypar.h:37
void interpolate(SimulationObject *const, const SimulationObject *const)
double * TotalBoundaryIntegral
Definition: hypar.h:386
int ghosts
Definition: hypar.h:52
Structure of MPI-related variables.
int VolumeIntegral(double *, double *, void *, void *)
double * VolumeIntegralInitial
Definition: hypar.h:380

◆ Solve()

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.

18 {
19  int tic = 0;
20 
21  /* write out iblank to file for visualization */
22  if (m_write_sg_solutions == 1) {
23  for (int ns = 0; ns < m_nsims_sg; ns++) {
24  if (m_sims_sg[ns].solver.flag_ib) {
25 
26  char fname_root[_MAX_STRING_SIZE_] = "iblank_sg";
27  if (m_nsims_sg > 1) {
28  char index[_MAX_STRING_SIZE_];
29  GetStringFromInteger(ns, index, (int)log10((m_nsims_sg)+1));
30  strcat(fname_root, "_");
31  strcat(fname_root, index);
32  }
33 
34  WriteArray( m_sims_sg[ns].solver.ndims,
35  1,
36  m_sims_sg[ns].solver.dim_global,
37  m_sims_sg[ns].solver.dim_local,
38  m_sims_sg[ns].solver.ghosts,
39  m_sims_sg[ns].solver.x,
40  m_sims_sg[ns].solver.iblank,
41  &(m_sims_sg[ns].solver),
42  &(m_sims_sg[ns].mpi),
43  fname_root );
44  }
45  }
46  }
47  if (m_sim_fg->solver.flag_ib) {
48 
49  char fname_root[_MAX_STRING_SIZE_] = "iblank";
51  1,
55  m_sim_fg->solver.x,
57  &(m_sim_fg->solver),
58  &(m_sim_fg->mpi),
59  fname_root );
60  }
61 
62  /* Define and initialize the time-integration object */
63  TimeIntegration TS;
64  if (!m_rank) printf("Setting up time integration.\n");
65  TimeInitialize((void*)m_sims_sg.data(), m_nsims_sg, m_rank, m_nproc, &TS);
66 
67  if (!m_rank) {
68  printf( "Solving in time (from %d to %d iterations)\n",
69  TS.restart_iter,TS.n_iter);
70  }
71 
72  for (TS.iter = TS.restart_iter; TS.iter < TS.n_iter; TS.iter++) {
73 
74  /* Write initial solution to file if this is the first iteration */
75  if (!TS.iter) OutputSolution(TS.waqt);
76 
77  /* Call pre-step function */
78  TimePreStep (&TS);
79 
80  /* Step in time */
81  TimeStep (&TS);
82 
83  /* Call post-step function */
84  TimePostStep (&TS);
85 
86  /* Print information to screen */
87  TimePrintStep(&TS);
88  tic++;
89 
90  /* Write intermediate solution to file */
91  if ((TS.iter+1)%m_sims_sg[0].solver.file_op_iter == 0) {
92  OutputSolution(TS.waqt);
93  tic = 0;
94  }
95 
96  }
97 
98  /* write a final solution file, if last iteration did not write one */
99  if (tic || (!TS.n_iter)) OutputSolution(TS.waqt);
100 
101  if (!m_rank) {
102  printf("Completed time integration (Final time: %f).\n",TS.waqt);
103  if (m_nsims_sg > 1) printf("\n");
104  }
105 
106  /* calculate error if exact solution has been provided */
107  CalculateError();
108 
109  TimeCleanup(&TS);
110 
111  return(0);
112 }
Structure of variables/parameters and function pointers for time integration.
int TimeInitialize(void *, int, int, int, void *)
int TimePreStep(void *)
Definition: TimePreStep.c:22
int flag_ib
Definition: hypar.h:441
double * iblank
Definition: hypar.h:436
int TimePrintStep(void *)
Definition: TimePrintStep.c:16
double * x
Definition: hypar.h:107
int ndims
Definition: hypar.h:26
int TimePostStep(void *)
Definition: TimePostStep.c:28
#define _MAX_STRING_SIZE_
Definition: basic.h:14
int WriteArray(int, int, int *, int *, int, double *, double *, void *, void *, char *)
Definition: WriteArray.c:27
std::vector< SimulationObject > m_sims_sg
void GetStringFromInteger(int, char *, int)
int TimeCleanup(void *)
Definition: TimeCleanup.c:18
int * dim_local
Definition: hypar.h:37
int TimeStep(void *)
Definition: TimeStep.c:13
int ghosts
Definition: hypar.h:52
int * dim_global
Definition: hypar.h:33

◆ WriteErrors()

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_runtimeMeasured runtime of solver
main_runtimeMeasured total runtime

Implements Simulation.

Definition at line 14 of file SparseGridsWriteErrors.cpp.

17 {
18  if (!m_rank) {
19 
20  /* Write sparse grids stuff, if asked for */
21  if (m_print_sg_errors == 1) {
22  for (int n = 0; n < m_nsims_sg; n++) {
23 
24  char err_fname[_MAX_STRING_SIZE_],
25  cons_fname[_MAX_STRING_SIZE_],
26  fc_fname[_MAX_STRING_SIZE_];
27 
28  strcpy(err_fname, "errors");
29  strcpy(cons_fname,"conservation");
30  strcpy(fc_fname, "function_counts");
31 
32  if (m_nsims_sg > 1) {
33 
34  strcat(err_fname,"_");
35  strcat(cons_fname,"_");
36  strcat(fc_fname,"_");
37 
38  char index[_MAX_STRING_SIZE_];
39  GetStringFromInteger(n, index, (int)log10(m_nsims_sg)+1);
40 
41  strcat(err_fname,index);
42  strcat(cons_fname,index);
43  strcat(fc_fname,index);
44  }
45 
46  strcat(err_fname,".dat");
47  strcat(cons_fname,".dat");
48  strcat(fc_fname,".dat");
49 
50  FILE *out;
51  /* write out solution errors and wall times to file */
52  int d;
53  out = fopen(err_fname,"w");
54  for (d=0; d<m_sims_sg[n].solver.ndims; d++) fprintf(out,"%4d ",m_sims_sg[n].solver.dim_global[d]);
55  for (d=0; d<m_sims_sg[n].solver.ndims; d++) fprintf(out,"%4d ",m_sims_sg[n].mpi.iproc[d]);
56  fprintf(out,"%1.16E ",m_sims_sg[n].solver.dt);
57  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]);
58  fprintf(out,"%1.16E %1.16E\n",solver_runtime,main_runtime);
59  fclose(out);
60  /* write out conservation errors to file */
61  out = fopen(cons_fname,"w");
62  for (d=0; d<m_sims_sg[n].solver.ndims; d++) fprintf(out,"%4d ",m_sims_sg[n].solver.dim_global[d]);
63  for (d=0; d<m_sims_sg[n].solver.ndims; d++) fprintf(out,"%4d ",m_sims_sg[n].mpi.iproc[d]);
64  fprintf(out,"%1.16E ",m_sims_sg[n].solver.dt);
65  for (d=0; d<m_sims_sg[n].solver.nvars; d++) fprintf(out,"%1.16E ",m_sims_sg[n].solver.ConservationError[d]);
66  fprintf(out,"\n");
67  fclose(out);
68  /* write out function call counts to file */
69  out = fopen(fc_fname,"w");
70  fprintf(out,"%d\n",m_sims_sg[n].solver.n_iter);
71  fprintf(out,"%d\n",m_sims_sg[n].solver.count_hyp);
72  fprintf(out,"%d\n",m_sims_sg[n].solver.count_par);
73  fprintf(out,"%d\n",m_sims_sg[n].solver.count_sou);
74  #ifdef with_petsc
75  fprintf(out,"%d\n",m_sims_sg[n].solver.count_RHSFunction);
76  fprintf(out,"%d\n",m_sims_sg[n].solver.count_IFunction);
77  fprintf(out,"%d\n",m_sims_sg[n].solver.count_IJacobian);
78  fprintf(out,"%d\n",m_sims_sg[n].solver.count_IJacFunction);
79  #endif
80  fclose(out);
81 
82  /* print solution errors, conservation errors, and wall times to screen */
83  if (m_sims_sg[n].solver.error[0] >= 0) {
84  printf("Computed errors for sparse grids domain %d:\n", n);
85  printf(" L1 Error: %1.16E\n",m_sims_sg[n].solver.error[0]);
86  printf(" L2 Error: %1.16E\n",m_sims_sg[n].solver.error[1]);
87  printf(" Linf Error: %1.16E\n",m_sims_sg[n].solver.error[2]);
88  }
89  if (!strcmp(m_sims_sg[n].solver.ConservationCheck,"yes")) {
90  printf("Conservation Errors:\n");
91  for (d=0; d<m_sims_sg[n].solver.nvars; d++) printf("\t%1.16E\n",m_sims_sg[n].solver.ConservationError[d]);
92  printf("\n");
93  }
94 
95  }
96  }
97 
98  /* First write stuff for the full grid solution */
99  {
100  char err_fname[_MAX_STRING_SIZE_];
101  strcpy(err_fname,"errors_fg");
102  strcat(err_fname,".dat");
103 
104  FILE *out;
105  /* write out solution errors and wall times to file */
106  int d;
107  out = fopen(err_fname,"w");
108  for (d=0; d<m_sim_fg->solver.ndims; d++) fprintf(out,"%4d ",m_sim_fg->solver.dim_global[d]);
109  for (d=0; d<m_sim_fg->solver.ndims; d++) fprintf(out,"%4d ",m_sim_fg->mpi.iproc[d]);
110  fprintf(out,"%1.16E ",m_sim_fg->solver.dt);
111  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]);
112  fprintf(out,"%1.16E %1.16E\n",solver_runtime,main_runtime);
113  fclose(out);
114 
115  /* print solution errors, conservation errors, and wall times to screen */
116  if (m_sim_fg->solver.error[0] >= 0) {
117  printf("Computed errors for full grid solution:\n");
118  printf(" L1 Error: %1.16E\n",m_sim_fg->solver.error[0]);
119  printf(" L2 Error: %1.16E\n",m_sim_fg->solver.error[1]);
120  printf(" Linf Error: %1.16E\n",m_sim_fg->solver.error[2]);
121  }
122  }
123 
124  /* report wall times */
125  printf("Solver runtime (in seconds): %1.16E\n",solver_runtime);
126  printf("Total runtime (in seconds): %1.16E\n",main_runtime);
127 
128  }
129 
130  return;
131 }
double error[3]
Definition: hypar.h:371
int ndims
Definition: hypar.h:26
#define _MAX_STRING_SIZE_
Definition: basic.h:14
std::vector< SimulationObject > m_sims_sg
void GetStringFromInteger(int, char *, int)
double dt
Definition: hypar.h:67
int * dim_global
Definition: hypar.h:33

◆ isDefined()

bool isDefined ( ) const
inlinevirtual

Return whether object is defined or not

Implements Simulation.

Definition at line 176 of file sparse_grids_simulation.h.

◆ mpiCommDup()

int mpiCommDup ( )
inlinevirtual

Create duplicate MPI communicators

Implements Simulation.

Definition at line 180 of file sparse_grids_simulation.h.

181  {
182  if (!m_sim_fg) {
183  fprintf(stderr, "Error in SparseGridsSimulation::mpiCommDup()\n");
184  fprintf(stderr, " m_sim_fg is not allocated on rank %d!\n", m_rank);
185  return 1;
186  }
187  MPI_Comm_dup(MPI_COMM_WORLD, &(m_sim_fg->mpi.world));
188  return 0;
189  }
MPI_Comm world

◆ usePetscTS()

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.

195  {
196  m_sim_fg->solver.use_petscTS = a_flag;
197  }
int use_petscTS
Definition: hypar.h:395

◆ SolvePETSc()

int SolvePETSc ( )
inlinevirtual

Run the simulation using PETSc time integrators

Implements Simulation.

Definition at line 200 of file sparse_grids_simulation.h.

201  {
202  int retval = ::SolvePETSc( (void*) m_sims_sg.data(),
203  m_nsims_sg,
204  m_rank,
205  m_nproc );
206  return retval;
207  }
std::vector< SimulationObject > m_sims_sg

◆ SanityChecks()

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.

19 {
20  int *dim_global_fg = m_sim_fg->solver.dim_global;
21  int *iproc_fg = m_sim_fg->mpi.iproc;
22 
23  /* check if grid sizes are same along all dimensions */
24  {
25  bool flag = true;
26  for (int d=0; d<m_ndims; d++) {
27  flag = flag && (dim_global_fg[d] == dim_global_fg[0]);
28  }
29  if (!flag) {
30  fprintf(stderr, "Error in SparseGridsSimulation::SanityChecks()\n");
31  fprintf(stderr, " full grid dimensions are not equal in all dimensions.\n");
32  return 1;
33  }
34  }
35 
36  /* check if grid size is a power of 2 */
37  {
38  bool flag = isPowerOfTwo(dim_global_fg[0]);
39  if (!flag) {
40  fprintf(stderr, "Error in SparseGridsSimulation::SanityChecks()\n");
41  fprintf(stderr, " full grid dimensions are not a power of 2.\n");
42  return 1;
43  }
44  }
45 
46  return 0;
47 }
int * dim_global
Definition: hypar.h:33

◆ ComputeSGDimsAndCoeffs()

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.

25 {
26  m_n_fg = log2(m_sim_fg->solver.dim_global[0]);
27  int d = m_ndims;
28 
29  double c1 = 1;
30  m_combination.clear();
31  for (int s = 1; s <= d; s++) {
32 
33  double coeff = c1 * cfunc(d-1,s-1);
34 
35  std::vector<GridDimensions> dims(0);
36  GetCTGridSizes((m_n_fg+(m_ndims-1)*(m_imin-1))+d-s, dims);
37  if (dims.size() == 0) {
38  fprintf(stderr, "Error in SparseGridsSimulation::ComputeSGDimsAndCoeffs()\n");
39  fprintf(stderr, " GetCTGridSize() returned empty vector!\n");
40  m_combination.clear();
41  return 1;
42  }
43  for (int i=0; i<dims.size(); i++) {
44  m_combination.push_back(std::pair<double,GridDimensions>(coeff,dims[i]));
45  }
46 
47  c1 *= -1;
48 
49  }
50 
51  return 0;
52 }
void GetCTGridSizes(const int, std::vector< GridDimensions > &)
std::vector< SGCTElem > m_combination
double cfunc(int a, int b)
int * dim_global
Definition: hypar.h:33

◆ ComputeProcessorDistribution()

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_iprocsProcessor distribution to compute
a_dimGrid dimensions

Definition at line 85 of file SparseGridsMiscFuncs.cpp.

87 {
88  a_iprocs.resize(m_ndims);
89  /* get the normal vector for the grid dimensions */
90  std::vector<double> dvec;
91  StdVecOps::createNormalVector(dvec, a_dim);
92 
93  /* calculate the maximum number of MPI ranks in each dimension */
94  int max_procs[m_ndims], min_procs[m_ndims];
95  for (int i=0; i<m_ndims; i++) {
96  max_procs[i] = a_dim[i]/_MIN_GRID_PTS_PER_PROC_;
97  min_procs[i] = 1;
98  }
99  int max_nproc; _ArrayProduct1D_(max_procs, m_ndims, max_nproc);
100  if (max_nproc < m_nproc) {
101  fprintf(stderr, "Error in SparseGridsSimulation::ComputeProcessorDistribution() - rank %d\n", m_rank);
102  fprintf(stderr, " Number of MPI ranks greater than the maximum number of MPI ranks that can be used.\n");
103  fprintf(stderr, " Please re-run with %d MPI ranks.\n", max_nproc);
104  return 1;
105  }
106 
107  /* find all the processor distributions that are okay, i.e., their product
108  * is the total number of MPI ranks #SparseGridsSimulation::m_nproc */
109  std::vector<ProcDistribution> iproc_candidates(0);
110  int iproc[m_ndims], ubound[m_ndims], done = 0;
111  for (int d=0; d<m_ndims; d++) ubound[d] = max_procs[d]+1;
112  _ArraySetValue_(iproc, m_ndims, 1);
113  while (!done) {
114  int prod; _ArrayProduct1D_(iproc, m_ndims, prod);
115  if (prod == m_nproc) {
116  ProcDistribution iproc_vec(m_ndims);
117  for (int d = 0; d < m_ndims; d++) iproc_vec[d] = iproc[d];
118  iproc_candidates.push_back(iproc_vec);
119  }
120  _ArrayIncrementIndexWithLBound_(m_ndims,ubound,min_procs,iproc,done);
121  }
122  if (iproc_candidates.size() == 0) {
123  fprintf(stderr, "Error in SparseGridsSimulation::ComputeProcessorDistribution() - rank %d\n", m_rank);
124  fprintf(stderr, " Unable to fine candidate iprocs!\n");
125  return 1;
126  }
127 
128  /* find the candidate that is closest to the normalized dim vector */
129  double min_norm = DBL_MAX;
130  for (int i = 0; i<iproc_candidates.size(); i++) {
131  std::vector<double> pvec;
132  StdVecOps::createNormalVector(pvec, iproc_candidates[i]);
133  double norm = StdVecOps::compute2Norm(pvec, dvec);
134  if (norm < min_norm) {
135  min_norm = norm;
136  a_iprocs = iproc_candidates[i];
137  }
138  }
139 
140  return 0;
141 }
void createNormalVector(std::vector< double > &a_normal_vec, const int *a_vec, const int a_size)
Definition: std_vec_ops.h:70
#define _MIN_GRID_PTS_PER_PROC_
Definition: basic.h:39
#define _ArrayIncrementIndexWithLBound_(N, imax, imin, i, done)
double compute2Norm(const std::vector< double > &a_a, const std::vector< double > &a_b)
Definition: std_vec_ops.h:102
#define _ArraySetValue_(x, size, value)
std::vector< int > ProcDistribution
#define _ArrayProduct1D_(x, size, p)

◆ GetCTGridSizes()

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_NDesired sum of the base2 logs of grid sizes
a_dimsVector of grid sizes

Definition at line 56 of file SparseGridsMiscFuncs.cpp.

59 {
60  a_dims.clear();
61 
62  int index[m_ndims], ubounds[m_ndims], lbounds[m_ndims];
64  _ArraySetValue_(ubounds, m_ndims, a_N);
65  _ArraySetValue_(lbounds, m_ndims, m_imin);
66 
67  int done = 0;
68  while (!done) {
69  int sum; _ArraySum1D_(index, sum, m_ndims);
70  if (sum == a_N) {
71  std::vector<int> tmp(m_ndims);
72  for (int d=0; d<m_ndims; d++) {
73  raiseto_int( tmp[d], 2,index[d] );
74  }
75  a_dims.push_back(tmp);
76  }
77  _ArrayIncrementIndexWithLBound_(m_ndims,ubounds,lbounds,index,done);
78  }
79 
80  return;
81 }
long sum(const std::vector< int > &a_iv)
Definition: std_vec_ops.h:18
#define _ArrayIncrementIndexWithLBound_(N, imax, imin, i, done)
#define _ArraySetValue_(x, size, value)
#define raiseto_int(y, x, a)
Definition: math_ops.h:42
#define _ArraySum1D_(x, a, size)

◆ InitializeBarebones()

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
simobjsimulation objects of type SimulationObject

Definition at line 194 of file SparseGridsMiscFuncs.cpp.

195 {
196  /* Set this to "true" since this is a barebones initialization */
197  simobj->is_barebones = 1;
198 
199  HyPar* solver = &(simobj->solver);
200  MPIVariables* mpi = &(simobj->mpi);
201 
202  /* allocations */
203  mpi->ip = (int*) calloc (m_ndims,sizeof(int));
204  mpi->is = (int*) calloc (m_ndims,sizeof(int));
205  mpi->ie = (int*) calloc (m_ndims,sizeof(int));
206  mpi->bcperiodic = (int*) calloc (m_ndims,sizeof(int));
207  solver->dim_local = (int*) calloc (m_ndims,sizeof(int));
208  solver->isPeriodic = (int*) calloc (m_ndims,sizeof(int));
209 
210 #ifndef serial
211  /* Domain partitioning */
212  int total_proc = 1;
213  for (int i=0; i<m_ndims; i++) total_proc *= mpi->iproc[i];
214 
215  /* calculate ndims-D rank of each process (ip[]) from rank in MPI_COMM_WORLD */
217  m_rank,
218  mpi->iproc,
219  mpi->ip ); CHECKERR(ierr);
220 
221  /* calculate local domain sizes along each dimension */
222  for (int i=0; i<m_ndims; i++) {
223  solver->dim_local[i] = MPIPartition1D( solver->dim_global[i],
224  mpi->iproc[i],
225  mpi->ip[i] );
226  }
227 
228  /* calculate local domain limits in terms of global domain */
230  m_rank,
231  &(simobj->mpi),
232  solver->dim_global,
233  mpi->is,
234  mpi->ie );
235  CHECKERR(ierr);
236 
237  /* create sub-communicators for parallel computations
238  * along grid lines in each dimension */
239  IERR MPICreateCommunicators(m_ndims,&(simobj->mpi)); CHECKERR(ierr);
240 
241  /* initialize periodic BC flags to zero */
242  for (int i=0; i<solver->ndims; i++) mpi->bcperiodic[i] = 0;
243 
244  /* create communication groups */
245  IERR MPICreateIOGroups(&(simobj->mpi)); CHECKERR(ierr);
246 
247 #else
248 
249  for (int i=0; i<m_ndims; i++) {
250  mpi->ip[i] = 0;
251  solver->dim_local[i] = solver->dim_global[i];
252  mpi->iproc[i] = 1;
253  mpi->is[i] = 0;
254  mpi->ie[i] = solver->dim_local[i];
255  mpi->bcperiodic[i] = 0;
256  }
257 
258 #endif
259 
260  solver->npoints_global
261  = solver->npoints_local
262  = solver->npoints_local_wghosts
263  = 1;
264  for (int i=0; i<m_ndims; i++) {
265  solver->npoints_global *= solver->dim_global[i];
266  solver->npoints_local *= solver->dim_local [i];
267  solver->npoints_local_wghosts *= (solver->dim_local[i]+2*solver->ghosts);
268  }
269 
270  /* Allocations */
271  if (!m_rank) printf("Allocating data arrays for full grid.\n");
272  solver->index = (int*) calloc (m_ndims,sizeof(int));
273  solver->stride_with_ghosts = (int*) calloc (solver->ndims,sizeof(int));
274  solver->stride_without_ghosts = (int*) calloc (solver->ndims,sizeof(int));
275  int accu1 = 1, accu2 = 1;
276  for (int i=0; i<solver->ndims; i++) {
277  solver->stride_with_ghosts[i] = accu1;
278  solver->stride_without_ghosts[i] = accu2;
279  accu1 *= (solver->dim_local[i]+2*solver->ghosts);
280  accu2 *= solver->dim_local[i];
281  }
282 
283 
284  int size;
285  /* state variables */
286  size = 1;
287  for (int i=0; i<m_ndims; i++) size *= (solver->dim_local[i]+2*solver->ghosts);
288  solver->u = (double*) calloc (solver->nvars*size,sizeof(double));
289 
290  /* grid */
291  size = 0;
292  for (int i=0; i<m_ndims; i++) size += (solver->dim_local[i]+2*solver->ghosts);
293  solver->x = (double*) calloc (size,sizeof(double));
294  solver->dxinv = (double*) calloc (size,sizeof(double));
295 
296  /* allocate MPI send/receive buffer arrays */
297  int bufdim[solver->ndims], maxbuf = 0;
298  for (int d = 0; d < solver->ndims; d++) {
299  bufdim[d] = 1;
300  for (int i = 0; i < solver->ndims; i++) {
301  if (i == d) bufdim[d] *= solver->ghosts;
302  else bufdim[d] *= solver->dim_local[i];
303  }
304  if (bufdim[d] > maxbuf) maxbuf = bufdim[d];
305  }
306  maxbuf *= solver->nvars;
307  mpi->maxbuf = maxbuf;
308  mpi->sendbuf = (double*) calloc (2*solver->ndims*maxbuf,sizeof(double));
309  mpi->recvbuf = (double*) calloc (2*solver->ndims*maxbuf,sizeof(double));
310 
311  solver->VolumeIntegral = (double*) calloc (solver->nvars,sizeof(double));
312  solver->VolumeIntegralInitial = (double*) calloc (solver->nvars,sizeof(double));
313  solver->TotalBoundaryIntegral = (double*) calloc (solver->nvars,sizeof(double));
314  solver->ConservationError = (double*) calloc (solver->nvars,sizeof(double));
315 
316  for (int i=0; i<solver->nvars; i++) solver->ConservationError[i] = -1;
317 
318  return(0);
319 }
double * recvbuf
int nvars
Definition: hypar.h:29
#define IERR
Definition: basic.h:16
#define CHECKERR(ierr)
Definition: basic.h:18
double * u
Definition: hypar.h:116
int npoints_local
Definition: hypar.h:42
double * ConservationError
Definition: hypar.h:374
double * x
Definition: hypar.h:107
int * stride_with_ghosts
Definition: hypar.h:414
int ndims
Definition: hypar.h:26
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
int * stride_without_ghosts
Definition: hypar.h:416
double * sendbuf
int npoints_global
Definition: hypar.h:40
int * index
Definition: hypar.h:102
int MPIPartition1D(int, int, int)
int * dim_local
Definition: hypar.h:37
double * TotalBoundaryIntegral
Definition: hypar.h:386
int MPILocalDomainLimits(int, int, void *, int *, int *, int *)
int ghosts
Definition: hypar.h:52
int MPICreateIOGroups(void *)
Definition: MPIIOGroups.c:37
Structure of MPI-related variables.
int npoints_local_wghosts
Definition: hypar.h:42
int MPIRanknD(int, int, int *, int *)
Definition: MPIRanknD.c:27
int * dim_global
Definition: hypar.h:33
int * isPeriodic
Definition: hypar.h:162
double * VolumeIntegral
Definition: hypar.h:378
int MPICreateCommunicators(int, void *)
double * VolumeIntegralInitial
Definition: hypar.h:380
double * dxinv
Definition: hypar.h:110

◆ InitializeSolversBarebones()

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
simsimulation objects of type SimulationObject

Definition at line 324 of file SparseGridsMiscFuncs.cpp.

325 {
326  if (sim->is_barebones != 1) {
327  fprintf(stderr, "Error in SparseGridsSimulation::InitializeSolversBarebones() - \n");
328  fprintf(stderr, " simulation object is not barebones type.\n");
329  }
330 
331  int ns;
332 
333  HyPar *solver = &(sim->solver);
334  MPIVariables *mpi = &(sim->mpi);
335 
336  solver->ParabolicFunction = NULL;
337  solver->SecondDerivativePar = NULL;
338  solver->FirstDerivativePar = NULL;
339  solver->InterpolateInterfacesPar = NULL;
340 
341  solver->interp = NULL;
342  solver->compact = NULL;
343  solver->lusolver = NULL;
344  solver->SetInterpLimiterVar = NULL;
345  solver->flag_nonlinearinterp = 0;
346  solver->time_integrator = NULL;
347  solver->msti = NULL;
348 
349  /* Solution output function */
350  solver->WriteOutput = NULL; /* default - no output */
351  solver->filename_index = NULL;
352  if (!strcmp(solver->output_mode,"serial")) {
353  solver->index_length = 5;
354  solver->filename_index = (char*) calloc (solver->index_length+1,sizeof(char));
355  int i; for (i=0; i<solver->index_length; i++) solver->filename_index[i] = '0';
356  solver->filename_index[solver->index_length] = (char) 0;
357  if (!strcmp(solver->op_file_format,"text")) {
358  solver->WriteOutput = WriteText;
359  strcpy(solver->solnfilename_extn,".dat");
360  } else if (!strcmp(solver->op_file_format,"tecplot2d")) {
361  solver->WriteOutput = WriteTecplot2D;
362  strcpy(solver->solnfilename_extn,".dat");
363  } else if (!strcmp(solver->op_file_format,"tecplot3d")) {
364  solver->WriteOutput = WriteTecplot3D;
365  strcpy(solver->solnfilename_extn,".dat");
366  } else if ((!strcmp(solver->op_file_format,"binary")) || (!strcmp(solver->op_file_format,"bin"))) {
367  solver->WriteOutput = WriteBinary;
368  strcpy(solver->solnfilename_extn,".bin");
369  } else if (!strcmp(solver->op_file_format,"none")) {
370  solver->WriteOutput = NULL;
371  } else {
372  fprintf(stderr,"Error (domain %d): %s is not a supported file format.\n",
373  ns, solver->op_file_format);
374  return(1);
375  }
376  if ((!strcmp(solver->op_overwrite,"no")) && solver->restart_iter) {
377  /* if it's a restart run, fast-forward the filename */
378  int t;
379  for (t=0; t<solver->restart_iter; t++)
380  if ((t+1)%solver->file_op_iter == 0) IncrementFilenameIndex(solver->filename_index,solver->index_length);
381  }
382  } else if (!strcmp(solver->output_mode,"parallel")) {
383  if (!strcmp(solver->op_file_format,"none")) solver->WriteOutput = NULL;
384  else {
385  /* only binary file writing supported in parallel mode */
386  /* use post-processing scripts to convert */
387  solver->WriteOutput = WriteBinary;
388  strcpy(solver->solnfilename_extn,".bin");
389  }
390  } else {
391 
392  fprintf(stderr,"Error (domain %d): %s is not a supported output mode.\n",
393  ns, solver->output_mode);
394  fprintf(stderr,"Should be \"serial\" or \"parallel\". \n");
395  return(1);
396 
397  }
398 
399  /* Solution plotting function */
400  strcpy(solver->plotfilename_extn,".png");
401 #ifdef with_python
402  solver->py_plt_func = NULL;
403  solver->py_plt_func_args = NULL;
404  {
405  char python_plotting_fname[_MAX_STRING_SIZE_] = "plotSolution";
406  PyObject* py_plot_name = PyUnicode_DecodeFSDefault(python_plotting_fname);
407  PyObject* py_plot_module = PyImport_Import(py_plot_name);
408  Py_DECREF(py_plot_name);
409  if (py_plot_module) {
410  solver->py_plt_func = PyObject_GetAttrString(py_plot_module, "plotSolution");
411  if (!solver->py_plt_func) {
412  if (!mpi->rank) {
413  printf("Unable to load plotSolution function from Python module.\n");
414  }
415  } else {
416  if (!mpi->rank) {
417  printf("Loaded Python module for plotting.\n");
418  printf("Loaded plotSolution function from Python module.\n");
419  }
420  }
421  } else {
422  if (!mpi->rank) {
423  printf("Unable to load Python module for plotting.\n");
424  }
425  }
426  }
427 #endif
428 
429  return 0;
430 }
int WriteBinary(int, int, int *, double *, double *, char *, int *)
Definition: WriteBinary.c:34
int(* SecondDerivativePar)(double *, double *, int, void *, void *)
Definition: hypar.h:247
int WriteTecplot3D(int, int, int *, double *, double *, char *, int *)
char * filename_index
Definition: hypar.h:197
int(* ParabolicFunction)(double *, double *, void *, void *, double)
Definition: hypar.h:256
void * time_integrator
Definition: hypar.h:165
char plotfilename_extn[_MAX_STRING_SIZE_]
Definition: hypar.h:203
char output_mode[_MAX_STRING_SIZE_]
Definition: hypar.h:183
void * py_plt_func
Definition: hypar.h:466
int restart_iter
Definition: hypar.h:58
int WriteText(int, int, int *, double *, double *, char *, int *)
Definition: WriteText.c:27
int(* SetInterpLimiterVar)(double *, double *, double *, int, void *, void *)
Definition: hypar.h:234
void IncrementFilenameIndex(char *, int)
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
#define _MAX_STRING_SIZE_
Definition: basic.h:14
char op_overwrite[_MAX_STRING_SIZE_]
Definition: hypar.h:191
int index_length
Definition: hypar.h:199
void * lusolver
Definition: hypar.h:368
char solnfilename_extn[_MAX_STRING_SIZE_]
Definition: hypar.h:201
void * interp
Definition: hypar.h:362
int WriteTecplot2D(int, int, int *, double *, double *, char *, int *)
int(* FirstDerivativePar)(double *, double *, int, int, void *, void *)
Definition: hypar.h:243
int(* InterpolateInterfacesPar)(double *, double *, int, void *, void *)
Definition: hypar.h:239
int flag_nonlinearinterp
Definition: hypar.h:411
int(* WriteOutput)(int, int, int *, double *, double *, char *, int *)
Definition: hypar.h:211
Structure of MPI-related variables.
void * msti
Definition: hypar.h:366
char op_file_format[_MAX_STRING_SIZE_]
Definition: hypar.h:186
void * py_plt_func_args
Definition: hypar.h:467
void * compact
Definition: hypar.h:364
int file_op_iter
Definition: hypar.h:171

◆ CleanupBarebones()

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
simsimulation objects of type SimulationObject

Definition at line 146 of file SparseGridsMiscFuncs.cpp.

147 {
148  if (sim->is_barebones == 0) {
149  fprintf(stderr, "Error in SparseGridsSimulation::CleanupBarebones()\n");
150  fprintf(stderr, " Simulation object is not a barebones one.\n");
151  return 1;
152  }
153 
154  HyPar* solver = &(sim->solver);
155  MPIVariables* mpi = &(sim->mpi);
156 
157  /* Free the communicators created */
158  IERR MPIFreeCommunicators(solver->ndims,mpi); CHECKERR(ierr);
159 
160  /* These variables are allocated in Initialize.c */
161  free(solver->dim_global);
162  free(solver->dim_global_ex);
163  free(solver->dim_local);
164  free(solver->index);
165  free(solver->isPeriodic);
166  free(solver->u);
167  free(solver->x);
168  free(solver->dxinv);
169  free(mpi->iproc);
170  free(mpi->ip);
171  free(mpi->is);
172  free(mpi->ie);
173  free(mpi->bcperiodic);
174  free(mpi->sendbuf);
175  free(mpi->recvbuf);
176  free(solver->VolumeIntegral);
177  free(solver->VolumeIntegralInitial);
178  free(solver->TotalBoundaryIntegral);
179  free(solver->ConservationError);
180  free(solver->stride_with_ghosts);
181  free(solver->stride_without_ghosts);
182  if (solver->filename_index) free(solver->filename_index);
183 
184  return(0);
185 }
double * recvbuf
#define IERR
Definition: basic.h:16
char * filename_index
Definition: hypar.h:197
#define CHECKERR(ierr)
Definition: basic.h:18
int * dim_global_ex
Definition: hypar.h:75
double * u
Definition: hypar.h:116
double * ConservationError
Definition: hypar.h:374
double * x
Definition: hypar.h:107
int * stride_with_ghosts
Definition: hypar.h:414
int ndims
Definition: hypar.h:26
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
int * stride_without_ghosts
Definition: hypar.h:416
double * sendbuf
int MPIFreeCommunicators(int, void *)
int * index
Definition: hypar.h:102
int * dim_local
Definition: hypar.h:37
double * TotalBoundaryIntegral
Definition: hypar.h:386
Structure of MPI-related variables.
int * dim_global
Definition: hypar.h:33
int * isPeriodic
Definition: hypar.h:162
double * VolumeIntegral
Definition: hypar.h:378
double * VolumeIntegralInitial
Definition: hypar.h:380
double * dxinv
Definition: hypar.h:110

◆ SetSolverParameters()

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_simSimulation object
a_dim_globalSpecified global grid sizes
a_iprocSpecified processor distibution
a_src_simSource simulation object
a_idxIndex or serial number
a_nsimsTotal number of simulations

Definition at line 11 of file SparseGridsSetSolverParameters.cpp.

18 {
19  a_dst_sim.solver.my_idx = a_idx;
20  a_dst_sim.solver.nsims = a_nsims;
21 
22  a_dst_sim.mpi.rank = m_rank;
23  a_dst_sim.mpi.nproc = m_nproc;
24 
25  a_dst_sim.solver.ndims = a_src_sim.solver.ndims;
26  a_dst_sim.solver.nvars = a_src_sim.solver.nvars;
27  a_dst_sim.solver.ghosts = a_src_sim.solver.ghosts;
28 
29  a_dst_sim.solver.dim_global = (int*) calloc (a_dst_sim.solver.ndims,sizeof(int));
30  a_dst_sim.mpi.iproc = (int*) calloc (a_dst_sim.solver.ndims,sizeof(int));
31  for (int d = 0; d < a_dst_sim.solver.ndims; d++) {
32  a_dst_sim.solver.dim_global[d] = a_dim_global[d];
33  a_dst_sim.mpi.iproc[d] = a_iproc[d];
34  }
35 
36  a_dst_sim.solver.n_iter = a_src_sim.solver.n_iter;
37  a_dst_sim.solver.restart_iter = a_src_sim.solver.restart_iter;
38 
39  strcpy(a_dst_sim.solver.time_scheme, a_src_sim.solver.time_scheme);
40  strcpy(a_dst_sim.solver.time_scheme_type, a_src_sim.solver.time_scheme_type);
41  strcpy(a_dst_sim.solver.spatial_scheme_hyp, a_src_sim.solver.spatial_scheme_hyp);
42  strcpy(a_dst_sim.solver.SplitHyperbolicFlux, a_src_sim.solver.SplitHyperbolicFlux);
43  strcpy(a_dst_sim.solver.interp_type, a_src_sim.solver.interp_type);
44  strcpy(a_dst_sim.solver.spatial_type_par, a_src_sim.solver.spatial_type_par);
45  strcpy(a_dst_sim.solver.spatial_scheme_par, a_src_sim.solver.spatial_scheme_par);
46 
47  a_dst_sim.solver.dt = a_src_sim.solver.dt;
48 
49  strcpy(a_dst_sim.solver.ConservationCheck, a_src_sim.solver.ConservationCheck);
50 
51  a_dst_sim.solver.screen_op_iter = a_src_sim.solver.screen_op_iter;
52  a_dst_sim.solver.file_op_iter = a_src_sim.solver.file_op_iter;
53 
54  strcpy(a_dst_sim.solver.op_file_format, a_src_sim.solver.op_file_format);
55  strcpy(a_dst_sim.solver.ip_file_type, a_src_sim.solver.ip_file_type);
56 
57  strcpy(a_dst_sim.solver.input_mode, a_src_sim.solver.input_mode);
58  strcpy(a_dst_sim.solver.output_mode, a_src_sim.solver.output_mode);
59  a_dst_sim.mpi.N_IORanks = a_src_sim.mpi.N_IORanks;
60 
61  strcpy(a_dst_sim.solver.op_overwrite, a_src_sim.solver.op_overwrite);
62  strcpy(a_dst_sim.solver.plot_solution, a_src_sim.solver.plot_solution);
63  strcpy(a_dst_sim.solver.model, a_src_sim.solver.model);
64  strcpy(a_dst_sim.solver.ib_filename, a_src_sim.solver.ib_filename);
65 
66  a_dst_sim.solver.flag_ib = a_src_sim.solver.flag_ib;
67 
68 #ifdef with_petsc
69  a_dst_sim.solver.use_petscTS = a_src_sim.solver.use_petscTS;
70 #endif
71 
72  return(0);
73 }
char plot_solution[_MAX_STRING_SIZE_]
Definition: hypar.h:194
char ib_filename[_MAX_STRING_SIZE_]
Definition: hypar.h:439
int nvars
Definition: hypar.h:29
char SplitHyperbolicFlux[_MAX_STRING_SIZE_]
Definition: hypar.h:92
int nsims
Definition: hypar.h:64
char ip_file_type[_MAX_STRING_SIZE_]
Definition: hypar.h:180
int n_iter
Definition: hypar.h:55
int flag_ib
Definition: hypar.h:441
char output_mode[_MAX_STRING_SIZE_]
Definition: hypar.h:183
char spatial_scheme_par[_MAX_STRING_SIZE_]
Definition: hypar.h:99
int restart_iter
Definition: hypar.h:58
int use_petscTS
Definition: hypar.h:395
char input_mode[_MAX_STRING_SIZE_]
Definition: hypar.h:177
int ndims
Definition: hypar.h:26
char ConservationCheck[_MAX_STRING_SIZE_]
Definition: hypar.h:376
char spatial_scheme_hyp[_MAX_STRING_SIZE_]
Definition: hypar.h:84
char time_scheme[_MAX_STRING_SIZE_]
Definition: hypar.h:78
char op_overwrite[_MAX_STRING_SIZE_]
Definition: hypar.h:191
char interp_type[_MAX_STRING_SIZE_]
Definition: hypar.h:88
char spatial_type_par[_MAX_STRING_SIZE_]
Definition: hypar.h:96
int my_idx
Definition: hypar.h:61
double dt
Definition: hypar.h:67
char model[_MAX_STRING_SIZE_]
Definition: hypar.h:263
int ghosts
Definition: hypar.h:52
char time_scheme_type[_MAX_STRING_SIZE_]
Definition: hypar.h:81
char op_file_format[_MAX_STRING_SIZE_]
Definition: hypar.h:186
int * dim_global
Definition: hypar.h:33
int file_op_iter
Definition: hypar.h:171
int screen_op_iter
Definition: hypar.h:168

◆ OutputSolution()

void OutputSolution ( double  a_time)
protected

Output solutions to file

Write solutions to file

Parameters
a_timesimulation time

Definition at line 11 of file SparseGridsOutputSolution.cpp.

12 {
13  /* if asked for, write individual sparse grids solutions */
14  if (m_write_sg_solutions == 1) {
15  for (int ns = 0; ns < m_nsims_sg; ns++) {
16  if (m_sims_sg[ns].solver.PhysicsOutput) {
17  m_sims_sg[ns].solver.PhysicsOutput( &(m_sims_sg[ns].solver),
18  &(m_sims_sg[ns].mpi),
19  a_time );
20  }
21  }
22  ::OutputSolution((void*)m_sims_sg.data(), m_nsims_sg, a_time);
23  }
24 
25  /* Combine the sparse grids solutions to full grid */
27 
28  /* Write the full grid solution */
29  ::OutputSolution((void*)m_sim_fg, 1, a_time);
30 
31  return;
32 }
void CombinationTechnique(SimulationObject *const)
std::vector< SimulationObject > m_sims_sg

◆ CombinationTechnique()

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_simtarget simulation object on which to combine

Definition at line 23 of file SparseGridsCombinationTechnique.cpp.

24 {
25  double** u_sg = (double**) calloc (m_nsims_sg, sizeof(double*));
26  std::vector<double> coeffs(m_nsims_sg, 0.0);
27  for (int n=0; n<m_nsims_sg; n++) {
28  u_sg[n] = m_sims_sg[n].solver.u;
29  coeffs[n] = m_combination[n]._coeff_;
30  }
31 
33  u_sg,
34  m_nsims_sg,
35  a_sim,
36  a_sim->solver.u,
37  coeffs.data() );
38 
39  /* done */
40  return;
41 }
void CombineSolutions(SimulationObject *, double *const *const, const int, SimulationObject *, double *const, const double *const)
double * u
Definition: hypar.h:116
std::vector< SGCTElem > m_combination
std::vector< SimulationObject > m_sims_sg

◆ interpolate() [1/2]

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_dstDestination object
a_srcSource object

Definition at line 15 of file SparseGridsInterpolate.cpp.

18 {
19  /* get the destination grid dimensions */
20  GridDimensions dim_dst;
21  StdVecOps::copyFrom(dim_dst, a_dst->solver.dim_global, m_ndims);
22 
23  /* get number of vector components */
24  int nvars = a_src->solver.nvars;
25  if (nvars != a_dst->solver.nvars) {
26  fprintf(stderr, "Error in SparseGridsSimulation::interpolate(): unequal nvars\n");
27  exit(1);
28  }
29 
30  double *ug_dst = NULL;
31  interpolate(dim_dst, &ug_dst, a_src);
32 
33  /* partition the destination */
35  (void*) &(a_dst->mpi),
36  (a_dst->mpi.rank ? NULL : ug_dst),
37  a_dst->solver.u,
38  a_dst->solver.dim_global,
39  a_dst->solver.dim_local,
40  a_dst->solver.ghosts,
41  nvars);
42 
43  if (!m_rank) {
44  free(ug_dst);
45  }
46 
47  return;
48 }
int MPIPartitionArraynDwGhosts(int, void *, double *, double *, int *, int *, int, int)
int nvars
Definition: hypar.h:29
std::vector< int > GridDimensions
double * u
Definition: hypar.h:116
int * dim_local
Definition: hypar.h:37
void interpolate(SimulationObject *const, const SimulationObject *const)
int ghosts
Definition: hypar.h:52
int * dim_global
Definition: hypar.h:33
void copyFrom(std::vector< int > &a_iv, const int *const a_iv_carr, int a_n)
Definition: std_vec_ops.h:36

◆ interpolate() [2/2]

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_dstgrid dimensions to interpolate to
a_u_dstpointer to array containing interpolated data
a_srcSource object

Definition at line 63 of file SparseGridsInterpolate.cpp.

67 {
68  if ((*a_u_dst) != NULL) {
69  fprintf(stderr, "Error in SparseGridsSimulation::interpolate() - \n");
70  fprintf(stderr, " a_u_dst is not NULL!\n");
71  exit(1);
72  }
73 
74  /* get the source grid dimensions */
75  GridDimensions dim_src;
76  StdVecOps::copyFrom(dim_src, a_src->solver.dim_global, m_ndims);
77 
78  /* get number of vector components and number of ghost points*/
79  int nvars = a_src->solver.nvars;
80  int ghosts = a_src->solver.ghosts;
81 
82  /* gather the source on rank 0 */
83  double *ug_src = NULL;
84  if (!m_rank) allocateDataArrays(dim_src, nvars, &ug_src, ghosts);
86  (void*) &(a_src->mpi),
87  ug_src,
88  a_src->solver.u,
89  a_src->solver.dim_global,
90  a_src->solver.dim_local,
91  ghosts,
92  nvars );
93 
94  std::vector<int> periodic_arr(m_ndims);
95  for (int i=0; i<m_ndims; i++) {
96  periodic_arr[i] = (m_is_periodic[i] ? 1 : 0);
97  }
98 
99  if (!m_rank) {
100  int ierr = ::InterpolateGlobalnDVar( a_dim_dst.data(),
101  a_u_dst,
102  dim_src.data(),
103  ug_src,
104  nvars,
105  ghosts,
106  m_ndims,
107  periodic_arr.data());
108  if (ierr) {
109  fprintf(stderr,"InterpolateGlobalnDVar() returned with error!\n");
110  exit(1);
111  }
112  }
113 
114  return;
115 }
int nvars
Definition: hypar.h:29
std::vector< int > GridDimensions
double * u
Definition: hypar.h:116
void allocateDataArrays(const GridDimensions &a_dim, const int a_nvars, double **const a_u, const int a_ngpt=0)
int MPIGatherArraynDwGhosts(int, void *, double *, double *, int *, int *, int, int)
int InterpolateGlobalnDVar(const int *const, double **const, const int *const, double *const, const int, const int, const int, const int *const)
std::vector< bool > m_is_periodic
int * dim_local
Definition: hypar.h:37
int ghosts
Definition: hypar.h:52
int * dim_global
Definition: hypar.h:33
void copyFrom(std::vector< int > &a_iv, const int *const a_iv_carr, int a_n)
Definition: std_vec_ops.h:36

◆ interpolateGrid()

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_dstDestination object
a_srcSource object

Definition at line 20 of file SparseGridsInterpolateGrid.cpp.

23 {
24  /* get the source and destination grid dimensions */
25  GridDimensions dim_src, dim_dst;
26  StdVecOps::copyFrom(dim_src, a_src->solver.dim_global, m_ndims);
27  StdVecOps::copyFrom(dim_dst, a_dst->solver.dim_global, m_ndims);
28 
29  /* gather the source on rank 0 */
30  double* xg_src = NULL;
31  if (!m_rank) allocateGridArrays(dim_src, &xg_src);
32  {
33  int offset_global, offset_local;
34  offset_global = offset_local = 0;
35  for (int d=0; d<m_ndims; d++) {
36  MPIGatherArray1D( (void*) &(a_src->mpi),
37  (a_src->mpi.rank ? NULL : &xg_src[offset_global]),
38  &(a_src->solver.x[offset_local+a_src->solver.ghosts]),
39  a_src->mpi.is[d],
40  a_src->mpi.ie[d],
41  a_src->solver.dim_local[d],
42  0 );
43  offset_global += a_src->solver.dim_global[d];
44  offset_local += a_src->solver.dim_local [d] + 2*a_src->solver.ghosts;
45  }
46  }
47 
48  /* now do the interpolation, dimension-by-dimension */
49  double *xg_dst = NULL;
50  if (!m_rank) {
51 
52  GridDimensions dim_to(m_ndims,0);
53  GridDimensions dim_from(m_ndims,0);
54 
55  double* x_from;
56  double* x_to;
57 
58  dim_to = dim_src;
59  x_to = xg_src;
60  x_from = NULL;
61 
62  for (int dir = 0; dir < m_ndims; dir++) {
63 
64  dim_from = dim_to;
65  dim_to[dir] = dim_dst[dir];
66 
67  if (dim_from[dir] == dim_to[dir]) continue;
68 
69  double fac = (dim_to[dir] > dim_from[dir] ?
70  (double)dim_to[dir]/(double)dim_from[dir]
71  : (double)dim_from[dir]/(double)dim_to[dir] );
72  if (!isPowerOfTwo((int)fac)) {
73  fprintf(stderr,"Error in SparseGridsSimulation::interpolate() - \n");
74  fprintf(stderr," refinement/coarsening factor not a power of 2!\n");
75  exit(1);
76  }
77 
78  if (x_from != NULL) free(x_from);
79  x_from = x_to;
80 
81  allocateGridArrays(dim_to, &x_to);
82  if (dim_to[dir] < dim_from[dir]) {
83  coarsenGrid1D(dim_from, dim_to, x_from, x_to, dir);
84  } else {
85  refineGrid1D(dim_from, dim_to, x_from, x_to, dir);
86  }
87 
88  }
89 
90  /* dim_to should be equal to dim_dst now */
91  for (int d = 0; d < m_ndims; d++) {
92  if (dim_to[d] != dim_dst[d]) {
93  fprintf(stderr,"Error in SparseGridsSimulation::interpolate() - \n");
94  fprintf(stderr," dim_to[%d] != dim_dst[%d]!\n", d, d);
95  exit(1);
96  }
97  }
98 
99  if (x_from != NULL) free(x_from);
100  xg_dst = x_to;
101 
102  }
103 
104  /* partition the destination */
105  int offset_global = 0, offset_local = 0;
106  for (int d=0; d<m_ndims; d++) {
107  MPIPartitionArray1D( (void*) &(a_dst->mpi),
108  (a_dst->mpi.rank ? NULL : &xg_dst[offset_global]),
109  &(a_dst->solver.x[offset_local+a_dst->solver.ghosts]),
110  a_dst->mpi.is[d],
111  a_dst->mpi.ie[d],
112  a_dst->solver.dim_local[d],
113  0);
114  offset_global += a_dst->solver.dim_global[d];
115  offset_local += a_dst->solver.dim_local [d] + 2*a_dst->solver.ghosts;
116  }
117 
118  if (!m_rank) {
119  free(xg_dst);
120  }
121 
122  /* exchange MPI-boundary values of x between processors */
123  {
124  int offset = 0;
125  for (int d = 0; d < m_ndims; d++) {
126  MPIExchangeBoundaries1D( (void*) &(a_dst->mpi),
127  &(a_dst->solver.x[offset]),
128  a_dst->solver.dim_local[d],
129  a_dst->solver.ghosts,
130  d,
131  m_ndims);
132  offset += (a_dst->solver.dim_local[d] + 2*a_dst->solver.ghosts);
133  }
134  }
135  /* fill in ghost values of x at physical boundaries by extrapolation */
136  {
137  int offset = 0;
138  for (int d = 0; d < m_ndims; d++) {
139  double* X = &(a_dst->solver.x[offset]);
140  int* dim = a_dst->solver.dim_local;
141  int ghosts = a_dst->solver.ghosts;
142  if (a_dst->mpi.ip[d] == 0) {
143  /* fill left boundary along this dimension */
144  for (int i = 0; i < ghosts; i++) {
145  int delta = ghosts - i;
146  X[i] = X[ghosts] + ((double) delta) * (X[ghosts]-X[ghosts+1]);
147  }
148  }
149  if (a_dst->mpi.ip[d] == a_dst->mpi.iproc[d]-1) {
150  /* fill right boundary along this dimension */
151  for (int i = dim[d]+ghosts; i < dim[d]+2*ghosts; i++) {
152  int delta = i - (dim[d]+ghosts-1);
153  X[i] = X[dim[d]+ghosts-1]
154  + ((double) delta) * (X[dim[d]+ghosts-1]-X[dim[d]+ghosts-2]);
155  }
156  }
157  offset += (dim[d] + 2*ghosts);
158  }
159  }
160  return;
161 }
void refineGrid1D(const GridDimensions &, const GridDimensions &, const double *const, double *const, int)
int MPIPartitionArray1D(void *, double *, double *, int, int, int, int)
std::vector< int > GridDimensions
double * x
Definition: hypar.h:107
int MPIExchangeBoundaries1D(void *, double *, int, int, int, int)
int MPIGatherArray1D(void *, double *, double *, int, int, int, int)
void allocateGridArrays(const GridDimensions &a_dim, double **const a_x, const int a_ngpt=0)
int * dim_local
Definition: hypar.h:37
int ghosts
Definition: hypar.h:52
void coarsenGrid1D(const GridDimensions &, const GridDimensions &, const double *const, double *const, int)
int * dim_global
Definition: hypar.h:33
void copyFrom(std::vector< int > &a_iv, const int *const a_iv_carr, int a_n)
Definition: std_vec_ops.h:36

◆ coarsenGrid1D()

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_srcGrid size of source data
a_dim_dstGrid size of destination data
a_x_srcSource grid
a_x_dstDestination grid
a_dirDimension along which to coarsen

Definition at line 168 of file SparseGridsInterpolateGrid.cpp.

174 {
175  for (int d = 0; d < m_ndims; d++) {
176  if ((d != a_dir) && (a_dim_src[d] != a_dim_dst[d])) {
177  fprintf(stderr, "Error in SparseGridsSimulation::coarsenGrid1D() -\n");
178  fprintf(stderr, " a_dim_src[%d] != a_dim_dst[%d]\n", d, d);
179  exit(1);
180  }
181  }
182 
183  int n_src = a_dim_src[a_dir];
184  int n_dst = a_dim_dst[a_dir];
185  if (n_dst > n_src) {
186  fprintf(stderr, "Error in SparseGridsSimulation::coarsenGrid1D() -\n");
187  fprintf(stderr, " destination grid is finer than source grid along a_dir!\n");
188  exit(1);
189  }
190 
191  double fac = ((double) n_src) / ((double) n_dst);
192  int stride = (int) fac;
193  if (std::abs(((double)stride)-fac) > _MACHINE_ZERO_) {
194  fprintf(stderr, "Error in SparseGridsSimulation::coarsenGrid1D() -\n");
195  fprintf(stderr, " non-integer coarsening factor!\n");
196  exit(1);
197  }
198 
199  const double* x_src = a_x_src;
200  double* x_dst = a_x_dst;
201 
202  for (int d = 0; d < m_ndims; d++) {
203 
204  if (d == a_dir) {
205 
206  for (int i = 0; i < a_dim_dst[d]; i++) {
207  double avg = 0;
208  for (int j=stride*i; j<stride*(i+1); j++) {
209  if (j >= n_src) {
210  fprintf(stderr, "Error in SparseGridsSimulation::coarsenGrid1D() -\n");
211  fprintf(stderr, " j >= n_src!\n");
212  exit(1);
213  }
214  avg += x_src[j];
215  }
216  avg /= (double) stride;
217  x_dst[i] = avg;
218  }
219 
220  } else {
221 
222  _ArrayCopy1D_(x_src, x_dst, a_dim_dst[d]);
223 
224  }
225 
226  x_src += a_dim_src[d];
227  x_dst += a_dim_dst[d];
228  }
229 
230  return;
231 }
#define _MACHINE_ZERO_
Definition: basic.h:26
#define _ArrayCopy1D_(x, y, size)

◆ refineGrid1D()

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_srcGrid size of source data
a_dim_dstGrid size of destination data
a_x_srcSource grid
a_x_dstDestination grid
a_dirDimension along which to refine

Definition at line 238 of file SparseGridsInterpolateGrid.cpp.

244 {
245  for (int d = 0; d < m_ndims; d++) {
246  if ((d != a_dir) && (a_dim_src[d] != a_dim_dst[d])) {
247  fprintf(stderr, "Error in SparseGridsSimulation::coarsenGrid1D() -\n");
248  fprintf(stderr, " a_dim_src[%d] != a_dim_dst[%d]\n", d, d);
249  exit(1);
250  }
251  }
252 
253  int n_src = a_dim_src[a_dir];
254  int n_dst = a_dim_dst[a_dir];
255  if (n_dst < n_src) {
256  fprintf(stderr, "Error in SparseGridsSimulation::refineGrid1D() -\n");
257  fprintf(stderr, " destination grid is coarser than source grid along a_dir!\n");
258  exit(1);
259  }
260 
261  double fac = ((double) n_dst) / ((double) n_src);
262  int stride = (int) fac;
263  if (std::abs(((double)stride)-fac) > _MACHINE_ZERO_) {
264  fprintf(stderr, "Error in SparseGridsSimulation::refineGrid1D() -\n");
265  fprintf(stderr, " non-integer refinement factor!\n");
266  exit(1);
267  }
268 
269  int offset = 0;
270  for (int d = 0; d < a_dir; d++) offset += a_dim_src[d];
271  const double* x_src = a_x_src + offset;
272  double* x_dst = a_x_dst + offset;
273 
274  fprintf(stderr, "Error in SparseGridsSimulation::refineGrid1D() -\n");
275  fprintf(stderr, " NOT YET IMPLEMENTED! Why do you need this?\n");
276  exit(1);
277 
278  return;
279 }
#define _MACHINE_ZERO_
Definition: basic.h:26

◆ fillGhostCells()

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_dimgrid dimensions of solution
a_ngptnumber of ghost cells
a_usolution array
a_nvarsnumber of vector components of the solution

Definition at line 18 of file SparseGridsFillGhostCells.cpp.

23 {
24  if (m_is_periodic.size() != m_ndims) {
25  fprintf(stderr, "Error in SparseGridsSimulation::fillGhostCells() -\n");
26  fprintf(stderr, "m_is_periodic.size() != m_ndims \n");
27  exit(1);
28  }
29 
30  for (int d = 0; d < m_ndims; d++) {
31 
32  int bounds[m_ndims];
33  _ArrayCopy1D_(a_dim, bounds, m_ndims);
34  bounds[d] = a_ngpt;
35 
36  int index[m_ndims];
37  _ArraySetValue_(index, m_ndims, 0);
38 
39  if (m_is_periodic[d]) {
40 
41  /* periodic case */
42 
43  int done = 0;
44  while (!done) {
45 
46  {
47  /* low end - face = 1 */
48 
49  int p_gpt = 0,
50  p_int = 0;
51 
52  int index_gpt[m_ndims];
53  _ArrayCopy1D_(index, index_gpt, m_ndims);
54  index_gpt[d] -= a_ngpt;
55  _ArrayIndex1D_(m_ndims, a_dim, index_gpt, a_ngpt, p_gpt);
56 
57  int index_int[m_ndims];
58  _ArrayCopy1D_(index, index_int, m_ndims);
59  index_int[d] += (a_dim[d]-a_ngpt);
60  _ArrayIndex1D_(m_ndims, a_dim, index_int, a_ngpt, p_int);
61 
62  _ArrayCopy1D_((a_u+a_nvars*p_int), (a_u+a_nvars*p_gpt), a_nvars);
63  }
64 
65  {
66  /* high end - face = -1 */
67 
68  int p_gpt = 0,
69  p_int = 0;
70 
71  int index_gpt[m_ndims];
72  _ArrayCopy1D_(index, index_gpt, m_ndims);
73  index_gpt[d] += a_dim[d];
74  _ArrayIndex1D_(m_ndims, a_dim, index_gpt, a_ngpt, p_gpt);
75 
76  int index_int[m_ndims];
77  _ArrayCopy1D_(index, index_int, m_ndims);
78  _ArrayIndex1D_(m_ndims, a_dim, index_int, a_ngpt, p_int);
79 
80  _ArrayCopy1D_((a_u+a_nvars*p_int), (a_u+a_nvars*p_gpt), a_nvars);
81  }
82 
83  _ArrayIncrementIndex_(m_ndims, bounds, index, done);
84 
85  }
86 
87  } else {
88 
89  /* not periodic - extrapolate */
90 
91  int done = 0;
92  while (!done) {
93 
94  {
95  /* low end - face = 1 */
96 
97  int p_gpt = 0,
98  p_int_0 = 0,
99  p_int_1 = 0,
100  p_int_2 = 0,
101  p_int_3 = 0;
102 
103  int index_gpt[m_ndims];
104  _ArrayCopy1D_(index, index_gpt, m_ndims);
105  index_gpt[d] -= a_ngpt;
106  _ArrayIndex1D_(m_ndims, a_dim, index_gpt, a_ngpt, p_gpt);
107 
108  int index_int[m_ndims];
109  _ArrayCopy1D_(index, index_int, m_ndims);
110 
111  index_int[d] = 0;
112  _ArrayIndex1D_(m_ndims, a_dim, index_int, a_ngpt, p_int_0);
113  index_int[d]++;
114  _ArrayIndex1D_(m_ndims, a_dim, index_int, a_ngpt, p_int_1);
115  index_int[d]++;
116  _ArrayIndex1D_(m_ndims, a_dim, index_int, a_ngpt, p_int_2);
117  index_int[d]++;
118  _ArrayIndex1D_(m_ndims, a_dim, index_int, a_ngpt, p_int_3);
119 
120  double alpha = - (double) (a_ngpt - index[d]);
121  double c0 = -((-2.0 + alpha)*(-1.0 + alpha)*alpha)/6.0;
122  double c1 = ((-2.0 + alpha)*(-1.0 + alpha)*(1.0 + alpha))/2.0;
123  double c2 = (alpha*(2.0 + alpha - alpha*alpha))/2.0;
124  double c3 = (alpha*(-1.0 + alpha*alpha))/6.0;
125 
126  for (int v = 0; v < a_nvars; v++) {
127 
128  a_u[p_gpt*a_nvars+v] = c0 * a_u[p_int_0*a_nvars+v]
129  + c1 * a_u[p_int_1*a_nvars+v]
130  + c2 * a_u[p_int_2*a_nvars+v]
131  + c3 * a_u[p_int_3*a_nvars+v];
132 
133  }
134 
135  }
136 
137  {
138  /* high end - face = -1 */
139 
140  int p_gpt = 0,
141  p_int_0 = 0,
142  p_int_1 = 0,
143  p_int_2 = 0,
144  p_int_3 = 0;
145 
146  int index_gpt[m_ndims];
147  _ArrayCopy1D_(index, index_gpt, m_ndims);
148  index_gpt[d] += a_dim[d];
149  _ArrayIndex1D_(m_ndims, a_dim, index_gpt, a_ngpt, p_gpt);
150 
151  int index_int[m_ndims];
152  _ArrayCopy1D_(index, index_int, m_ndims);
153 
154  index_int[d] = a_dim[d]-1;
155  _ArrayIndex1D_(m_ndims, a_dim, index, a_ngpt, p_int_0);
156  index_int[d]--;
157  _ArrayIndex1D_(m_ndims, a_dim, index, a_ngpt, p_int_1);
158  index_int[d]--;
159  _ArrayIndex1D_(m_ndims, a_dim, index, a_ngpt, p_int_2);
160  index_int[d]--;
161  _ArrayIndex1D_(m_ndims, a_dim, index, a_ngpt, p_int_3);
162 
163  double alpha = - (double) (index[d]+1);
164  double c0 = -((-2.0 + alpha)*(-1.0 + alpha)*alpha)/6.0;
165  double c1 = ((-2.0 + alpha)*(-1.0 + alpha)*(1.0 + alpha))/2.0;
166  double c2 = (alpha*(2.0 + alpha - alpha*alpha))/2.0;
167  double c3 = (alpha*(-1.0 + alpha*alpha))/6.0;
168 
169  for (int v = 0; v < a_nvars; v++) {
170 
171  a_u[p_gpt*a_nvars+v] = c0 * a_u[p_int_0*a_nvars+v]
172  + c1 * a_u[p_int_1*a_nvars+v]
173  + c2 * a_u[p_int_2*a_nvars+v]
174  + c3 * a_u[p_int_3*a_nvars+v];
175 
176  }
177 
178  }
179 
180  _ArrayIncrementIndex_(m_ndims, bounds, index, done);
181 
182  }
183 
184  }
185 
186  }
187 
188  return;
189 }
#define _ArrayIndex1D_(N, imax, i, ghost, index)
std::vector< bool > m_is_periodic
#define _ArraySetValue_(x, size, value)
#define _ArrayIncrementIndex_(N, imax, i, done)
#define _ArrayCopy1D_(x, y, size)

◆ CalculateError()

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.

18 {
19  int exact_flag;
20  double *uex = NULL;
21 
22  /* read in full grid exact solution, if available */
23  {
24  HyPar* solver = &(m_sim_fg->solver);
25  MPIVariables* mpi = &(m_sim_fg->mpi);
26  long size = solver->nvars * solver->npoints_local_wghosts;
27  uex = (double*) calloc (size, sizeof(double));
28 
29  char fname_root[_MAX_STRING_SIZE_] = "exact";
30  IERR ExactSolution( solver,
31  mpi,
32  uex,
33  fname_root,
34  &exact_flag ); CHECKERR(ierr);
35  }
36 
37  for (int n=0; n<m_nsims_sg; n++) {
38 // TimeError( &(m_sims_sg[n].solver),
39 // &(m_sims_sg[n].mpi),
40 // uex );
41  }
42 
43  if (!exact_flag) {
44 
45  /* No exact solution available */
46  for (int n=0; n<m_nsims_sg; n++) {
47  m_sims_sg[n].solver.error[0]
48  = m_sims_sg[n].solver.error[1]
49  = m_sims_sg[n].solver.error[2]
50  = -1;
51  }
52  m_sim_fg->solver.error[0]
53  = m_sim_fg->solver.error[1]
54  = m_sim_fg->solver.error[2]
55  = -1;
56 
57  } else {
58 
59  std::vector<int> periodic_arr(m_ndims);
60  for (int i=0; i<m_ndims; i++) {
61  periodic_arr[i] = (m_is_periodic[i] ? 1 : 0);
62  }
63 
64  double *uex2 = NULL;
65 
66  if (m_print_sg_errors == 1) {
67  long size = m_sim_fg->solver.nvars
69  uex2 = (double*) calloc(size, sizeof(double));
70  _ArrayCopy1D_(uex, uex2, size);
71  }
72 
73  /* calculate error for full grid */
74  computeError( *m_sim_fg, uex );
75 
76  /* calculate error for sparse grids */
77  if (m_print_sg_errors == 1) {
78  for (int n = 0; n < m_nsims_sg; n++) {
79 
80  GridDimensions dim_fg(m_ndims,0);
82 
83  GridDimensions dim_sg(m_ndims,0);
84  StdVecOps::copyFrom(dim_sg, m_sims_sg[n].solver.dim_global, m_ndims);
85 
86  /* assemble the global exact solution on full grid */
87  double *uex_global_fg = NULL;
88  if (!m_rank) {
89  allocateDataArrays( dim_fg,
91  &uex_global_fg,
93  }
95  (void*) &(m_sim_fg->mpi),
96  uex_global_fg,
97  uex2,
101  m_sim_fg->solver.nvars );
102 
103  /* interpolate to sparse grid -
104  * this will delete the full grid array*/
105  double *uex_global_sg = NULL;
106  if (!m_rank) {
107  int ierr = ::InterpolateGlobalnDVar( dim_sg.data(),
108  &uex_global_sg,
109  dim_fg.data(),
110  uex_global_fg,
111  m_sims_sg[n].solver.nvars,
112  m_sims_sg[n].solver.ghosts,
113  m_ndims,
114  periodic_arr.data() );
115  if (ierr) {
116  fprintf(stderr,"InterpolateGlobalnDVar() returned with error!\n");
117  exit(1);
118  }
119  }
120 
121  /* allocate local exact solution on this sparse grid */
122  long size = m_sims_sg[n].solver.nvars
123  * m_sims_sg[n].solver.npoints_local_wghosts;
124  double* uex_sg = (double*) calloc(size, sizeof(double));
125 
126  /* partition the global exact solution to local on this sparse grid */
128  (void*) &(m_sims_sg[n].mpi),
129  (m_rank ? NULL : uex_global_sg),
130  uex_sg,
131  m_sims_sg[n].solver.dim_global,
132  m_sims_sg[n].solver.dim_local,
133  m_sims_sg[n].solver.ghosts,
134  m_sims_sg[n].solver.nvars );
135 
136  /* delete the global exact solution array */
137  if (!m_rank) free(uex_global_sg);
138 
139  /* compute the error */
140  computeError( m_sims_sg[n], uex_sg);
141 
142  /* delete the exact solution array */
143  free(uex_sg);
144  }
145  free(uex2);
146  }
147 
148  }
149 
150  free(uex);
151  return;
152 }
int MPIPartitionArraynDwGhosts(int, void *, double *, double *, int *, int *, int, int)
int nvars
Definition: hypar.h:29
#define IERR
Definition: basic.h:16
#define CHECKERR(ierr)
Definition: basic.h:18
std::vector< int > GridDimensions
double error[3]
Definition: hypar.h:371
void allocateDataArrays(const GridDimensions &a_dim, const int a_nvars, double **const a_u, const int a_ngpt=0)
int MPIGatherArraynDwGhosts(int, void *, double *, double *, int *, int *, int, int)
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
#define _MAX_STRING_SIZE_
Definition: basic.h:14
int InterpolateGlobalnDVar(const int *const, double **const, const int *const, double *const, const int, const int, const int, const int *const)
std::vector< SimulationObject > m_sims_sg
std::vector< bool > m_is_periodic
void computeError(SimulationObject &, double *const)
int * dim_local
Definition: hypar.h:37
int ExactSolution(void *, void *, double *, char *, int *)
Definition: ExactSolution.c:16
int ghosts
Definition: hypar.h:52
Structure of MPI-related variables.
int npoints_local_wghosts
Definition: hypar.h:42
#define _ArrayCopy1D_(x, y, size)
int * dim_global
Definition: hypar.h:33
void copyFrom(std::vector< int > &a_iv, const int *const a_iv_carr, int a_n)
Definition: std_vec_ops.h:36

◆ computeError()

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_simSimulation object
a_uexExact solution

Definition at line 155 of file SparseGridsCalculateError.cpp.

157 {
158  static const double tolerance = 1e-15;
159 
160  HyPar* solver = &(a_sim.solver);
161  MPIVariables* mpi = &(a_sim.mpi);
162 
163  if (a_uex == NULL) {
164  fprintf(stderr, "Error in SparseGrids::computeError() -\n");
165  fprintf(stderr, " exact solution pointer in NULL.\n");
166  exit(1);
167  }
168 
169  /* calculate solution norms (for relative error) */
170  double sum, global_sum;
171  double solution_norm[3] = {0.0,0.0,0.0};
172  /* L1 */
173  sum = ArraySumAbsnD (solver->nvars,solver->ndims,solver->dim_local,
174  solver->ghosts,solver->index,a_uex);
175  global_sum = 0; MPISum_double(&global_sum,&sum,1,&mpi->world);
176  solution_norm[0] = global_sum/((double)solver->npoints_global);
177  /* L2 */
178  sum = ArraySumSquarenD(solver->nvars,solver->ndims,solver->dim_local,
179  solver->ghosts,solver->index,a_uex);
180  global_sum = 0; MPISum_double(&global_sum,&sum,1,&mpi->world);
181  solution_norm[1] = sqrt(global_sum/((double)solver->npoints_global));
182  /* Linf */
183  sum = ArrayMaxnD (solver->nvars,solver->ndims,solver->dim_local,
184  solver->ghosts,solver->index,a_uex);
185  global_sum = 0; MPIMax_double(&global_sum,&sum,1,&mpi->world);
186  solution_norm[2] = global_sum;
187 
188  /* compute error = difference between exact and numerical solution */
189  long size = solver->nvars*solver->npoints_local_wghosts;
190  _ArrayAXPY_(solver->u,-1.0,a_uex,size);
191 
192  /* calculate L1 norm of error */
193  sum = ArraySumAbsnD (solver->nvars,solver->ndims,solver->dim_local,
194  solver->ghosts,solver->index,a_uex);
195  global_sum = 0; MPISum_double(&global_sum,&sum,1,&mpi->world);
196  solver->error[0] = global_sum/((double)solver->npoints_global);
197 
198  /* calculate L2 norm of error */
199  sum = ArraySumSquarenD(solver->nvars,solver->ndims,solver->dim_local,
200  solver->ghosts,solver->index,a_uex);
201  global_sum = 0; MPISum_double(&global_sum,&sum,1,&mpi->world);
202  solver->error[1] = sqrt(global_sum/((double)solver->npoints_global));
203 
204  /* calculate Linf norm of error */
205  sum = ArrayMaxnD (solver->nvars,solver->ndims,solver->dim_local,
206  solver->ghosts,solver->index,a_uex);
207  global_sum = 0; MPIMax_double(&global_sum,&sum,1,&mpi->world);
208  solver->error[2] = global_sum;
209 
210  /*
211  decide whether to normalize and report relative errors,
212  or report absolute errors.
213  */
214  if ( (solution_norm[0] > tolerance)
215  && (solution_norm[1] > tolerance)
216  && (solution_norm[2] > tolerance) ) {
217  solver->error[0] /= solution_norm[0];
218  solver->error[1] /= solution_norm[1];
219  solver->error[2] /= solution_norm[2];
220  }
221 
222  return;
223 }
int nvars
Definition: hypar.h:29
int MPISum_double(double *, double *, int, void *)
Definition: MPISum.c:39
double error[3]
Definition: hypar.h:371
long sum(const std::vector< int > &a_iv)
Definition: std_vec_ops.h:18
double * u
Definition: hypar.h:116
int ndims
Definition: hypar.h:26
INLINE double ArraySumAbsnD(int, int, int *, int, int *, double *)
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
int MPIMax_double(double *, double *, int, void *)
Definition: MPIMax.c:38
#define _ArrayAXPY_(x, a, y, size)
int npoints_global
Definition: hypar.h:40
MPI_Comm world
int * index
Definition: hypar.h:102
int * dim_local
Definition: hypar.h:37
INLINE double ArrayMaxnD(int, int, int *, int, int *, double *)
int ghosts
Definition: hypar.h:52
Structure of MPI-related variables.
int npoints_local_wghosts
Definition: hypar.h:42
INLINE double ArraySumSquarenD(int, int, int *, int, int *, double *)

◆ isPowerOfTwo()

bool isPowerOfTwo ( int  x)
inlineprotected

Checks if an integer is a power of 2

Definition at line 317 of file sparse_grids_simulation.h.

318  {
319  if (x == 0) return false;
320 
321  while (x > 1) {
322  if (x%2 != 0) return false;
323  x /= 2;
324  }
325  return true;
326  }

◆ cfunc()

double cfunc ( int  a,
int  b 
)
inlineprotected

The combination/choose function from probability stuff

Definition at line 329 of file sparse_grids_simulation.h.

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

◆ factorial()

int factorial ( int  a)
inlineprotected

Factorial function

Definition at line 335 of file sparse_grids_simulation.h.

336  {
337  int retval = 1.0;
338  for (int i=1; i<=a; i++) retval *= i;
339  return retval;
340  }

◆ allocateGridArrays()

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.

346  {
347  GridDimensions dim_wghosts = a_dim;
348  for (int i=0; i<dim_wghosts.size(); i++) {
349  dim_wghosts[i] += (2*a_ngpt);
350  }
351  long size_x = StdVecOps::sum(dim_wghosts);
352  (*a_x) = (double*) calloc (size_x, sizeof(double));
353  for (int i=0; i<size_x; i++) (*a_x)[i] = 0.0;
354  return;
355  }
std::vector< int > GridDimensions
long sum(const std::vector< int > &a_iv)
Definition: std_vec_ops.h:18

◆ allocateDataArrays()

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.

362  {
363  GridDimensions dim_wghosts = a_dim;
364  for (int i=0; i<dim_wghosts.size(); i++) {
365  dim_wghosts[i] += (2*a_ngpt);
366  }
367  long size_u = a_nvars*StdVecOps::product(dim_wghosts);
368  (*a_u) = (double*) calloc (size_u, sizeof(double));
369  for (int i=0; i<size_u; i++) (*a_u)[i] = 0.0;
370  return;
371  }
std::vector< int > GridDimensions
long product(const std::vector< int > &a_iv)
Definition: std_vec_ops.h:27

Field Documentation

◆ m_is_defined

bool m_is_defined
protected

Boolean to show if this object is defined

Definition at line 212 of file sparse_grids_simulation.h.

◆ m_nsims_sg

int m_nsims_sg
protected

Number of sparse grids simulations

Definition at line 214 of file sparse_grids_simulation.h.

◆ m_ndims

int m_ndims
protected

Number of spatial dimensions

Definition at line 215 of file sparse_grids_simulation.h.

◆ m_rank

int m_rank
protected

MPI rank of this process

Definition at line 216 of file sparse_grids_simulation.h.

◆ m_nproc

int m_nproc
protected

Total number of MPI ranks

Definition at line 216 of file sparse_grids_simulation.h.

◆ m_interp_order

int m_interp_order
protected

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

Definition at line 220 of file sparse_grids_simulation.h.

◆ m_is_periodic

std::vector<bool> m_is_periodic
protected

Periodicity along each dimension

Definition at line 222 of file sparse_grids_simulation.h.

◆ m_write_sg_solutions

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.

◆ m_print_sg_errors

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.

◆ m_sim_fg

SimulationObject* m_sim_fg
protected

full grid simulation object

Definition at line 230 of file sparse_grids_simulation.h.

◆ m_sims_sg

std::vector<SimulationObject> m_sims_sg
protected

vector of sparse grids simulation objects

Definition at line 231 of file sparse_grids_simulation.h.

◆ m_n_fg

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.

◆ m_imin

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.

◆ m_combination

std::vector<SGCTElem> m_combination
protected

Coefficients and grid dimensions for the combination technique

Definition at line 236 of file sparse_grids_simulation.h.

◆ m_iprocs

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: