LinearADRInitialize.c File Reference

Initialize the linear advection-diffusion-reaction module. More...

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

Go to the source code of this file.

Functions
int	LinearADRAdvectionField (void *, void *, int, int, int *)
double	LinearADRComputeCFL (void *, void *, double, double)
double	LinearADRComputeDiffNumber (void *, void *, double, double)
int	LinearADRAdvection (double *, double *, int, void *, double)
int	LinearADRDiffusionG (double *, double *, int, void *, double)
int	LinearADRDiffusionH (double *, double *, int, int, void *, double)
int	LinearADRReaction ()
int	LinearADRUpwind (double *, double *, double *, double *, double *, double *, int, void *, double)
int	LinearADRCenteredFlux (double *, double *, double *, double *, double *, double *, int, void *, double)
int	LinearADRWriteAdvField (void *, void *, double)
int	LinearADRAdvectionJacobian (double *, double *, void *, int, int, int)
int	LinearADRDiffusionJacobian (double *, double *, void *, int, int)
int	LinearADRInitialize (void *s, void *m)

Detailed Description

Initialize the linear advection-diffusion-reaction module.

Author

Debojyoti Ghosh

Definition in file LinearADRInitialize.c.

Function Documentation

int LinearADRAdvectionField	(	void *	s,
		void *	m,
		int	idx,
		int	nsims,
		int *	dim_adv
	)

Set the advection field over the domain - reads the spatially-varying advection data from a file, if available. The array to store the field must already be allocated.

For a single simulation, it reads in the data from a file. The data must have the same grid dimensions as the solution.

For an ensemble simulation, it reads in the advection field from the file corresponding to the idx (index of this simulation). The data in this file must have the same grid dimensions as the domain it is being read in for.

For a sparse grids simulation, it reads in the advection field from the file with data that has the grid dimensions as the full grid. The field on the current grid is obtained by interpolation.

Parameters

s	Solver object of type HyPar
m	MPI object of type MPIVariables
idx	Index of this simulation
nsims	Total number of simulations
dim_adv	Grid dimensions of the advection field stored in file

Definition at line 33 of file LinearADRAdvectionField.c.

{
  HyPar         *solver = (HyPar*) s;
  MPIVariables  *mpi    = (MPIVariables*) m;
  LinearADR     *param  = (LinearADR*) solver->physics;
  double        *adv    = param->a;

  /* sanity checks */
  if (param->constant_advection == 1) {
    fprintf(stderr,"Error in LinearADRAdvectionField(): param->constant_advection == 1!\n");
    return(1);
  }
  if (!strcmp(param->adv_filename,"none")) {
    fprintf(stderr,"Error in LinearADRAdvectionField(): param->adv_filename is \"none\"!\n");
    return(1);
  }

  int d, done, flag;
  int *dim = solver->dim_local;
  int ghosts = solver->ghosts;

  /* number of advection field components at a grid point is
   * number of spatial dimensions times number of solution 
   * components */
  int adv_nvar = solver->ndims * solver->nvars;
  if (param->adv_arr_size != adv_nvar*solver->npoints_local_wghosts) {
    fprintf(stderr,"Error in LinearADRAdvectionField(): Incorrect adv_arr_size\n");
    return(1);
  }

  char fname_root[_MAX_STRING_SIZE_];
  strcpy(fname_root, param->adv_filename);
  if (idx >= 0) {
    if (nsims > 1) {
      char index[_MAX_STRING_SIZE_];
      GetStringFromInteger(idx, index, (int)log10(nsims)+1);
      strcat(fname_root, "_");
      strcat(fname_root, index);
    }
  }

  /* read spatially-varying advection field from provided file */
  if (dim_adv == NULL) {
    int ierr = ReadArray( solver->ndims,
                          adv_nvar,
                          solver->dim_global,
                          solver->dim_local,
                          solver->ghosts,
                          solver,
                          mpi,
                          NULL,
                          adv,
                          fname_root,
                          &flag);
    if (ierr) {
      fprintf(stderr,"Error in LinearADRAdvectionField()\n");
      fprintf(stderr,"  ReadArray() returned error!\n");
      return ierr;
    }
  } else {
    int ierr = ReadArraywInterp(solver->ndims,
                                adv_nvar,
                                solver->dim_global,
                                solver->dim_local,
                                dim_adv,
                                solver->ghosts,
                                solver,
                                mpi,
                                NULL,
                                adv,
                                fname_root,
                                &flag);
    if (ierr) {
      fprintf(stderr,"Error in LinearADRAdvectionField()\n");
      fprintf(stderr,"  ReadArraywInterp() returned error!\n");
      return ierr;
    }
  }

  if (!flag) {
    /* if advection file not available, set it to zero */
    _ArraySetValue_(adv,solver->npoints_local_wghosts*adv_nvar,0.0);
  }

  /* if parallel, exchange MPI-boundary ghost point data */
  IERR MPIExchangeBoundariesnD( solver->ndims,
                                adv_nvar,
                                solver->dim_local,
                                solver->ghosts,
                                mpi,
                                adv );

  /* Along each dimension:
     - If boundaries are periodic, then set ghost points accordingly if
       number of MPI ranks along that dimension is 1 (if more than 1, 
       MPIExchangeBoundariesnD() called above has taken care of it).
     - Else, extrapolate at the boundaries. */
  int indexb[solver->ndims], 
      indexi[solver->ndims], 
      bounds[solver->ndims], 
      offset[solver->ndims];
  for (d = 0; d < solver->ndims; d++) {
    if (solver->isPeriodic[d] && (mpi->iproc[d] == 1)) {
      _ArrayCopy1D_(dim,bounds,solver->ndims); bounds[d] = ghosts;
      /* left boundary */
      done = 0; _ArraySetValue_(indexb,solver->ndims,0);
      _ArraySetValue_(offset,solver->ndims,0); offset[d] = -ghosts;
      while (!done) {
        _ArrayCopy1D_(indexb,indexi,solver->ndims); indexi[d] = indexb[d] + dim[d] - ghosts;
        int p1; _ArrayIndex1DWO_(solver->ndims,dim,indexb,offset,ghosts,p1);
        int p2; _ArrayIndex1D_  (solver->ndims,dim,indexi,ghosts,p2);
        _ArrayCopy1D_((adv+p2*adv_nvar), (adv+p1*adv_nvar), adv_nvar);
        _ArrayIncrementIndex_(solver->ndims,bounds,indexb,done);
      }
      /* right boundary */
      done = 0; _ArraySetValue_(indexb,solver->ndims,0);
      _ArraySetValue_(offset,solver->ndims,0); offset[d] = dim[d];
      while (!done) {
        _ArrayCopy1D_(indexb,indexi,solver->ndims);
        int p1; _ArrayIndex1DWO_(solver->ndims,dim,indexb,offset,ghosts,p1);
        int p2; _ArrayIndex1D_  (solver->ndims,dim,indexi,ghosts,p2);
        _ArrayCopy1D_((adv+p2*adv_nvar), (adv+p1*adv_nvar), adv_nvar);
        _ArrayIncrementIndex_(solver->ndims,bounds,indexb,done);
      }
    } else {
      /* left boundary */
      if (!mpi->ip[d]) {
        _ArrayCopy1D_(dim,bounds,solver->ndims); bounds[d] = ghosts;
        _ArraySetValue_(offset,solver->ndims,0); offset[d] = -ghosts;
        done = 0; _ArraySetValue_(indexb,solver->ndims,0);
        while (!done) {
          _ArrayCopy1D_(indexb,indexi,solver->ndims); indexi[d] = ghosts-1-indexb[d];
          int p1; _ArrayIndex1DWO_(solver->ndims,dim,indexb,offset,ghosts,p1);
          int p2; _ArrayIndex1D_  (solver->ndims,dim,indexi,ghosts,p2);
          _ArrayCopy1D_((adv+p2*adv_nvar), (adv+p1*adv_nvar), adv_nvar);
          _ArrayIncrementIndex_(solver->ndims,bounds,indexb,done);
        }
      }
      /* right boundary */
      if (mpi->ip[d] == mpi->iproc[d]-1) {
        _ArrayCopy1D_(dim,bounds,solver->ndims); bounds[d] = ghosts;
        _ArraySetValue_(offset,solver->ndims,0); offset[d] = dim[d];
        done = 0; _ArraySetValue_(indexb,solver->ndims,0);
        while (!done) {
          _ArrayCopy1D_(indexb,indexi,solver->ndims); indexi[d] = dim[d]-1-indexb[d];
          int p1; _ArrayIndex1DWO_(solver->ndims,dim,indexb,offset,ghosts,p1);
          int p2; _ArrayIndex1D_  (solver->ndims,dim,indexi,ghosts,p2);
          _ArrayCopy1D_((adv+p2*adv_nvar), (adv+p1*adv_nvar), adv_nvar);
          _ArrayIncrementIndex_(solver->ndims,bounds,indexb,done);
        }
      }
    }
  }

  return(0);
}

double LinearADRComputeCFL	(	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 15 of file LinearADRComputeCFL.c.

{
  HyPar         *solver = (HyPar*)        s;
  LinearADR     *params = (LinearADR*)    solver->physics;

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

  double  max_cfl = 0;
  if (params->constant_advection == 1) {

    int d, i, v;
    for (d = 0; d < ndims; d++) {
      for (i = 0; i < dim[d]; i++) {
        for (v = 0; v < nvars; v++) {
          double dxinv; _GetCoordinate_(d,i,dim,ghosts,solver->dxinv,dxinv);
          double local_cfl = params->a[nvars*d+v]*dt*dxinv;
          if (local_cfl > max_cfl) max_cfl = local_cfl;
        }
      }
    }

  } else if (params->constant_advection == 0) {

    int d;
    for (d = 0; d < ndims; d++) {
      int index[ndims];
      int done = 0; _ArraySetValue_(index,ndims,0);
      while (!done) {
        int p; _ArrayIndex1D_(ndims,dim,index,ghosts,p);
        double dxinv; _GetCoordinate_(0,index[0],dim,ghosts,solver->dxinv,dxinv);
        int v;
        for (v = 0; v < nvars; v++) {
          double a = params->a[nvars*ndims*p+nvars*d+v];
          double local_cfl = a*dt*dxinv; 
          if (local_cfl > max_cfl) max_cfl = local_cfl;
        }
        _ArrayIncrementIndex_(ndims,dim,index,done);
      }
    }

  }

  return(max_cfl);
}

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

Computes the maximum diffusion number over the domain. Note that the diffusion number 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 14 of file LinearADRComputeDiffNumber.c.

{
  HyPar         *solver = (HyPar*)        s;
  LinearADR     *params = (LinearADR*)    solver->physics;
  int           d, i, v;

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

  double  max_diffno = 0;
  for (d = 0; d < ndims; d++) {
    for (i = 0; i < dim[d]; i++) {
      for (v = 0; v < nvars; v++) {
        double dxinv;  _GetCoordinate_(d,i,dim,ghosts,solver->dxinv,dxinv);
        double local_diffno =   params->d[nvars*d+v] * dt * dxinv * dxinv;
        if (local_diffno > max_diffno) max_diffno = local_diffno;
      }
    }
  }

  return(max_diffno);
}

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

Evaluate the advection term in the linear advection-diffusion-reaction model:
Compute

\begin{equation} a_d u \end{equation}

given \(u\) and \(d\) in the hyperbolic term

\begin{equation} \sum_d \frac {\partial} {\partial x_d} \left( a_d u \right) \end{equation}

• If the advection field is spatially varying, then the advection speed array (LinearADR::a) is expected to have the size (LinearADR::adv_arr_size): (number of grid points with ghosts - HyPar::npoints_local_wghosts) X (number of spatial dimensions - HyPar::ndims) X (number of solution components - HyPar::nvars), with the innermost loop being the solution components and outermost loop being the grid points.
• If the advection field is constant, then the advection speed array (LinearADR::a) is expected to have the size (LinearADR::adv_arr_size): (number of spatial dimensions - HyPar::ndims) X (number of solution components - HyPar::nvars), with the inner loop being the solution components.

Parameters

f	Array to hold the computed flux (same size and layout as u)
u	Array containing the solution
dir	Spatial dimension \(d\)
s	Solver object of type HyPar
t	Current time

Definition at line 37 of file LinearADRAdvection.c.

{
  HyPar     *solver = (HyPar*)     s;
  LinearADR *param  = (LinearADR*) solver->physics;
  double    *adv    = param->a;
  int       i, v;

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

  int index[ndims], 
      bounds[ndims], 
      offset[ndims];

  /* set bounds for array index to include ghost points */
  _ArrayCopy1D_(dim,bounds,ndims);
  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);
  if (param->constant_advection == 1) {
    while (!done) {
      int p; _ArrayIndex1DWO_(ndims,dim,index,offset,ghosts,p);
      for (v = 0; v < nvars; v++) {
        f[nvars*p+v] = param->a[nvars*dir+v] * u[nvars*p+v];
      }
      _ArrayIncrementIndex_(ndims,bounds,index,done);
    }
  } else if (param->constant_advection == 0) {
    while (!done) {
      int p; _ArrayIndex1DWO_(ndims,dim,index,offset,ghosts,p);
      for (v = 0; v < nvars; v++) {
        double a = adv[nvars*ndims*p+nvars*dir+v];
        f[nvars*p+v] = a * u[nvars*p+v];
      }
      _ArrayIncrementIndex_(ndims,bounds,index,done);
    }
  } else {
    while (!done) {
      int p; _ArrayIndex1DWO_(ndims,dim,index,offset,ghosts,p);
      for (v = 0; v < nvars; v++) {
        f[nvars*p+v] = 0.0;
      }
      _ArrayIncrementIndex_(ndims,bounds,index,done);
    }
  }

  return(0);
}

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

Evaluate the diffusion term in the linear advection-diffusion-reaction model for a "pure Laplacian" type operator (no cross derivatives):
Compute

\begin{equation} \nu_d u \end{equation}

given \(u\) and \(d\) in the parabolic term

\begin{equation} \sum_d \frac {\partial^2} {\partial x_d^2} \left( \nu_d u \right) \end{equation}

Parameters

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

Definition at line 26 of file LinearADRDiffusion.c.

{
  HyPar     *solver = (HyPar*)     s;
  LinearADR *param  = (LinearADR*) solver->physics;
  int       i, v;

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

  /* set bounds for array index to include ghost points */
  _ArrayCopy1D_(dim,bounds,ndims);
  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);
    for (v = 0; v < nvars; v++) f[nvars*p+v] = param->d[nvars*dir+v] * u[nvars*p+v];
    _ArrayIncrementIndex_(ndims,bounds,index,done);
  }

  return(0);
}

int LinearADRDiffusionH	(	double *	f,
		double *	u,
		int	dir1,
		int	dir2,
		void *	s,
		double	t
	)

Evaluate the diffusion term in the linear advection-diffusion-reaction model for a parabolic operator with no cross derivatives:
Compute

\begin{equation} \nu_d u \end{equation}

given \(u\) and \(d_1,d_2\) in the parabolic term

\begin{equation} \sum_{d_1}\sum_{d_2} \frac {\partial^2} {\partial x_{d_1} \partial x_{d_2}} \left( \nu_d u \right) \end{equation}

Note: it's not correctly implemented. Will implement when necessary.

Parameters

f	Array to hold the computed diffusion term (same size and layout as u)
u	Array containing the solution
dir1	First spatial dimension of the double derivative \(d_1\)
dir2	Second spatial dimension of the double derivative \(d_2\)
s	Solver object of type HyPar
t	Current time

Definition at line 74 of file LinearADRDiffusion.c.

{
  HyPar     *solver = (HyPar*)     s;
  LinearADR *param  = (LinearADR*) solver->physics;
  int       i, v;

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

  if (dir1 == dir2) {

    int dir = dir1;

    /* set bounds for array index to include ghost points */
    _ArrayCopy1D_(dim,bounds,ndims);
    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);
      for (v = 0; v < nvars; v++) f[nvars*p+v] = param->d[nvars*dir+v] * u[nvars*p+v];
      _ArrayIncrementIndex_(ndims,bounds,index,done);
    }

  } else _ArraySetValue_(f,solver->npoints_local_wghosts*nvars,0.0);


  return(0);
}

int LinearADRReaction

(

)

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

Upwinding scheme for linear advection

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
s	Solver object of type HyPar
t	Current solution time

Definition at line 16 of file LinearADRUpwind.c.

{
  HyPar     *solver = (HyPar*)      s;
  LinearADR *param  = (LinearADR*)  solver->physics;
  int       done,v;

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

  double *a = param->a;

  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;

  if (param->constant_advection == 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);
        for (v = 0; v < nvars; v++) {
          fI[nvars*p+v] = (a[nvars*dir+v] > 0 ? fL[nvars*p+v] : fR[nvars*p+v] );
        }
      }
      _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
    }

  } else if (param->constant_advection == 0) {

    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);
        for (v = 0; v < nvars; v++) {
          double eigL = a[nvars*ndims*pL+nvars*dir+v],
                 eigR = a[nvars*ndims*pR+nvars*dir+v];
          if ((eigL > 0) && (eigR > 0)) {
             fI[nvars*p+v] = fL[nvars*p+v];
          } else if ((eigL < 0) && (eigR < 0)) {
             fI[nvars*p+v] = fR[nvars*p+v];
          } else { 
             double alpha = max(absolute(eigL), absolute(eigR));
             fI[nvars*p+v] = 0.5 * (fL[nvars*p+v] + fR[nvars*p+v] - alpha * (uR[nvars*p+v] - uL[nvars*p+v]));
          }
        }
      }
      _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
    }

  } else {

    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);
        for (v = 0; v < nvars; v++) {
          fI[nvars*p+v] = 0.0;
        }
      }
      _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
    }

  }

  return(0);
}

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

Centered scheme for linear advection

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
s	Solver object of type HyPar
t	Current solution time

Definition at line 103 of file LinearADRUpwind.c.

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

  int ndims = solver->ndims;
  int nvars = solver->nvars;
  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;
    
  int 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);
      for (int v = 0; v < nvars; v++) {
        fI[nvars*p+v] = 0.5 * (fL[nvars*p+v] + fR[nvars*p+v]);
      }
    }
    _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
  }

  return(0);
}

int LinearADRWriteAdvField	(	void *	s,
		void *	m,
		double	a_t
	)

Write out the advection field to file

Parameters

s	Solver object of type HyPar
m	MPI object of type MPIVariables
a_t	Current simulation time

Definition at line 18 of file LinearADRWriteAdvField.c.

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

  if (param->constant_advection == 0) {

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

    int adv_nvar = solver->ndims * solver->nvars;
    WriteArray(  solver->ndims,
                 adv_nvar,
                 solver->dim_global,
                 solver->dim_local,
                 solver->ghosts,
                 solver->x,
                 param->a,
                 solver,mpi,
                 fname_root );

    if (!strcmp(solver->plot_solution, "yes")) {
      PlotArray( solver->ndims,
                 adv_nvar,
                 solver->dim_global,
                 solver->dim_local,
                 solver->ghosts,
                 solver->x,
                 param->a,
                 a_t,
                 solver,mpi,
                 fname_root );
    }

  }

  return(0);
}

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

Function to compute the flux Jacobian for the hyperbolic (advection) part of the linear-advection-diffusion-reaction model.

Parameters

Jac	Jacobian matrix of size 1 (nvar = 1)
u	solution at a grid point
p	object containing physics-related parameters
dir	dimension (x/y/z)
nvars	number of 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 11 of file LinearADRJacobian.c.

{
  LinearADR *param = (LinearADR*) p;

  if (param->a) {
    *Jac =    (1-absolute(upw))*absolute(param->a[dir])
           +  absolute(upw) * (1+upw) * max(0,param->a[dir]) * 0.5
           -  absolute(upw) * (1-upw) * min(0,param->a[dir]) * 0.5 ;
  } else {
    /* no advection term */
    *Jac = 0.0;
  }

  return 0;
}

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

Function to compute the Jacobian for the parabolic (diffusion) part of the linear-advection-diffusion-reaction model.

Parameters

Jac	Jacobian matrix of size 1 (nvar = 1)
u	solution at a grid point
p	object containing physics-related parameters
dir	dimension (x/y/z)
nvars	number of components

Definition at line 36 of file LinearADRJacobian.c.

{
  LinearADR *param = (LinearADR*) p;

  int v;
  for (v = 0; v < nvars; v++) {
    Jac[nvars*v+v] = -param->d[nvars*dir+v];
  }

  return 0;
}

int LinearADRInitialize	(	void *	s,
		void *	m
	)

Initialize the linear advection-diffusion-reaction physics module - allocate and set physics-related parameters, read physics-related inputs from file, 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
advection_filename	char[]	LinearADR::adv_filename	"none"
advection	double[]	LinearADR::a	0
diffusion	double[]	LinearADR::d	0
centered_flux	char[]	LinearADR::centered_flux	"no"

Note:

• "physics.inp" is optional; if absent, default values will be used.
• Please do not specify both "advection" and "advection_filename"!

Parameters

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

Definition at line 59 of file LinearADRInitialize.c.

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

  static int count = 0;

  /* default values */
  physics->constant_advection = -1;
  strcpy(physics->adv_filename,"none");
  physics->a = NULL;
  physics->adv_arr_size = -1;
  strcpy(physics->centered_flux,"no");

  physics->d = (double*) calloc (solver->ndims*solver->nvars,sizeof(double));
  _ArraySetValue_(physics->d,solver->ndims*solver->nvars,0.0);

  /* reading physical model specific inputs - all processes */
  if (!mpi->rank) {

    FILE *in;
    if (!count) printf("Reading physical model inputs from file \"physics.inp\".\n");
    in = fopen("physics.inp","r");

    if (!in) {

      fprintf(stderr,"Error: File \"physics.inp\" not found.\n");
      return(1);

    } 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, "advection_filename")) {
            if (physics->constant_advection != -1) {
              fprintf(stderr,"Error in LinearADRInitialize():\n");
              fprintf(stderr,"Maybe advection_filename and advection are both specified.\n");
              return 1;
            } else {
              /* read filename for spatially-varying advection field */
              physics->constant_advection = 0;
              physics->adv_arr_size = solver->npoints_local_wghosts*solver->ndims*solver->nvars;
              physics->a = (double*) calloc ( physics->adv_arr_size, sizeof(double) );
              _ArraySetValue_(physics->a, physics->adv_arr_size, 0.0);
              ferr = fscanf(in,"%s",physics->adv_filename); if (ferr != 1) return(ferr);
            }
          } else if (!strcmp(word, "advection")) {
            if (physics->constant_advection != -1) {
              fprintf(stderr,"Error in LinearADRInitialize():\n");
              fprintf(stderr,"Maybe advection_filename and advection are both specified.\n");
              return 1;
            } else {
              /* read advection coefficients */
              physics->constant_advection = 1;
              physics->adv_arr_size = solver->ndims*solver->nvars;
              physics->a = (double*) calloc ( physics->adv_arr_size, sizeof(double) );
              for (i=0; i<solver->ndims*solver->nvars; i++) ferr = fscanf(in,"%lf",&physics->a[i]);
              if (ferr != 1) return(1);
            }
          } else if (!strcmp(word, "diffusion")) {
            /* read diffusion coefficients */
            for (i=0; i<solver->ndims*solver->nvars; i++) ferr = fscanf(in,"%lf",&physics->d[i]);
            if (ferr != 1) return(1);
          } else if (!strcmp(word, "centered_flux")) {
            ferr = fscanf(in, "%s", physics->centered_flux); 
            if (ferr != 1) return(ferr);
          } 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
  MPIBroadcast_integer(&physics->constant_advection,1,0,&mpi->world);
  MPIBroadcast_integer(&physics->adv_arr_size,1,0,&mpi->world);
#endif

  if (mpi->rank) {
    physics->a = (double*) calloc (physics->adv_arr_size,sizeof(double));
    _ArraySetValue_(physics->a, physics->adv_arr_size, 0.0);
  }

  if (physics->constant_advection == 1) {
#ifndef serial
    MPIBroadcast_double(physics->a,physics->adv_arr_size,0,&mpi->world);
#endif
  } else if (physics->constant_advection == 0) {
#ifndef serial
    MPIBroadcast_character(physics->adv_filename, _MAX_STRING_SIZE_,0,&mpi->world);
#endif
  }

#ifndef serial
  MPIBroadcast_double(physics->d,solver->ndims*solver->nvars,0,&mpi->world);
  MPIBroadcast_character(physics->centered_flux, _MAX_STRING_SIZE_,0,&mpi->world);
#endif

  if (!strcmp(solver->SplitHyperbolicFlux,"yes")) {
    if (!mpi->rank) {
      fprintf(stderr,"Error in LinearADRInitialize: This physical model does not have a splitting ");
      fprintf(stderr,"of the hyperbolic term defined.\n");
    }
    return(1);
  }

  /* initializing physical model-specific functions */
  solver->ComputeCFL         = LinearADRComputeCFL;
  solver->ComputeDiffNumber  = LinearADRComputeDiffNumber;
  solver->FFunction          = LinearADRAdvection;
  solver->GFunction          = LinearADRDiffusionG;
  solver->HFunction          = LinearADRDiffusionH;
  solver->SFunction          = LinearADRReaction;
  solver->JFunction          = LinearADRAdvectionJacobian;
  solver->KFunction          = LinearADRDiffusionJacobian;

  if (!strcmp(physics->centered_flux,"no")) {
    solver->Upwind = LinearADRUpwind;
  } else {
    solver->Upwind = LinearADRCenteredFlux;
  }

  if (physics->constant_advection == 0) {
    solver->PhysicsInput = LinearADRAdvectionField;
    solver->PhysicsOutput = LinearADRWriteAdvField;
  } else {
    solver->PhysicsInput = NULL;
    solver->PhysicsOutput = NULL;
  }

  count++;
  return(0);
}