HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
SparseGridsCalculateError.cpp
Go to the documentation of this file.
1 
6 #include <arrayfunctions.h>
7 #include <mathfunctions_cpp.h>
8 #include <std_vec_ops.h>
10 
11 extern "C" int ExactSolution(void*,void*,double*,char*,int*);
12 
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);
81  StdVecOps::copyFrom(dim_fg, m_sim_fg->solver.dim_global, m_ndims);
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  }
94  MPIGatherArraynDwGhosts( m_ndims,
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 }
153 
156  double* const a_uex )
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 }
224 
std::vector< SimulationObject > m_sims_sg
int npoints_local_wghosts
Definition: hypar.h:42
Contains some vector ops.
int * dim_global
Definition: hypar.h:33
double * u
Definition: hypar.h:116
void copyFrom(std::vector< int > &a_iv, const int *const a_iv_carr, int a_n)
Definition: std_vec_ops.h:36
int MPIPartitionArraynDwGhosts(int, void *, double *, double *, int *, int *, int, int)
int npoints_global
Definition: hypar.h:40
INLINE double ArraySumAbsnD(int, int, int *, int, int *, double *)
Contains function definitions for common mathematical functions for C++ code.
int * dim_local
Definition: hypar.h:37
int MPISum_double(double *, double *, int, void *)
Definition: MPISum.c:39
int MPIMax_double(double *, double *, int, void *)
Definition: MPIMax.c:38
std::vector< bool > m_is_periodic
int ghosts
Definition: hypar.h:52
MPI_Comm world
#define _MAX_STRING_SIZE_
Definition: basic.h:14
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)
INLINE double ArraySumSquarenD(int, int, int *, int, int *, double *)
int InterpolateGlobalnDVar(const int *const, double **const, const int *const, double *const, const int, const int, const int, const int *const)
#define _ArrayCopy1D_(x, y, size)
int nvars
Definition: hypar.h:29
#define CHECKERR(ierr)
Definition: basic.h:18
long sum(const std::vector< int > &a_iv)
Definition: std_vec_ops.h:18
int ExactSolution(void *, void *, double *, char *, int *)
Definition: ExactSolution.c:16
void computeError(SimulationObject &, double *const)
Structure defining a simulation.
int ndims
Definition: hypar.h:26
int * index
Definition: hypar.h:102
Contains macros and function definitions for common array operations.
#define IERR
Definition: basic.h:16
Structure of MPI-related variables.
INLINE double ArrayMaxnD(int, int, int *, int, int *, double *)
#define _ArrayAXPY_(x, a, y, size)
std::vector< int > GridDimensions
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
double error[3]
Definition: hypar.h:371