Euler1DInitialize.c File Reference

Initialize the 1D Euler equations module. More...

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <basic.h>
#include <arrayfunctions.h>
#include <physicalmodels/euler1d.h>
#include <mpivars.h>
#include <hypar.h>

Go to the source code of this file.

Functions
double	Euler1DComputeCFL (void *, void *, double, double)
int	Euler1DFlux (double *, double *, int, void *, double)
int	Euler1DStiffFlux (double *, double *, int, void *, double)
int	Euler1DSource (double *, double *, void *, void *, double)
int	Euler1DUpwindRoe (double *, double *, double *, double *, double *, double *, int, void *, double)
int	Euler1DUpwinddFRoe (double *, double *, double *, double *, double *, double *, int, void *, double)
int	Euler1DUpwindRF (double *, double *, double *, double *, double *, double *, int, void *, double)
int	Euler1DUpwinddFRF (double *, double *, double *, double *, double *, double *, int, void *, double)
int	Euler1DUpwindLLF (double *, double *, double *, double *, double *, double *, int, void *, double)
int	Euler1DUpwinddFLLF (double *, double *, double *, double *, double *, double *, int, void *, double)
int	Euler1DUpwindSWFS (double *, double *, double *, double *, double *, double *, int, void *, double)
int	Euler1DUpwindRusanov (double *, double *, double *, double *, double *, double *, int, void *, double)
int	Euler1DJacobian (double *, double *, void *, int, int, int)
int	Euler1DStiffJacobian (double *, double *, void *, int, int, int)
int	Euler1DRoeAverage (double *, double *, double *, void *)
int	Euler1DLeftEigenvectors (double *, double *, void *, int)
int	Euler1DRightEigenvectors (double *, double *, void *, int)
int	Euler1DGravityField (void *, void *)
int	Euler1DSourceUpwindLLF (double *, double *, double *, double *, int, void *, double)
int	Euler1DSourceUpwindRoe (double *, double *, double *, double *, int, void *, double)
int	Euler1DModifiedSolution (double *, double *, int, void *, void *, double)
int	Euler1DPreStep (double *, void *, void *, double)
int	Euler1DInitialize (void *s, void *m)

Detailed Description

Initialize the 1D Euler equations module.

Author

Debojyoti Ghosh

Definition in file Euler1DInitialize.c.

Function Documentation

Euler1DComputeCFL()

double Euler1DComputeCFL	(	void *	s,
		void *	m,
		double	dt,
		double	t
	)

Computes the maximum CFL number over the domain. Note that the CFL is computed over the local domain on this processor only.

Parameters

s	Solver object of type HyPar
m	MPI object of type MPIVariables
dt	Time step size for which to compute the CFL
t	Time

Definition at line 17 of file Euler1DComputeCFL.c.

{
  HyPar     *solver = (HyPar*)   s;
  Euler1D   *param  = (Euler1D*) solver->physics;

  int *dim    = solver->dim_local;
  int ghosts  = solver->ghosts;
  int ndims   = solver->ndims;
  int index[ndims];
  double *u   = solver->u;

  double max_cfl = 0;
  int done = 0; _ArraySetValue_(index,ndims,0);
  while (!done) {
    int p; _ArrayIndex1D_(ndims,dim,index,ghosts,p);
    double rho, v, e, P, c, dxinv, local_cfl;
    _Euler1DGetFlowVar_((u+_MODEL_NVARS_*p),rho,v,e,P,param);
    _GetCoordinate_(0,index[0],dim,ghosts,solver->dxinv,dxinv); /* 1/dx */
    c         = sqrt(param->gamma*P/rho); /* speed of sound */
    local_cfl = (absolute(v)+c)*dt*dxinv; /* local cfl for this grid point */
    if (local_cfl > max_cfl) max_cfl = local_cfl;
    _ArrayIncrementIndex_(ndims,dim,index,done);
  }

  return(max_cfl);
}

Euler1DFlux()

int Euler1DFlux	(	double *	f,
		double *	u,
		int	dir,
		void *	s,
		double	t
	)

Compute the hyperbolic flux over the local domain.

\begin{equation} {\bf F}\left({\bf u}\right) = \left[\begin{array}{c} \rho u \\ \rho u^2 + p \\ (e+p) u \end{array}\right] \end{equation}

Parameters

f	Array to hold the computed flux (same size and layout as u)
u	Array containing the conserved solution
dir	Spatial dimension (unused since this is a 1D system)
s	Solver object of type HyPar
t	Current time

Definition at line 16 of file Euler1DFlux.c.

{
  HyPar             *solver = (HyPar*)   s;
  Euler1D           *param  = (Euler1D*) solver->physics;
  int               *dim    = solver->dim_local;
  int               ghosts  = solver->ghosts;
  static const int  ndims   = _MODEL_NDIMS_;
  static const int  nvars   = _MODEL_NVARS_;
  static int        index[_MODEL_NDIMS_], bounds[_MODEL_NDIMS_], offset[_MODEL_NDIMS_];

  /* set bounds for array index to include ghost points */
  _ArrayAddCopy1D_(dim,(2*ghosts),bounds,ndims);
  /* set offset such that index is compatible with ghost point arrangement */
  _ArraySetValue_(offset,ndims,-ghosts);

  int done = 0; _ArraySetValue_(index,ndims,0);
  while (!done) {
    int p; _ArrayIndex1DWO_(ndims,dim,index,offset,ghosts,p);
    double rho, v, e, P;
    _Euler1DGetFlowVar_((u+nvars*p),rho,v,e,P,param);
    _Euler1DSetFlux_((f+nvars*p),rho,v,e,P);
    _ArrayIncrementIndex_(ndims,bounds,index,done);
  }

  return(0);
}

Euler1DStiffFlux()

int Euler1DStiffFlux	(	double *	f,
		double *	u,
		int	dir,
		void *	s,
		double	t
	)

Compute the stiff component of the hyperbolic flux over the local domain.

\begin{equation} {\bf F}_F\left({\bf u}\right) = A_f\left({\bf u}\right){\bf u} \end{equation}

where \(A_f\left({\bf u}\right) \) is the fast Jacobian representing the acoustic waves only. A linearized formulation is used where the fast Jacobian \(A_f\) is computed for the solution at the beginning of each time step in Euler1DPreStep.

See also

_Euler1DSetStiffFlux_, _Euler1DSetLinearizedStiffFlux_, _Euler1DSetStiffJac_

Reference:

• Ghosh, D., Constantinescu, E. M., Semi-Implicit Time Integration of Atmospheric Flows with Characteristic-Based Flux Partitioning, SIAM Journal on Scientific Computing, 38 (3), 2016, A1848-A1875, http://dx.doi.org/10.1137/15M1044369

Parameters

f	Array to hold the computed flux (same size and layout as u)
u	Array containing the conserved solution
dir	Spatial dimension (unused since this is a 1D system)
s	Solver object of type HyPar
t	Current time

Definition at line 64 of file Euler1DFlux.c.

{
  HyPar             *solver = (HyPar*)   s;
  Euler1D           *param  = (Euler1D*) solver->physics;
  int               *dim    = solver->dim_local;
  int               ghosts  = solver->ghosts;
  static const int  ndims   = _MODEL_NDIMS_;
  static const int  nvars   = _MODEL_NVARS_;
  static int        index[_MODEL_NDIMS_], bounds[_MODEL_NDIMS_], offset[_MODEL_NDIMS_];

  /* set bounds for array index to include ghost points */
  _ArrayAddCopy1D_(dim,(2*ghosts),bounds,ndims);
  /* set offset such that index is compatible with ghost point arrangement */
  _ArraySetValue_(offset,ndims,-ghosts);

  int done = 0; _ArraySetValue_(index,ndims,0);
  while (!done) {
    int p; _ArrayIndex1DWO_(ndims,dim,index,offset,ghosts,p);
    _Euler1DSetLinearizedStiffFlux_((f+nvars*p),(u+nvars*p),(param->fast_jac+nvars*nvars*p));
    _ArrayIncrementIndex_(ndims,bounds,index,done);
  }

  return(0);
}

Euler1DSource()

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

Compute the gravitational source terms for the 1D Euler equations. The source term is computed according to the balanced formulation introduced in the reference below. The source term is reformulated and "discretized" in a similar fashion as the hyperbolic flux to ensure that the hydrostatic balance is maintained to machine precision.

• Xing, Shu, "High Order Well-Balanced WENO Scheme for the Gas Dynamics Equations Under Gravitational Fields", J. Sci. Comput., 54, 2013, pp. 645–662, http://dx.doi.org/10.1007/s10915-012-9585-8.

Parameters

source	Computed source terms (array size & layout same as u)
u	Solution (conserved variables)
s	Solver object of type HyPar
m	MPI object of type MPIVariables
t	Current solution time

Definition at line 23 of file Euler1DSource.c.

{
  HyPar         *solver = (HyPar* ) s;
  MPIVariables  *mpi = (MPIVariables*) m;
  Euler1D       *param  = (Euler1D*) solver->physics;

  if (param->grav == 0.0)  return(0); /* no gravitational forces */

  int     v, done, p, p1, p2;
  double  *SourceI = solver->fluxI; /* interace source term       */
  double  *SourceC = solver->fluxC; /* cell-centered source term  */
  double  *SourceL = solver->fL;
  double  *SourceR = solver->fR;

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

  /* set interface dimensions */
  _ArrayCopy1D_(dim,dim_interface,ndims); dim_interface[_XDIR_]++;
  /* calculate the split source function exp(-phi/RT) */
  IERR Euler1DSourceFunction(SourceC,u,x,solver,mpi,t); CHECKERR(ierr);
  /* calculate the left and right interface source terms */
  IERR solver->InterpolateInterfacesHyp(SourceL,SourceC,u,x, 1,_XDIR_,solver,mpi,0); CHECKERR(ierr);
  IERR solver->InterpolateInterfacesHyp(SourceR,SourceC,u,x,-1,_XDIR_,solver,mpi,0); CHECKERR(ierr);
  /* calculate the final interface source term */
  IERR param->SourceUpwind(SourceI,SourceL,SourceR,u,_XDIR_,solver,t);
  /* calculate the final cell-centered source term */
  done = 0; _ArraySetValue_(index,ndims,0);
  while (!done) {
    _ArrayCopy1D_(index,index1,ndims);
    _ArrayCopy1D_(index,index2,ndims); index2[_XDIR_]++;
    _ArrayIndex1D_(ndims,dim          ,index ,ghosts,p );
    _ArrayIndex1D_(ndims,dim_interface,index1,0     ,p1);
    _ArrayIndex1D_(ndims,dim_interface,index2,0     ,p2);
    double dx_inverse;   _GetCoordinate_(_XDIR_,index[_XDIR_],dim,ghosts,dxinv,dx_inverse);
    double rho, vel, e, P; _Euler1DGetFlowVar_((u+_MODEL_NVARS_*p),rho,vel,e,P,param);
    double term[_MODEL_NVARS_] = {0.0, rho, rho*vel};
    for (v=0; v<_MODEL_NVARS_; v++) {
      source[_MODEL_NVARS_*p+v] += (  (term[v]*(1.0/param->grav_field[p])) 
                                    * (SourceI[_MODEL_NVARS_*p2+v]-SourceI[_MODEL_NVARS_*p1+v])*dx_inverse );
    }
    vel = P; /* useless statement to avoid compiler warning */
    _ArrayIncrementIndex_(ndims,dim,index,done);
  }

  return(0);
}

Euler1DUpwindRoe()

int Euler1DUpwindRoe	(	double *	fI,
		double *	fL,
		double *	fR,
		double *	uL,
		double *	uR,
		double *	u,
		int	dir,
		void *	s,
		double	t
	)

Roe's upwinding scheme.

\begin{equation} {\bf f}_{j+1/2} = \frac{1}{2}\left[ {\bf f}_{j+1/2}^L + {\bf f}_{j+1/2}^R - \left| A\left({\bf u}_{j+1/2}^L,{\bf u}_{j+1/2}^R\right) \right| \left( {\bf u}_{j+1/2}^R - {\bf u}_{j+1/2}^L \right)\right] \end{equation}

• Roe, P. L., “Approximate Riemann solvers, parameter vectors, and difference schemes,” Journal of Computational Physics, Vol. 43, No. 2, 1981, pp. 357–372, http://dx.doi.org/10.1016/0021-9991(81)90128-5.

This upwinding scheme is modified for the balanced discretization of the 1D Euler equations when there is a non-zero gravitational force. See the reference below. For flows without any gravitational forces, it reduces to its original form.

• Xing, Shu, "High Order Well-Balanced WENO Scheme for the Gas Dynamics Equations Under Gravitational Fields", J. Sci. Comput., 54, 2013, pp. 645–662, http://dx.doi.org/10.1007/s10915-012-9585-8.

Parameters

fI	Computed upwind interface flux
fL	Left-biased reconstructed interface flux
fR	Right-biased reconstructed interface flux
uL	Left-biased reconstructed interface solution
uR	Right-biased reconstructed interface solution
u	Cell-centered solution
dir	Spatial dimension (unused since this is a 1D system)
s	Solver object of type HyPar
t	Current solution time

Definition at line 30 of file Euler1DUpwind.c.

{
  HyPar     *solver = (HyPar*)    s;
  Euler1D   *param  = (Euler1D*)  solver->physics;
  int       done,k;
  _DECLARE_IERR_;

  int ndims = solver->ndims;
  int ghosts= solver->ghosts;
  int *dim  = solver->dim_local;

  int index_outer[ndims], index_inter[ndims], bounds_outer[ndims], bounds_inter[ndims];
  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] =  1;
  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_], L[_MODEL_NVARS_*_MODEL_NVARS_], 
                DL[_MODEL_NVARS_*_MODEL_NVARS_], modA[_MODEL_NVARS_*_MODEL_NVARS_];

  done = 0; _ArraySetValue_(index_outer,ndims,0);
  while (!done) {
    _ArrayCopy1D_(index_outer,index_inter,ndims);
    for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
      int p; _ArrayIndex1D_(ndims,bounds_inter,index_inter,0,p);
      int indexL[ndims]; _ArrayCopy1D_(index_inter,indexL,ndims); indexL[dir]--;
      int indexR[ndims]; _ArrayCopy1D_(index_inter,indexR,ndims);
      int pL; _ArrayIndex1D_(ndims,dim,indexL,ghosts,pL);
      int pR; _ArrayIndex1D_(ndims,dim,indexR,ghosts,pR);
      double udiff[_MODEL_NVARS_], uavg[_MODEL_NVARS_],udiss[_MODEL_NVARS_];

      /* Roe's upwinding scheme */

      udiff[0] = 0.5 * (uR[_MODEL_NVARS_*p+0] - uL[_MODEL_NVARS_*p+0]);
      udiff[1] = 0.5 * (uR[_MODEL_NVARS_*p+1] - uL[_MODEL_NVARS_*p+1]);
      udiff[2] = 0.5 * (uR[_MODEL_NVARS_*p+2] - uL[_MODEL_NVARS_*p+2]);

      _Euler1DRoeAverage_         (uavg,(u+_MODEL_NVARS_*pL),(u+_MODEL_NVARS_*pR),param); 
      _Euler1DEigenvalues_        (uavg,D,param,0);
      _Euler1DLeftEigenvectors_   (uavg,L,param,0);
      _Euler1DRightEigenvectors_  (uavg,R,param,0);

      double kappa = max(param->grav_field[pL],param->grav_field[pR]);
      k = 0; D[k] = kappa*absolute(D[k]);
      k = 4; D[k] = kappa*absolute(D[k]);
      k = 8; D[k] = kappa*absolute(D[k]);

      MatMult3(3,DL,D,L);
      MatMult3(3,modA,R,DL);

      udiss[0] = modA[0*_MODEL_NVARS_+0]*udiff[0] + modA[0*_MODEL_NVARS_+1]*udiff[1] + modA[0*_MODEL_NVARS_+2]*udiff[2];
      udiss[1] = modA[1*_MODEL_NVARS_+0]*udiff[0] + modA[1*_MODEL_NVARS_+1]*udiff[1] + modA[1*_MODEL_NVARS_+2]*udiff[2];
      udiss[2] = modA[2*_MODEL_NVARS_+0]*udiff[0] + modA[2*_MODEL_NVARS_+1]*udiff[1] + modA[2*_MODEL_NVARS_+2]*udiff[2];

      fI[_MODEL_NVARS_*p+0] = 0.5 * (fL[_MODEL_NVARS_*p+0]+fR[_MODEL_NVARS_*p+0]) - udiss[0];
      fI[_MODEL_NVARS_*p+1] = 0.5 * (fL[_MODEL_NVARS_*p+1]+fR[_MODEL_NVARS_*p+1]) - udiss[1];
      fI[_MODEL_NVARS_*p+2] = 0.5 * (fL[_MODEL_NVARS_*p+2]+fR[_MODEL_NVARS_*p+2]) - udiss[2];
    }
    _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
  }

  return(0);
}

Euler1DUpwinddFRoe()

int Euler1DUpwinddFRoe	(	double *	fI,
		double *	fL,
		double *	fR,
		double *	uL,
		double *	uR,
		double *	u,
		int	dir,
		void *	s,
		double	t
	)

The Roe upwinding scheme (Euler1DUpwindRoe) for the partitioned hyperbolic flux that comprises of the acoustic waves only (see Euler1DStiffFlux, _Euler1DSetStiffFlux_). Thus, only the characteristic fields / eigen-modes corresponding to \( u\pm a\) are used.

Parameters

fI	Computed upwind interface flux
fL	Left-biased reconstructed interface flux
fR	Right-biased reconstructed interface flux
uL	Left-biased reconstructed interface solution
uR	Right-biased reconstructed interface solution
u	Cell-centered solution
dir	Spatial dimension (unused since this is a 1D system)
s	Solver object of type HyPar
t	Current solution time

Definition at line 453 of file Euler1DUpwind.c.

{
  HyPar     *solver = (HyPar*)    s;
  Euler1D   *param  = (Euler1D*)  solver->physics;
  int       done,k;
  _DECLARE_IERR_;

  int     ndims = solver->ndims;
  int     ghosts= solver->ghosts;
  int     *dim  = solver->dim_local;
  double  *uref = param->solution;

  int index_outer[ndims], index_inter[ndims], bounds_outer[ndims], bounds_inter[ndims];
  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] =  1;
  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_], L[_MODEL_NVARS_*_MODEL_NVARS_], 
                DL[_MODEL_NVARS_*_MODEL_NVARS_], modA[_MODEL_NVARS_*_MODEL_NVARS_];

  done = 0; _ArraySetValue_(index_outer,ndims,0);
  while (!done) {
    _ArrayCopy1D_(index_outer,index_inter,ndims);
    for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
      int p; _ArrayIndex1D_(ndims,bounds_inter,index_inter,0,p);
      int indexL[ndims]; _ArrayCopy1D_(index_inter,indexL,ndims); indexL[dir]--;
      int indexR[ndims]; _ArrayCopy1D_(index_inter,indexR,ndims);
      int pL; _ArrayIndex1D_(ndims,dim,indexL,ghosts,pL);
      int pR; _ArrayIndex1D_(ndims,dim,indexR,ghosts,pR);
      double udiff[_MODEL_NVARS_], uavg[_MODEL_NVARS_],udiss[_MODEL_NVARS_];

      /* Roe's upwinding scheme */

      udiff[0] = 0.5 * (uR[_MODEL_NVARS_*p+0] - uL[_MODEL_NVARS_*p+0]);
      udiff[1] = 0.5 * (uR[_MODEL_NVARS_*p+1] - uL[_MODEL_NVARS_*p+1]);
      udiff[2] = 0.5 * (uR[_MODEL_NVARS_*p+2] - uL[_MODEL_NVARS_*p+2]);

      _Euler1DRoeAverage_         (uavg,(uref+_MODEL_NVARS_*pL),(uref+_MODEL_NVARS_*pR),param); 
      _Euler1DEigenvalues_        (uavg,D,param,0);
      _Euler1DLeftEigenvectors_   (uavg,L,param,0);
      _Euler1DRightEigenvectors_  (uavg,R,param,0);

      double kappa = max(param->grav_field[pL],param->grav_field[pR]);
      k = 0; D[k] = 0.0; 
      k = 4; D[k] = kappa*absolute(D[k]);
      k = 8; D[k] = kappa*absolute(D[k]);

      MatMult3(3,DL,D,L);
      MatMult3(3,modA,R,DL);

      udiss[0] = modA[0*_MODEL_NVARS_+0]*udiff[0] + modA[0*_MODEL_NVARS_+1]*udiff[1] + modA[0*_MODEL_NVARS_+2]*udiff[2];
      udiss[1] = modA[1*_MODEL_NVARS_+0]*udiff[0] + modA[1*_MODEL_NVARS_+1]*udiff[1] + modA[1*_MODEL_NVARS_+2]*udiff[2];
      udiss[2] = modA[2*_MODEL_NVARS_+0]*udiff[0] + modA[2*_MODEL_NVARS_+1]*udiff[1] + modA[2*_MODEL_NVARS_+2]*udiff[2];

      fI[_MODEL_NVARS_*p+0] = 0.5 * (fL[_MODEL_NVARS_*p+0]+fR[_MODEL_NVARS_*p+0]) - udiss[0];
      fI[_MODEL_NVARS_*p+1] = 0.5 * (fL[_MODEL_NVARS_*p+1]+fR[_MODEL_NVARS_*p+1]) - udiss[1];
      fI[_MODEL_NVARS_*p+2] = 0.5 * (fL[_MODEL_NVARS_*p+2]+fR[_MODEL_NVARS_*p+2]) - udiss[2];
    }
    _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
  }

  return(0);
}

Euler1DUpwindRF()

int Euler1DUpwindRF	(	double *	fI,
		double *	fL,
		double *	fR,
		double *	uL,
		double *	uR,
		double *	u,
		int	dir,
		void *	s,
		double	t
	)

Characteristic-based Roe-fixed upwinding scheme.

\begin{align} \alpha_{j+1/2}^{k,L} &= \sum_{k=1}^3 {\bf l}_{j+1/2}^k \cdot {\bf f}_{j+1/2}^{k,L}, \\ \alpha_{j+1/2}^{k,R} &= \sum_{k=1}^3 {\bf l}_{j+1/2}^k \cdot {\bf f}_{j+1/2}^{k,R}, \\ v_{j+1/2}^{k,L} &= \sum_{k=1}^3 {\bf l}_{j+1/2}^k \cdot {\bf u}_{j+1/2}^{k,L}, \\ v_{j+1/2}^{k,R} &= \sum_{k=1}^3 {\bf l}_{j+1/2}^k \cdot {\bf u}_{j+1/2}^{k,R}, \\ \alpha_{j+1/2}^k &= \left\{ \begin{array}{cc} \alpha_{j+1/2}^{k,L} & {\rm if}\ \lambda_{j,j+1/2,j+1} > 0 \\ \alpha_{j+1/2}^{k,R} & {\rm if}\ \lambda_{j,j+1/2,j+1} < 0 \\ \frac{1}{2}\left[ \alpha_{j+1/2}^{k,L} + \alpha_{j+1/2}^{k,R} - \left(\max_{\left[j,j+1\right]} \lambda\right) \left( v_{j+1/2}^{k,R} - v_{j+1/2}^{k,L} \right) \right] & {\rm otherwise} \end{array}\right., \\ {\bf f}_{j+1/2} &= \sum_{k=1}^3 \alpha_{j+1/2}^k {\bf r}_{j+1/2}^k \end{align}

where \({\bf l}\), \({\bf r}\), and \(\lambda\) are the left-eigenvectors, right-eigenvectors and eigenvalues. The subscripts denote the grid locations.

• C.-W. Shu, and S. Osher, "Efficient implementation of essentially non-oscillatory schemes, II", J. Comput. Phys., 83 (1989), pp. 32–78, http://dx.doi.org/10.1016/0021-9991(89)90222-2.

Note that this upwinding scheme cannot be used for solving flows with non-zero gravitational forces.

Parameters

fI	Computed upwind interface flux
fL	Left-biased reconstructed interface flux
fR	Right-biased reconstructed interface flux
uL	Left-biased reconstructed interface solution
uR	Right-biased reconstructed interface solution
u	Cell-centered solution
dir	Spatial dimension (unused since this is a 1D system)
s	Solver object of type HyPar
t	Current solution time

Definition at line 115 of file Euler1DUpwind.c.

{
  HyPar     *solver = (HyPar*)    s;
  Euler1D   *param  = (Euler1D*)  solver->physics;
  int       done,k;
  _DECLARE_IERR_;

  int ndims   = solver->ndims;
  int *dim    = solver->dim_local;
  int ghosts  = solver->ghosts;

  int index_outer[ndims], index_inter[ndims], bounds_outer[ndims], bounds_inter[ndims];
  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] =  1;
  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_], L[_MODEL_NVARS_*_MODEL_NVARS_];

  done = 0; _ArraySetValue_(index_outer,ndims,0);
  while (!done) {
    _ArrayCopy1D_(index_outer,index_inter,ndims);
    for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
      int p; _ArrayIndex1D_(ndims,bounds_inter,index_inter,0,p);
      int indexL[ndims]; _ArrayCopy1D_(index_inter,indexL,ndims); indexL[dir]--;
      int indexR[ndims]; _ArrayCopy1D_(index_inter,indexR,ndims);
      int pL; _ArrayIndex1D_(ndims,dim,indexL,ghosts,pL);
      int pR; _ArrayIndex1D_(ndims,dim,indexR,ghosts,pR);
      double uavg[_MODEL_NVARS_], fcL[_MODEL_NVARS_], fcR[_MODEL_NVARS_], 
             ucL[_MODEL_NVARS_], ucR[_MODEL_NVARS_], fc[_MODEL_NVARS_];
      double kappa = max(param->grav_field[pL],param->grav_field[pR]);

      /* Roe-Fixed upwinding scheme */

      _Euler1DRoeAverage_       (uavg,(u+_MODEL_NVARS_*pL),(u+_MODEL_NVARS_*pR),param);
      _Euler1DEigenvalues_      (uavg,D,param,0);
      _Euler1DLeftEigenvectors_ (uavg,L,param,0);
      _Euler1DRightEigenvectors_(uavg,R,param,0);

      /* calculate characteristic fluxes and variables */
      MatVecMult3(_MODEL_NVARS_,ucL,L,(uL+_MODEL_NVARS_*p));
      MatVecMult3(_MODEL_NVARS_,ucR,L,(uR+_MODEL_NVARS_*p));
      MatVecMult3(_MODEL_NVARS_,fcL,L,(fL+_MODEL_NVARS_*p));
      MatVecMult3(_MODEL_NVARS_,fcR,L,(fR+_MODEL_NVARS_*p));

      for (k = 0; k < _MODEL_NVARS_; k++) {
        double eigL,eigC,eigR;
        _Euler1DEigenvalues_((u+_MODEL_NVARS_*pL),D,param,0);
        eigL = D[k*_MODEL_NVARS_+k];
        _Euler1DEigenvalues_((u+_MODEL_NVARS_*pR),D,param,0);
        eigR = D[k*_MODEL_NVARS_+k];
        _Euler1DEigenvalues_(uavg,D,param,0);
        eigC = D[k*_MODEL_NVARS_+k];

        if ((eigL > 0) && (eigC > 0) && (eigR > 0))       fc[k] = fcL[k];
        else if ((eigL < 0) && (eigC < 0) && (eigR < 0))  fc[k] = fcR[k];
        else {
          double alpha = kappa * max3(absolute(eigL),absolute(eigC),absolute(eigR));
          fc[k] = 0.5 * (fcL[k] + fcR[k] + alpha * (ucL[k]-ucR[k]));
        }

      }

      /* calculate the interface flux from the characteristic flux */
      MatVecMult3(_MODEL_NVARS_,(fI+_MODEL_NVARS_*p),R,fc);
    }
    _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
  }

  return(0);
}

Euler1DUpwinddFRF()

int Euler1DUpwinddFRF	(	double *	fI,
		double *	fL,
		double *	fR,
		double *	uL,
		double *	uR,
		double *	u,
		int	dir,
		void *	s,
		double	t
	)

The characteristic-based Roe-fixed upwinding scheme (Euler1DUpwindRF) for the partitioned hyperbolic flux that comprises of the acoustic waves only (see Euler1DStiffFlux, _Euler1DSetStiffFlux_). Thus, only the characteristic fields / eigen-modes corresponding to \( u\pm a\) are used.

Reference:

• Ghosh, D., Constantinescu, E. M., Semi-Implicit Time Integration of Atmospheric Flows with Characteristic-Based Flux Partitioning, SIAM Journal on Scientific Computing, 38 (3), 2016, A1848-A1875, http://dx.doi.org/10.1137/15M1044369

Parameters

fI	Computed upwind interface flux
fL	Left-biased reconstructed interface flux
fR	Right-biased reconstructed interface flux
uL	Left-biased reconstructed interface solution
uR	Right-biased reconstructed interface solution
u	Cell-centered solution
dir	Spatial dimension (unused since this is a 1D system)
s	Solver object of type HyPar
t	Current solution time

Definition at line 534 of file Euler1DUpwind.c.

{
  HyPar     *solver = (HyPar*)    s;
  Euler1D   *param  = (Euler1D*)  solver->physics;
  int       done,k;
  _DECLARE_IERR_;

  int     ndims   = solver->ndims;
  int     *dim    = solver->dim_local;
  int     ghosts  = solver->ghosts;
  double  *uref   = param->solution;

  int index_outer[ndims], index_inter[ndims], bounds_outer[ndims], bounds_inter[ndims];
  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] =  1;
  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_], L[_MODEL_NVARS_*_MODEL_NVARS_];

  done = 0; _ArraySetValue_(index_outer,ndims,0);
  while (!done) {
    _ArrayCopy1D_(index_outer,index_inter,ndims);
    for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
      int p; _ArrayIndex1D_(ndims,bounds_inter,index_inter,0,p);
      int indexL[ndims]; _ArrayCopy1D_(index_inter,indexL,ndims); indexL[dir]--;
      int indexR[ndims]; _ArrayCopy1D_(index_inter,indexR,ndims);
      int pL; _ArrayIndex1D_(ndims,dim,indexL,ghosts,pL);
      int pR; _ArrayIndex1D_(ndims,dim,indexR,ghosts,pR);
      double uavg[_MODEL_NVARS_], fcL[_MODEL_NVARS_], fcR[_MODEL_NVARS_], 
             ucL[_MODEL_NVARS_], ucR[_MODEL_NVARS_], fc[_MODEL_NVARS_];
      double kappa = max(param->grav_field[pL],param->grav_field[pR]);

      /* Roe-Fixed upwinding scheme */

      _Euler1DRoeAverage_       (uavg,(uref+_MODEL_NVARS_*pL),(uref+_MODEL_NVARS_*pR),param);
      _Euler1DEigenvalues_      (uavg,D,param,0);
      _Euler1DLeftEigenvectors_ (uavg,L,param,0);
      _Euler1DRightEigenvectors_(uavg,R,param,0);

      /* calculate characteristic fluxes and variables */
      MatVecMult3(_MODEL_NVARS_,ucL,L,(uL+_MODEL_NVARS_*p));
      MatVecMult3(_MODEL_NVARS_,ucR,L,(uR+_MODEL_NVARS_*p));
      MatVecMult3(_MODEL_NVARS_,fcL,L,(fL+_MODEL_NVARS_*p));
      MatVecMult3(_MODEL_NVARS_,fcR,L,(fR+_MODEL_NVARS_*p));

      for (k = 0; k < _MODEL_NVARS_; k++) {
        double eigL,eigC,eigR;
        _Euler1DEigenvalues_((u+_MODEL_NVARS_*pL),D,param,0);
        eigL = (k == 0? 0.0 : D[k*_MODEL_NVARS_+k]);
        _Euler1DEigenvalues_((u+_MODEL_NVARS_*pR),D,param,0);
        eigR = (k == 0? 0.0 : D[k*_MODEL_NVARS_+k]);
        _Euler1DEigenvalues_(uavg,D,param,0);
        eigC = (k == 0? 0.0 : D[k*_MODEL_NVARS_+k]);

        if ((eigL > 0) && (eigC > 0) && (eigR > 0))       fc[k] = fcL[k];
        else if ((eigL < 0) && (eigC < 0) && (eigR < 0))  fc[k] = fcR[k];
        else {
          double alpha = kappa * max3(absolute(eigL),absolute(eigC),absolute(eigR));
          fc[k] = 0.5 * (fcL[k] + fcR[k] + alpha * (ucL[k]-ucR[k]));
        }

      }

      /* calculate the interface flux from the characteristic flux */
      MatVecMult3(_MODEL_NVARS_,(fI+_MODEL_NVARS_*p),R,fc);
    }
    _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
  }

  return(0);
}

Euler1DUpwindLLF()

int Euler1DUpwindLLF	(	double *	fI,
		double *	fL,
		double *	fR,
		double *	uL,
		double *	uR,
		double *	u,
		int	dir,
		void *	s,
		double	t
	)

Characteristic-based local Lax-Friedrich upwinding scheme.

\begin{align} \alpha_{j+1/2}^{k,L} &= \sum_{k=1}^3 {\bf l}_{j+1/2}^k \cdot {\bf f}_{j+1/2}^{k,L}, \\ \alpha_{j+1/2}^{k,R} &= \sum_{k=1}^3 {\bf l}_{j+1/2}^k \cdot {\bf f}_{j+1/2}^{k,R}, \\ v_{j+1/2}^{k,L} &= \sum_{k=1}^3 {\bf l}_{j+1/2}^k \cdot {\bf u}_{j+1/2}^{k,L}, \\ v_{j+1/2}^{k,R} &= \sum_{k=1}^3 {\bf l}_{j+1/2}^k \cdot {\bf u}_{j+1/2}^{k,R}, \\ \alpha_{j+1/2}^k &= \frac{1}{2}\left[ \alpha_{j+1/2}^{k,L} + \alpha_{j+1/2}^{k,R} - \left(\max_{\left[j,j+1\right]} \lambda\right) \left( v_{j+1/2}^{k,R} - v_{j+1/2}^{k,L} \right) \right], \\ {\bf f}_{j+1/2} &= \sum_{k=1}^3 \alpha_{j+1/2}^k {\bf r}_{j+1/2}^k \end{align}

where \({\bf l}\), \({\bf r}\), and \(\lambda\) are the left-eigenvectors, right-eigenvectors and eigenvalues. The subscripts denote the grid locations.

• C.-W. Shu, and S. Osher, "Efficient implementation of essentially non-oscillatory schemes, II", J. Comput. Phys., 83 (1989), pp. 32–78, http://dx.doi.org/10.1016/0021-9991(89)90222-2.

This upwinding scheme is modified for the balanced discretization of the 1D Euler equations when there is a non-zero gravitational force. See the reference below. For flows without any gravitational forces, it reduces to its original form.

• Xing, Shu, "High Order Well-Balanced WENO Scheme for the Gas Dynamics Equations Under Gravitational Fields", J. Sci. Comput., 54, 2013, pp. 645–662, http://dx.doi.org/10.1007/s10915-012-9585-8.

Parameters

fI	Computed upwind interface flux
fL	Left-biased reconstructed interface flux
fR	Right-biased reconstructed interface flux
uL	Left-biased reconstructed interface solution
uR	Right-biased reconstructed interface solution
u	Cell-centered solution
dir	Spatial dimension (unused since this is a 1D system)
s	Solver object of type HyPar
t	Current solution time

Definition at line 212 of file Euler1DUpwind.c.

{
  HyPar     *solver = (HyPar*)    s;
  Euler1D   *param  = (Euler1D*)  solver->physics;
  int       done,k;
  _DECLARE_IERR_;

  int ndims   = solver->ndims;
  int *dim    = solver->dim_local;
  int ghosts  = solver->ghosts;

  int index_outer[ndims], index_inter[ndims], bounds_outer[ndims], bounds_inter[ndims];
  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] =  1;
  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_], L[_MODEL_NVARS_*_MODEL_NVARS_];

  done = 0; _ArraySetValue_(index_outer,ndims,0);
  while (!done) {
    _ArrayCopy1D_(index_outer,index_inter,ndims);
    for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
      int p; _ArrayIndex1D_(ndims,bounds_inter,index_inter,0,p);
      int indexL[ndims]; _ArrayCopy1D_(index_inter,indexL,ndims); indexL[dir]--;
      int indexR[ndims]; _ArrayCopy1D_(index_inter,indexR,ndims);
      int pL; _ArrayIndex1D_(ndims,dim,indexL,ghosts,pL);
      int pR; _ArrayIndex1D_(ndims,dim,indexR,ghosts,pR);
      double uavg[_MODEL_NVARS_], fcL[_MODEL_NVARS_], fcR[_MODEL_NVARS_], 
             ucL[_MODEL_NVARS_], ucR[_MODEL_NVARS_], fc[_MODEL_NVARS_];
      double kappa = max(param->grav_field[pL],param->grav_field[pR]);

      /* Local Lax-Friedrich upwinding scheme */

      _Euler1DRoeAverage_       (uavg,(u+_MODEL_NVARS_*pL),(u+_MODEL_NVARS_*pR),param);
      _Euler1DEigenvalues_      (uavg,D,param,0);
      _Euler1DLeftEigenvectors_ (uavg,L,param,0);
      _Euler1DRightEigenvectors_(uavg,R,param,0);

      /* calculate characteristic fluxes and variables */
      MatVecMult3(_MODEL_NVARS_,ucL,L,(uL+_MODEL_NVARS_*p));
      MatVecMult3(_MODEL_NVARS_,ucR,L,(uR+_MODEL_NVARS_*p));
      MatVecMult3(_MODEL_NVARS_,fcL,L,(fL+_MODEL_NVARS_*p));
      MatVecMult3(_MODEL_NVARS_,fcR,L,(fR+_MODEL_NVARS_*p));

      for (k = 0; k < _MODEL_NVARS_; k++) {
        double eigL,eigC,eigR;
        _Euler1DEigenvalues_((u+_MODEL_NVARS_*pL),D,param,0);
        eigL = D[k*_MODEL_NVARS_+k];
        _Euler1DEigenvalues_((u+_MODEL_NVARS_*pR),D,param,0);
        eigR = D[k*_MODEL_NVARS_+k];
        _Euler1DEigenvalues_(uavg,D,param,0);
        eigC = D[k*_MODEL_NVARS_+k];

        double alpha = kappa * max3(absolute(eigL),absolute(eigC),absolute(eigR));
        fc[k] = 0.5 * (fcL[k] + fcR[k] + alpha * (ucL[k]-ucR[k]));
      }

      /* calculate the interface flux from the characteristic flux */
      MatVecMult3(_MODEL_NVARS_,(fI+_MODEL_NVARS_*p),R,fc);
    }
    _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
  }

  return(0);
}

Euler1DUpwinddFLLF()

int Euler1DUpwinddFLLF	(	double *	fI,
		double *	fL,
		double *	fR,
		double *	uL,
		double *	uR,
		double *	u,
		int	dir,
		void *	s,
		double	t
	)

The characteristic-based local Lax-Friedrich upwinding scheme (Euler1DUpwindLLF) for the partitioned hyperbolic flux that comprises of the acoustic waves only (see Euler1DStiffFlux, _Euler1DSetStiffFlux_). Thus, only the characteristic fields / eigen-modes corresponding to \( u\pm a\) are used.

Reference:

• Ghosh, D., Constantinescu, E. M., Semi-Implicit Time Integration of Atmospheric Flows with Characteristic-Based Flux Partitioning, SIAM Journal on Scientific Computing, 38 (3), 2016, A1848-A1875, http://dx.doi.org/10.1137/15M1044369

Parameters

fI	Computed upwind interface flux
fL	Left-biased reconstructed interface flux
fR	Right-biased reconstructed interface flux
uL	Left-biased reconstructed interface solution
uR	Right-biased reconstructed interface solution
u	Cell-centered solution
dir	Spatial dimension (unused since this is a 1D system)
s	Solver object of type HyPar
t	Current solution time

Definition at line 623 of file Euler1DUpwind.c.

{
  HyPar     *solver = (HyPar*)    s;
  Euler1D   *param  = (Euler1D*)  solver->physics;
  int       done,k;
  _DECLARE_IERR_;

  int     ndims   = solver->ndims;
  int     *dim    = solver->dim_local;
  int     ghosts  = solver->ghosts;
  double  *uref   = param->solution;

  int index_outer[ndims], index_inter[ndims], bounds_outer[ndims], bounds_inter[ndims];
  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] =  1;
  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_], L[_MODEL_NVARS_*_MODEL_NVARS_];

  done = 0; _ArraySetValue_(index_outer,ndims,0);
  while (!done) {
    _ArrayCopy1D_(index_outer,index_inter,ndims);
    for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
      int p; _ArrayIndex1D_(ndims,bounds_inter,index_inter,0,p);
      int indexL[ndims]; _ArrayCopy1D_(index_inter,indexL,ndims); indexL[dir]--;
      int indexR[ndims]; _ArrayCopy1D_(index_inter,indexR,ndims);
      int pL; _ArrayIndex1D_(ndims,dim,indexL,ghosts,pL);
      int pR; _ArrayIndex1D_(ndims,dim,indexR,ghosts,pR);
      double uavg[_MODEL_NVARS_], fcL[_MODEL_NVARS_], fcR[_MODEL_NVARS_], 
             ucL[_MODEL_NVARS_], ucR[_MODEL_NVARS_], fc[_MODEL_NVARS_];
      double kappa = max(param->grav_field[pL],param->grav_field[pR]);

      /* Local Lax-Friedrich upwinding scheme */

      _Euler1DRoeAverage_       (uavg,(uref+_MODEL_NVARS_*pL),(uref+_MODEL_NVARS_*pR),param);
      _Euler1DEigenvalues_      (uavg,D,param,0);
      _Euler1DLeftEigenvectors_ (uavg,L,param,0);
      _Euler1DRightEigenvectors_(uavg,R,param,0);

      /* calculate characteristic fluxes and variables */
      MatVecMult3(_MODEL_NVARS_,ucL,L,(uL+_MODEL_NVARS_*p));
      MatVecMult3(_MODEL_NVARS_,ucR,L,(uR+_MODEL_NVARS_*p));
      MatVecMult3(_MODEL_NVARS_,fcL,L,(fL+_MODEL_NVARS_*p));
      MatVecMult3(_MODEL_NVARS_,fcR,L,(fR+_MODEL_NVARS_*p));

      for (k = 0; k < _MODEL_NVARS_; k++) {
        double eigL,eigC,eigR;
        _Euler1DEigenvalues_((u+_MODEL_NVARS_*pL),D,param,0);
        eigL = (k == 0? 0.0 : D[k*_MODEL_NVARS_+k]);
        _Euler1DEigenvalues_((u+_MODEL_NVARS_*pR),D,param,0);
        eigR = (k == 0? 0.0 : D[k*_MODEL_NVARS_+k]);
        _Euler1DEigenvalues_(uavg,D,param,0);
        eigC = (k == 0? 0.0 : D[k*_MODEL_NVARS_+k]);

        double alpha = kappa * max3(absolute(eigL),absolute(eigC),absolute(eigR));
        fc[k] = 0.5 * (fcL[k] + fcR[k] + alpha * (ucL[k]-ucR[k]));
      }

      /* calculate the interface flux from the characteristic flux */
      MatVecMult3(_MODEL_NVARS_,(fI+_MODEL_NVARS_*p),R,fc);
    }
    _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
  }

  return(0);
}

Euler1DUpwindSWFS()

int Euler1DUpwindSWFS	(	double *	fI,
		double *	fL,
		double *	fR,
		double *	uL,
		double *	uR,
		double *	u,
		int	dir,
		void *	s,
		double	t
	)

Steger-Warming Flux-Splitting scheme

• Steger, J.L., Warming, R.F., "Flux vector splitting of the inviscid gasdynamic equations with application to finite-difference methods", J. Comput. Phys., 40(2), 1981, pp. 263-293, http://dx.doi.org/10.1016/0021-9991(81)90210-2.

Note that this method cannot be used for flows with non-zero gravitational forces.

Parameters

fI	Computed upwind interface flux
fL	Left-biased reconstructed interface flux
fR	Right-biased reconstructed interface flux
uL	Left-biased reconstructed interface solution
uR	Right-biased reconstructed interface solution
u	Cell-centered solution
dir	Spatial dimension (unused since this is a 1D system)
s	Solver object of type HyPar
t	Current solution time

Definition at line 293 of file Euler1DUpwind.c.

{
  HyPar     *solver = (HyPar*)    s;
  Euler1D   *param  = (Euler1D*)  solver->physics;
  int       done;
  _DECLARE_IERR_;

  int ndims = solver->ndims;
  int *dim  = solver->dim_local;

  int index_outer[ndims], index_inter[ndims], bounds_outer[ndims], bounds_inter[ndims];
  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] =  1;
  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
  static double fp[_MODEL_NVARS_], fm[_MODEL_NVARS_],uavg[_MODEL_NVARS_];

  done = 0; _ArraySetValue_(index_outer,ndims,0);
  while (!done) {
    _ArrayCopy1D_(index_outer,index_inter,ndims);
    for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
      int p; _ArrayIndex1D_(ndims,bounds_inter,index_inter,0,p);
      double rho,v,e,P,c,gamma=param->gamma,term,Mach;

      /* Steger Warming flux splitting */
      _Euler1DRoeAverage_(uavg,(uL+_MODEL_NVARS_*p),(uR+_MODEL_NVARS_*p),param); 
      _Euler1DGetFlowVar_(uavg,rho,v,e,P,param);
      Mach = v/sqrt(gamma*P/rho);

      if (Mach < -1.0) {

        _ArrayCopy1D3_((fR+_MODEL_NVARS_*p),(fI+_MODEL_NVARS_*p),_MODEL_NVARS_);

      } else if (Mach < 1.0) {

        _Euler1DGetFlowVar_((uL+_MODEL_NVARS_*p),rho,v,e,P,param);
        c = sqrt(gamma*P/rho);
        term = rho/(2.0*gamma);

        fp[0] = term * (2*gamma*v + c - v);
        fp[1] = term * (2*(gamma-1.0)*v*v + (v+c)*(v+c));
        fp[2] = term * ((gamma-1.0)*v*v*v + 0.5*(v+c)*(v+c)*(v+c) + ((3.0-gamma)*(v+c)*c*c)/(2.0*(gamma-1.0)));

        _Euler1DGetFlowVar_((uR+_MODEL_NVARS_*p),rho,v,e,P,param);
        c = sqrt(gamma*P/rho);
        term = rho/(2.0*gamma);

        fm[0] = term * (v - c);
        fm[1] = term * (v-c) * (v-c);
        fm[2] = term * (0.5*(v-c)*(v-c)*(v-c) + ((3.0-gamma)*(v-c)*c*c)/(2.0*(gamma-1.0)));

        _ArrayAdd1D_((fI+_MODEL_NVARS_*p),fp,fm,_MODEL_NVARS_);

      } else {

        _ArrayCopy1D3_((fL+_MODEL_NVARS_*p),(fI+_MODEL_NVARS_*p),_MODEL_NVARS_);

      }

    }
    _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
  }

  return(0);
}

Euler1DUpwindRusanov()

int Euler1DUpwindRusanov	(	double *	fI,
		double *	fL,
		double *	fR,
		double *	uL,
		double *	uR,
		double *	u,
		int	dir,
		void *	s,
		double	t
	)

Rusanov's upwinding scheme.

\begin{equation} {\bf f}_{j+1/2} = \frac{1}{2}\left[ {\bf f}_{j+1/2}^L + {\bf f}_{j+1/2}^R - \max_{j,j+1} \nu_j \left( {\bf u}_{j+1/2}^R - {\bf u}_{j+1/2}^L \right)\right] \end{equation}

where \(\nu = c + \left|u\right|\).

• Rusanov, V. V., "The calculation of the interaction of non-stationary shock waves and obstacles," USSR Computational Mathematics and Mathematical Physics, Vol. 1, No. 2, 1962, pp. 304–320

Parameters

fI	Computed upwind interface flux
fL	Left-biased reconstructed interface flux
fR	Right-biased reconstructed interface flux
uL	Left-biased reconstructed interface solution
uR	Right-biased reconstructed interface solution
u	Cell-centered solution
dir	Spatial dimension (unused since this is a 1D system)
s	Solver object of type HyPar
t	Current solution time

Definition at line 376 of file Euler1DUpwind.c.

{
  HyPar     *solver = (HyPar*)    s;
  Euler1D   *param  = (Euler1D*)  solver->physics;
  int       done,k;
  _DECLARE_IERR_;

  int ndims = solver->ndims;
  int ghosts= solver->ghosts;
  int *dim  = solver->dim_local;

  int index_outer[ndims], index_inter[ndims], bounds_outer[ndims], bounds_inter[ndims];
  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] =  1;
  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;

  static double udiff[_MODEL_NVARS_], uavg[_MODEL_NVARS_];

  done = 0; _ArraySetValue_(index_outer,ndims,0);
  while (!done) {

    _ArrayCopy1D_(index_outer,index_inter,ndims);

    for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {

      int p; _ArrayIndex1D_(ndims,bounds_inter,index_inter,0,p);
      int indexL[ndims]; _ArrayCopy1D_(index_inter,indexL,ndims); indexL[dir]--;
      int indexR[ndims]; _ArrayCopy1D_(index_inter,indexR,ndims);
      int pL; _ArrayIndex1D_(ndims,dim,indexL,ghosts,pL);
      int pR; _ArrayIndex1D_(ndims,dim,indexR,ghosts,pR);

      _Euler1DRoeAverage_(uavg,(u+_MODEL_NVARS_*pL),(u+_MODEL_NVARS_*pR),param);
      for (k = 0; k < _MODEL_NVARS_; k++) udiff[k] = 0.5 * (uR[_MODEL_NVARS_*p+k] - uL[_MODEL_NVARS_*p+k]);

      double rho, uvel, E, P, c;

      _Euler1DGetFlowVar_((u+_MODEL_NVARS_*pL),rho,uvel,E,P,param);
      c = param->gamma*P/rho;
      double alphaL = c + absolute(uvel);

      _Euler1DGetFlowVar_((u+_MODEL_NVARS_*pR),rho,uvel,E,P,param);
      c = param->gamma*P/rho;
      double alphaR = c + absolute(uvel);

      _Euler1DGetFlowVar_(uavg,rho,uvel,E,P,param);
      c = param->gamma*P/rho;
      double alphaavg = c + absolute(uvel);

      double kappa = max(param->grav_field[pL],param->grav_field[pR]);
      double alpha = kappa*max3(alphaL,alphaR,alphaavg);

      for (k = 0; k < _MODEL_NVARS_; k++) {
        fI[_MODEL_NVARS_*p+k] = 0.5 * (fL[_MODEL_NVARS_*p+k]+fR[_MODEL_NVARS_*p+k]) - alpha*udiff[k];
      }

    }

    _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);

  }

  return(0);
}

Euler1DJacobian()

int Euler1DJacobian	(	double *	Jac,
		double *	u,
		void *	p,
		int	dir,
		int	nvars,
		int	upw
	)

Function to compute the flux Jacobian of the 1D Euler equations, given the solution at a grid point. The Jacobian is square matrix of size nvar=3, and is returned as a 1D array (double) of 9 elements in row-major format.

Parameters

Jac	Jacobian matrix: 1D array of size nvar^2 = 9
u	solution at a grid point (array of size nvar = 3)
p	object containing the physics-related parameters
dir	dimension (x/y/z) (not used, since this is 1D system)
nvars	number of vector components
upw	0 -> send back complete Jacobian, 1 -> send back Jacobian of right(+)-moving flux, -1 -> send back Jacobian of left(-)-moving flux

Definition at line 14 of file Euler1DJacobian.c.

{
  Euler1D       *param = (Euler1D*) p;
  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_], 
                L[_MODEL_NVARS_*_MODEL_NVARS_], DL[_MODEL_NVARS_*_MODEL_NVARS_];

  /* get the eigenvalues and left,right eigenvectors */
  _Euler1DEigenvalues_      (u,D,param,0);
  _Euler1DLeftEigenvectors_ (u,L,param,0);
  _Euler1DRightEigenvectors_(u,R,param,0);

  int aupw = absolute(upw), k;
  k = 0; D[k] = absolute( (1-aupw)*D[k] + 0.5*aupw*(1+upw)*max(0,D[k]) + 0.5*aupw*(1-upw)*min(0,D[k]) );
  k = 4; D[k] = absolute( (1-aupw)*D[k] + 0.5*aupw*(1+upw)*max(0,D[k]) + 0.5*aupw*(1-upw)*min(0,D[k]) );
  k = 8; D[k] = absolute( (1-aupw)*D[k] + 0.5*aupw*(1+upw)*max(0,D[k]) + 0.5*aupw*(1-upw)*min(0,D[k]) );

  MatMult3(3,DL,D,L);
  MatMult3(3,Jac,R,DL);

  return(0);
}

Euler1DStiffJacobian()

int Euler1DStiffJacobian	(	double *	Jac,
		double *	u,
		void *	p,
		int	dir,
		int	nvars,
		int	upw
	)

Function to compute the Jacobian of the fast flux (representing the acoustic waves) of the 1D Euler equations, given the solution at a grid point (see _Euler1DSetStiffFlux_, _Euler1DSetStiffJac_). The Jacobian is square matrix of size nvar=3, and is returned as a 1D array (double) of 9 elements in row-major format.

Parameters

Jac	Jacobian matrix: 1D array of size nvar^2 = 9
u	solution at a grid point (array of size nvar = 3)
p	object containing the physics-related parameters
dir	dimension (x/y/z) (not used, since this is 1D system)
nvars	number of vector components
upw	0 -> send back complete Jacobian, 1 -> send back Jacobian of right(+)-moving flux, -1 -> send back Jacobian of left(-)-moving flux

Definition at line 48 of file Euler1DJacobian.c.

{
  Euler1D       *param = (Euler1D*) p;
  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_], 
                L[_MODEL_NVARS_*_MODEL_NVARS_], DL[_MODEL_NVARS_*_MODEL_NVARS_];

  /* get the eigenvalues and left,right eigenvectors */
  _Euler1DEigenvalues_      (u,D,param,0);
  _Euler1DLeftEigenvectors_ (u,L,param,0);
  _Euler1DRightEigenvectors_(u,R,param,0);

  int aupw = absolute(upw), k;
  k = 0; D[k] = 0.0; /* remove the entropy eigenmode */
  k = 4; D[k] = absolute( (1-aupw)*D[k] + 0.5*aupw*(1+upw)*max(0,D[k]) + 0.5*aupw*(1-upw)*min(0,D[k]) );
  k = 8; D[k] = absolute( (1-aupw)*D[k] + 0.5*aupw*(1+upw)*max(0,D[k]) + 0.5*aupw*(1-upw)*min(0,D[k]) );

  MatMult3(3,DL,D,L);
  MatMult3(3,Jac,R,DL);

  return(0);
}

Euler1DRoeAverage()

int Euler1DRoeAverage	(	double *	uavg,
		double *	uL,
		double *	uR,
		void *	p
	)

Compute the Roe-averaged state for the 1D Euler equations. This function just calls the macro _Euler1DRoeAverage_ and is not used by any functions within the 1D Euler module. However, it's necessary to define it and provide it to the the solver object (HyPar) so that it can then send it to interpolation functions for a characteristic-based reconstruction.

Parameters

uavg	The computed Roe-averaged state
uL	Left state (conserved variables)
uR	Right state (conserved variables)
p	Object of type Euler1D with physics-related variables

Definition at line 16 of file Euler1DFunctions.c.

{
  Euler1D *param  = (Euler1D*) p;
  _Euler1DRoeAverage_(uavg,uL,uR,param);
  return(0);
}

Euler1DLeftEigenvectors()

int Euler1DLeftEigenvectors	(	double *	u,
		double *	L,
		void *	p,
		int	dir
	)

Compute the left eigenvections for the 1D Euler equations. This function just calls the macro _Euler1DLeftEigenvectors_ and is not used by any functions within the 1D Euler module. However, it's necessary to define it and provide it to the the solver object (HyPar) so that it can then send it to interpolation functions for a characteristic-based reconstruction.

Parameters

u	Conserved solution at a grid point
L	Array of size nvar^2 = 3^2 to save the matrix of left eigenvectors in (row-major format).
p	Object of type Euler1D with physics-related variables
dir	Spatial dimension (not used, since this is a 1D system)

Definition at line 19 of file Euler1DEigen.c.

{
  Euler1D *param  = (Euler1D*)  p;
  _Euler1DLeftEigenvectors_(u,L,param,dir);
  return(0);
}

Euler1DRightEigenvectors()

int Euler1DRightEigenvectors	(	double *	u,
		double *	R,
		void *	p,
		int	dir
	)

Compute the right eigenvections for the 1D Euler equations. This function just calls the macro _Euler1DRightEigenvectors_ and is not used by any functions within the 1D Euler module. However, it's necessary to define it and provide it to the the solver object (HyPar) so that it can then send it to interpolation functions for a characteristic-based reconstruction.

Parameters

u	Conserved solution at a grid point
R	Array of size nvar^2 = 3^2 to save the matrix of right eigenvectors in (row-major format).
p	Object of type Euler1D with physics-related variables
dir	Spatial dimension (not used, since this is a 1D system)

Definition at line 38 of file Euler1DEigen.c.

{
  Euler1D *param  = (Euler1D*)  p;
  _Euler1DRightEigenvectors_(u,R,param,dir);
  return(0);
}

Euler1DGravityField()

int Euler1DGravityField	(	void *	s,
		void *	m
	)

Compute the gravitational field over the domain. See

• Xing, Y., Shu, C.-W., "High Order Well-Balanced WENO Scheme for the Gas Dynamics Equations Under Gravitational Fields", Journal of Scientific Computing, 54, 2013, pp. 645-662. http://dx.doi.org/10.1007/s10915-012-9585-8
  for details on the treatment of the gravitational source term. This function computes the graviational potential function. If there is no gravity, the gravitational field is 1.0 all over the domain.

Parameters

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

Definition at line 23 of file Euler1DGravityField.c.

{
  HyPar         *solver = (HyPar*)        s;
  MPIVariables  *mpi    = (MPIVariables*) m;
  Euler1D       *param  = (Euler1D*)      solver->physics;

  double  *S      = param->grav_field;
  int     *dim    = solver->dim_local;
  int     ghosts  = solver->ghosts;
  int     index[_MODEL_NDIMS_], bounds[_MODEL_NDIMS_],
          offset[_MODEL_NDIMS_], d, done;

  /* set bounds for array index to include ghost points */
  _ArrayCopy1D_(dim,bounds,_MODEL_NDIMS_);
  for (d=0; d<_MODEL_NDIMS_; d++) bounds[d] += 2*ghosts;
  /* set offset such that index is compatible with ghost point arrangement */
  _ArraySetValue_(offset,_MODEL_NDIMS_,-ghosts);

  /* set the value of the gravity field */
  done = 0; _ArraySetValue_(index,_MODEL_NDIMS_,0);
  while (!done) {
    int p;                _ArrayIndex1DWO_(_MODEL_NDIMS_,dim,index,offset,ghosts,p);
    double xcoord;        _GetCoordinate_(_XDIR_,index[_XDIR_]-ghosts,dim,ghosts,solver->x,xcoord);
    if (param->grav_type == 0) {
      S[p] = exp(-param->grav*xcoord);
    } else if (param->grav_type == 1) {
      double pi = 4.0 * atan(1.0);
      double phi = -sin(2*pi*xcoord)/(2*pi);
      S[p] = exp(-phi);
    }
    _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds,index,done);
  }

  if (param->grav_type != 1) {
    /* a sensible simulation will not specify peridic boundary conditions 
    * along a direction in which gravity acts (unless the gravitational field
    * itself is sinusoidal), so extrapolate the gravity field at the boundaries */
    int indexb[_MODEL_NDIMS_], indexi[_MODEL_NDIMS_];
    for (d = 0; d < _MODEL_NDIMS_; d++) {
      /* left boundary */
      if (!mpi->ip[d]) {
        _ArrayCopy1D_(dim,bounds,_MODEL_NDIMS_); bounds[d] = ghosts;
        _ArraySetValue_(offset,_MODEL_NDIMS_,0); offset[d] = -ghosts;
        done = 0; _ArraySetValue_(indexb,_MODEL_NDIMS_,0);
        while (!done) {
          _ArrayCopy1D_(indexb,indexi,_MODEL_NDIMS_); indexi[d] = ghosts-1-indexb[d];
          int p1; _ArrayIndex1DWO_(_MODEL_NDIMS_,dim,indexb,offset,ghosts,p1);
          int p2; _ArrayIndex1D_  (_MODEL_NDIMS_,dim,indexi,ghosts,p2);
          S[p1] = S[p2];
          _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds,indexb,done);
        }
      }
      /* right boundary */
      if (mpi->ip[d] == mpi->iproc[d]-1) {
        _ArrayCopy1D_(dim,bounds,_MODEL_NDIMS_); bounds[d] = ghosts;
        _ArraySetValue_(offset,_MODEL_NDIMS_,0); offset[d] = dim[d];
        done = 0; _ArraySetValue_(indexb,_MODEL_NDIMS_,0);
        while (!done) {
          _ArrayCopy1D_(indexb,indexi,_MODEL_NDIMS_); indexi[d] = dim[d]-1-indexb[d];
          int p1; _ArrayIndex1DWO_(_MODEL_NDIMS_,dim,indexb,offset,ghosts,p1);
          int p2; _ArrayIndex1D_  (_MODEL_NDIMS_,dim,indexi,ghosts,p2);
          S[p1] = S[p2];
          _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds,indexb,done);
        }
      }
    }
  }

  return(0);
}

Euler1DSourceUpwindLLF()

int Euler1DSourceUpwindLLF	(	double *	fI,
		double *	fL,
		double *	fR,
		double *	u,
		int	dir,
		void *	s,
		double	t
	)

Compute the "upwind" source term in the balanced formulation introduced in the reference below. The "upwind" state is just the arithmetic average of the left and right states.

• Xing, Shu, "High Order Well-Balanced WENO Scheme for the Gas Dynamics Equations Under Gravitational Fields", J. Sci. Comput., 54, 2013, pp. 645–662, http://dx.doi.org/10.1007/s10915-012-9585-8.

Parameters

fI	Computed interface source term ("upwinded")
fL	Left-biased interface source term
fR	Right-biased interface source term
u	Solution (conserved variables)
dir	Spatial dimension (unused since this is a 1D case)
s	Solver object of type HyPar
t	Current solution time

Definition at line 21 of file Euler1DSourceUpwind.c.

{
  HyPar     *solver = (HyPar*)    s;
  int       done,k;
  _DECLARE_IERR_;

  int ndims = solver->ndims;
  int *dim  = solver->dim_local;

  int index_outer[ndims], index_inter[ndims], bounds_outer[ndims], bounds_inter[ndims];
  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] =  1;
  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;

  done = 0; _ArraySetValue_(index_outer,ndims,0);
  while (!done) {
    _ArrayCopy1D_(index_outer,index_inter,ndims);
    for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
      int p;  _ArrayIndex1D_(ndims,bounds_inter,index_inter,0,p);
      /* Local Lax-Friedrich upwinding scheme */
      for (k = 0; k < _MODEL_NVARS_; k++) 
        (fI+_MODEL_NVARS_*p)[k] = 0.5 * ((fL+_MODEL_NVARS_*p)[k] + (fR+_MODEL_NVARS_*p)[k]);
    }
    _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
  }

  return(0);
}

Euler1DSourceUpwindRoe()

int Euler1DSourceUpwindRoe	(	double *	fI,
		double *	fL,
		double *	fR,
		double *	u,
		int	dir,
		void *	s,
		double	t
	)

Compute the "upwind" source term in the balanced formulation introduced in the reference below. The "upwind" state is just the arithmetic average of the left and right states.

• Xing, Shu, "High Order Well-Balanced WENO Scheme for the Gas Dynamics Equations Under Gravitational Fields", J. Sci. Comput., 54, 2013, pp. 645–662, http://dx.doi.org/10.1007/s10915-012-9585-8.

Parameters

fI	Computed interface source term ("upwinded")
fL	Left-biased interface source term
fR	Right-biased interface source term
u	Solution (conserved variables)
dir	Spatial dimension (unused since this is a 1D case)
s	Solver object of type HyPar
t	Current solution time

Definition at line 64 of file Euler1DSourceUpwind.c.

{
  HyPar     *solver = (HyPar*)    s;
  int       done,k;
  _DECLARE_IERR_;

  int ndims = solver->ndims;
  int *dim  = solver->dim_local;

  int index_outer[ndims], index_inter[ndims], bounds_outer[ndims], bounds_inter[ndims];
  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] =  1;
  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;

  done = 0; _ArraySetValue_(index_outer,ndims,0);
  while (!done) {
    _ArrayCopy1D_(index_outer,index_inter,ndims);
    for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
      int p;  _ArrayIndex1D_(ndims,bounds_inter,index_inter,0,p);
      /* Local Lax-Friedrich upwinding scheme */
      for (k = 0; k < _MODEL_NVARS_; k++) 
        (fI+_MODEL_NVARS_*p)[k] = 0.5 * ((fL+_MODEL_NVARS_*p)[k] + (fR+_MODEL_NVARS_*p)[k]);
    }
    _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
  }

  return(0);
}

Euler1DModifiedSolution()

int Euler1DModifiedSolution	(	double *	uC,
		double *	u,
		int	d,
		void *	s,
		void *	m,
		double	waqt
	)

Compute the modified solution for the upwinding step in a balanced conservative finite-difference algorithm for the 1D Euler equations with gravitational sources. If no gravitational forces exist, the modified solution is identical to the solution.

Refer to

• Xing, Y., Shu, C.-W., "High Order Well-Balanced WENO Scheme for the Gas Dynamics Equations Under Gravitational Fields", Journal of Scientific Computing, 54, 2013, pp. 645-662, http://dx.doi.org/10.1007/s10915-012-9585-8. See Eq. (3.8) on Page 651 on why this modification is needed.

Parameters

uC	The modified solution (same array size and layout as u)
u	The solution (conserved variables)
d	Spatial dimension (unused since this is a 1D system)
s	Solver object of type HyPar
m	MPI object of type MPIVariables
waqt	Current solution time

Definition at line 23 of file Euler1DModifiedSolution.c.

{
  HyPar         *solver = (HyPar*)         s;
  Euler1D       *param  = (Euler1D*)       solver->physics;

  int     ghosts  = solver->ghosts;
  int     *dim    = solver->dim_local;
  int     ndims   = solver->ndims;
  int     index[ndims], bounds[ndims], offset[ndims];

  /* set bounds for array index to include ghost points */
  _ArrayCopy1D_(dim,bounds,ndims);
  int i; for (i=0; i<ndims; i++) bounds[i] += 2*ghosts;

  /* set offset such that index is compatible with ghost point arrangement */
  _ArraySetValue_(offset,ndims,-ghosts);

  int done = 0; _ArraySetValue_(index,ndims,0);
  while (!done) {
    int p; _ArrayIndex1DWO_(ndims,dim,index,offset,ghosts,p); p *= _MODEL_NVARS_;
    _ArrayScaleCopy1D_((u+p),(1.0/param->grav_field[p/_MODEL_NVARS_]),(uC+p),_MODEL_NVARS_);
    _ArrayIncrementIndex_(ndims,bounds,index,done);
  }
  return(0);
}

Euler1DPreStep()

int Euler1DPreStep	(	double *	u,
		void *	s,
		void *	m,
		double	waqt
	)

1D Euler-specific function called at the beginning of each time-step: For a simulation that splits the hyperbolic flux into its acoustic and entropy components for implicit- explicit time-integration, this function computes the fast Jacobian at the beginning of a time step for the linearized formulation.

See also

_Euler1DSetLinearizedStiffFlux_, _Euler1DSetStiffJac_, Euler1DStiffFlux

Parameters

u	Solution (conserved variables)
s	Solver object of type HyPar
m	MPI object of type MPIVariables
waqt	Current solution time

Definition at line 16 of file Euler1DPreStep.c.

{
  HyPar             *solver = (HyPar*)   s;
  Euler1D           *param  = (Euler1D*) solver->physics;
  int               *dim    = solver->dim_local;
  int               ghosts  = solver->ghosts;
  static const int  ndims   = _MODEL_NDIMS_;
  static const int  JacSize = _MODEL_NVARS_*_MODEL_NVARS_;
  static int        index[_MODEL_NDIMS_], bounds[_MODEL_NDIMS_], offset[_MODEL_NDIMS_];

  /* set bounds for array index to include ghost points */
  _ArrayAddCopy1D_(dim,(2*ghosts),bounds,ndims);
  /* set offset such that index is compatible with ghost point arrangement */
  _ArraySetValue_(offset,ndims,-ghosts);
  /* copy the solution to act as a reference for linearization */
  _ArrayCopy1D_(u,param->solution,(solver->npoints_local_wghosts*_MODEL_NVARS_));

  int done = 0; _ArraySetValue_(index,ndims,0);
  while (!done) {
    int p; _ArrayIndex1DWO_(ndims,dim,index,offset,ghosts,p);
    double rho, v, e, P;
    _Euler1DGetFlowVar_((u+_MODEL_NVARS_*p),rho,v,e,P,param);
    _Euler1DSetStiffJac_((param->fast_jac+JacSize*p),rho,v,e,P,param->gamma);
    _ArrayIncrementIndex_(ndims,bounds,index,done);
  }

  return(0);
}

Euler1DInitialize()

int Euler1DInitialize	(	void *	s,
		void *	m
	)

Function to initialize the 1D inviscid Euler equations (Euler1D) module: Sets the default parameters, read in and set physics-related parameters, and set the physics-related function pointers in HyPar.

This file reads the file "physics.inp" that must have the following format:

begin
    <keyword>   <value>
    <keyword>   <value>
    <keyword>   <value>
    ...
    <keyword>   <value>
end

where the list of keywords are:

Keyword name	Type	Variable	Default value -------—
gamma	double	Euler1D::gamma	1.4
gravity	double	Euler1D::grav	0.0
grav_type	int	Euler1D::grav_type	0
upwinding	char[]	Euler1D::upw_choice	"roe" (_ROE_)

Note: "physics.inp" is optional; if absent, default values will be used.

Parameters

s	Solver object of type HyPar
m	Object of type MPIVariables containing MPI-related info

Definition at line 71 of file Euler1DInitialize.c.

{
  HyPar         *solver  = (HyPar*)         s;
  MPIVariables  *mpi     = (MPIVariables*)  m; 
  Euler1D       *physics = (Euler1D*)       solver->physics;
  int           ferr, d;

  static int count = 0;

  if (solver->nvars != _MODEL_NVARS_) {
    fprintf(stderr,"Error in Euler1DInitialize(): nvars has to be %d.\n",_MODEL_NVARS_);
    return(1);
  }
  if (solver->ndims != _MODEL_NDIMS_) {
    fprintf(stderr,"Error in Euler1DInitialize(): ndims has to be %d.\n",_MODEL_NDIMS_);
    return(1);
  }

  /* default values */
  physics->gamma      = 1.4; 
  physics->grav       = 0.0;
  physics->grav_type  = 0;
  strcpy(physics->upw_choice,"roe");

  /* reading physical model specific inputs */
  if (!mpi->rank) {
    FILE *in;
    if (!count) printf("Reading physical model inputs from file \"physics.inp\".\n");
    in = fopen("physics.inp","r");
    if (!in) printf("Warning: File \"physics.inp\" not found. Using default values.\n");
    else {
      char word[_MAX_STRING_SIZE_];
      ferr = fscanf(in,"%s",word); if (ferr != 1) return(1);
      if (!strcmp(word, "begin")){
          while (strcmp(word, "end")){
              ferr = fscanf(in,"%s",word); if (ferr != 1) return(1);
          if (!strcmp(word, "gamma")) { 
            ferr = fscanf(in,"%lf",&physics->gamma); 
            if (ferr != 1) return(1);
          } else if (!strcmp(word, "gravity")) { 
            ferr = fscanf(in,"%lf",&physics->grav); 
            if (ferr != 1) return(1);
          } else if (!strcmp(word, "gravity_type")) { 
            ferr = fscanf(in,"%d",&physics->grav_type); 
            if (ferr != 1) return(1);
          } else if (!strcmp(word,"upwinding")) {
            ferr = fscanf(in,"%s",physics->upw_choice); 
            if (ferr != 1) return(1);
          } else if (strcmp(word,"end")) {
            char useless[_MAX_STRING_SIZE_];
            ferr = fscanf(in,"%s",useless); if (ferr != 1) return(ferr);
            printf("Warning: keyword %s in file \"physics.inp\" with value %s not ",word,useless);
            printf("recognized or extraneous. Ignoring.\n");
          }
        }
        } else {
          fprintf(stderr,"Error: Illegal format in file \"physics.inp\".\n");
        return(1);
        }
    }
    fclose(in);
  }

#ifndef serial
  IERR MPIBroadcast_double    (&physics->gamma    ,1,0,&mpi->world);                  CHECKERR(ierr);
  IERR MPIBroadcast_double    (&physics->grav     ,1,0,&mpi->world);                  CHECKERR(ierr);
  IERR MPIBroadcast_integer   (&physics->grav_type,1,0,&mpi->world);                  CHECKERR(ierr);
  IERR MPIBroadcast_character (physics->upw_choice,_MAX_STRING_SIZE_,0,&mpi->world);  CHECKERR(ierr);
#endif

  if ((physics->grav != 0.0) && (strcmp(physics->upw_choice,_LLF_)) && (strcmp(physics->upw_choice,_ROE_))) {
    if (!mpi->rank) {
      fprintf(stderr,"Error in Euler1DInitialize: %s or %s upwinding is needed for flows ",_LLF_,_ROE_);
      fprintf(stderr,"with gravitational forces.\n");
    }
    return(1);
  }

  /* initializing physical model-specific functions */
  solver->PreStep            = Euler1DPreStep;
  solver->ComputeCFL         = Euler1DComputeCFL;
  solver->FFunction          = Euler1DFlux;
  solver->SFunction          = Euler1DSource;
  solver->UFunction          = Euler1DModifiedSolution;
  if      (!strcmp(physics->upw_choice,_ROE_    )) solver->Upwind = Euler1DUpwindRoe;
  else if (!strcmp(physics->upw_choice,_RF_     )) solver->Upwind = Euler1DUpwindRF;
  else if (!strcmp(physics->upw_choice,_LLF_    )) solver->Upwind = Euler1DUpwindLLF;
  else if (!strcmp(physics->upw_choice,_SWFS_   )) solver->Upwind = Euler1DUpwindSWFS;
  else if (!strcmp(physics->upw_choice,_RUSANOV_)) solver->Upwind = Euler1DUpwindRusanov;
  else {
    if (!mpi->rank) fprintf(stderr,"Error in Euler1DInitialize(): %s is not a valid upwinding scheme.\n",
                            physics->upw_choice);
    return(1);
  }
  if (!strcmp(solver->SplitHyperbolicFlux,"yes")) {
    solver->dFFunction = Euler1DStiffFlux;
    solver->JFunction  = Euler1DStiffJacobian;
    if      (!strcmp(physics->upw_choice,_ROE_ )) solver->UpwinddF = Euler1DUpwinddFRoe;
    else if (!strcmp(physics->upw_choice,_RF_  )) solver->UpwinddF = Euler1DUpwinddFRF;
    else if (!strcmp(physics->upw_choice,_LLF_ )) solver->UpwinddF = Euler1DUpwinddFLLF;
    else {
      if (!mpi->rank) {
        fprintf(stderr,"Error in Euler1DInitialize(): %s is not a valid upwinding scheme ",
                physics->upw_choice);
        fprintf(stderr,"when split form of the hyperbolic flux is used. Use %s, %s or %s.\n",
                _ROE_,_RF_,_LLF_);
      }
      return(1);
    }
  } else {
    solver->dFFunction = NULL;
    solver->UpwinddF   = NULL;
    solver->JFunction  = Euler1DJacobian;
  }
  solver->AveragingFunction     = Euler1DRoeAverage;
  solver->GetLeftEigenvectors   = Euler1DLeftEigenvectors;
  solver->GetRightEigenvectors  = Euler1DRightEigenvectors;
   
  if      (!strcmp(physics->upw_choice,_LLF_ )) physics->SourceUpwind = Euler1DSourceUpwindLLF;
  else if (!strcmp(physics->upw_choice,_ROE_ )) physics->SourceUpwind = Euler1DSourceUpwindRoe;

  /* allocate array to hold the gravity field */
  int *dim    = solver->dim_local;
  int ghosts  = solver->ghosts;
  int size = 1; for (d=0; d<_MODEL_NDIMS_; d++) size *= (dim[d] + 2*ghosts);
  physics->grav_field = (double*) calloc (size, sizeof(double));
  physics->fast_jac   = (double*) calloc (size*_MODEL_NVARS_*_MODEL_NVARS_, sizeof(double));
  physics->solution   = (double*) calloc (size*_MODEL_NVARS_, sizeof(double));
  IERR Euler1DGravityField(solver,mpi); CHECKERR(ierr);

  count++;
  return(0);
}