timeintegration_cpp.h File Reference

Contains C++ function declarations for time integration. More...

#include <timeintegration_struct.h>

Go to the source code of this file.

Functions
int	TimeInitialize (void *, int, int, int, void *)
int	TimeCleanup (void *)
int	TimePreStep (void *)
int	TimeStep (void *)
int	TimePostStep (void *)
int	TimePrintStep (void *)
int	TimeError (void *, void *, double *)
int	TimeGetAuxSolutions (int *, double **, void *, int, int)

Detailed Description

Contains C++ function declarations for time integration.

Author

Debojyoti Ghosh

Definition in file timeintegration_cpp.h.

Function Documentation

int TimeInitialize	(	void *	s,
		int	nsims,
		int	rank,
		int	nproc,
		void *	ts
	)

Initialize the time integration

Initialize time integration: This function initializes all that is required for time integration.

• It sets the number of iterations, time step size, simulation time, etc.
• It allocates solution, right-hand-side, and stage solution arrays needed by specific time integration methods.
• It calls the method-specific initialization functions.

Parameters

s	Array of simulation objects of type SimulationObject
nsims	number of simulation objects
rank	MPI rank of this process
nproc	number of MPI processes
ts	Time integration object of type TimeIntegration

Definition at line 28 of file TimeInitialize.c.

{
  TimeIntegration*  TS  = (TimeIntegration*) ts;
  SimulationObject* sim = (SimulationObject*) s;
  int ns, d, i;

  if (!sim) return(1);


  TS->simulation    = sim;
  TS->nsims         = nsims;

  TS->rank          = rank;
  TS->nproc         = nproc;

  TS->n_iter        = sim[0].solver.n_iter;
  TS->restart_iter  = sim[0].solver.restart_iter;
  TS->dt            = sim[0].solver.dt;

  TS->waqt          = (double) TS->restart_iter * TS->dt;
  TS->max_cfl       = 0.0;
  TS->norm          = 0.0;
  TS->TimeIntegrate = sim[0].solver.TimeIntegrate;
  TS->iter_wctime_total = 0.0;

  TS->u_offsets = (long*) calloc (nsims, sizeof(long));
  TS->u_sizes = (long*) calloc (nsims, sizeof(long));

  TS->u_size_total = 0;
  for (ns = 0; ns < nsims; ns++) {
    TS->u_offsets[ns] = TS->u_size_total;
    TS->u_sizes[ns] = sim[ns].solver.npoints_local_wghosts * sim[ns].solver.nvars ;
    TS->u_size_total += TS->u_sizes[ns];
  }

  TS->u   = (double*) calloc (TS->u_size_total,sizeof(double));
  TS->rhs = (double*) calloc (TS->u_size_total,sizeof(double));
  _ArraySetValue_(TS->u  ,TS->u_size_total,0.0);
  _ArraySetValue_(TS->rhs,TS->u_size_total,0.0);

  /* initialize arrays to NULL, then allocate as necessary */
  TS->U             = NULL;
  TS->Udot          = NULL;
  TS->BoundaryFlux  = NULL;

  TS->bf_offsets = (int*) calloc (nsims, sizeof(int));
  TS->bf_sizes = (int*) calloc (nsims, sizeof(int));
  TS->bf_size_total = 0;
  for (ns = 0; ns < nsims; ns++) {
    TS->bf_offsets[ns] = TS->bf_size_total;
    TS->bf_sizes[ns] =  2 * sim[ns].solver.ndims * sim[ns].solver.nvars;
    TS->bf_size_total += TS->bf_sizes[ns];
  }

#if defined(HAVE_CUDA)
  if (sim[0].solver.use_gpu) {

    if (!strcmp(sim[0].solver.time_scheme,_RK_)) {

      /* explicit Runge-Kutta methods */
      ExplicitRKParameters  *params = (ExplicitRKParameters*)  sim[0].solver.msti;
      int nstages = params->nstages;

      gpuMalloc((void**)&TS->gpu_U, TS->u_size_total*nstages*sizeof(double));
      gpuMalloc((void**)&TS->gpu_Udot, TS->u_size_total*nstages*sizeof(double));
      gpuMalloc((void**)&TS->gpu_BoundaryFlux, TS->bf_size_total*nstages*sizeof(double));
      gpuMemset(TS->gpu_U, 0, TS->u_size_total*nstages*sizeof(double));
      gpuMemset(TS->gpu_Udot, 0, TS->u_size_total*nstages*sizeof(double));
      gpuMemset(TS->gpu_BoundaryFlux, 0, TS->bf_size_total*nstages*sizeof(double));

    } else {

      fprintf(stderr,"ERROR in TimeInitialize(): %s is not yet implemented on GPUs.\n", 
              sim[0].solver.time_scheme );
      return 1;

    }

  } else {
#endif

    if (!strcmp(sim[0].solver.time_scheme,_RK_)) {
  
      /* explicit Runge-Kutta methods */
      ExplicitRKParameters  *params = (ExplicitRKParameters*)  sim[0].solver.msti;
      int nstages = params->nstages;
      TS->U     = (double**) calloc (nstages,sizeof(double*));
      TS->Udot  = (double**) calloc (nstages,sizeof(double*));
      for (i = 0; i < nstages; i++) {
        TS->U[i]    = (double*) calloc (TS->u_size_total,sizeof(double));
        TS->Udot[i] = (double*) calloc (TS->u_size_total,sizeof(double));
      }
  
      TS->BoundaryFlux = (double**) calloc (nstages,sizeof(double*));
      for (i=0; i<nstages; i++) {
        TS->BoundaryFlux[i] = (double*) calloc (TS->bf_size_total,sizeof(double));
      }
  
    } else if (!strcmp(sim[0].solver.time_scheme,_FORWARD_EULER_)) {
  
      int nstages = 1;
      TS->BoundaryFlux = (double**) calloc (nstages,sizeof(double*));
      for (i=0; i<nstages; i++) {
        TS->BoundaryFlux[i] = (double*) calloc (TS->bf_size_total,sizeof(double));
      }
  
    } else if (!strcmp(sim[0].solver.time_scheme,_GLM_GEE_)) {
  
      /* General Linear Methods with Global Error Estimate */
      GLMGEEParameters *params = (GLMGEEParameters*) sim[0].solver.msti;
      int nstages = params->nstages;
      int r       = params->r;
      TS->U     = (double**) calloc (2*r-1  ,sizeof(double*));
      TS->Udot  = (double**) calloc (nstages,sizeof(double*));
      for (i=0; i<2*r-1; i++)   TS->U[i]    = (double*) calloc (TS->u_size_total,sizeof(double));
      for (i=0; i<nstages; i++) TS->Udot[i] = (double*) calloc (TS->u_size_total,sizeof(double));
  
      TS->BoundaryFlux = (double**) calloc (nstages,sizeof(double*));
      for (i=0; i<nstages; i++) {
        TS->BoundaryFlux[i] = (double*) calloc (TS->bf_size_total,sizeof(double));
      }
  
  
      if (!strcmp(params->ee_mode,_GLM_GEE_YYT_)) {
        for (ns = 0; ns < nsims; ns++) {
          for (i=0; i<r-1; i++) {
            _ArrayCopy1D_(  (sim[ns].solver.u),
                            (TS->U[r+i] + TS->u_offsets[ns]),
                            (TS->u_sizes[ns])  );
          }
        }
      } else {
        for (i=0; i<r-1; i++) {
          _ArraySetValue_(TS->U[r+i],TS->u_size_total,0.0);
        }
      }
    }
#if defined(HAVE_CUDA)
  }
#endif

  /* set right-hand side function pointer */
  TS->RHSFunction = TimeRHSFunctionExplicit;

  /* open files for writing */
  if (!rank) {
    if (sim[0].solver.write_residual) TS->ResidualFile = (void*) fopen("residual.out","w");
    else                              TS->ResidualFile = NULL;
  } else                              TS->ResidualFile = NULL;

  for (ns = 0; ns < nsims; ns++) {
    sim[ns].solver.time_integrator = TS;
  }

  return 0;
}

int TimeCleanup

(

void *

ts

)

Clean up variables related to time integration

Clean up all allocations related to time integration

Parameters

ts	Object of type TimeIntegration

Definition at line 18 of file TimeCleanup.c.

{
  TimeIntegration* TS = (TimeIntegration*) ts;
  SimulationObject* sim = (SimulationObject*) TS->simulation;
  int ns, nsims = TS->nsims;

  /* close files opened for writing */
  if (!TS->rank) if (sim[0].solver.write_residual) fclose((FILE*)TS->ResidualFile);

#if defined(HAVE_CUDA)
  if (sim[0].solver.use_gpu) {
    if (!strcmp(sim[0].solver.time_scheme,_RK_)) {
      gpuFree(TS->gpu_U);
      gpuFree(TS->gpu_Udot);
      gpuFree(TS->gpu_BoundaryFlux);
    }
  } else {
#endif
    if (!strcmp(sim[0].solver.time_scheme,_RK_)) {
      int i;
      ExplicitRKParameters  *params = (ExplicitRKParameters*)  sim[0].solver.msti;
      for (i=0; i<params->nstages; i++) free(TS->U[i]);            free(TS->U);
      for (i=0; i<params->nstages; i++) free(TS->Udot[i]);         free(TS->Udot);
      for (i=0; i<params->nstages; i++) free(TS->BoundaryFlux[i]); free(TS->BoundaryFlux);
    } else if (!strcmp(sim[0].solver.time_scheme,_FORWARD_EULER_)) {
      int nstages = 1, i;
      for (i=0; i<nstages; i++) free(TS->BoundaryFlux[i]); free(TS->BoundaryFlux);
    } else if (!strcmp(sim[0].solver.time_scheme,_GLM_GEE_)) {
      int i;
      GLMGEEParameters  *params = (GLMGEEParameters*)  sim[0].solver.msti;
      for (i=0; i<2*params->r-1  ; i++) free(TS->U[i]);            free(TS->U);
      for (i=0; i<params->nstages; i++) free(TS->Udot[i]);         free(TS->Udot);
      for (i=0; i<params->nstages; i++) free(TS->BoundaryFlux[i]); free(TS->BoundaryFlux);
    }
#if defined(HAVE_CUDA)
  }
#endif

  /* deallocate arrays */
  free(TS->u_offsets);
  free(TS->u_sizes);
  free(TS->bf_offsets);
  free(TS->bf_sizes);
  free(TS->u  );
  free(TS->rhs);
  for (ns = 0; ns < nsims; ns++) {
    sim[ns].solver.time_integrator = NULL;
  }
  return(0);
}

int TimePreStep

(

void *

ts

)

Function called at the beginning of a time step

Pre-time-step function: This function is called before each time step. Some notable things this does are:

• Computes CFL and diffusion numbers.
• Call the physics-specific pre-time-step function, if defined.

Parameters

ts	Object of type TimeIntegration

Definition at line 22 of file TimePreStep.c.

{
  TimeIntegration*  TS  = (TimeIntegration*) ts;
  _DECLARE_IERR_;

  SimulationObject* sim = (SimulationObject*) TS->simulation;
  int ns, nsims = TS->nsims;

  gettimeofday(&TS->iter_start_time,NULL);

  for (ns = 0; ns < nsims; ns++) {

    HyPar*        solver = &(sim[ns].solver);
    MPIVariables* mpi    = &(sim[ns].mpi);

    double *u = NULL;
#if defined(HAVE_CUDA)
    if (solver->use_gpu) {
      u = solver->gpu_u;
    } else{ 
#endif
      u = solver->u;
#if defined(HAVE_CUDA)
    }
#endif

    /* apply boundary conditions and exchange data over MPI interfaces */

    solver->ApplyBoundaryConditions( solver,
                                     mpi,
                                     u,
                                     NULL,
                                     TS->waqt );

    solver->ApplyIBConditions( solver,
                               mpi,
                               u,
                               TS->waqt );

#if defined(HAVE_CUDA)
    if (solver->use_gpu) {
      gpuMPIExchangeBoundariesnD( solver->ndims,
                                  solver->nvars,
                                  solver->gpu_dim_local,
                                  solver->ghosts,
                                  mpi,
                                  u );
    } else {
#endif
      MPIExchangeBoundariesnD( solver->ndims,
                               solver->nvars,
                               solver->dim_local,
                               solver->ghosts,
                               mpi,
                               u );
#if defined(HAVE_CUDA)
    }
#endif

    if ((TS->iter+1)%solver->screen_op_iter == 0) {

#if defined(HAVE_CUDA)
      if (!solver->use_gpu) {
#endif
        _ArrayCopy1D_(  solver->u,
                        (TS->u + TS->u_offsets[ns]),
                        (solver->npoints_local_wghosts*solver->nvars) );
#if defined(HAVE_CUDA)
      }
#endif

      /* compute max CFL and diffusion number over the domain */
      if (solver->ComputeCFL) {
        double local_max_cfl  = -1.0;
        local_max_cfl  = solver->ComputeCFL (solver,mpi,TS->dt,TS->waqt);
        MPIMax_double(&TS->max_cfl ,&local_max_cfl ,1,&mpi->world);
      } else {
        TS->max_cfl = -1;
      }
      if (solver->ComputeDiffNumber) {
        double local_max_diff = -1.0;
        local_max_diff = solver->ComputeDiffNumber (solver,mpi,TS->dt,TS->waqt);
        MPIMax_double(&TS->max_diff,&local_max_diff,1,&mpi->world);
      } else {
        TS->max_diff = -1;
      }

    }

    /* set the step boundary flux integral value to zero */
#if defined(HAVE_CUDA)
    if (solver->use_gpu) {
      gpuArraySetValue(solver->StepBoundaryIntegral,2*solver->ndims*solver->nvars,0.0);
    } else {
#endif
      _ArraySetValue_(solver->StepBoundaryIntegral,2*solver->ndims*solver->nvars,0.0);
#if defined(HAVE_CUDA)
    }
#endif

    if (solver->PreStep) {
      solver->PreStep(u,solver,mpi,TS->waqt);
    }

  }

  return 0;
}

int TimeStep

(

void *

ts

)

Take one step in time

Advance one time step.

Parameters

ts	Object of type TimeIntegration

Definition at line 13 of file TimeStep.c.

{
  TimeIntegration *TS  = (TimeIntegration*) ts;
  SimulationObject* sim = (SimulationObject*) TS->simulation;

  if (TS->TimeIntegrate) { TS->TimeIntegrate(TS); }

  return(0);
}

int TimePostStep

(

void *

ts

)

Function called at the end of a time step

Post-time-step function: this function is called at the end of each time step.

• It updates the current simulation time.
• It calls functions to print information and to write transient solution to file.
• It will also call any physics-specific post-time-step function, if defined.

Parameters

ts	Object of type TimeIntegration

Definition at line 28 of file TimePostStep.c.

{
  TimeIntegration* TS = (TimeIntegration*) ts;
  SimulationObject* sim = (SimulationObject*) TS->simulation;
  int ns, nsims = TS->nsims;

  /* update current time */
  TS->waqt += TS->dt;

  if ((TS->iter+1)%sim[0].solver.screen_op_iter == 0) {

#if defined(HAVE_CUDA)
    if (sim[0].solver.use_gpu) {
      TS->norm = -1;
    } else {
#endif
      /* Calculate norm for this time step */
      double sum = 0.0;
      double npts = 0;
      for (ns = 0; ns < nsims; ns++) {
        _ArrayAXPY_(sim[ns].solver.u,-1.0,(TS->u+TS->u_offsets[ns]),TS->u_sizes[ns]);
        sum += ArraySumSquarenD( sim[ns].solver.nvars,
                                 sim[ns].solver.ndims,
                                 sim[ns].solver.dim_local,
                                 sim[ns].solver.ghosts,
                                 sim[ns].solver.index,
                                 (TS->u+TS->u_offsets[ns]) );
        npts += (double)sim[ns].solver.npoints_global;
      }

      double global_sum = 0;
      MPISum_double(  &global_sum,
                      &sum,1,
                      &(sim[0].mpi.world) );

      TS->norm = sqrt(global_sum/npts);
#if defined(HAVE_CUDA)
    }
#endif

    /* write to file */
    if (TS->ResidualFile) {
      fprintf((FILE*)TS->ResidualFile,"%10d\t%E\t%E\n",TS->iter+1,TS->waqt,TS->norm);
    }

  }


#if defined(HAVE_CUDA)
  if (!sim[0].solver.use_gpu) {
#endif
    for (ns = 0; ns < nsims; ns++) {
  
      if (!strcmp(sim[ns].solver.ConservationCheck,"yes")) {
        /* calculate volume integral of the solution at this time step */
        IERR sim[ns].solver.VolumeIntegralFunction( sim[ns].solver.VolumeIntegral,
                                                    sim[ns].solver.u,
                                                    &(sim[ns].solver),
                                                    &(sim[ns].mpi) ); CHECKERR(ierr);
        /* calculate surface integral of the flux at this time step */
        IERR sim[ns].solver.BoundaryIntegralFunction( &(sim[ns].solver),
                                                      &(sim[ns].mpi)); CHECKERR(ierr);
        /* calculate the conservation error at this time step       */
        IERR sim[ns].solver.CalculateConservationError( &(sim[ns].solver),
                                                        &(sim[ns].mpi)); CHECKERR(ierr);
      }
  
      if (sim[ns].solver.PostStep) { 
        sim[ns].solver.PostStep(sim[ns].solver.u,&(sim[ns].solver),&(sim[ns].mpi),TS->waqt,TS->iter);
      }
  
    }
#if defined(HAVE_CUDA)
  }
#endif

  gettimeofday(&TS->iter_end_time,NULL);
  long long walltime;
  walltime = (  (TS->iter_end_time.tv_sec * 1000000 + TS->iter_end_time.tv_usec)
              - (TS->iter_start_time.tv_sec * 1000000 + TS->iter_start_time.tv_usec));
  TS->iter_wctime = (double) walltime / 1000000.0;
  TS->iter_wctime_total += TS->iter_wctime;

  double global_total = 0, global_wctime = 0, global_mpi_total = 0, global_mpi_wctime = 0;

  MPIMax_double(&global_wctime, &TS->iter_wctime, 1, &(sim[0].mpi.world));
  MPIMax_double(&global_total, &TS->iter_wctime_total, 1, &(sim[0].mpi.world));

#if defined(HAVE_CUDA)
  MPIMax_double(&global_mpi_wctime, &sim[0].mpi.wctime, 1, &(sim[0].mpi.world));
  MPIMax_double(&global_mpi_total, &sim[0].mpi.wctime_total, 1, &(sim[0].mpi.world));
#endif

  return(0);
}

int TimePrintStep

(

void *

ts

)

Print time integration related information

Print information to screen (also calls any physics-specific printing function, if defined).

Parameters

ts	Object of type TimeIntegration

Definition at line 16 of file TimePrintStep.c.

{
  TimeIntegration* TS = (TimeIntegration*) ts;
  SimulationObject* sim = (SimulationObject*) TS->simulation;
  int ns, nsims = TS->nsims;

  if ((!TS->rank) && ((TS->iter+1)%sim[0].solver.screen_op_iter == 0)) {
    if (nsims > 1) {
      printf("--\n");
      printf("iter=%7d,  t=%1.3e\n", TS->iter+1, TS->waqt);
      if (TS->max_cfl >= 0) printf("  CFL=%1.3E\n", TS->max_cfl);
      if (TS->norm >= 0) printf("  norm=%1.4E\n", TS->norm);
      printf("  wctime=%1.1E (s)\n",TS->iter_wctime);
    } else {
      printf("iter=%7d  ",TS->iter+1 );
      printf("t=%1.3E  ",TS->waqt );
      if (TS->max_cfl >= 0) printf("CFL=%1.3E  ",TS->max_cfl );
      if (TS->norm >= 0) printf("norm=%1.4E  ",TS->norm );
      printf("wctime: %1.1E (s)  ",TS->iter_wctime);
    }

    /* calculate and print conservation error */
    if (!strcmp(sim[0].solver.ConservationCheck,"yes")) {

      double error = 0;
      for  (ns = 0; ns < nsims; ns++) {
        int v;
        for (v=0; v<sim[ns].solver.nvars; v++) {
          error += (sim[ns].solver.ConservationError[v] * sim[ns].solver.ConservationError[v]);
        }
      }
      error = sqrt(error);
      printf("  cons_err=%1.4E\n", error);

    } else {

      if (nsims == 1) printf("\n");

    }

    /* print physics-specific info, if available */
    for (ns = 0; ns < nsims; ns++) {
      if (sim[ns].solver.PrintStep) {
        if (nsims > 1) printf("Physics-specific output for domain %d:\n", ns);
        sim[ns].solver.PrintStep( &(sim[ns].solver),
                                  &(sim[ns].mpi),
                                  TS->waqt );
      }
    }
    if (nsims > 1) {
      printf("--\n");
      printf("\n");
    }
  }

  return(0);
}

int TimeError	(	void *	s,
		void *	m,
		double *	uex
	)

Compute/estimate error in solution

Computes the time integration error in the solution: If the time integration method chosen has a mechanism to compute the error in the numerical integration, it is computed in this function. In addition, if an exact solution is available (see function argument uex, the error of the computed solution (stored in HyPar::u) with respect to this exact solution is also computed. These errors are written to a text file whose name depends on the time integration method being used.

Time integration method with error estimation currently implemented:

• _GLM_GEE_ (General Linear Methods with Global Error Estimation) (see TimeGLMGEE(), TimeGLMGEEInitialize() )

Parameters

s	Solver object of type HyPar
m	MPI object of type MPIVariables
uex	Exact solution (stored with the same array layout as HyPar::u (may be NULL)

Definition at line 27 of file TimeError.c.

{
  HyPar               *solver = (HyPar*)           s;
  MPIVariables        *mpi    = (MPIVariables*)    m;
  TimeIntegration     *TS     = (TimeIntegration*) solver->time_integrator;
  int                 size    = solver->npoints_local_wghosts * solver->nvars;
  double              sum     = 0.0, global_sum = 0.0;
  static const double tolerance = 1e-15;
  if (!TS) return(0);

  if (!strcmp(solver->time_scheme,_GLM_GEE_)) {
    /* For GLM-GEE methods, calculate the norm of the estimated global error */
    GLMGEEParameters *params = (GLMGEEParameters*) solver->msti;
    double error[6] = {0,0,0,0,0,0}, *Uerr;
    if (!strcmp(params->ee_mode,_GLM_GEE_YEPS_)) Uerr = TS->U[params->r];
    else {
      Uerr = TS->U[0];
      _ArraySubtract1D_(Uerr,solver->u,TS->U[params->r],size);
      _ArrayScale1D_(Uerr,(1.0/(1.0-params->gamma)),size);
    }

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

    /* calculate L1 norm of error */
    sum = ArraySumAbsnD   (solver->nvars,solver->ndims,solver->dim_local,
                           solver->ghosts,solver->index,Uerr);
    global_sum = 0; MPISum_double(&global_sum,&sum,1,&mpi->world);
    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,Uerr);
    global_sum = 0; MPISum_double(&global_sum,&sum,1,&mpi->world);
    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,Uerr);
    global_sum = 0; MPIMax_double(&global_sum,&sum,1,&mpi->world);
    error[2] = global_sum;

    if (uex) {
      _ArrayAXBY_(TS->Udot[0],1.0,solver->u,-1.0,uex,size);
      _ArrayAXPY_(Uerr,-1.0,TS->Udot[0],size);
      /* calculate L1 norm of error */
      sum = ArraySumAbsnD   (solver->nvars,solver->ndims,solver->dim_local,
                             solver->ghosts,solver->index,TS->Udot[0]);
      global_sum = 0; MPISum_double(&global_sum,&sum,1,&mpi->world);
      error[3] = global_sum/((double)solver->npoints_global);
      /* calculate L2 norm of error */
      sum = ArraySumSquarenD(solver->nvars,solver->ndims,solver->dim_local,
                             solver->ghosts,solver->index,TS->Udot[0]);
      global_sum = 0; MPISum_double(&global_sum,&sum,1,&mpi->world);
      error[4] = 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,TS->Udot[0]);
      global_sum = 0; MPIMax_double(&global_sum,&sum,1,&mpi->world);
      error[5] = global_sum;
    } else error[3] = error[4] = error[5] = -1;

    if (   (sol_norm[0] > tolerance)
        && (sol_norm[1] > tolerance)
        && (sol_norm[2] > tolerance) ) {
      error[0] /= sol_norm[0];
      error[1] /= sol_norm[1];
      error[2] /= sol_norm[2];
      if (uex) {
        error[3] /= sol_norm[0];
        error[4] /= sol_norm[1];
        error[5] /= sol_norm[2];
      }
    }

    /* write to file */
    if (!mpi->rank) {
      FILE *out;
      out = fopen("glm_err.dat","w");
      fprintf(out,"%1.16E  %1.16E  %1.16E  %1.16E  ",TS->dt,error[0],error[1],error[2]);
      fprintf(out,"%1.16E  %1.16E  %1.16E\n",error[3],error[4],error[5]);
      fclose(out);
      printf("Estimated time integration errors (GLM-GEE time-integration):-\n");
      printf("  L1         Error           : %1.16E\n",error[0]);
      printf("  L2         Error           : %1.16E\n",error[1]);
      printf("  Linfinity  Error           : %1.16E\n",error[2]);
    }
  }

  return(0);
}

int TimeGetAuxSolutions	(	int *	N,
		double **	uaux,
		void *	s,
		int	nu,
		int	ns
	)

Function to get auxiliary solutions if available (for example, in GLM-GEE methods)

Return auxiliary solution: Some time integrators may have the concept of auxiliary solutions that they evolve along with the main solution HyPar::u (these may be used for error estimation, for example). This function returns a pointer to such an auxiliary solution. Note that the auxiliary solution has the same dimensions and array layout as the main solution.

• Call with the final argument n less than 0 to get the total number of auxiliary solutions (this value will be stored in the argument N at exit).
• Call with the final argument n less than or equal to zero to get the n-th auxiliary solution (the argument uaux will point to this at exit).

Note that auxiliary solutions are numbered in the C convention: 0,1,...,N-1.

Time integration methods which use auxiliary solutions currently implemented:

• General Linear Methods with Global Error Estimators (GLM-GEE - _GLM_GEE_) (see TimeGLMGEE(), TimeGLMGEEInitialize() ).

Parameters

N	Number of auxiliary solutions
uaux	Pointer to the array holding the auxiliary solution
s	Solver object of type HyPar
nu	Index of the auxiliary solution to return
ns	Index of the simulation domain of which the auxiliary solution to return

Definition at line 29 of file TimeGetAuxSolutions.c.

{
  HyPar           *solver = (HyPar*) s;
  TimeIntegration *TS     = (TimeIntegration*) solver->time_integrator;

  if (nu >= 0) {
    if (!strcmp(solver->time_scheme,_GLM_GEE_)) {
      GLMGEEParameters *params = (GLMGEEParameters*) solver->msti;
      *uaux = (TS->U[params->r+nu] + TS->u_offsets[ns]);
    }
  } else {
    if (!TS) *N = 0;
    else {
      if (!strcmp(solver->time_scheme,_GLM_GEE_)) {
        GLMGEEParameters *params = (GLMGEEParameters*) solver->msti;
        *N = params->r - 1;
      } else *N = 0;
    }
  }

  return(0);
}