HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
TimeError.c
Go to the documentation of this file.
1 
6 #include <math.h>
7 #include <string.h>
8 #include <arrayfunctions.h>
9 #include <mpivars.h>
10 #include <hypar.h>
11 #include <timeintegration.h>
12 
28  void *s,
29  void *m,
30  double *uex
33  )
34 {
35  HyPar *solver = (HyPar*) s;
36  MPIVariables *mpi = (MPIVariables*) m;
38  int size = solver->npoints_local_wghosts * solver->nvars;
39  double sum = 0.0, global_sum = 0.0;
40  static const double tolerance = 1e-15;
41  if (!TS) return(0);
42 
43  if (!strcmp(solver->time_scheme,_GLM_GEE_)) {
44  /* For GLM-GEE methods, calculate the norm of the estimated global error */
45  GLMGEEParameters *params = (GLMGEEParameters*) solver->msti;
46  double error[6] = {0,0,0,0,0,0}, *Uerr;
47  if (!strcmp(params->ee_mode,_GLM_GEE_YEPS_)) Uerr = TS->U[params->r];
48  else {
49  Uerr = TS->U[0];
50  _ArraySubtract1D_(Uerr,solver->u,TS->U[params->r],size);
51  _ArrayScale1D_(Uerr,(1.0/(1.0-params->gamma)),size);
52  }
53 
54  /* calculate solution norm for relative errors */
55  double sol_norm[3] = {0.0,0.0,0.0};
56  /* L1 */
57  sum = ArraySumAbsnD (solver->nvars,solver->ndims,solver->dim_local,
58  solver->ghosts,solver->index,solver->u);
59  global_sum = 0; MPISum_double(&global_sum,&sum,1,&mpi->world);
60  sol_norm[0] = global_sum/((double)solver->npoints_global);
61  /* L2 */
62  sum = ArraySumSquarenD(solver->nvars,solver->ndims,solver->dim_local,
63  solver->ghosts,solver->index,solver->u);
64  global_sum = 0; MPISum_double(&global_sum,&sum,1,&mpi->world);
65  sol_norm[1] = sqrt(global_sum/((double)solver->npoints_global));
66  /* Linf */
67  sum = ArrayMaxnD (solver->nvars,solver->ndims,solver->dim_local,
68  solver->ghosts,solver->index,solver->u);
69  global_sum = 0; MPIMax_double(&global_sum,&sum,1,&mpi->world);
70  sol_norm[2] = global_sum;
71 
72  /* calculate L1 norm of error */
73  sum = ArraySumAbsnD (solver->nvars,solver->ndims,solver->dim_local,
74  solver->ghosts,solver->index,Uerr);
75  global_sum = 0; MPISum_double(&global_sum,&sum,1,&mpi->world);
76  error[0] = global_sum/((double)solver->npoints_global);
77  /* calculate L2 norm of error */
78  sum = ArraySumSquarenD(solver->nvars,solver->ndims,solver->dim_local,
79  solver->ghosts,solver->index,Uerr);
80  global_sum = 0; MPISum_double(&global_sum,&sum,1,&mpi->world);
81  error[1] = sqrt(global_sum/((double)solver->npoints_global));
82  /* calculate Linf norm of error */
83  sum = ArrayMaxnD (solver->nvars,solver->ndims,solver->dim_local,
84  solver->ghosts,solver->index,Uerr);
85  global_sum = 0; MPIMax_double(&global_sum,&sum,1,&mpi->world);
86  error[2] = global_sum;
87 
88  if (uex) {
89  _ArrayAXBY_(TS->Udot[0],1.0,solver->u,-1.0,uex,size);
90  _ArrayAXPY_(Uerr,-1.0,TS->Udot[0],size);
91  /* calculate L1 norm of error */
92  sum = ArraySumAbsnD (solver->nvars,solver->ndims,solver->dim_local,
93  solver->ghosts,solver->index,TS->Udot[0]);
94  global_sum = 0; MPISum_double(&global_sum,&sum,1,&mpi->world);
95  error[3] = global_sum/((double)solver->npoints_global);
96  /* calculate L2 norm of error */
97  sum = ArraySumSquarenD(solver->nvars,solver->ndims,solver->dim_local,
98  solver->ghosts,solver->index,TS->Udot[0]);
99  global_sum = 0; MPISum_double(&global_sum,&sum,1,&mpi->world);
100  error[4] = sqrt(global_sum/((double)solver->npoints_global));
101  /* calculate Linf norm of error */
102  sum = ArrayMaxnD (solver->nvars,solver->ndims,solver->dim_local,
103  solver->ghosts,solver->index,TS->Udot[0]);
104  global_sum = 0; MPIMax_double(&global_sum,&sum,1,&mpi->world);
105  error[5] = global_sum;
106  } else error[3] = error[4] = error[5] = -1;
107 
108  if ( (sol_norm[0] > tolerance)
109  && (sol_norm[1] > tolerance)
110  && (sol_norm[2] > tolerance) ) {
111  error[0] /= sol_norm[0];
112  error[1] /= sol_norm[1];
113  error[2] /= sol_norm[2];
114  if (uex) {
115  error[3] /= sol_norm[0];
116  error[4] /= sol_norm[1];
117  error[5] /= sol_norm[2];
118  }
119  }
120 
121  /* write to file */
122  if (!mpi->rank) {
123  FILE *out;
124  out = fopen("glm_err.dat","w");
125  fprintf(out,"%1.16E %1.16E %1.16E %1.16E ",TS->dt,error[0],error[1],error[2]);
126  fprintf(out,"%1.16E %1.16E %1.16E\n",error[3],error[4],error[5]);
127  fclose(out);
128  printf("Estimated time integration errors (GLM-GEE time-integration):-\n");
129  printf(" L1 Error : %1.16E\n",error[0]);
130  printf(" L2 Error : %1.16E\n",error[1]);
131  printf(" Linfinity Error : %1.16E\n",error[2]);
132  }
133  }
134 
135  return(0);
136 }
int npoints_local_wghosts
Definition: hypar.h:42
double * u
Definition: hypar.h:116
#define _ArrayAXBY_(z, a, x, b, y, size)
int TimeError(void *, void *, double *)
Definition: TimeError.c:27
INLINE double ArraySumAbsnD(int, int, int *, int, int *, double *)
Structure of variables/parameters and function pointers for time integration.
MPI related function definitions.
int MPISum_double(double *, double *, int, void *)
Definition: MPISum.c:39
#define _ArrayScale1D_(x, a, size)
int MPIMax_double(double *, double *, int, void *)
Definition: MPIMax.c:38
#define _ArraySubtract1D_(x, a, b, size)
#define _GLM_GEE_
Structure containing the parameters for the General Linear Methods with Global Error Estimators (GLM-...
INLINE double ArraySumSquarenD(int, int, int *, int, int *, double *)
char time_scheme[_MAX_STRING_SIZE_]
Definition: hypar.h:78
#define _GLM_GEE_YEPS_
int nvars
Definition: hypar.h:29
long sum(const std::vector< int > &a_iv)
Definition: std_vec_ops.h:18
Contains structure definition for hypar.
Contains function declarations for time integration.
void * msti
Definition: hypar.h:366
Contains macros and function definitions for common array operations.
Structure of MPI-related variables.
INLINE double ArrayMaxnD(int, int, int *, int, int *, double *)
#define _ArrayAXPY_(x, a, y, size)
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
void * time_integrator
Definition: hypar.h:165