CalculateError.c File Reference

Computes the error in the solution. More...

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <basic.h>
#include <common.h>
#include <arrayfunctions.h>
#include <timeintegration.h>
#include <mpivars.h>
#include <hypar.h>

Go to the source code of this file.

Functions
int	ExactSolution (void *, void *, double *, char *, int *)
int	CalculateError (void *s, void *m)

Detailed Description

Computes the error in the solution.

Author

Debojyoti Ghosh

Definition in file CalculateError.c.

Function Documentation

int ExactSolution	(	void *	s,
		void *	m,
		double *	uex,
		char *	fname,
		int *	flag
	)

Read in the exact solution, if available.

Parameters

s	Solver object of type HyPar
m	MPI object of type MPIVariables
uex	Array to hold the exact solution, if available
fname	Filename root from which to read exact solution
flag	Flag to indicate if exact solution was available

Definition at line 16 of file ExactSolution.c.

{
  HyPar  *solver = (HyPar*) s;

  int same_size = 1;
  for (int i=0; i<solver->ndims; i++) {
    if (solver->dim_global[i] != solver->dim_global_ex[i]) same_size = 0;
  }

  if (same_size) {
    ReadArray( solver->ndims,
               solver->nvars,
               solver->dim_global,
               solver->dim_local,
               solver->ghosts,
               solver,
               m,
               NULL,
               uex,
               fname,
               flag);
  } else {
    ReadArraywInterp( solver->ndims,
                      solver->nvars,
                      solver->dim_global,
                      solver->dim_local,
                      solver->dim_global_ex,
                      solver->ghosts,
                      solver,
                      m,
                      NULL,
                      uex,
                      fname,
                      flag);
  }

  return(0);
}

int CalculateError	(	void *	s,
		void *	m
	)

Calculates the error in the solution if the exact solution is available. If the exact solution is not available, the errors are reported as zero. The exact solution should be provided in the file "exact.inp" in the same format as the initial solution.

Parameters

s	Solver object of type HyPar
m	MPI object of type MPIVariables

Definition at line 25 of file CalculateError.c.

{
  HyPar         *solver     = (HyPar*)        s;
  MPIVariables  *mpi        = (MPIVariables*) m;
  int           exact_flag  = 0, i, size;
  double        sum         = 0, global_sum = 0;
  double        *uex        = NULL;
  _DECLARE_IERR_;

  size = solver->nvars;
  for (i = 0; i < solver->ndims; i++) 
    size *= (solver->dim_local[i]+2*solver->ghosts);
  uex = (double*) calloc (size, sizeof(double));

  char fname_root[_MAX_STRING_SIZE_] = "exact";
  if (solver->nsims > 1) {
    char index[_MAX_STRING_SIZE_];
    GetStringFromInteger(solver->my_idx, index, (int)log10(solver->nsims)+1);
    strcat(fname_root, "_");
    strcat(fname_root, index);
  }

  static const double tolerance = 1e-15;
  IERR ExactSolution( solver,
                      mpi,
                      uex,
                      fname_root,
                      &exact_flag ); CHECKERR(ierr);

  if (!exact_flag) {

    /* No exact solution */
    IERR TimeError(solver,mpi,NULL); CHECKERR(ierr);
    solver->error[0] = solver->error[1] = solver->error[2] = -1;

  } else {

    IERR TimeError(solver,mpi,uex); CHECKERR(ierr);

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

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

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

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

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

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

  free(uex);
  return(0);
}