InitializeSolvers.c File Reference

Initialize all solvers. More...

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <io.h>
#include <tridiagLU.h>
#include <timeintegration.h>
#include <interpolation.h>
#include <firstderivative.h>
#include <secondderivative.h>
#include <mpivars.h>
#include <simulation_object.h>
#include <Python.h>

Go to the source code of this file.

Functions
int	ApplyBoundaryConditions (void *, void *, double *, double *, double)
	Applies the boundary conditions specified for each boundary zone. More...
int	ApplyIBConditions (void *, void *, double *, double)
int	HyperbolicFunction (double *, double *, void *, void *, double, int, int(*)(double *, double *, int, void *, double), int(*)(double *, double *, double *, double *, double *, double *, int, void *, double))
int	ParabolicFunctionNC1Stage (double *, double *, void *, void *, double)
int	ParabolicFunctionNC2Stage (double *, double *, void *, void *, double)
int	ParabolicFunctionNC1_5Stage (double *, double *, void *, void *, double)
int	ParabolicFunctionCons1Stage (double *, double *, void *, void *, double)
int	SourceFunction (double *, double *, void *, void *, double)
int	VolumeIntegral (double *, double *, void *, void *)
int	BoundaryIntegral (void *, void *)
int	CalculateConservationError (void *, void *)
void	IncrementFilenameIndex (char *, int)
int	NonLinearInterpolation (double *, void *, void *, double, int(*)(double *, double *, int, void *, double))
int	gpuHyperbolicFunction (double *, double *, void *, void *, double, int, int(*)(double *, double *, int, void *, double), int(*)(double *, double *, double *, double *, double *, double *, int, void *, double))
int	InitializeSolvers (void *s, int nsims)

Detailed Description

Initialize all solvers.

Author

Debojyoti Ghosh

Definition in file InitializeSolvers.c.

Function Documentation

int ApplyBoundaryConditions	(	void *	s,
		void *	m,
		double *	x,
		double *	xref,
		double	waqt
	)

Applies the boundary conditions specified for each boundary zone.

The solver object (of type HyPar) contains an oject of type DomainBoundary that contains all the boundary information (dimension, extent, face, type, etc). This function iterates through each of the boundary zones (HyPar::boundary[HyPar::nBoundaryZones]) and calls the corresponding boundary condition function.

The variable flag indicates if the array x is the solution, or a delta-solution (from implicit time-integration methods).

Parameters

s	Object of type HyPar containing solver-related variables
m	Object of type MPIVariables containing MPI-related variables
x	The solution vector on which the boundary conditions are to be applied
xref	Reference solution vector, if needed
waqt	Current simulation time

Definition at line 28 of file ApplyBoundaryConditions.c.

{
  HyPar           *solver   = (HyPar*)          s;
  DomainBoundary  *boundary = (DomainBoundary*) solver->boundary;
  MPIVariables    *mpi      = (MPIVariables*)   m;
  int             nb        = solver->nBoundaryZones;

  int* dim_local;
#if defined(HAVE_CUDA)
  if (solver->use_gpu) {
    dim_local = solver->gpu_dim_local;
  } else {
#endif
    dim_local = solver->dim_local;
#if defined(HAVE_CUDA)
  }
#endif

  /* Apply domain boundary conditions to x */
  int n;
  for (n = 0; n < nb; n++) {
    boundary[n].BCFunctionU(&boundary[n],mpi,solver->ndims,solver->nvars,
                            dim_local,solver->ghosts,x,waqt);
  }

  return(0);
}

int ApplyIBConditions	(	void *	s,
		void *	m,
		double *	x,
		double	waqt
	)

Call the physics-specific function that applies the immersed boundary conditions on the immersed boundary points.

Parameters

s	Object of type HyPar containing solver-related variables
m	Object of type MPIVariables containing MPI-related variables
x	The solution vector to apply immersed BCs on
waqt	Current simulation time

Definition at line 14 of file ApplyIBConditions.c.

{
  HyPar           *solver   = (HyPar*)          s;
  MPIVariables    *mpi      = (MPIVariables*)   m;

  /* Apply immersed boundary conditions, if applicable */
#if defined(HAVE_CUDA)
  if (solver->use_gpu) {
    if (solver->flag_ib) {
      fprintf(stderr, "ERROR: immersed boundaries have not yet been implemented on GPU.\n");
    }
  } else {
#endif
    if (solver->flag_ib) solver->IBFunction(solver,mpi,x,waqt);
#if defined(HAVE_CUDA)
  }
#endif

  return 0;
}

int HyperbolicFunction	(	double *	hyp,
		double *	u,
		void *	s,
		void *	m,
		double	t,
		int	LimFlag,
		int(*)(double *, double *, int, void *, double)	FluxFunction,
		int(*)(double *, double *, double *, double *, double *, double *, int, void *, double)	UpwindFunction
	)

This function computes the hyperbolic term of the governing equations, expressed as follows:-

\begin{equation} {\bf F} \left({\bf u}\right) = \sum_{d=0}^{D-1} \frac {\partial {\bf f}_d\left({\bf u}\right)} {\partial x_d} \end{equation}

using a conservative finite-difference discretization is space:

\begin{equation} \hat{\bf F} \left({\bf u}\right) = \sum_{d=0}^{D-1} \frac {1} {\Delta x_d} \left[ \hat{\bf f}_{d,j+1/2} - \hat{\bf f}_{d,j-1/2} \right] \end{equation}

where \(d\) denotes the spatial dimension, \(D\) denotes the total number of spatial dimensions, the hat denotes the discretized quantity, and \(j\) is the grid coordinate along dimension \(d\). The approximation to the flux function \({\bf f}_d\) at the interfaces \(j\pm1/2\), denoted by \(\hat{\bf f}_{d,j\pm 1/2}\), are computed using the function ReconstructHyperbolic().

Parameters

hyp	Array to hold the computed hyperbolic term (shares the same layout as u
u	Solution array
s	Solver object of type HyPar
m	MPI object of type MPIVariables
t	Current simulation time
LimFlag	Flag to indicate if the nonlinear coefficients for solution-dependent interpolation method should be recomputed (see ReconstructHyperbolic() for an explanation on why this is needed)
FluxFunction	Function pointer to the flux function for the hyperbolic term
UpwindFunction	Function pointer to the upwinding function for the hyperbolic term

Definition at line 31 of file HyperbolicFunction.c.

{
  HyPar         *solver = (HyPar*)        s;
  MPIVariables  *mpi    = (MPIVariables*) m;
  int           d, v, i, done;
  double        *FluxI  = solver->fluxI; /* interface flux     */
  double        *FluxC  = solver->fluxC; /* cell centered flux */
  _DECLARE_IERR_;

  int     ndims  = solver->ndims;
  int     nvars  = solver->nvars;
  int     ghosts = solver->ghosts;
  int     *dim   = solver->dim_local;
  int     size   = solver->npoints_local_wghosts;
  double  *x     = solver->x;
  double  *dxinv = solver->dxinv;
  int     index[ndims], index1[ndims], index2[ndims], dim_interface[ndims];

  LimFlag = (LimFlag && solver->flag_nonlinearinterp && solver->SetInterpLimiterVar);

  _ArraySetValue_(hyp,size*nvars,0.0);
  _ArraySetValue_(solver->StageBoundaryIntegral,2*ndims*nvars,0.0);
  if (!FluxFunction) return(0); /* zero hyperbolic term */
  solver->count_hyp++;

#if defined(CPU_STAT)
  double cpu_time = 0.0;
  clock_t cpu_start, cpu_end;
#endif

  int offset = 0;
  for (d = 0; d < ndims; d++) {
    _ArrayCopy1D_(dim,dim_interface,ndims); dim_interface[d]++;
    int size_cellcenter = 1; for (i = 0; i < ndims; i++) size_cellcenter *= (dim[i] + 2*ghosts);
    int size_interface = 1; for (i = 0; i < ndims; i++) size_interface *= dim_interface[i];

    /* evaluate cell-centered flux */
    IERR FluxFunction(FluxC,u,d,solver,t); CHECKERR(ierr);
    /* compute interface fluxes */
    IERR ReconstructHyperbolic(FluxI,FluxC,u,x+offset,d,solver,mpi,t,LimFlag,UpwindFunction);
    CHECKERR(ierr);

    /* calculate the first derivative */
    done = 0; _ArraySetValue_(index,ndims,0);
    int p, p1, p2;

#if defined(CPU_STAT)
    cpu_start = clock();
#endif

    while (!done) {
      _ArrayCopy1D_(index,index1,ndims);
      _ArrayCopy1D_(index,index2,ndims); index2[d]++;
      _ArrayIndex1D_(ndims,dim          ,index ,ghosts,p);
      _ArrayIndex1D_(ndims,dim_interface,index1,0     ,p1);
      _ArrayIndex1D_(ndims,dim_interface,index2,0     ,p2);
      for (v=0; v<nvars; v++) hyp[nvars*p+v] += dxinv[offset+ghosts+index[d]]
                                              * (FluxI[nvars*p2+v]-FluxI[nvars*p1+v]);
      /* boundary flux integral */
      if (index[d] == 0)
        for (v=0; v<nvars; v++) solver->StageBoundaryIntegral[(2*d+0)*nvars+v] -= FluxI[nvars*p1+v];
      if (index[d] == dim[d]-1)
        for (v=0; v<nvars; v++) solver->StageBoundaryIntegral[(2*d+1)*nvars+v] += FluxI[nvars*p2+v];

      _ArrayIncrementIndex_(ndims,dim,index,done);
    }

#if defined(CPU_STAT)
    cpu_end = clock();
    cpu_time += (double)(cpu_end - cpu_start) / CLOCKS_PER_SEC;
#endif

    offset += dim[d] + 2*ghosts;
  }

#if defined(CPU_STAT)
  printf("HyperbolicFunction CPU time = %8.6lf\n", cpu_time);
#endif

  if (solver->flag_ib) _ArrayBlockMultiply_(hyp,solver->iblank,size,nvars);

  return(0);
}

int ParabolicFunctionNC1Stage	(	double *	par,
		double *	u,
		void *	s,
		void *	m,
		double	t
	)

Evaluate the parabolic term using a conservative finite-difference spatial discretization: The parabolic term is assumed to be of the form:

\begin{equation} {\bf P}\left({\bf u}\right) = \sum_{d=0}^{D-1} \frac {\partial^2 {\bf g}_d\left(\bf u\right)} {\partial x_d^2}, \end{equation}

which is discretized at a grid point as:

\begin{equation} \left.{\bf P}\left({\bf u}\right)\right|_j = \sum_{d=0}^{D-1} \frac { \mathcal{L}_d \left[ {\bf g}_d \right] } {\Delta x_d^2}, \end{equation}

where \(d\) is the spatial dimension index, \(D\) is the total number of spatial dimensions (HyPar::ndims), and \(j\) is the grid index along \(d\). \(\mathcal{L}_d\) represents the finite-difference approximation to the Laplacian along the \(d\) spatial dimension, computed using HyPar::SecondDerivativePar.

Note: this form of the parabolic term does not allow for cross-derivatives.

To use this form of the parabolic term:

• specify "par_space_type" in solver.inp as "nonconservative-1stage" (HyPar::spatial_type_par).
• the physical model must specify \({\bf g}_d\left({\bf u}\right)\) through HyPar::GFunction.

Parameters

par	array to hold the computed parabolic term
u	solution
s	Solver object of type HyPar
m	MPI object of type MPIVariables
t	Current simulation time

Definition at line 31 of file ParabolicFunctionNC1Stage.c.

{
  HyPar         *solver = (HyPar*)        s;
  MPIVariables  *mpi    = (MPIVariables*) m;
  double        *Func   = solver->fluxC;
  double        *Deriv2 = solver->Deriv2;
  int           d, v, i, done;
  _DECLARE_IERR_;

  int     ndims  = solver->ndims;
  int     nvars  = solver->nvars;
  int     ghosts = solver->ghosts;
  int     *dim   = solver->dim_local;
  double  *dxinv = solver->dxinv;
  int     size   = solver->npoints_local_wghosts;

  if (!solver->GFunction) return(0); /* zero parabolic terms */
  solver->count_par++;

  int index[ndims];
  _ArraySetValue_(par,size*nvars,0.0);

  int offset = 0;
  for (d = 0; d < ndims; d++) {
    int size_deriv = 1; for (i=0; i<ndims; i++) size_deriv *= dim[i];

    /* calculate the diffusion function */
    IERR solver->GFunction(Func,u,d,solver,t); CHECKERR(ierr);
    IERR solver->SecondDerivativePar(Deriv2,Func,d,solver,mpi); CHECKERR(ierr);

    /* calculate the final term - second derivative of the diffusion function */
    done = 0; _ArraySetValue_(index,ndims,0);
    int p;
    while (!done) {
      _ArrayIndex1D_(ndims,dim,index,ghosts,p);
      for (v=0; v<nvars; v++)
        par[nvars*p+v] += (   dxinv[offset+ghosts+index[d]]*dxinv[offset+ghosts+index[d]] 
                            * Deriv2[nvars*p+v] );
      _ArrayIncrementIndex_(ndims,dim,index,done);
    }

    offset += dim[d] + 2*ghosts;
  }

  if (solver->flag_ib) _ArrayBlockMultiply_(par,solver->iblank,size,nvars);
  return(0);
}

int ParabolicFunctionNC2Stage	(	double *	par,
		double *	u,
		void *	s,
		void *	m,
		double	t
	)

Evaluate the parabolic term using a "1.5"-stage finite-difference spatial discretization: The parabolic term is assumed to be of the form:

\begin{equation} {\bf P}\left({\bf u}\right) = \sum_{d1=0}^{D-1}\sum_{d2=0}^{D-1} \frac {\partial^2 h_{d1,d2}\left(\bf u\right)} {\partial x_{d1} \partial x_{d2}}, \end{equation}

where \(d1\) and \(d2\) are spatial dimension indices, and \(D\) is the total number of spatial dimensions (HyPar::ndims). This term is discretized at a grid point as:

\begin{equation} \left.{\bf P}\left({\bf u}\right)\right|_j = \sum_{d1=0}^{D-1} \sum_{d2=0}^{D-1} \frac { \mathcal{D}_{d1}\mathcal{D}_{d2} \left[ {\bf h}_{d1,d2} \right] } {\Delta x_{d1} \Delta x_{d2}}, \end{equation}

where \(\mathcal{D}\) denotes the finite-difference approximation to the first derivative. Each of the first derivative approximations are \(\mathcal{D}_{d1}\) and \(\mathcal{D}_{d2}\) are computed separately, and thus the cross-derivative is evaluated in two steps using HyPar::FirstDerivativePar.

Notes:

• This form of the parabolic term does allow for cross-derivatives ( \( d1 \ne d2 \)).
• A \(n\)-th order central approximation to the second derivative can be expressed as a conjugation of two \((n-1)\)-th order approximations to the first derivative, one forward and one backward. Computing it this way avoids odd-even decoupling. Thus, where possible HyPar::FirstDerivativePar should point to the function computing \((n-1)\)-th order first derivative where \(n\) is the desired order. Currently, this is implemented only for \(n=2\). For other values of \(n\), the first derivative is also computed with a \(n\)-th order approximation.

To use this form of the parabolic term:

• specify "par_space_type" in solver.inp as "nonconservative-2stage" (HyPar::spatial_type_par).
• the physical model must specify \({\bf h}_{d1,d2}\left({\bf u}\right)\) through HyPar::HFunction.

See Also

ParabolicFunctionNC1_5Stage()

Parameters

par	array to hold the computed parabolic term
u	solution
s	Solver object of type HyPar
m	MPI object of type MPIVariables
t	Current simulation time

Definition at line 40 of file ParabolicFunctionNC2Stage.c.

{
  HyPar         *solver = (HyPar*)        s;
  MPIVariables  *mpi    = (MPIVariables*) m;
  double        *Func   = solver->fluxC;
  double        *Deriv1 = solver->Deriv1;
  double        *Deriv2 = solver->Deriv2;
  int           d, d1, d2, v, p, done;
  double        dxinv1, dxinv2;
  _DECLARE_IERR_;

  int     ndims  = solver->ndims;
  int     nvars  = solver->nvars;
  int     ghosts = solver->ghosts;
  int     *dim   = solver->dim_local;
  double  *dxinv = solver->dxinv;
  int     size   = solver->npoints_local_wghosts;

  printf("HFunction is defined? = %p\n", solver->HFunction);

  if (!solver->HFunction) return(0); /* zero parabolic terms */
  solver->count_par++;

  int index[ndims];
  _ArraySetValue_(par,size*nvars,0.0);

  for (d1 = 0; d1 < ndims; d1++) {
    for (d2 = 0; d2 < ndims; d2++) {

      /* calculate the diffusion function */
      IERR solver->HFunction(Func,u,d1,d2,solver,t);                    CHECKERR(ierr);
      IERR solver->FirstDerivativePar(Deriv1,Func  ,d1, 1,solver,mpi);  CHECKERR(ierr);
      IERR MPIExchangeBoundariesnD(ndims,nvars,dim,ghosts,mpi,Deriv1);  CHECKERR(ierr);
      IERR solver->FirstDerivativePar(Deriv2,Deriv1,d2,-1,solver,mpi);  CHECKERR(ierr);

      /* calculate the final term - second derivative of the diffusion function */
      done = 0; _ArraySetValue_(index,ndims,0);
      while (!done) {
        _ArrayIndex1D_(ndims,dim,index,ghosts,p);
        _GetCoordinate_(d1,index[d1],dim,ghosts,dxinv,dxinv1);
        _GetCoordinate_(d2,index[d2],dim,ghosts,dxinv,dxinv2);
        for (v=0; v<nvars; v++) par[nvars*p+v] += (dxinv1*dxinv2 * Deriv2[nvars*p+v]);
        _ArrayIncrementIndex_(ndims,dim,index,done);
      }

    }
  }

  if (solver->flag_ib) _ArrayBlockMultiply_(par,solver->iblank,size,nvars);
  return(0);
}

int ParabolicFunctionNC1_5Stage	(	double *	par,
		double *	u,
		void *	s,
		void *	m,
		double	t
	)

Evaluate the parabolic term using a "1.5"-stage finite-difference spatial discretization: The parabolic term is assumed to be of the form:

\begin{equation} {\bf P}\left({\bf u}\right) = \sum_{d1=0}^{D-1}\sum_{d2=0}^{D-1} \frac {\partial^2 h_{d1,d2}\left(\bf u\right)} {\partial x_{d1} \partial x_{d2}}, \end{equation}

where \(d1\) and \(d2\) are spatial dimension indices, and \(D\) is the total number of spatial dimensions (HyPar::ndims). This term is discretized at a grid point as:

\begin{equation} \left.{\bf P}\left({\bf u}\right)\right|_j = \sum_{d1=0}^{D-1} \sum_{d2=0}^{D-1} \frac { \mathcal{D}_{d1}\mathcal{D}_{d2} \left[ {\bf h}_{d1,d2} \right] } {\Delta x_{d1} \Delta x_{d2}}, \end{equation}

where \(\mathcal{D}\) denotes the finite-difference approximation to the first derivative.

• When \(d = d1 = d2\), the operator \(\mathcal{D}^2_d\left[\cdot\right] = \mathcal{D}_{d}\mathcal{D}_{d}\left[\cdot\right]\) is computed in one step as the finite-difference approximation to the second derivative (Laplacian) \(\mathcal{L}_d\left[\cdot\right]\) using HyPar::SecondDerivativePar.
• When \(d1 \ne d2 \), each of the first derivative approximations are \(\mathcal{D}_{d1}\) and \(\mathcal{D}_{d2}\) are computed separately, and thus the cross-derivative is evaluated in two steps using HyPar::FirstDerivativePar.

Note: this form of the parabolic term does allow for cross-derivatives ( \( d1 \ne d2 \)).

To use this form of the parabolic term:

• specify "par_space_type" in solver.inp as "nonconservative-1.5stage" (HyPar::spatial_type_par).
• the physical model must specify \({\bf h}_{d1,d2}\left({\bf u}\right)\) through HyPar::HFunction.

See Also

ParabolicFunctionNC2Stage()

Parameters

par	array to hold the computed parabolic term
u	solution
s	Solver object of type HyPar
m	MPI object of type MPIVariables
t	Current simulation time

Definition at line 35 of file ParabolicFunctionNC1.5Stage.c.

{
  HyPar         *solver = (HyPar*)        s;
  MPIVariables  *mpi    = (MPIVariables*) m;
  double        *Func   = solver->fluxC;
  double        *Deriv1 = solver->Deriv1;
  double        *Deriv2 = solver->Deriv2;
  int           d, d1, d2, v, p, done;
  double        dxinv1, dxinv2;
  _DECLARE_IERR_;

  int     ndims   = solver->ndims;
  int     nvars   = solver->nvars;
  int     ghosts  = solver->ghosts;
  int     *dim    = solver->dim_local;
  double  *dxinv  = solver->dxinv;
  int     size   = solver->npoints_local_wghosts;

  if (!solver->HFunction) return(0); /* zero parabolic terms */
  solver->count_par++;

  int index[ndims];
  _ArraySetValue_(par,size*nvars,0.0);

  for (d1 = 0; d1 < ndims; d1++) {
    for (d2 = 0; d2 < ndims; d2++) {

      /* calculate the diffusion function */
      IERR solver->HFunction(Func,u,d1,d2,solver,t);                      CHECKERR(ierr);
      if (d1 == d2) {
        /* second derivative with respect to the same independent variable */
        IERR solver->SecondDerivativePar(Deriv2,Func,d1,solver,mpi);      CHECKERR(ierr);
      } else {
        /* cross derivative */
        IERR solver->FirstDerivativePar(Deriv1,Func  ,d1, 1,solver,mpi);  CHECKERR(ierr);
        IERR MPIExchangeBoundariesnD(ndims,nvars,dim,ghosts,mpi,Deriv1);  CHECKERR(ierr);
        IERR solver->FirstDerivativePar(Deriv2,Deriv1,d2,-1,solver,mpi);  CHECKERR(ierr);
      }

      /* calculate the final term - second derivative of the diffusion function */
      done = 0; _ArraySetValue_(index,ndims,0);
      while (!done) {
        _ArrayIndex1D_(ndims,dim,index,ghosts,p);
        _GetCoordinate_(d1,index[d1],dim,ghosts,dxinv,dxinv1);
        _GetCoordinate_(d2,index[d2],dim,ghosts,dxinv,dxinv2);
        for (v=0; v<nvars; v++) par[nvars*p+v] += (dxinv1*dxinv2 * Deriv2[nvars*p+v]);
        _ArrayIncrementIndex_(ndims,dim,index,done);
      }

    }
  }

  if (solver->flag_ib) _ArrayBlockMultiply_(par,solver->iblank,size,nvars);
  return(0);
}

int ParabolicFunctionCons1Stage	(	double *	par,
		double *	u,
		void *	s,
		void *	m,
		double	t
	)

Evaluate the parabolic term using a conservative finite-difference spatial discretization: The parabolic term is assumed to be of the form:

\begin{equation} {\bf P}\left({\bf u}\right) = \sum_{d=0}^{D-1} \frac {\partial^2 {\bf g}_d\left(\bf u\right)} {\partial x_d^2}, \end{equation}

which is discretized at a grid point as:

\begin{equation} \left.{\bf P}\left({\bf u}\right)\right|_j = \sum_{d=0}^{D-1} \frac { \hat{\bf g}_{j+1/2} - \hat{\bf g}_{j-1/2} } {\Delta x_d^2}, \end{equation}

where \(d\) is the spatial dimension index, \(D\) is the total number of spatial dimensions (HyPar::ndims), and \(j\) is the grid index along \(d\). \(\hat{\bf g}_d\) is the numerical approximation to the second primitive of \({\bf g}_d\left({\bf u}\right)\), computed using HyPar::InterpolateInterfacesPar.

Note: this form of the parabolic term does not allow for cross-derivatives.

To use this form of the parabolic term:

• specify "par_space_type" in solver.inp as "conservative-1stage" (HyPar::spatial_type_par).
• the physical model must specify \({\bf g}_d\left({\bf u}\right)\) through HyPar::GFunction.

Reference: Liu, Y., Shu, C.-W., Zhang, M., "High order finite difference WENO schemes for nonlinear degenerate parabolic equations", SIAM J. Sci. Comput., 33 (2), 2011, pp. 939-965, http://dx.doi.org/10.1137/100791002

Parameters

par	array to hold the computed parabolic term
u	solution
s	Solver object of type HyPar
m	MPI object of type MPIVariables
t	Current simulation time

Definition at line 34 of file ParabolicFunctionCons1Stage.c.

{
  HyPar         *solver = (HyPar*)        s;
  MPIVariables  *mpi    = (MPIVariables*) m;
  double        *FluxI  = solver->fluxI; /* interface flux     array */
  double        *Func   = solver->fluxC; /* diffusion function array */
  int           d, v, i, done;
  _DECLARE_IERR_;

  int     ndims  = solver->ndims;
  int     nvars  = solver->nvars;
  int     ghosts = solver->ghosts;
  int     *dim   = solver->dim_local;
  double  *dxinv = solver->dxinv;
  int     size   = solver->npoints_local_wghosts;

  int index[ndims], index1[ndims], index2[ndims], dim_interface[ndims];

  _ArraySetValue_(par,size*nvars,0.0);
  if (!solver->GFunction) return(0); /* zero parabolic term */
  solver->count_par++;

  int offset = 0;
  for (d = 0; d < ndims; d++) {
    _ArrayCopy1D_(dim,dim_interface,ndims); dim_interface[d]++;
    int size_cellcenter = 1; for (i = 0; i < ndims; i++) size_cellcenter *= (dim[i] + 2*ghosts);
    int size_interface = 1; for (i = 0; i < ndims; i++) size_interface *= dim_interface[i];

    /* evaluate cell-centered flux */
    IERR solver->GFunction(Func,u,d,solver,t); CHECKERR(ierr);
    /* compute interface fluxes */
    IERR solver->InterpolateInterfacesPar(FluxI,Func,d,solver,mpi); CHECKERR(ierr);

    /* calculate the second derivative */
    done = 0; _ArraySetValue_(index,ndims,0);
    int p, p1, p2;
    while (!done) {
      _ArrayCopy1D_(index,index1,ndims);
      _ArrayCopy1D_(index,index2,ndims); index2[d]++;
      _ArrayIndex1D_(ndims,dim          ,index ,ghosts,p);
      _ArrayIndex1D_(ndims,dim_interface,index1,0     ,p1);
      _ArrayIndex1D_(ndims,dim_interface,index2,0     ,p2);
      for (v=0; v<nvars; v++) 
        par[nvars*p+v] +=  ((dxinv[offset+ghosts+index[d]] * dxinv[offset+ghosts+index[d]]) 
                          * (FluxI[nvars*p2+v] - FluxI[nvars*p1+v]));
      _ArrayIncrementIndex_(ndims,dim,index,done);
    }

    offset += dim[d] + 2*ghosts;
  }

  if (solver->flag_ib) _ArrayBlockMultiply_(par,solver->iblank,size,nvars);
  return(0);
}

int SourceFunction	(	double *	source,
		double *	u,
		void *	s,
		void *	m,
		double	t
	)

Evaluate the source term \({\bf S}\left({\bf u}\right)\) in the governing equation, if the physical model specifies one. In addition, if the simulation requires a sponge boundary treatment, the sponge BC function is called.

Parameters

source	the computed source term
u	solution
s	solver object of type HyPar
m	MPI object of type MPIVariables
t	Current simulation time

Definition at line 22 of file SourceFunction.c.

{
  HyPar           *solver   = (HyPar*)        s;
  MPIVariables    *mpi      = (MPIVariables*) m;

  /* initialize to zero */
  int size = solver->ndof_cells_wghosts;
#if defined(HAVE_CUDA)
  if (solver->use_gpu) {
    gpuArraySetValue(source,size, 0.0);
  } else {
#endif
    _ArraySetValue_(source,size,0.0);
#if defined(HAVE_CUDA)
  }
#endif

  /* call the source function of the physics model, if available */
  if (solver->SFunction) {
    solver->SFunction(source,u,solver,mpi,t);
    solver->count_sou++;
  }

  /* Apart from other source terms, implement sponge BC as a source */
  DomainBoundary* boundary = (DomainBoundary*) solver->boundary;
  int n;
  int nb = solver->nBoundaryZones;
#if defined(HAVE_CUDA)
  if (solver->use_gpu) {
    for (n = 0; n < nb; n++) {
      if (!strcmp(boundary[n].bctype,_SPONGE_)) {
        fprintf(stderr,"ERROR: Sponge BC not yet implemented on GPU.\n");
      }
    }
  } else {
#endif
    for (n = 0; n < nb; n++) {
      if (!strcmp(boundary[n].bctype,_SPONGE_)) {
        BCSpongeSource( &boundary[n],
                        solver->ndims,
                        solver->nvars,
                        solver->ghosts,
                        solver->dim_local,
                        solver->x,
                        u,
                        source  ); 
      }
    }
#if defined(HAVE_CUDA)
  }
#endif


  return(0);
}

int VolumeIntegral	(	double *	VolumeIntegral,
		double *	u,
		void *	s,
		void *	m
	)

Compute the volume integral of the solution.

Parameters

VolumeIntegral	The computed volume integral
u	Solution
s	Solver object of type HyPar
m	MPI object of type MPIVariables

Definition at line 14 of file VolumeIntegral.c.

{
  HyPar         *solver = (HyPar*)        s;
  MPIVariables  *mpi    = (MPIVariables*) m;

  int ndims   = solver->ndims;
  int nvars   = solver->nvars;
  int *dim    = solver->dim_local;
  int ghosts  = solver->ghosts;
  int d,v;

  /* calculate local volume integral of the solution */
  int index[ndims], done = 0;
  double *local_integral = (double*) calloc (nvars,sizeof(double));
  _ArraySetValue_(local_integral,nvars,0.0);
  _ArraySetValue_(index,ndims,0);
  while (!done) {
    int p; _ArrayIndex1D_(ndims,dim,index,ghosts,p);
    double dxinv[ndims];
    for (d=0; d<ndims; d++) { _GetCoordinate_(d,index[d],dim,ghosts,solver->dxinv,dxinv[d]); }
    double dV = 1.0; for (d=0; d<ndims; d++) dV *= (1.0/dxinv[d]);
    for (v=0; v<nvars; v++) local_integral[v] += (u[p*nvars+v]*dV);
    _ArrayIncrementIndex_(ndims,dim,index,done);
  }
  /* sum over all processors to get global integral of the solution */
  IERR MPISum_double(VolumeIntegral,local_integral,nvars,&mpi->world); CHECKERR(ierr);
  free(local_integral);

  return(0);
}

int BoundaryIntegral	(	void *	s,
		void *	m
	)

Computes the flux integral over the boundary. The local flux integral (on this processor) is computed for physical as well as MPI boundaries. The global boundary integral is computed by summing the local integrals over all the processors, since the contributions from the MPI boundaries cancel out.

Parameters

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

Definition at line 17 of file BoundaryIntegral.c.

{
  HyPar         *solver = (HyPar*)        s;
  MPIVariables  *mpi    = (MPIVariables*) m;

  int ndims   = solver->ndims;
  int nvars   = solver->nvars;
  int *dim    = solver->dim_local;
  int ghosts  = solver->ghosts;
  int d,v,k;

  double *local_integral  = (double*) calloc (nvars,sizeof(double));
  double *global_integral = (double*) calloc (nvars,sizeof(double));

  /* calculate the local boundary integral on each process */
  _ArraySetValue_(local_integral,nvars,0.0);
  for (d=0; d<ndims; d++) {
    for (v=0; v<nvars; v++) {
      double dxinv[ndims], dS = 1.0; 
      for (k=0; k<ndims; k++) { 
        /* uniform grid assumed */
        _GetCoordinate_(k,dim[k]/2,dim,ghosts,solver->dxinv,dxinv[k]); 
      }
      for (k=0; k<ndims; k++) if (k!=d) dS *= (1.0/dxinv[k]);
      local_integral[v] += (solver->StepBoundaryIntegral[(2*d+0)*nvars+v])*dS;
      local_integral[v] += (solver->StepBoundaryIntegral[(2*d+1)*nvars+v])*dS;
    }
  }

  /* add across process to calculate global boundary integral 
   * (internal (MPI) boundaries must cancel out)
   */
  IERR MPISum_double(global_integral,local_integral,nvars,&mpi->world); CHECKERR(ierr);

  /* add to the total boundary integral */
  _ArrayAXPY_(global_integral,1.0,solver->TotalBoundaryIntegral,nvars);

  free(local_integral);
  free(global_integral);
  return(0);
}

int CalculateConservationError	(	void *	s,
		void *	m
	)

Calculates the error (L2) in conservation by computing the difference between the initial volume integral of the solution, and the sum of the current volume integral and the time integral of the boundary flux from the start of the simulation to the current simulation time.

Parameters

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

Definition at line 16 of file CalculateConservationError.c.

{
  HyPar         *solver = (HyPar*) s;
  int           v,nvars = solver->nvars;
  double        error;

  double base[nvars];
  for (v=0; v<nvars; v++) {
    if (absolute(solver->VolumeIntegralInitial[v]) > 1.0)
      base[v] = absolute(solver->VolumeIntegralInitial[v]);
    else base[v] = 1.0;
  }
  
  for (v=0; v<nvars; v++) {
    error =  (solver->VolumeIntegral[v]+solver->TotalBoundaryIntegral[v]-solver->VolumeIntegralInitial[v]) 
           * (solver->VolumeIntegral[v]+solver->TotalBoundaryIntegral[v]-solver->VolumeIntegralInitial[v]);
    solver->ConservationError[v] = sqrt(error)/base[v];
  }
  
  return(0);
}

void IncrementFilenameIndex	(	char *	f,
		int	len
	)

  Increment the output filename of the form op_nnnnn.dat by 1,

i.e., the number represented by the string nnnnn is incremented by 1.

Increment a character string representing an integer by 1. For example: "00002" -> "00003"; "3421934" -> "3421935"; "999" -> "000". The string can be of arbitrary length.

Parameters

f	Character string representing the integer
len	Length of the string

Definition at line 46 of file IncrementFilename.c.

{
  int i;
  for (i=len-1; i>=0; i--) {
    if (f[i] == '9') {
      f[i] = '0';
      if (!i) fprintf(stderr,"Warning: file increment hit max limit. Resetting to zero.\n");
    } else {
      f[i]++;
      break;
    }
  }
}

int NonLinearInterpolation	(	double *	u,
		void *	s,
		void *	m,
		double	t,
		int(*)(double *, double *, int, void *, double)	FluxFunction
	)

Compute the interpolation coefficients of a nonlinear interpolation method, based on the smoothness of the flux function of the solution passed as an argument. This function is provided separately so that these coefficients can be pre-computed and stored for future use.

In the implementation of non-linear methods (such as WENO5), the calculation of the non-linear coefficients, and the actual evaluation of the interpolant are separated into different functions. This provides flexibility on when the nonlinear coefficients are evaluated. Some scenarios are as follows:

• For explicit time integration, it will be computed every time the hyperbolic flux term is being computed.
• For implicit time integration, consistency or linearization may require that the coefficients be computed and "frozen" at the beginning of each stage. Thus, this function can be called at the beginning of each time integration stage, while HyperbolicFunction() is called with the argument LimFlag = 0.

Parameters

u	Solution array
s	Solver object of type HyPar
m	MPI object of type MPIVariables
t	Current solution time
FluxFunction	The flux function \({\bf f}_d\left({\bf u}\right)\), whose properties (typically smoothness) is used to evaluate the nonlinear coefficients

Definition at line 28 of file NonLinearInterpolation.c.

{
  HyPar         *solver = (HyPar*)        s;
  MPIVariables  *mpi    = (MPIVariables*) m;
  double        *FluxC  = solver->fluxC; /* cell centered flux */
  _DECLARE_IERR_;

  int flag = (FluxFunction && solver->flag_nonlinearinterp && solver->SetInterpLimiterVar);
  if (flag) {;
    int     ndims  = solver->ndims;
    int     nvars  = solver->nvars;
    int     ghosts = solver->ghosts;
    int     *dim   = solver->dim_local;
    double  *x     = solver->x;

    int size = 1, d;
    for (d=0; d<ndims; d++) size *= (dim[d] + 2*ghosts);

    /* apply boundary conditions and exchange data over MPI interfaces */
    IERR solver->ApplyBoundaryConditions(solver,mpi,u,NULL,t);      CHECKERR(ierr);
    IERR solver->ApplyIBConditions(solver,mpi,u,t);                 CHECKERR(ierr);
    IERR MPIExchangeBoundariesnD(ndims,nvars,dim,ghosts,mpi,u);     CHECKERR(ierr);

    int offset = 0;
    for (d = 0; d < ndims; d++) {
      /* evaluate cell-centered flux */
      IERR FluxFunction(FluxC,u,d,solver,t);                          CHECKERR(ierr);
      /* calculate non-linear interpolation coefficients */
      IERR solver->SetInterpLimiterVar(FluxC,u,x+offset,d,solver,mpi);CHECKERR(ierr);
      offset += dim[d] + 2*ghosts;
    }
  }

  return(0);
}

int gpuHyperbolicFunction	(	double *	hyp,
		double *	u,
		void *	s,
		void *	m,
		double	t,
		int	LimFlag,
		int(*)(double *, double *, int, void *, double)	FluxFunction,
		int(*)(double *, double *, double *, double *, double *, double *, int, void *, double)	UpwindFunction
	)

This function computes the hyperbolic term of the governing equations, expressed as follows:-

\begin{equation} {\bf F} \left({\bf u}\right) = \sum_{d=0}^{D-1} \frac {\partial {\bf f}_d\left({\bf u}\right)} {\partial x_d} \end{equation}

using a conservative finite-difference discretization is space:

\begin{equation} \hat{\bf F} \left({\bf u}\right) = \sum_{d=0}^{D-1} \frac {1} {\Delta x_d} \left[ \hat{\bf f}_{d,j+1/2} - \hat{\bf f}_{d,j-1/2} \right] \end{equation}

where \(d\) denotes the spatial dimension, \(D\) denotes the total number of spatial dimensions, the hat denotes the discretized quantity, and \(j\) is the grid coordinate along dimension \(d\). The approximation to the flux function \({\bf f}_d\) at the interfaces \(j\pm1/2\), denoted by \(\hat{\bf f}_{d,j\pm 1/2}\), are computed using the function ReconstructHyperbolic().

Parameters

hyp	Array to hold the computed hyperbolic term (shares the same layout as u
u	Solution array
s	Solver object of type HyPar
m	MPI object of type MPIVariables
t	Current simulation time
LimFlag	Flag to indicate if the nonlinear coefficients for solution-dependent interpolation method should be recomputed (see ReconstructHyperbolic() for an explanation on why this is needed)
FluxFunction	Function pointer to the flux function for the hyperbolic term
UpwindFunction	Function pointer to the upwinding function for the hyperbolic term

Definition at line 473 of file HyperbolicFunction_GPU.cu.

{
    HyPar         *solver = (HyPar*)        s;
    MPIVariables  *mpi    = (MPIVariables*) m;
    int           d;
    double        *FluxI  = solver->fluxI; /* interface flux     */
    double        *FluxC  = solver->fluxC; /* cell centered flux */

    int     ndims     = solver->ndims;
    int     nvars     = solver->nvars;
    int     ghosts    = solver->ghosts;
    int     *dim      = solver->dim_local;
    int     size      = solver->npoints_local_wghosts;
    double  *x        = solver->gpu_x;
    double  *dxinv    = solver->gpu_dxinv;
    double  *FluxI_p1, *FluxI_p2;

    LimFlag = (LimFlag && solver->flag_nonlinearinterp && solver->SetInterpLimiterVar);

    gpuMemset(hyp, 0, size*nvars*sizeof(double));
    gpuMemset(solver->StageBoundaryBuffer, 0, solver->StageBoundaryBuffer_size*sizeof(double));
    gpuMemset(solver->StageBoundaryIntegral, 0, 2*ndims*nvars*sizeof(double));

    if (!FluxFunction) return(0); /* zero hyperbolic term */
    /*solver->count_hyp++;*/

    int npoints_grid = 1; for (d = 0; d < ndims; d++) npoints_grid *= dim[d];
    int nblocks = (npoints_grid - 1) / GPU_THREADS_PER_BLOCK + 1;
    int offset = 0;

#if defined(GPU_STAT)
    cudaEvent_t start, stop;
    float milliseconds = 0;

    checkCuda( cudaEventCreate(&start) );
    checkCuda( cudaEventCreate(&stop) );
#endif

    for (d = 0; d < ndims; d++) {
        int npoints_dim_interface = 1; for (int i=0; i<ndims; i++) npoints_dim_interface *= (i==d) ? (dim[i]+1) : dim[i];
        /* evaluate cell-centered flux */
        FluxFunction(FluxC,u,d,solver,t);
        /* compute interface fluxes */
        ReconstructHyperbolic(FluxI,FluxC,u,x+offset,d,solver,mpi,t,LimFlag,UpwindFunction);

        FluxI_p1 = solver->StageBoundaryBuffer+(solver->gpu_npoints_boundary_offset[d]*nvars);
        FluxI_p2 = solver->StageBoundaryBuffer+(solver->gpu_npoints_boundary_offset[d]+solver->gpu_npoints_boundary[d])*nvars;

#if defined(GPU_STAT)
        checkCuda( cudaEventRecord(start, 0) );
#endif
        if (ndims == 3) {
          HyperbolicFunction_dim3_kernel<<<nblocks, GPU_THREADS_PER_BLOCK>>>(
              npoints_grid, size, npoints_dim_interface, d, nvars, ghosts, offset,
              solver->gpu_dim_local, FluxI, dxinv,
              hyp, FluxI_p1, FluxI_p2
          );
        } else {
          fprintf(stderr,"gpuHyperbolicFunction(): Not implemented for ndims = %d!\n", ndims);
          exit(1);
        }
        cudaDeviceSynchronize();

#if defined(GPU_STAT)
        checkCuda( cudaEventRecord(stop, 0) );
        checkCuda( cudaEventSynchronize(stop) );
        checkCuda( cudaEventElapsedTime(&milliseconds, start, stop) );
        printf("%-50s GPU time = %.6f dir = %d bandwidth (GB/s) = %.2f\n",
                "HyperbolicFunction", milliseconds*1e-3, d, (1e-6*2*npoints_grid*nvars*sizeof(double))/milliseconds);

        checkCuda( cudaEventRecord(start, 0) );
#endif

        StageBoundaryIntegral_kernel<GPU_THREADS_PER_BLOCK><<<1, GPU_THREADS_PER_BLOCK, GPU_THREADS_PER_BLOCK*nvars*sizeof(double)>>>(
          solver->gpu_npoints_boundary[d], nvars, -1, FluxI_p1,
          solver->StageBoundaryIntegral + 2*d*nvars
        );
        StageBoundaryIntegral_kernel<GPU_THREADS_PER_BLOCK><<<1, GPU_THREADS_PER_BLOCK, GPU_THREADS_PER_BLOCK*nvars*sizeof(double)>>>(
          solver->gpu_npoints_boundary[d], nvars, 1, FluxI_p2,
          solver->StageBoundaryIntegral + (2*d+1)*nvars
        );
        cudaDeviceSynchronize();

#if defined(GPU_STAT)
        checkCuda( cudaEventRecord(stop, 0) );
        checkCuda( cudaEventSynchronize(stop) );
        checkCuda( cudaEventElapsedTime(&milliseconds, start, stop) );
        printf("%-50s GPU time = %.6f dir = %d bandwidth (GB/s) = %.2f\n",
                "StageBoundaryIntegral", milliseconds*1e-3, d, (1e-6*2*solver->gpu_npoints_boundary[d]*nvars*sizeof(double))/milliseconds);
#endif
        offset += dim[d] + 2*ghosts;
    }

    if (solver->flag_ib) gpuArrayBlockMultiply(hyp, solver->gpu_iblank, size, nvars);

#if defined(GPU_STAT)
    checkCuda(cudaEventDestroy(start));
    checkCuda(cudaEventDestroy(stop));
#endif

    return(0);
}

int InitializeSolvers	(	void *	s,
		int	nsims
	)

This function initializes all solvers-specific function pointers depending on user input. The specific functions used for spatial discretization, time integration, and solution output are set here.

Parameters

s	Array of simulation objects of type SimulationObject
nsims	number of simulation objects

Definition at line 52 of file InitializeSolvers.c.

{
  SimulationObject *sim = (SimulationObject*) s;
  int ns;
  _DECLARE_IERR_;

  if (nsims == 0) return 0;

  if (!sim[0].mpi.rank) {
    printf("Initializing solvers.\n");
  }

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

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

    solver->ApplyBoundaryConditions = ApplyBoundaryConditions;
    solver->ApplyIBConditions = ApplyIBConditions;
    solver->SourceFunction = SourceFunction;
#if defined(HAVE_CUDA)
    if (solver->use_gpu) {
      solver->HyperbolicFunction = gpuHyperbolicFunction;
    } else {
#endif
      solver->HyperbolicFunction = HyperbolicFunction;
#if defined(HAVE_CUDA)
    }
#endif
    solver->VolumeIntegralFunction      = VolumeIntegral;
    solver->BoundaryIntegralFunction    = BoundaryIntegral;
    solver->CalculateConservationError  = CalculateConservationError;
    solver->NonlinearInterp             = NonLinearInterpolation;

    /* choose the type of parabolic discretization */
    solver->ParabolicFunction         = NULL;
    solver->SecondDerivativePar       = NULL;
    solver->FirstDerivativePar        = NULL;
    solver->InterpolateInterfacesPar  = NULL;

#if defined(HAVE_CUDA)
    if (solver->use_gpu) {

      if (!strcmp(solver->spatial_type_par,_NC_2STAGE_)) {
  
        solver->ParabolicFunction = ParabolicFunctionNC2Stage;

        if (!strcmp(solver->spatial_scheme_par,_FOURTH_ORDER_CENTRAL_)) {
          solver->FirstDerivativePar = gpuFirstDerivativeFourthOrderCentral;
        } else {
          fprintf(stderr,"ERROR (domain %d): scheme %s is not supported on GPU!",
                  ns, solver->spatial_scheme_par);
          return 1;
        }
  
      }

    } else {
#endif

      if (!strcmp(solver->spatial_type_par,_NC_1STAGE_)) {
  
        solver->ParabolicFunction = ParabolicFunctionNC1Stage;
        if (!strcmp(solver->spatial_scheme_par,_SECOND_ORDER_CENTRAL_)) {
          solver->SecondDerivativePar      = SecondDerivativeSecondOrderCentral;
        } else if (!strcmp(solver->spatial_scheme_par,_FOURTH_ORDER_CENTRAL_)) {
          solver->SecondDerivativePar      = SecondDerivativeFourthOrderCentral;
        } else {
          fprintf(stderr,"Error (domain %d): %s is not a supported ",
                  ns, solver->spatial_scheme_par);
          fprintf(stderr,"spatial scheme of type %s for the parabolic terms.\n",
                  solver->spatial_type_par);
        }
  
      } else if (!strcmp(solver->spatial_type_par,_NC_2STAGE_)) {
  
        solver->ParabolicFunction = ParabolicFunctionNC2Stage;
        if (!strcmp(solver->spatial_scheme_par,_SECOND_ORDER_CENTRAL_)) {
          solver->FirstDerivativePar       = FirstDerivativeFirstOrder;
          /* why first order? see ParabolicFunctionNC2Stage.c. 2nd order central
             approximation to the 2nd derivative can be expressed as a conjugation
             of 1st order approximations to the 1st derivative (one forward and
             one backward) -- this prevents odd-even decoupling */
        } else if (!strcmp(solver->spatial_scheme_par,_FOURTH_ORDER_CENTRAL_)) {
          solver->FirstDerivativePar       = FirstDerivativeFourthOrderCentral;
          /* why 4th order? I could not derive the decomposition of the
             4th order central approximation to the 2nd derivative! Some problems
             may show odd-even decoupling */
        } else {
          fprintf(stderr,"Error (domain %d): %s is not a supported ",
                  ns, solver->spatial_scheme_par);
          fprintf(stderr,"spatial scheme of type %s for the parabolic terms.\n",
                solver->spatial_type_par);
        }
  
      } else if (!strcmp(solver->spatial_type_par,_NC_1_5STAGE_)) {
  
        solver->ParabolicFunction = ParabolicFunctionNC1_5Stage;
        if (!strcmp(solver->spatial_scheme_par,_SECOND_ORDER_CENTRAL_)) {
          solver->FirstDerivativePar       = FirstDerivativeSecondOrderCentral;
          solver->SecondDerivativePar      = SecondDerivativeSecondOrderCentral;
        } else if (!strcmp(solver->spatial_scheme_par,_FOURTH_ORDER_CENTRAL_)) {
          solver->FirstDerivativePar       = FirstDerivativeFourthOrderCentral;
          solver->SecondDerivativePar      = SecondDerivativeFourthOrderCentral;
        } else {
          fprintf(stderr,"Error (domain %d): %s is not a supported ",
                  ns, solver->spatial_scheme_par);
          fprintf(stderr,"spatial scheme of type %s for the parabolic terms.\n",
                solver->spatial_type_par);
        }
  
      } else if (!strcmp(solver->spatial_type_par,_CONS_1STAGE_)) {
  
        solver->ParabolicFunction = ParabolicFunctionCons1Stage;
        if (!strcmp(solver->spatial_scheme_par,_SECOND_ORDER_CENTRAL_)) {
          solver->InterpolateInterfacesPar = Interp2PrimSecondOrder;
        } else {
          fprintf(stderr,"Error (domain %d): %s is not a supported ",
                  ns, solver->spatial_scheme_par);
          fprintf(stderr,"spatial scheme of type %s for the parabolic terms.\n",
                solver->spatial_type_par);
        }
  
      } else {
  
        fprintf(stderr,"Error (domain %d): %s is not a supported ",
                ns, solver->spatial_type_par);
        fprintf(stderr,"spatial discretization type for the parabolic terms.\n");
        return(1);
  
      }

#if defined(HAVE_CUDA)
    }
#endif

    /* Spatial interpolation for hyperbolic term */
    solver->interp                = NULL;
    solver->compact               = NULL;
    solver->lusolver              = NULL;
    solver->SetInterpLimiterVar   = NULL;
    solver->flag_nonlinearinterp  = 1;
    if (strcmp(solver->interp_type,_CHARACTERISTIC_) && strcmp(solver->interp_type,_COMPONENTS_)) {
      fprintf(stderr,"Error in InitializeSolvers() (domain %d): %s is not a ",
              ns, solver->interp_type);
      fprintf(stderr,"supported interpolation type.\n");
      return(1);
    }

#if defined(HAVE_CUDA)
    if (solver->use_gpu) {

      if (!strcmp(solver->spatial_scheme_hyp,_FIFTH_ORDER_WENO_)) {
  
        /* Fifth order WENO scheme */
        if ((solver->nvars > 1) && (!strcmp(solver->interp_type,_CHARACTERISTIC_))) {
          fprintf(stderr,
                  "Error (domain %d): characteristic-based WENO5 is not yet implemented on GPUs.\n",
                  ns );
          return 1;
        } else {
          solver->InterpolateInterfacesHyp = gpuInterp1PrimFifthOrderWENO;
        }
        solver->interp = (WENOParameters*) calloc(1,sizeof(WENOParameters));
        IERR WENOInitialize(solver,mpi,solver->spatial_scheme_hyp,solver->interp_type); CHECKERR(ierr);
        solver->flag_nonlinearinterp = !(((WENOParameters*)solver->interp)->no_limiting);
  
      } else {

        fprintf(stderr,
                "Error (domain %d): %s is a not a supported spatial interpolation scheme on GPUs.\n",
                ns, solver->spatial_scheme_hyp);
        return 1;
      }

    } else {
#endif

      if (!strcmp(solver->spatial_scheme_hyp,_FIRST_ORDER_UPWIND_)) {
  
        /* First order upwind scheme */
        if ((solver->nvars > 1) && (!strcmp(solver->interp_type,_CHARACTERISTIC_))) {
          solver->InterpolateInterfacesHyp = Interp1PrimFirstOrderUpwindChar;
        } else {
          solver->InterpolateInterfacesHyp = Interp1PrimFirstOrderUpwind;
        }
  
      } else if (!strcmp(solver->spatial_scheme_hyp,_SECOND_ORDER_CENTRAL_)) {
  
        /* Second order central scheme */
        if ((solver->nvars > 1) && (!strcmp(solver->interp_type,_CHARACTERISTIC_))) {
          solver->InterpolateInterfacesHyp = Interp1PrimSecondOrderCentralChar;
        } else {
          solver->InterpolateInterfacesHyp = Interp1PrimSecondOrderCentral;
        }
  
      } else if (!strcmp(solver->spatial_scheme_hyp,_SECOND_ORDER_MUSCL_)) {
  
        /* Second order MUSCL scheme */
        if ((solver->nvars > 1) && (!strcmp(solver->interp_type,_CHARACTERISTIC_))) {
          solver->InterpolateInterfacesHyp = Interp1PrimSecondOrderMUSCLChar;
        } else {
          solver->InterpolateInterfacesHyp = Interp1PrimSecondOrderMUSCL;
        }
        solver->interp = (MUSCLParameters*) calloc(1,sizeof(MUSCLParameters));
        IERR MUSCLInitialize(solver,mpi); CHECKERR(ierr);
  
      } else if (!strcmp(solver->spatial_scheme_hyp,_THIRD_ORDER_MUSCL_)) {
  
        /* Third order MUSCL scheme */
        if ((solver->nvars > 1) && (!strcmp(solver->interp_type,_CHARACTERISTIC_))) {
          solver->InterpolateInterfacesHyp = Interp1PrimThirdOrderMUSCLChar;
        } else {
          solver->InterpolateInterfacesHyp = Interp1PrimThirdOrderMUSCL;
        }
        solver->interp = (MUSCLParameters*) calloc(1,sizeof(MUSCLParameters));
        IERR MUSCLInitialize(solver,mpi); CHECKERR(ierr);
  
      } else if (!strcmp(solver->spatial_scheme_hyp,_FOURTH_ORDER_CENTRAL_)) {
  
        /* Fourth order central scheme */
        if ((solver->nvars > 1) && (!strcmp(solver->interp_type,_CHARACTERISTIC_))) {
          solver->InterpolateInterfacesHyp = Interp1PrimFourthOrderCentralChar;
        } else {
          solver->InterpolateInterfacesHyp = Interp1PrimFourthOrderCentral;
        }
  
      } else if (!strcmp(solver->spatial_scheme_hyp,_FIFTH_ORDER_UPWIND_)) {
  
        /* Fifth order upwind scheme */
        if ((solver->nvars > 1) && (!strcmp(solver->interp_type,_CHARACTERISTIC_))) {
          solver->InterpolateInterfacesHyp = Interp1PrimFifthOrderUpwindChar;
        } else {
          solver->InterpolateInterfacesHyp = Interp1PrimFifthOrderUpwind;
        }
  
      } else if (!strcmp(solver->spatial_scheme_hyp,_FIFTH_ORDER_COMPACT_UPWIND_)) {
  
        /* Fifth order compact upwind scheme */
        if ((solver->nvars > 1) && (!strcmp(solver->interp_type,_CHARACTERISTIC_))) {
          solver->InterpolateInterfacesHyp = Interp1PrimFifthOrderCompactUpwindChar;
        } else {
          solver->InterpolateInterfacesHyp = Interp1PrimFifthOrderCompactUpwind;
        }
        solver->compact = (CompactScheme*) calloc(1,sizeof(CompactScheme));
        IERR CompactSchemeInitialize(solver,mpi,solver->interp_type);
        solver->lusolver = (TridiagLU*) calloc (1,sizeof(TridiagLU));
        IERR tridiagLUInit(solver->lusolver,&mpi->world);CHECKERR(ierr);
  
      } else if (!strcmp(solver->spatial_scheme_hyp,_FIFTH_ORDER_WENO_)) {
  
        /* Fifth order WENO scheme */
        if ((solver->nvars > 1) && (!strcmp(solver->interp_type,_CHARACTERISTIC_))) {
          solver->InterpolateInterfacesHyp = Interp1PrimFifthOrderWENOChar;
        } else {
          solver->InterpolateInterfacesHyp = Interp1PrimFifthOrderWENO;
        }
        solver->interp = (WENOParameters*) calloc(1,sizeof(WENOParameters));
        IERR WENOInitialize(solver,mpi,solver->spatial_scheme_hyp,solver->interp_type); CHECKERR(ierr);
        solver->flag_nonlinearinterp = !(((WENOParameters*)solver->interp)->no_limiting);
  
      } else if (!strcmp(solver->spatial_scheme_hyp,_FIFTH_ORDER_CRWENO_)) {
  
        /* Fifth order CRWENO scheme */
        if ((solver->nvars > 1) && (!strcmp(solver->interp_type,_CHARACTERISTIC_))) {
          solver->InterpolateInterfacesHyp = Interp1PrimFifthOrderCRWENOChar;
        } else {
          solver->InterpolateInterfacesHyp = Interp1PrimFifthOrderCRWENO;
        }
        solver->interp = (WENOParameters*) calloc(1,sizeof(WENOParameters));
        IERR WENOInitialize(solver,mpi,solver->spatial_scheme_hyp,solver->interp_type); CHECKERR(ierr);
        solver->flag_nonlinearinterp = !(((WENOParameters*)solver->interp)->no_limiting);
        solver->compact = (CompactScheme*) calloc(1,sizeof(CompactScheme));
        IERR CompactSchemeInitialize(solver,mpi,solver->interp_type);
        solver->lusolver = (TridiagLU*) calloc (1,sizeof(TridiagLU));
        IERR tridiagLUInit(solver->lusolver,&mpi->world);CHECKERR(ierr);
  
      } else if (!strcmp(solver->spatial_scheme_hyp,_FIFTH_ORDER_HCWENO_)) {
  
        /* Fifth order HCWENO scheme */
        if ((solver->nvars > 1) && (!strcmp(solver->interp_type,_CHARACTERISTIC_))) {
          solver->InterpolateInterfacesHyp = Interp1PrimFifthOrderHCWENOChar;
        } else {
          solver->InterpolateInterfacesHyp = Interp1PrimFifthOrderHCWENO;
        }
        solver->interp = (WENOParameters*) calloc(1,sizeof(WENOParameters));
        IERR WENOInitialize(solver,mpi,solver->spatial_scheme_hyp,solver->interp_type); CHECKERR(ierr);
        solver->flag_nonlinearinterp = !(((WENOParameters*)solver->interp)->no_limiting);
        solver->compact = (CompactScheme*) calloc(1,sizeof(CompactScheme));
        IERR CompactSchemeInitialize(solver,mpi,solver->interp_type);
        solver->lusolver = (TridiagLU*) calloc (1,sizeof(TridiagLU));
        IERR tridiagLUInit(solver->lusolver,&mpi->world);CHECKERR(ierr);
  
      } else {

        fprintf(stderr,"Error (domain %d): %s is a not a supported spatial interpolation scheme.\n",
                ns, solver->spatial_scheme_hyp);
        return(1);
      }

#if defined(HAVE_CUDA)
    }
#endif

    /* Time integration */
    solver->time_integrator = NULL;
#ifdef with_petsc
    if (solver->use_petscTS) {
      /* dummy -- not used */
      solver->TimeIntegrate = TimeForwardEuler;
      solver->msti = NULL;
    } else {
      if (!strcmp(solver->time_scheme,_FORWARD_EULER_)) {
        solver->TimeIntegrate = TimeForwardEuler;
        solver->msti = NULL;
      } else if (!strcmp(solver->time_scheme,_RK_)) {
        solver->TimeIntegrate = TimeRK;
        solver->msti = (ExplicitRKParameters*) calloc (1,sizeof(ExplicitRKParameters));
        IERR TimeExplicitRKInitialize(solver->time_scheme,solver->time_scheme_type,
                                      solver->msti,mpi); CHECKERR(ierr);
      } else if (!strcmp(solver->time_scheme,_GLM_GEE_)) {
        solver->TimeIntegrate = TimeGLMGEE;
        solver->msti = (GLMGEEParameters*) calloc (1,sizeof(GLMGEEParameters));
        IERR TimeGLMGEEInitialize(solver->time_scheme,solver->time_scheme_type,
                                  solver->msti,mpi); CHECKERR(ierr);
      } else {
        fprintf(stderr,"Error (domain %d): %s is a not a supported time-integration scheme.\n",
                ns, solver->time_scheme);
        return(1);
      }
    }
#else
    if (!strcmp(solver->time_scheme,_FORWARD_EULER_)) {
      solver->TimeIntegrate = TimeForwardEuler;
      solver->msti = NULL;
    } else if (!strcmp(solver->time_scheme,_RK_)) {
      solver->TimeIntegrate = TimeRK;
      solver->msti = (ExplicitRKParameters*) calloc (1,sizeof(ExplicitRKParameters));
      IERR TimeExplicitRKInitialize(solver->time_scheme,solver->time_scheme_type,
                                    solver->msti,mpi); CHECKERR(ierr);
    } else if (!strcmp(solver->time_scheme,_GLM_GEE_)) {
      solver->TimeIntegrate = TimeGLMGEE;
      solver->msti = (GLMGEEParameters*) calloc (1,sizeof(GLMGEEParameters));
      IERR TimeGLMGEEInitialize(solver->time_scheme,solver->time_scheme_type,
                                solver->msti,mpi); CHECKERR(ierr);
    } else {
      fprintf(stderr,"Error (domain %d): %s is a not a supported time-integration scheme.\n",
              ns, solver->time_scheme);
      return(1);
    }
#endif

    /* Solution output function */
    solver->WriteOutput    = NULL; /* default - no output */
    solver->filename_index = NULL;
    strcpy(solver->op_fname_root, "op");
#ifdef with_librom
    strcpy(solver->op_rom_fname_root, "op_rom");
#endif
    strcpy(solver->aux_op_fname_root, "ts0");
    if (!strcmp(solver->output_mode,"serial")) {
      solver->index_length = 5;
      solver->filename_index = (char*) calloc (solver->index_length+1,sizeof(char));
      int i; for (i=0; i<solver->index_length; i++) solver->filename_index[i] = '0';
      solver->filename_index[solver->index_length] = (char) 0;
      if (!strcmp(solver->op_file_format,"text")) {
        solver->WriteOutput = WriteText;
        strcpy(solver->solnfilename_extn,".dat");
      } else if (!strcmp(solver->op_file_format,"tecplot2d")) {
        solver->WriteOutput = WriteTecplot2D;
        strcpy(solver->solnfilename_extn,".dat");
      } else if (!strcmp(solver->op_file_format,"tecplot3d")) {
        solver->WriteOutput = WriteTecplot3D;
        strcpy(solver->solnfilename_extn,".dat");
      } else if ((!strcmp(solver->op_file_format,"binary")) || (!strcmp(solver->op_file_format,"bin"))) {
        solver->WriteOutput = WriteBinary;
        strcpy(solver->solnfilename_extn,".bin");
      } else if (!strcmp(solver->op_file_format,"none")) {
        solver->WriteOutput = NULL;
      } else {
        fprintf(stderr,"Error (domain %d): %s is not a supported file format.\n",
                ns, solver->op_file_format);
        return(1);
      }
      if ((!strcmp(solver->op_overwrite,"no")) && solver->restart_iter) {
        /* if it's a restart run, fast-forward the filename */
        int t;
        for (t=0; t<solver->restart_iter; t++)
          if ((t+1)%solver->file_op_iter == 0) IncrementFilenameIndex(solver->filename_index,solver->index_length);
      }
    } else if (!strcmp(solver->output_mode,"parallel")) {
      if (!strcmp(solver->op_file_format,"none")) solver->WriteOutput = NULL;
      else {
        /* only binary file writing supported in parallel mode */
        /* use post-processing scripts to convert              */
        solver->WriteOutput = WriteBinary;
        strcpy(solver->solnfilename_extn,".bin");
      }
    } else {
      fprintf(stderr,"Error (domain %d): %s is not a supported output mode.\n",
              ns, solver->output_mode);
      fprintf(stderr,"Should be \"serial\" or \"parallel\".    \n");
      return(1);
    }

    /* Solution plotting function */
    strcpy(solver->plotfilename_extn,".png");
#ifdef with_python
    solver->py_plt_func = NULL;
    solver->py_plt_func_args = NULL;
    {
      char python_plotting_fname[_MAX_STRING_SIZE_] = "plotSolution";
      PyObject* py_plot_name = PyUnicode_DecodeFSDefault(python_plotting_fname);
      PyObject* py_plot_module = PyImport_Import(py_plot_name);
      Py_DECREF(py_plot_name);
      if (py_plot_module) {
        solver->py_plt_func = PyObject_GetAttrString(py_plot_module, "plotSolution");
        if (!solver->py_plt_func) {
          if (!mpi->rank) {
            printf("Unable to load plotSolution function from Python module.\n");
          }
        } else {
          if (!mpi->rank) {
            printf("Loaded Python module for plotting.\n");
            printf("Loaded plotSolution function from Python module.\n");
          }
        }
      } else {
        if (!mpi->rank) {
          printf("Unable to load Python module for plotting.\n");
        }
      }
    }
#endif

  }

  return(0);
}