boundaryconditions.h File Reference

Containts the structures and definitions for boundary condition implementation. More...

#include <basic.h>

Go to the source code of this file.

Data Structures
struct	DomainBoundary
	Structure containing the variables and function pointers defining a boundary. More...

Macros
#define	_PERIODIC_   "periodic"
#define	_EXTRAPOLATE_   "extrapolate"
#define	_DIRICHLET_   "dirichlet"
#define	_REFLECT_   "reflect"
#define	_SPONGE_   "sponge"
#define	_NOSLIP_WALL_   "noslip-wall"
#define	_THERMAL_NOSLIP_WALL_   "thermal-noslip-wall"
#define	_SLIP_WALL_   "slip-wall"
#define	_THERMAL_SLIP_WALL_   "thermal-slip-wall"
#define	_SUBSONIC_INFLOW_   "subsonic-inflow"
#define	_SUBSONIC_OUTFLOW_   "subsonic-outflow"
#define	_SUBSONIC_AMBIVALENT_   "subsonic-ambivalent"
#define	_SUPERSONIC_INFLOW_   "supersonic-inflow"
#define	_SUPERSONIC_OUTFLOW_   "supersonic-outflow"
#define	_TURBULENT_SUPERSONIC_INFLOW_   "turbulent-supersonic-inflow"
#define	_NO_FLUX_BC_   "numa-nfbc"
#define	_SW_SLIP_WALL_   "shallow-water-slip-wall"
#define	_SW_NOSLIP_WALL_   "shallow-water-noslip-wall"

Functions
int	BCInitialize (void *, int)
int	BCCleanup (void *, int)
int	BCPeriodicU (void *, void *, int, int, int *, int, double *, double)
int	BCExtrapolateU (void *, void *, int, int, int *, int, double *, double)
int	BCDirichletU (void *, void *, int, int, int *, int, double *, double)
int	BCReflectU (void *, void *, int, int, int *, int, double *, double)
int	BCNoslipWallU (void *, void *, int, int, int *, int, double *, double)
int	BCThermalNoslipWallU (void *, void *, int, int, int *, int, double *, double)
int	BCSlipWallU (void *, void *, int, int, int *, int, double *, double)
int	BCThermalSlipWallU (void *, void *, int, int, int *, int, double *, double)
int	BCSubsonicInflowU (void *, void *, int, int, int *, int, double *, double)
int	BCSubsonicOutflowU (void *, void *, int, int, int *, int, double *, double)
int	BCSubsonicAmbivalentU (void *, void *, int, int, int *, int, double *, double)
int	BCSupersonicInflowU (void *, void *, int, int, int *, int, double *, double)
int	BCSupersonicOutflowU (void *, void *, int, int, int *, int, double *, double)
int	BCTurbulentSupersonicInflowU (void *, void *, int, int, int *, int, double *, double)
int	BCNoFluxU (void *, void *, int, int, int *, int, double *, double)
int	BCSWSlipWallU (void *, void *, int, int, int *, int, double *, double)
int	BCSpongeSource (void *, int, int, int, int *, double *, double *, double *)
int	BCSpongeUDummy (void *, void *, int, int, int *, int, double *, double)
int	BCReadTemperatureData (void *, void *, int, int, int *)
int	BCReadTurbulentInflowData (void *, void *, int, int, int *)
int	gpuBCPeriodicU (void *, void *, int, int, int *, int, double *, double)
int	gpuBCSlipWallU (void *, void *, int, int, int *, int, double *, double)

Detailed Description

Containts the structures and definitions for boundary condition implementation.

Author

Debojyoti Ghosh

Definition in file boundaryconditions.h.

Macro Definition Documentation

_PERIODIC_

#define _PERIODIC_   "periodic"

Periodic boundary conditions

See also

BCPeriodicU

Definition at line 12 of file boundaryconditions.h.

_EXTRAPOLATE_

#define _EXTRAPOLATE_   "extrapolate"

Extrapolative boundary conditions (values at ghost cells are copied from interior)

See also

BCExtrapolateU

Definition at line 14 of file boundaryconditions.h.

_DIRICHLET_

#define _DIRICHLET_   "dirichlet"

Dirichlet boundary conditions (values at ghost cells specified through input)

See also

BCDirichletU

Definition at line 16 of file boundaryconditions.h.

_REFLECT_

#define _REFLECT_   "reflect"

Reflective boundary conditions (values at ghost cells negative of interior)

See also

BCReflectU

Definition at line 18 of file boundaryconditions.h.

_SPONGE_

#define _SPONGE_   "sponge"

Sponge boundary conditions

See also

BCSpongeSource, BCSpongeUDummy

Definition at line 20 of file boundaryconditions.h.

_NOSLIP_WALL_

#define _NOSLIP_WALL_   "noslip-wall"

Viscous wall boundary condition (specific to Navier-Stokes)

See also

BCNoslipWallU

Definition at line 24 of file boundaryconditions.h.

_THERMAL_NOSLIP_WALL_

#define _THERMAL_NOSLIP_WALL_   "thermal-noslip-wall"

Viscous thermal wall boundary condition where wall temperature is specified (specific to Navier-Stokes)

See also

BCThermalNoslipWallU

Definition at line 26 of file boundaryconditions.h.

_SLIP_WALL_

#define _SLIP_WALL_   "slip-wall"

Inviscid wall boundary condition (specific to Euler/Navier-Stokes)

See also

BCSlipWallU

Definition at line 28 of file boundaryconditions.h.

_THERMAL_SLIP_WALL_

#define _THERMAL_SLIP_WALL_   "thermal-slip-wall"

Inviscid thermal wall boundary condition where wall temperature is specified (specific to Euler/Navier-Stokes)

See also

BCThermalSlipWallU

Definition at line 30 of file boundaryconditions.h.

_SUBSONIC_INFLOW_

#define _SUBSONIC_INFLOW_   "subsonic-inflow"

Subsonic inflow boundary condition: density and velocity are specified in the input, pressure is extrapolated from the interior (specific to Euler/Navier-Stokes)

See also

BCSubsonicInflowU

Definition at line 32 of file boundaryconditions.h.

_SUBSONIC_OUTFLOW_

#define _SUBSONIC_OUTFLOW_   "subsonic-outflow"

Subsonic outflow boundary condition: pressure is specified in the input, density and velocity are extrapolated from the interior (specific to Euler/Navier-Stokes)

See also

BCSubsonicOutflowU

Definition at line 34 of file boundaryconditions.h.

_SUBSONIC_AMBIVALENT_

#define _SUBSONIC_AMBIVALENT_   "subsonic-ambivalent"

Subsonic "ambivalent" boundary condition: (specific to Euler/Navier-Stokes) the velocity at the boundary is extrapolated from the interior, and based on that, either subsonic outflow or inflow boundary conditions are applied. Specify all flow quantities (density, velocity, and pressure) in the input; depending on whether it is inflow or outflow, the appropriate quantities will be used.

See also

BCSubsonicAmbivalentU

Definition at line 39 of file boundaryconditions.h.

_SUPERSONIC_INFLOW_

#define _SUPERSONIC_INFLOW_   "supersonic-inflow"

Supersonic inflow boundary condition: density, velocity, and pressure are specified in the input (specific to Euler/Navier-Stokes)

See also

BCSupersonicInflowU

Definition at line 41 of file boundaryconditions.h.

_SUPERSONIC_OUTFLOW_

#define _SUPERSONIC_OUTFLOW_   "supersonic-outflow"

Supersonic outflow boundary condition: all flow quantities are extrapolated from the interior (specific to Euler/Navier-Stokes)

See also

BCSupersonicOutflowU

Definition at line 43 of file boundaryconditions.h.

_TURBULENT_SUPERSONIC_INFLOW_

#define _TURBULENT_SUPERSONIC_INFLOW_   "turbulent-supersonic-inflow"

Turbulent, supersonic inflow boundary condition: density, velocity, and pressure are specified in the input, along with turbulent fluctuations (specific to Euler/Navier-Stokes)

See also

BCTurbulentSupersonicInflowU, BCReadTurbulentInflowData

Definition at line 45 of file boundaryconditions.h.

_NO_FLUX_BC_

#define _NO_FLUX_BC_   "numa-nfbc"

No-flux boundary condition (specific to NUMA)

See also

BCNoFluxU

Definition at line 49 of file boundaryconditions.h.

_SW_SLIP_WALL_

#define _SW_SLIP_WALL_   "shallow-water-slip-wall"

Slip boundary condition (specific to shallow water equations)

See also

BCSWSlipWallU

Definition at line 53 of file boundaryconditions.h.

_SW_NOSLIP_WALL_

#define _SW_NOSLIP_WALL_   "shallow-water-noslip-wall"

Viscous wall boundary condition (specific to shallow water equations) (not implemented yet)

Definition at line 55 of file boundaryconditions.h.

Function Documentation

BCInitialize()

int BCInitialize	(	void *	b,
		int	flag_gpu
	)

Function to initialize the boundary conditions

Assign the function pointers for boundary condition application depending on the boundary type, for a given boundary object

Parameters

b	Boundary object of type DomainBoundary
flag_gpu	Flag to indicate if GPU is being used

Definition at line 12 of file BCInitialize.c.

{
  DomainBoundary *boundary = (DomainBoundary*) b;

#if defined(HAVE_CUDA)
  if (flag_gpu) {

    if      (!strcmp(boundary->bctype,_PERIODIC_)) {
      boundary->BCFunctionU = gpuBCPeriodicU;
    } else if (!strcmp(boundary->bctype,_SLIP_WALL_)) {
      boundary->BCFunctionU = gpuBCSlipWallU;
    } else {
      fprintf(stderr,"[GPU] Error in BCInitialize(): \"%s\" is not a supported boundary condition.\n",
              boundary->bctype);
      return 1;
    }

  } else {
#endif

    if      (!strcmp(boundary->bctype,_PERIODIC_                    )) boundary->BCFunctionU = BCPeriodicU;
    else if (!strcmp(boundary->bctype,_EXTRAPOLATE_                 )) boundary->BCFunctionU = BCExtrapolateU;
    else if (!strcmp(boundary->bctype,_DIRICHLET_                   )) boundary->BCFunctionU = BCDirichletU;  
    else if (!strcmp(boundary->bctype,_REFLECT_                     )) boundary->BCFunctionU = BCReflectU;    
    else if (!strcmp(boundary->bctype,_SPONGE_                      )) boundary->BCFunctionU = BCSpongeUDummy;    
    else if (!strcmp(boundary->bctype,_NOSLIP_WALL_                 )) boundary->BCFunctionU = BCNoslipWallU;    
    else if (!strcmp(boundary->bctype,_SLIP_WALL_                   )) boundary->BCFunctionU = BCSlipWallU;    
    else if (!strcmp(boundary->bctype,_THERMAL_SLIP_WALL_           )) boundary->BCFunctionU = BCThermalSlipWallU;    
    else if (!strcmp(boundary->bctype,_THERMAL_NOSLIP_WALL_         )) boundary->BCFunctionU = BCThermalNoslipWallU;    
    else if (!strcmp(boundary->bctype,_SW_SLIP_WALL_                )) boundary->BCFunctionU = BCSWSlipWallU;    
    else if (!strcmp(boundary->bctype,_SUBSONIC_OUTFLOW_            )) boundary->BCFunctionU = BCSubsonicOutflowU;    
    else if (!strcmp(boundary->bctype,_SUBSONIC_INFLOW_             )) boundary->BCFunctionU = BCSubsonicInflowU;    
    else if (!strcmp(boundary->bctype,_SUBSONIC_AMBIVALENT_         )) boundary->BCFunctionU = BCSubsonicAmbivalentU;    
    else if (!strcmp(boundary->bctype,_SUPERSONIC_OUTFLOW_          )) boundary->BCFunctionU = BCSupersonicOutflowU;    
    else if (!strcmp(boundary->bctype,_SUPERSONIC_INFLOW_           )) boundary->BCFunctionU = BCSupersonicInflowU;    
    else if (!strcmp(boundary->bctype,_TURBULENT_SUPERSONIC_INFLOW_ )) boundary->BCFunctionU = BCTurbulentSupersonicInflowU;    
    else if (!strcmp(boundary->bctype,_NO_FLUX_BC_                  )) boundary->BCFunctionU = BCNoFluxU;    
    else {
      fprintf(stderr,"Error in BCInitialize(): \"%s\" is not a supported boundary condition.\n",
              boundary->bctype);
      return(1);
    }

#if defined(HAVE_CUDA)
  }
#endif

  return 0;
}

BCCleanup()

int BCCleanup	(	void *	b,
		int	flag_gpu
	)

Function to clean up boundary conditions-related variables and arrays

Cleans up a boundary object of type DomainBoundary

Parameters

b	Boundary object of type DomainBoundary
flag_gpu	Flag indicating if GPU is being used

Definition at line 12 of file BCCleanup.c.

{
  DomainBoundary *boundary = (DomainBoundary*) b;
  free(boundary->xmin);
  free(boundary->xmax);
  free(boundary->is);
  free(boundary->ie);
  if (boundary->DirichletValue) free(boundary->DirichletValue);
  if (boundary->SpongeValue   ) free(boundary->SpongeValue   );
  if (boundary->FlowVelocity  ) free(boundary->FlowVelocity  );

  if (boundary->UnsteadyDirichletSize) free(boundary->UnsteadyDirichletSize);
  if (boundary->UnsteadyDirichletData) free(boundary->UnsteadyDirichletData);

  if (boundary->UnsteadyTemperatureSize)  free(boundary->UnsteadyTemperatureSize);
  if (boundary->UnsteadyTimeLevels)       free(boundary->UnsteadyTimeLevels);
  if (boundary->UnsteadyTemperatureData)  free(boundary->UnsteadyTemperatureData);

#if defined(HAVE_CUDA)
  if (flag_gpu) {
    gpuFree(boundary->gpu_bounds);
    gpuFree(boundary->gpu_is);
    gpuFree(boundary->gpu_ie);
    if (boundary->FlowVelocity) gpuFree(boundary->gpu_FlowVelocity  );
  }
#endif

  return 0;
}

BCPeriodicU()

int BCPeriodicU	(	void *	b,
		void *	m,
		int	ndims,
		int	nvars,
		int *	size,
		int	ghosts,
		double *	phi,
		double	waqt
	)

Periodic boundary conditions for the solution vector U

Applies periodic boundary conditions: Implemented by copying the solution from the other end of the domain into the physical boundary ghost points.

Note**: This function only acts if the the number of processors is 1 along the spatial dimension this boundary corresponds to. If there are more than 1 processors along this dimension, periodicity is handled by MPIExchangeBoundariesnD() to minimize communication.

Parameters

b	Boundary object of type DomainBoundary
m	MPI object of type MPIVariables
ndims	Number of spatial dimensions
nvars	Number of variables/DoFs per grid point
size	Integer array with the number of grid points in each spatial dimension
ghosts	Number of ghost points
phi	The solution array on which to apply the boundary condition
waqt	Current solution time

Definition at line 19 of file BCPeriodic.c.

{
  DomainBoundary *boundary = (DomainBoundary*) b;
  MPIVariables   *mpi      = (MPIVariables*)   m;

  int dim   = boundary->dim;
  int face  = boundary->face;

  if ((boundary->on_this_proc) && (mpi->iproc[dim] == 1)) {
    int bounds[ndims], index1[ndims], index2[ndims];
    _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
    _ArraySetValue_(index1,ndims,0);
    _ArraySetValue_(index2,ndims,0);
    int done = 0;
    while (!done) {
      int p1 = 0, p2 = 0;
      _ArrayCopy1D_(index1,index2,ndims);
      if (face == 1) {
        index2[dim] = index1[dim] + size[dim]-ghosts;
        _ArrayIndex1DWO_(ndims,size,index1,boundary->is,ghosts,p1);
        _ArrayIndex1D_(ndims,size,index2,ghosts,p2);
      } else if (face == -1) {
        _ArrayIndex1DWO_(ndims,size,index1,boundary->is,ghosts,p1);
        _ArrayIndex1D_(ndims,size,index1,ghosts,p2);
      }
      _ArrayCopy1D_((phi+nvars*p2),(phi+nvars*p1),nvars);
      _ArrayIncrementIndex_(ndims,bounds,index1,done);
    }
  }

  return(0);
}

BCExtrapolateU()

int BCExtrapolateU	(	void *	b,
		void *	m,
		int	ndims,
		int	nvars,
		int *	size,
		int	ghosts,
		double *	phi,
		double	waqt
	)

extrapolate boundary conditions for the solution vector U

Apply the extrapolative boundary condition: Values at the physical boundary ghost points are extrapolated from the interior points adjacent to the boundary

Parameters

b	Boundary object of type DomainBoundary
m	MPI object of type MPIVariables
ndims	Number of spatial dimensions
nvars	Number of variables/DoFs per grid point
size	Integer array with the number of grid points in each spatial dimension
ghosts	Number of ghost points
phi	The solution array on which to apply the boundary condition
waqt	Current solution time

Definition at line 13 of file BCExtrapolate.c.

{
  DomainBoundary *boundary = (DomainBoundary*) b;

  int dim   = boundary->dim;
  int face  = boundary->face;

  if (boundary->on_this_proc) {
    int bounds[ndims], indexb[ndims], indexi[ndims];
    _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
    _ArraySetValue_(indexb,ndims,0);
    int done = 0;
    while (!done) {
      _ArrayCopy1D_(indexb,indexi,ndims);
      _ArrayAdd1D_(indexi,indexi,boundary->is,ndims);
      if (face == 1)        indexi[dim] = ghosts-1-indexb[dim];
      else if (face == -1)  indexi[dim] = size[dim]-indexb[dim]-1;
      else                  return(1);
      int p1,p2;
      _ArrayIndex1DWO_(ndims,size,indexb,boundary->is,ghosts,p1);
      _ArrayIndex1D_(ndims,size,indexi,ghosts,p2);
      _ArrayCopy1D_((phi+nvars*p2),(phi+nvars*p1),nvars);
      _ArrayIncrementIndex_(ndims,bounds,indexb,done);
    }
  }
  return(0);
}

BCDirichletU()

int BCDirichletU	(	void *	b,
		void *	m,
		int	ndims,
		int	nvars,
		int *	size,
		int	ghosts,
		double *	phi,
		double	waqt
	)

Dirichlet boundary conditions for the solution vector U

Applies (steady) Dirichlet boundary conditions for the solution: the ghost points at the physical boundaries are set to specified values

Parameters

b	Boundary object of type DomainBoundary
m	MPI object of type MPIVariables
ndims	Number of spatial dimensions
nvars	Number of variables/DoFs per grid point
size	Integer array with the number of grid points in each spatial dimension
ghosts	Number of ghost points
phi	The solution array on which to apply the boundary condition
waqt	Current solution time

Definition at line 13 of file BCDirichlet.c.

{
  DomainBoundary *boundary = (DomainBoundary*) b;

  if (boundary->on_this_proc) {
    int bounds[ndims], indexb[ndims];
    _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
    _ArraySetValue_(indexb,ndims,0); 
    int done = 0;
    while (!done) {
      int p; _ArrayIndex1DWO_(ndims,size  ,indexb,boundary->is,ghosts,p);
      _ArrayCopy1D_((boundary->DirichletValue),(phi+nvars*p),nvars);
      _ArrayIncrementIndex_(ndims,bounds,indexb,done);
    }
  }
  return(0);
}

BCReflectU()

int BCReflectU	(	void *	b,
		void *	m,
		int	ndims,
		int	nvars,
		int *	size,
		int	ghosts,
		double *	phi,
		double	waqt
	)

Reflection boundary conditions for the solution vector U

Applies the reflection boundary condition: The values at the physical boundary ghost points are set to the negative of the interior values adjacent to the boundary.

Parameters

b	Boundary object of type DomainBoundary
m	MPI object of type MPIVariables
ndims	Number of spatial dimensions
nvars	Number of variables/DoFs per grid point
size	Integer array with the number of grid points in each spatial dimension
ghosts	Number of ghost points
phi	The solution array on which to apply the boundary condition
waqt	Current solution time

Definition at line 14 of file BCReflect.c.

{
  DomainBoundary *boundary = (DomainBoundary*) b;

  int dim   = boundary->dim;
  int face  = boundary->face;

  if (boundary->on_this_proc) {
    int bounds[ndims], indexb[ndims], indexi[ndims];
    _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
    _ArraySetValue_(indexb,ndims,0);
    int done = 0;
    while (!done) {
      int p1, p2;
      _ArrayCopy1D_(indexb,indexi,ndims);
      _ArrayAdd1D_(indexi,indexi,boundary->is,ndims);
      if      (face ==  1) indexi[dim] = ghosts-1-indexb[dim];
      else if (face == -1) indexi[dim] = size[dim]-indexb[dim]-1;
      else return(1);
      _ArrayIndex1DWO_(ndims,size,indexb,boundary->is,ghosts,p1);
      _ArrayIndex1D_(ndims,size,indexi,ghosts,p2);
      _ArrayScaleCopy1D_((phi+nvars*p2),(-1.0),(phi+nvars*p1),nvars);
      _ArrayIncrementIndex_(ndims,bounds,indexb,done);
    }
  }
  return(0);
}

BCNoslipWallU()

int BCNoslipWallU	(	void *	b,
		void *	m,
		int	ndims,
		int	nvars,
		int *	size,
		int	ghosts,
		double *	phi,
		double	waqt
	)

No-slip wall (viscous) boundary conditions for the solution vector U

Applies the no-slip wall boundary conditions: Used to simulate viscous walls. The density and pressure at the physical boundary ghost points are extrapolated from the interior, while the velocities are set such that the interpolated velocity at the boundary face is the specified wall velocity. This boundary condition is specific to the two and three dimensional Navier-Stokes systems (NavierStokes2D, NavierStokes3D).

Parameters

b	Boundary object of type DomainBoundary
m	MPI object of type MPIVariables
ndims	Number of spatial dimensions
nvars	Number of variables/DoFs per grid point
size	Integer array with the number of grid points in each spatial dimension
ghosts	Number of ghost points
phi	The solution array on which to apply the boundary condition
waqt	Current solution time

Definition at line 20 of file BCNoslipWall.c.

{
  DomainBoundary *boundary = (DomainBoundary*) b;

  int dim   = boundary->dim;
  int face  = boundary->face;

  if (ndims == 2) {

    /* create a fake physics object */
    double gamma;
    gamma = boundary->gamma;
    double inv_gamma_m1 = 1.0/(gamma-1.0);

    if (boundary->on_this_proc) {
      int bounds[ndims], indexb[ndims], indexi[ndims];
      _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
      _ArraySetValue_(indexb,ndims,0);
      int done = 0;
      while (!done) {
        int p1, p2;
        _ArrayCopy1D_(indexb,indexi,ndims);
        _ArrayAdd1D_(indexi,indexi,boundary->is,ndims);
        if      (face ==  1) indexi[dim] = ghosts-1-indexb[dim];
        else if (face == -1) indexi[dim] = size[dim]-indexb[dim]-1;
        else return(1);
        _ArrayIndex1DWO_(ndims,size,indexb,boundary->is,ghosts,p1);
        _ArrayIndex1D_(ndims,size,indexi,ghosts,p2);
        
        /* flow variables in the interior */
        double rho, uvel, vvel, energy, pressure;
        double rho_gpt, uvel_gpt, vvel_gpt, energy_gpt, pressure_gpt;
        _NavierStokes2DGetFlowVar_((phi+nvars*p2),rho,uvel,vvel,energy,pressure,gamma);
        /* set the ghost point values */
        rho_gpt = rho;
        pressure_gpt = pressure;
        uvel_gpt = 2.0*boundary->FlowVelocity[0] - uvel;
        vvel_gpt = 2.0*boundary->FlowVelocity[1] - vvel;
        energy_gpt = inv_gamma_m1*pressure_gpt
                    + 0.5 * rho_gpt * (uvel_gpt*uvel_gpt + vvel_gpt*vvel_gpt);

        phi[nvars*p1+0] = rho_gpt;
        phi[nvars*p1+1] = rho_gpt * uvel_gpt;
        phi[nvars*p1+2] = rho_gpt * vvel_gpt;
        phi[nvars*p1+3] = energy_gpt;

        _ArrayIncrementIndex_(ndims,bounds,indexb,done);
      }
    }

  } else if (ndims == 3) {

    double gamma;
    gamma = boundary->gamma;
    double inv_gamma_m1 = 1.0/(gamma-1.0);

    if (boundary->on_this_proc) {
      int bounds[ndims], indexb[ndims], indexi[ndims];
      _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
      _ArraySetValue_(indexb,ndims,0);
      int done = 0;
      while (!done) {
        int p1, p2;
        _ArrayCopy1D_(indexb,indexi,ndims);
        _ArrayAdd1D_(indexi,indexi,boundary->is,ndims);
        if      (face ==  1) indexi[dim] = ghosts-1-indexb[dim];
        else if (face == -1) indexi[dim] = size[dim]-indexb[dim]-1;
        else return(1);
        _ArrayIndex1DWO_(ndims,size,indexb,boundary->is,ghosts,p1);
        _ArrayIndex1D_(ndims,size,indexi,ghosts,p2);
        
        /* flow variables in the interior */
        double rho, uvel, vvel, wvel, energy, pressure;
        double rho_gpt, uvel_gpt, vvel_gpt, wvel_gpt, energy_gpt, pressure_gpt;
        _NavierStokes3DGetFlowVar_((phi+nvars*p2),_NavierStokes3D_stride_,rho,uvel,vvel,wvel,energy,pressure,gamma);
        /* set the ghost point values */
        rho_gpt = rho;
        pressure_gpt = pressure;
        uvel_gpt = 2.0*boundary->FlowVelocity[0] - uvel;
        vvel_gpt = 2.0*boundary->FlowVelocity[1] - vvel;
        wvel_gpt = 2.0*boundary->FlowVelocity[2] - wvel;
        energy_gpt = inv_gamma_m1*pressure_gpt
                    + 0.5 * rho_gpt 
                    * (uvel_gpt*uvel_gpt + vvel_gpt*vvel_gpt + wvel_gpt*wvel_gpt);

        phi[nvars*p1+0] = rho_gpt;
        phi[nvars*p1+1] = rho_gpt * uvel_gpt;
        phi[nvars*p1+2] = rho_gpt * vvel_gpt;
        phi[nvars*p1+3] = rho_gpt * wvel_gpt;
        phi[nvars*p1+4] = energy_gpt;

        _ArrayIncrementIndex_(ndims,bounds,indexb,done);
      }
    }

  }
  return(0);
}

BCThermalNoslipWallU()

int BCThermalNoslipWallU	(	void *	b,
		void *	m,
		int	ndims,
		int	nvars,
		int *	size,
		int	ghosts,
		double *	phi,
		double	waqt
	)

No-slip thermal wall (viscous) boundary conditions for the solution vector U

Applies the thermal no-slip-wall boundary condition: This is specific to the 3D Navier-Stokes system (NavierStokes3D). It is used for simulating walls boundaries, where the temperature is specified. The density at the ghost points is extrapolated from the interior. The velocity at the ghost points is set such that the interpolated velocity at the boundary face is the specified wall velocity. The pressure at the ghost points is set by multiplying the extrapolated density by the specified temperature.

Note: It is assumed that the temperature already contains the gas constant factor, i.e., \( T = P/\rho\).

Parameters

b	Boundary object of type DomainBoundary
m	MPI object of type MPIVariables
ndims	Number of spatial dimensions
nvars	Number of variables/DoFs per grid point
size	Integer array with the number of grid points in each spatial dimension
ghosts	Number of ghost points
phi	The solution array on which to apply the boundary condition
waqt	Current solution time

Definition at line 27 of file BCThermalNoslipWall.c.

{
  DomainBoundary *boundary = (DomainBoundary*) b;

  int dim   = boundary->dim;
  int face  = boundary->face;

  if (ndims == 3) {

    /* create a fake physics object */
    double gamma; 
    gamma = boundary->gamma;
    double inv_gamma_m1 = 1.0/(gamma-1.0);

    if (boundary->on_this_proc) {

      int     *temperature_field_size = boundary->UnsteadyTemperatureSize;
      int     n_time_levels           = temperature_field_size[dim];
      double  *time_levels            = boundary->UnsteadyTimeLevels;
      double  *temperature_data       = boundary->UnsteadyTemperatureData;
    
      int it = n_time_levels - 1;
      while ((time_levels[it] > waqt) && (it > 0))  it--;

      int bounds[ndims], indexb[ndims], indexi[ndims];
      _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
      _ArraySetValue_(indexb,ndims,0);

      int done = 0;
      while (!done) {

        int p1, p2;
        _ArrayCopy1D_(indexb,indexi,ndims);
        _ArrayAdd1D_(indexi,indexi,boundary->is,ndims);
        if      (face ==  1) indexi[dim] = ghosts-1-indexb[dim];
        else if (face == -1) indexi[dim] = size[dim]-indexb[dim]-1;
        else return(1);
        _ArrayIndex1DWO_(ndims,size,indexb,boundary->is,ghosts,p1);
        _ArrayIndex1D_(ndims,size,indexi,ghosts,p2);

        /* get the specified temperature */
        int index1[ndims]; _ArrayCopy1D_(indexb,index1,ndims);
        index1[dim] = it;
        int q; _ArrayIndex1D_(ndims,temperature_field_size,index1,0,q);
        double temperature_b = temperature_data[q];
        
        /* flow variables in the interior */
        double rho, uvel, vvel, wvel, energy, pressure;
        _NavierStokes3DGetFlowVar_((phi+nvars*p2),_NavierStokes3D_stride_,rho,uvel,vvel,wvel,energy,pressure,gamma);
        /* set the ghost point values */
        double rho_gpt, uvel_gpt, vvel_gpt, wvel_gpt, energy_gpt, pressure_gpt;
        rho_gpt = rho;
        uvel_gpt = 2.0*boundary->FlowVelocity[0] - uvel;
        vvel_gpt = 2.0*boundary->FlowVelocity[1] - vvel;
        wvel_gpt = 2.0*boundary->FlowVelocity[2] - wvel;
        pressure_gpt = rho_gpt * temperature_b;
        energy_gpt = inv_gamma_m1*pressure_gpt 
                    + 0.5 * rho_gpt 
                    * (uvel_gpt*uvel_gpt + vvel_gpt*vvel_gpt + wvel_gpt*wvel_gpt);

        phi[nvars*p1+0] = rho_gpt;
        phi[nvars*p1+1] = rho_gpt * uvel_gpt;
        phi[nvars*p1+2] = rho_gpt * vvel_gpt;
        phi[nvars*p1+3] = rho_gpt * wvel_gpt;
        phi[nvars*p1+4] = energy_gpt;

        _ArrayIncrementIndex_(ndims,bounds,indexb,done);
      }
    }

  }
  return(0);
}

BCSlipWallU()

int BCSlipWallU	(	void *	b,
		void *	m,
		int	ndims,
		int	nvars,
		int *	size,
		int	ghosts,
		double *	phi,
		double	waqt
	)

Slip (inviscid) wall boundary conditions for the solution vector U

Applies the slip-wall boundary condition: This is specific to the two and three dimensional Euler and Navier-Stokes systems (Euler2D, NavierStokes2D, NavierStokes3D). It is used for simulating inviscid walls or symmetric boundaries. The pressure, density, and tangential velocity at the ghost points are extrapolated from the interior, while the normal velocity at the ghost points is set such that the interpolated value at the boundary face is equal to the specified wall velocity.

Parameters

b	Boundary object of type DomainBoundary
m	MPI object of type MPIVariables
ndims	Number of spatial dimensions
nvars	Number of variables/DoFs per grid point
size	Integer array with the number of grid points in each spatial dimension
ghosts	Number of ghost points
phi	The solution array on which to apply the boundary condition
waqt	Current solution time

Definition at line 22 of file BCSlipWall.c.

{
  DomainBoundary *boundary = (DomainBoundary*) b;

  int dim   = boundary->dim;
  int face  = boundary->face;

  if (ndims == 1) {

    /* create a fake physics object */
    Euler1D physics; 
    double gamma; 
    gamma = physics.gamma = boundary->gamma;
    double inv_gamma_m1 = 1.0/(gamma-1.0);

    if (boundary->on_this_proc) {
      int bounds[ndims], indexb[ndims], indexi[ndims];
      _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
      _ArraySetValue_(indexb,ndims,0);
      int done = 0;
      while (!done) {
        int p1, p2;
        _ArrayCopy1D_(indexb,indexi,ndims);
        _ArrayAdd1D_(indexi,indexi,boundary->is,ndims);
        if      (face ==  1) indexi[dim] = ghosts-1-indexb[dim];
        else if (face == -1) indexi[dim] = size[dim]-indexb[dim]-1;
        else return(1);
        _ArrayIndex1DWO_(ndims,size,indexb,boundary->is,ghosts,p1);
        _ArrayIndex1D_(ndims,size,indexi,ghosts,p2);
        
        /* flow variables in the interior */
        double rho, uvel, energy, pressure;
        double rho_gpt, uvel_gpt, energy_gpt, pressure_gpt;
        _Euler1DGetFlowVar_((phi+nvars*p2),rho,uvel,energy,pressure,(&physics));
        /* set the ghost point values */
        rho_gpt = rho;
        pressure_gpt = pressure;
        uvel_gpt = 2.0*boundary->FlowVelocity[_XDIR_] - uvel;
        energy_gpt = inv_gamma_m1*pressure_gpt 
                    + 0.5 * rho_gpt * uvel_gpt*uvel_gpt;

        phi[nvars*p1+0] = rho_gpt;
        phi[nvars*p1+1] = rho_gpt * uvel_gpt;
        phi[nvars*p1+2] = energy_gpt;

        _ArrayIncrementIndex_(ndims,bounds,indexb,done);
      }
    }

  } else if (ndims == 2) {

    /* create a fake physics object */
    Euler2D physics; 
    double gamma; 
    gamma = physics.gamma = boundary->gamma;
    double inv_gamma_m1 = 1.0/(gamma-1.0);

    if (boundary->on_this_proc) {
      int bounds[ndims], indexb[ndims], indexi[ndims];
      _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
      _ArraySetValue_(indexb,ndims,0);
      int done = 0;
      while (!done) {
        int p1, p2;
        _ArrayCopy1D_(indexb,indexi,ndims);
        _ArrayAdd1D_(indexi,indexi,boundary->is,ndims);
        if      (face ==  1) indexi[dim] = ghosts-1-indexb[dim];
        else if (face == -1) indexi[dim] = size[dim]-indexb[dim]-1;
        else return(1);
        _ArrayIndex1DWO_(ndims,size,indexb,boundary->is,ghosts,p1);
        _ArrayIndex1D_(ndims,size,indexi,ghosts,p2);
        
        /* flow variables in the interior */
        double rho, uvel, vvel, energy, pressure;
        double rho_gpt, uvel_gpt, vvel_gpt, energy_gpt, pressure_gpt;
        _Euler2DGetFlowVar_((phi+nvars*p2),rho,uvel,vvel,energy,pressure,(&physics));
        /* set the ghost point values */
        rho_gpt = rho;
        pressure_gpt = pressure;
        if (dim == _XDIR_) {
          uvel_gpt = 2.0*boundary->FlowVelocity[_XDIR_] - uvel;
          vvel_gpt = vvel;
        } else if (dim == _YDIR_) {
          uvel_gpt = uvel;
          vvel_gpt = 2.0*boundary->FlowVelocity[_YDIR_] - vvel;
        } else {
          uvel_gpt = 0.0;
          vvel_gpt = 0.0;
        }
        energy_gpt = inv_gamma_m1*pressure_gpt 
                    + 0.5 * rho_gpt * (uvel_gpt*uvel_gpt + vvel_gpt*vvel_gpt);

        phi[nvars*p1+0] = rho_gpt;
        phi[nvars*p1+1] = rho_gpt * uvel_gpt;
        phi[nvars*p1+2] = rho_gpt * vvel_gpt;
        phi[nvars*p1+3] = energy_gpt;

        _ArrayIncrementIndex_(ndims,bounds,indexb,done);
      }
    }

  } else if (ndims == 3) {

    /* create a fake physics object */
    double gamma; 
    gamma = boundary->gamma;
    double inv_gamma_m1 = 1.0/(gamma-1.0);

    if (boundary->on_this_proc) {
      int bounds[ndims], indexb[ndims], indexi[ndims];
      _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
      _ArraySetValue_(indexb,ndims,0);
      int done = 0;
      while (!done) {
        int p1, p2;
        _ArrayCopy1D_(indexb,indexi,ndims);
        _ArrayAdd1D_(indexi,indexi,boundary->is,ndims);
        if      (face ==  1) indexi[dim] = ghosts-1-indexb[dim];
        else if (face == -1) indexi[dim] = size[dim]-indexb[dim]-1;
        else return(1);
        _ArrayIndex1DWO_(ndims,size,indexb,boundary->is,ghosts,p1);
        _ArrayIndex1D_(ndims,size,indexi,ghosts,p2);
        
        /* flow variables in the interior */
        double rho, uvel, vvel, wvel, energy, pressure;
        double rho_gpt, uvel_gpt, vvel_gpt, wvel_gpt, energy_gpt, pressure_gpt;
        _NavierStokes3DGetFlowVar_((phi+nvars*p2),_NavierStokes3D_stride_,rho,uvel,vvel,wvel,energy,pressure,gamma);
        /* set the ghost point values */
        rho_gpt = rho;
        pressure_gpt = pressure;
        if (dim == _XDIR_) {
          uvel_gpt = 2.0*boundary->FlowVelocity[_XDIR_] - uvel;
          vvel_gpt = vvel;
          wvel_gpt = wvel;
        } else if (dim == _YDIR_) {
          uvel_gpt = uvel;
          vvel_gpt = 2.0*boundary->FlowVelocity[_YDIR_] - vvel;
          wvel_gpt = wvel;
        } else if (dim == _ZDIR_) {
          uvel_gpt = uvel;
          vvel_gpt = vvel;
          wvel_gpt = 2.0*boundary->FlowVelocity[_ZDIR_] - wvel;
        } else {
          uvel_gpt = 0.0;
          vvel_gpt = 0.0;
          wvel_gpt = 0.0;
        }
        energy_gpt = inv_gamma_m1*pressure_gpt 
                    + 0.5 * rho_gpt 
                    * (uvel_gpt*uvel_gpt + vvel_gpt*vvel_gpt + wvel_gpt*wvel_gpt);

        phi[nvars*p1+0] = rho_gpt;
        phi[nvars*p1+1] = rho_gpt * uvel_gpt;
        phi[nvars*p1+2] = rho_gpt * vvel_gpt;
        phi[nvars*p1+3] = rho_gpt * wvel_gpt;
        phi[nvars*p1+4] = energy_gpt;

        _ArrayIncrementIndex_(ndims,bounds,indexb,done);
      }
    }

  }
  return(0);
}

BCThermalSlipWallU()

int BCThermalSlipWallU	(	void *	b,
		void *	m,
		int	ndims,
		int	nvars,
		int *	size,
		int	ghosts,
		double *	phi,
		double	waqt
	)

Slip (inviscid) thermal wall boundary conditions for the solution vector U

Applies the thermal slip-wall boundary condition: This is specific to the 3D Navier-Stokes system (NavierStokes3D). It is used for simulating inviscid walls or symmetric boundaries, where the temperature is specified. The density and the tangential velocity at the ghost points are extrapolated for the interior. The normal velocity at the ghost points is set such that the interpolated velocity at the boundary face is the specified wall velocity. The pressure at the ghost points is set by multiplying the extrapolated density by the specified temperature.

Note: It is assumed that the temperature already contains the gas constant factor, i.e., \( T = P/\rho\).

Parameters

b	Boundary object of type DomainBoundary
m	MPI object of type MPIVariables
ndims	Number of spatial dimensions
nvars	Number of variables/DoFs per grid point
size	Integer array with the number of grid points in each spatial dimension
ghosts	Number of ghost points
phi	The solution array on which to apply the boundary condition
waqt	Current solution time

Definition at line 27 of file BCThermalSlipWall.c.

{
  DomainBoundary *boundary = (DomainBoundary*) b;

  int dim   = boundary->dim;
  int face  = boundary->face;

  if (ndims == 3) {

    /* create a fake physics object */
    double gamma; 
    gamma = boundary->gamma;
    double inv_gamma_m1 = 1.0/(gamma-1.0);

    if (boundary->on_this_proc) {

      int     *temperature_field_size = boundary->UnsteadyTemperatureSize;
      int     n_time_levels           = temperature_field_size[dim];
      double  *time_levels            = boundary->UnsteadyTimeLevels;
      double  *temperature_data       = boundary->UnsteadyTemperatureData;
    
      int it = n_time_levels - 1;
      while ((time_levels[it] > waqt) && (it > 0))  it--;

      int bounds[ndims], indexb[ndims], indexi[ndims];
      _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
      _ArraySetValue_(indexb,ndims,0);

      int done = 0;
      while (!done) {

        int p1, p2;
        _ArrayCopy1D_(indexb,indexi,ndims);
        _ArrayAdd1D_(indexi,indexi,boundary->is,ndims);
        if      (face ==  1) indexi[dim] = ghosts-1-indexb[dim];
        else if (face == -1) indexi[dim] = size[dim]-indexb[dim]-1;
        else return(1);
        _ArrayIndex1DWO_(ndims,size,indexb,boundary->is,ghosts,p1);
        _ArrayIndex1D_(ndims,size,indexi,ghosts,p2);

        /* get the specified temperature */
        int index1[ndims]; _ArrayCopy1D_(indexb,index1,ndims);
        index1[dim] = it;
        int q; _ArrayIndex1D_(ndims,temperature_field_size,index1,0,q);
        double temperature_b = temperature_data[q];
        
        /* flow variables in the interior */
        double rho, uvel, vvel, wvel, energy, pressure;
        _NavierStokes3DGetFlowVar_((phi+nvars*p2),_NavierStokes3D_stride_,rho,uvel,vvel,wvel,energy,pressure,gamma);
        /* set the ghost point values */
        double rho_gpt, uvel_gpt, vvel_gpt, wvel_gpt, energy_gpt, pressure_gpt;
        rho_gpt = rho;
        if (dim == _XDIR_) {
          uvel_gpt = 2.0*boundary->FlowVelocity[_XDIR_] - uvel;
          vvel_gpt = vvel;
          wvel_gpt = wvel;
        } else if (dim == _YDIR_) {
          uvel_gpt = uvel;
          vvel_gpt = 2.0*boundary->FlowVelocity[_YDIR_] - vvel;
          wvel_gpt = wvel;
        } else if (dim == _ZDIR_) {
          uvel_gpt = uvel;
          vvel_gpt = vvel;
          wvel_gpt = 2.0*boundary->FlowVelocity[_ZDIR_] - wvel;
        } else {
          uvel_gpt = 0.0;
          vvel_gpt = 0.0;
          wvel_gpt = 0.0;
        }
        pressure_gpt = rho_gpt * temperature_b;
        energy_gpt = inv_gamma_m1*pressure_gpt 
                    + 0.5 * rho_gpt 
                    * (uvel_gpt*uvel_gpt + vvel_gpt*vvel_gpt + wvel_gpt*wvel_gpt);

        phi[nvars*p1+0] = rho_gpt;
        phi[nvars*p1+1] = rho_gpt * uvel_gpt;
        phi[nvars*p1+2] = rho_gpt * vvel_gpt;
        phi[nvars*p1+3] = rho_gpt * wvel_gpt;
        phi[nvars*p1+4] = energy_gpt;

        _ArrayIncrementIndex_(ndims,bounds,indexb,done);
      }
    }

  }
  return(0);
}

BCSubsonicInflowU()

int BCSubsonicInflowU	(	void *	b,
		void *	m,
		int	ndims,
		int	nvars,
		int *	size,
		int	ghosts,
		double *	phi,
		double	waqt
	)

Subsonic inflow boundary conditions for the solution vector U

Applies the subsonic inflow boundary condition: The density and velocity at the physical boundary ghost points are specified, while the pressure is extrapolated from the interior of the domain. This boundary condition is specific to two and three dimension Euler and Navier-Stokes systems (Euler2D, NavierStokes2D, NavierStokes3D).

Parameters

b	Boundary object of type DomainBoundary
m	MPI object of type MPIVariables
ndims	Number of spatial dimensions
nvars	Number of variables/DoFs per grid point
size	Integer array with the number of grid points in each spatial dimension
ghosts	Number of ghost points
phi	The solution array on which to apply the boundary condition
waqt	Current solution time

Definition at line 20 of file BCSubsonicInflow.c.

{
  DomainBoundary *boundary = (DomainBoundary*) b;

  int dim   = boundary->dim;
  int face  = boundary->face;

  if (ndims == 2) {

    /* create a fake physics object */
    Euler2D physics; 
    double gamma;
    gamma = physics.gamma = boundary->gamma;
    double inv_gamma_m1 = 1.0/(gamma-1.0);

    if (boundary->on_this_proc) {
      int bounds[ndims], indexb[ndims], indexi[ndims];
      _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
      _ArraySetValue_(indexb,ndims,0);
      int done = 0;
      while (!done) {
        int p1, p2;
        _ArrayCopy1D_(indexb,indexi,ndims);
        _ArrayAdd1D_(indexi,indexi,boundary->is,ndims);
        if      (face ==  1) indexi[dim] = ghosts-1-indexb[dim];
        else if (face == -1) indexi[dim] = size[dim]-indexb[dim]-1;
        else return(1);
        _ArrayIndex1DWO_(ndims,size,indexb,boundary->is,ghosts,p1);
        _ArrayIndex1D_(ndims,size,indexi,ghosts,p2);
        
        /* flow variables in the interior */
        double rho, uvel, vvel, energy, pressure;
        double rho_gpt, uvel_gpt, vvel_gpt, energy_gpt, pressure_gpt;
        _Euler2DGetFlowVar_((phi+nvars*p2),rho,uvel,vvel,energy,pressure,(&physics));
        /* set the ghost point values */
        rho_gpt = boundary->FlowDensity;
        pressure_gpt = pressure;
        uvel_gpt = boundary->FlowVelocity[0];
        vvel_gpt = boundary->FlowVelocity[1];
        energy_gpt = inv_gamma_m1*pressure_gpt
                    + 0.5 * rho_gpt * (uvel_gpt*uvel_gpt + vvel_gpt*vvel_gpt);

        phi[nvars*p1+0] = rho_gpt;
        phi[nvars*p1+1] = rho_gpt * uvel_gpt;
        phi[nvars*p1+2] = rho_gpt * vvel_gpt;
        phi[nvars*p1+3] = energy_gpt;

        _ArrayIncrementIndex_(ndims,bounds,indexb,done);
      }
    }

  } else if (ndims == 3) {

    /* create a fake physics object */
    double gamma;
    gamma = boundary->gamma;
    double inv_gamma_m1 = 1.0/(gamma-1.0);

    if (boundary->on_this_proc) {
      int bounds[ndims], indexb[ndims], indexi[ndims];
      _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
      _ArraySetValue_(indexb,ndims,0);
      int done = 0;
      while (!done) {
        int p1, p2;
        _ArrayCopy1D_(indexb,indexi,ndims);
        _ArrayAdd1D_(indexi,indexi,boundary->is,ndims);
        if      (face ==  1) indexi[dim] = ghosts-1-indexb[dim];
        else if (face == -1) indexi[dim] = size[dim]-indexb[dim]-1;
        else return(1);
        _ArrayIndex1DWO_(ndims,size,indexb,boundary->is,ghosts,p1);
        _ArrayIndex1D_(ndims,size,indexi,ghosts,p2);
        
        /* flow variables in the interior */
        double rho, uvel, vvel, wvel, energy, pressure;
        double rho_gpt, uvel_gpt, vvel_gpt, wvel_gpt, energy_gpt, pressure_gpt;
        _NavierStokes3DGetFlowVar_((phi+nvars*p2),_NavierStokes3D_stride_,rho,uvel,vvel,wvel,energy,pressure,gamma);
        /* set the ghost point values */
        rho_gpt = boundary->FlowDensity;
        pressure_gpt = pressure;
        uvel_gpt = boundary->FlowVelocity[0];
        vvel_gpt = boundary->FlowVelocity[1];
        wvel_gpt = boundary->FlowVelocity[2];
        energy_gpt = inv_gamma_m1*pressure_gpt
                    + 0.5 * rho_gpt 
                    * (uvel_gpt*uvel_gpt + vvel_gpt*vvel_gpt + wvel_gpt*wvel_gpt);

        phi[nvars*p1+0] = rho_gpt;
        phi[nvars*p1+1] = rho_gpt * uvel_gpt;
        phi[nvars*p1+2] = rho_gpt * vvel_gpt;
        phi[nvars*p1+3] = rho_gpt * wvel_gpt;
        phi[nvars*p1+4] = energy_gpt;

        _ArrayIncrementIndex_(ndims,bounds,indexb,done);
      }
    }

  }
  return(0);
}

BCSubsonicOutflowU()

int BCSubsonicOutflowU	(	void *	b,
		void *	m,
		int	ndims,
		int	nvars,
		int *	size,
		int	ghosts,
		double *	phi,
		double	waqt
	)

Subsonic outflow boundary conditions for the solution vector U

Applies the subsonic outflow boundary condition: The pressure at the physical boundary ghost points is specified, while the density and velocity are extrapolated from the interior. This boundary condition is specific to two and three dimensional Euler/ Navier-Stokes systems (Euler2D, NavierStokes2D, NavierStokes3D).

Parameters

b	Boundary object of type DomainBoundary
m	MPI object of type MPIVariables
ndims	Number of spatial dimensions
nvars	Number of variables/DoFs per grid point
size	Integer array with the number of grid points in each spatial dimension
ghosts	Number of ghost points
phi	The solution array on which to apply the boundary condition
waqt	Current solution time

Definition at line 19 of file BCSubsonicOutflow.c.

{
  DomainBoundary *boundary = (DomainBoundary*) b;

  int dim   = boundary->dim;
  int face  = boundary->face;

  if (ndims == 2) {

    /* create a fake physics object */
    Euler2D physics; 
    double gamma;
    gamma = physics.gamma = boundary->gamma;
    double inv_gamma_m1 = 1.0/(gamma-1.0);

    if (boundary->on_this_proc) {
      int bounds[ndims], indexb[ndims], indexi[ndims];
      _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
      _ArraySetValue_(indexb,ndims,0);
      int done = 0;
      while (!done) {
        int p1, p2;
        _ArrayCopy1D_(indexb,indexi,ndims);
        _ArrayAdd1D_(indexi,indexi,boundary->is,ndims);
        if      (face ==  1) indexi[dim] = ghosts-1-indexb[dim];
        else if (face == -1) indexi[dim] = size[dim]-indexb[dim]-1;
        else return(1);
        _ArrayIndex1DWO_(ndims,size,indexb,boundary->is,ghosts,p1);
        _ArrayIndex1D_(ndims,size,indexi,ghosts,p2);
        
        /* flow variables in the interior */
        double rho, uvel, vvel, energy, pressure;
        double rho_gpt, uvel_gpt, vvel_gpt, energy_gpt, pressure_gpt;
        _Euler2DGetFlowVar_((phi+nvars*p2),rho,uvel,vvel,energy,pressure,(&physics));
        /* set the ghost point values */
        rho_gpt = rho;
        pressure_gpt = pressure; /* useless statement to avoid compiler warning */
        pressure_gpt = boundary->FlowPressure;
        uvel_gpt = uvel;
        vvel_gpt = vvel;
        energy_gpt = inv_gamma_m1*pressure_gpt
                    + 0.5 * rho_gpt * (uvel_gpt*uvel_gpt + vvel_gpt*vvel_gpt);

        phi[nvars*p1+0] = rho_gpt;
        phi[nvars*p1+1] = rho_gpt * uvel_gpt;
        phi[nvars*p1+2] = rho_gpt * vvel_gpt;
        phi[nvars*p1+3] = energy_gpt;

        _ArrayIncrementIndex_(ndims,bounds,indexb,done);
      }
    }

  } else if (ndims == 3) {

    /* create a fake physics object */
    double gamma;
    gamma = boundary->gamma;
    double inv_gamma_m1 = 1.0/(gamma-1.0);

    if (boundary->on_this_proc) {
      int bounds[ndims], indexb[ndims], indexi[ndims];
      _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
      _ArraySetValue_(indexb,ndims,0);
      int done = 0;
      while (!done) {
        int p1, p2;
        _ArrayCopy1D_(indexb,indexi,ndims);
        _ArrayAdd1D_(indexi,indexi,boundary->is,ndims);
        if      (face ==  1) indexi[dim] = ghosts-1-indexb[dim];
        else if (face == -1) indexi[dim] = size[dim]-indexb[dim]-1;
        else return(1);
        _ArrayIndex1DWO_(ndims,size,indexb,boundary->is,ghosts,p1);
        _ArrayIndex1D_(ndims,size,indexi,ghosts,p2);
        
        /* flow variables in the interior */
        double rho, uvel, vvel, wvel, energy, pressure;
        double rho_gpt, uvel_gpt, vvel_gpt, wvel_gpt, energy_gpt, pressure_gpt;
        _NavierStokes3DGetFlowVar_((phi+nvars*p2),_NavierStokes3D_stride_,rho,uvel,vvel,wvel,energy,pressure,gamma);
        /* set the ghost point values */
        rho_gpt = rho;
        pressure_gpt = pressure; /* useless statement to avoid compiler warning */
        pressure_gpt = boundary->FlowPressure;
        uvel_gpt = uvel;
        vvel_gpt = vvel;
        wvel_gpt = wvel;
        energy_gpt = inv_gamma_m1*pressure_gpt
                    + 0.5 * rho_gpt 
                    * (uvel_gpt*uvel_gpt + vvel_gpt*vvel_gpt + wvel_gpt*wvel_gpt);

        phi[nvars*p1+0] = rho_gpt;
        phi[nvars*p1+1] = rho_gpt * uvel_gpt;
        phi[nvars*p1+2] = rho_gpt * vvel_gpt;
        phi[nvars*p1+3] = rho_gpt * wvel_gpt;
        phi[nvars*p1+4] = energy_gpt;

        _ArrayIncrementIndex_(ndims,bounds,indexb,done);
      }
    }

  }
  return(0);
}

BCSubsonicAmbivalentU()

int BCSubsonicAmbivalentU	(	void *	b,
		void *	m,
		int	ndims,
		int	nvars,
		int *	size,
		int	ghosts,
		double *	phi,
		double	waqt
	)

Subsonic "ambivalent" boundary conditions for the solution vector U

Applies the subsonic "ambivalent" boundary condition: The velocity at the boundary face is extrapolated from the interior and its dot product with the boundary normal (pointing into the domain) is computed.

• If it is positive, subsonic inflow boundary conditions are applied: the density and velocity at the physical boundary ghost points are specified, while the pressure is extrapolated from the interior of the domain.
• If it is negative, subsonic outflow boundary conditions are applied: the pressure at the physical boundary ghost points is specified, while the density and velocity are extrapolated from the interior.

This boundary condition is specific to two and three dimension Euler and Navier-Stokes systems (Euler2D, NavierStokes2D, NavierStokes3D).

Parameters

b	Boundary object of type DomainBoundary
m	MPI object of type MPIVariables
ndims	Number of spatial dimensions
nvars	Number of variables/DoFs per grid point
size	Integer array with the number of grid points in each spatial dimension
ghosts	Number of ghost points
phi	The solution array on which to apply the boundary condition
waqt	Current solution time

Definition at line 28 of file BCSubsonicAmbivalent.c.

{
  DomainBoundary *boundary = (DomainBoundary*) b;

  int dim   = boundary->dim;
  int face  = boundary->face;

  if (ndims == 2) {

    /* create a fake physics object */
    Euler2D physics; 
    double gamma;
    gamma = physics.gamma = boundary->gamma;
    double inv_gamma_m1 = 1.0/(gamma-1.0);

    /* boundary normal (pointing into the domain) */
    double nx, ny;
    if (dim == 0) {
      nx = 1.0;
      ny = 0.0;
    } else if (dim == 1) {
      nx = 0.0;
      ny = 1.0;
    }
    nx *= (double) face;
    ny *= (double) face;

    if (boundary->on_this_proc) {
      int bounds[ndims], indexb[ndims], indexi[ndims], indexj[ndims];
      _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
      _ArraySetValue_(indexb,ndims,0);
      int done = 0;
      while (!done) {
        int p1, p2;
        double rho, uvel, vvel, energy, pressure;

        /* compute boundary face velocity  - 2nd order */
        _ArrayCopy1D_(indexb,indexi,ndims);
        _ArrayAdd1D_(indexi,indexi,boundary->is,ndims);
        _ArrayCopy1D_(indexi,indexj,ndims);
        if (face ==  1) {
          indexi[dim] = 0;
          indexj[dim] = indexi[dim] + 1;
        } else if (face == -1) {
          indexi[dim] = size[dim]-1;
          indexj[dim] = indexi[dim] - 1;
        }
        _ArrayIndex1D_(ndims,size,indexi,ghosts,p1);
        _ArrayIndex1D_(ndims,size,indexj,ghosts,p2);
        double uvel1, uvel2, uvelb,
               vvel1, vvel2, vvelb;
        _Euler2DGetFlowVar_((phi+nvars*p1),rho,uvel1,vvel1,energy,pressure,(&physics));
        _Euler2DGetFlowVar_((phi+nvars*p2),rho,uvel2,vvel2,energy,pressure,(&physics));
        uvelb = 1.5*uvel1 - 0.5*uvel2;
        vvelb = 1.5*vvel1 - 0.5*vvel2;
        double vel_normal = uvelb*nx + vvelb*ny;

        _ArrayCopy1D_(indexb,indexi,ndims);
        _ArrayAdd1D_(indexi,indexi,boundary->is,ndims);
        if      (face ==  1) indexi[dim] = ghosts-1-indexb[dim];
        else if (face == -1) indexi[dim] = size[dim]-indexb[dim]-1;
        else return(1);
        _ArrayIndex1DWO_(ndims,size,indexb,boundary->is,ghosts,p1);
        _ArrayIndex1D_(ndims,size,indexi,ghosts,p2);
        
        /* flow variables in the interior */
        _Euler2DGetFlowVar_((phi+nvars*p2),rho,uvel,vvel,energy,pressure,(&physics));

        /* set the ghost point values */
        double rho_gpt, uvel_gpt, vvel_gpt, energy_gpt, pressure_gpt;
        if (vel_normal > 0) {
          /* inflow */
          rho_gpt = boundary->FlowDensity;
          pressure_gpt = pressure;
          uvel_gpt = boundary->FlowVelocity[0];
          vvel_gpt = boundary->FlowVelocity[1];
        } else {
          /* outflow */
          rho_gpt = rho;
          pressure_gpt = boundary->FlowPressure;
          uvel_gpt = uvel;
          vvel_gpt = vvel;
        }
        energy_gpt = inv_gamma_m1*pressure_gpt
                    + 0.5 * rho_gpt * (uvel_gpt*uvel_gpt + vvel_gpt*vvel_gpt);

        phi[nvars*p1+0] = rho_gpt;
        phi[nvars*p1+1] = rho_gpt * uvel_gpt;
        phi[nvars*p1+2] = rho_gpt * vvel_gpt;
        phi[nvars*p1+3] = energy_gpt;

        _ArrayIncrementIndex_(ndims,bounds,indexb,done);
      }
    }

  } else if (ndims == 3) {

    /* create a fake physics object */
    double gamma;
    gamma = boundary->gamma;
    double inv_gamma_m1 = 1.0/(gamma-1.0);

    /* boundary normal (pointing into the domain) */
    double nx, ny, nz;
    if (dim == 0) {
      nx = 1.0;
      ny = 0.0;
      nz = 0.0;
    } else if (dim == 1) {
      nx = 0.0;
      ny = 1.0;
      nz = 0.0;
    } else if (dim == 2) {
      nx = 0.0;
      ny = 0.0;
      nz = 1.0;
    }
    nx *= (double) face;
    ny *= (double) face;
    nz *= (double) face;

    if (boundary->on_this_proc) {
      int bounds[ndims], indexb[ndims], indexi[ndims], indexj[ndims];
      _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
      _ArraySetValue_(indexb,ndims,0);
      int done = 0;
      while (!done) {
        int p1, p2;
        double rho, uvel, vvel, wvel, energy, pressure;

        /* compute boundary face velocity  - 2nd order */
        _ArrayCopy1D_(indexb,indexi,ndims);
        _ArrayAdd1D_(indexi,indexi,boundary->is,ndims);
        _ArrayCopy1D_(indexi,indexj,ndims);
        if (face ==  1) {
          indexi[dim] = 0;
          indexj[dim] = indexi[dim] + 1;
        } else if (face == -1) {
          indexi[dim] = size[dim]-1;
          indexj[dim] = indexi[dim] - 1;
        }
        _ArrayIndex1D_(ndims,size,indexi,ghosts,p1);
        _ArrayIndex1D_(ndims,size,indexj,ghosts,p2);
        double uvel1, uvel2, uvelb,
               vvel1, vvel2, vvelb,
               wvel1, wvel2, wvelb;
        _NavierStokes3DGetFlowVar_((phi+nvars*p1),_NavierStokes3D_stride_,rho,uvel1,vvel1,wvel1,energy,pressure,gamma);
        _NavierStokes3DGetFlowVar_((phi+nvars*p2),_NavierStokes3D_stride_,rho,uvel2,vvel2,wvel2,energy,pressure,gamma);
        uvelb = 1.5*uvel1 - 0.5*uvel2;
        vvelb = 1.5*vvel1 - 0.5*vvel2;
        wvelb = 1.5*wvel1 - 0.5*wvel2;
        double vel_normal = uvelb*nx + vvelb*ny + wvelb*nz;

        _ArrayCopy1D_(indexb,indexi,ndims);
        _ArrayAdd1D_(indexi,indexi,boundary->is,ndims);
        if      (face ==  1) indexi[dim] = ghosts-1-indexb[dim];
        else if (face == -1) indexi[dim] = size[dim]-indexb[dim]-1;
        else return(1);
        _ArrayIndex1DWO_(ndims,size,indexb,boundary->is,ghosts,p1);
        _ArrayIndex1D_(ndims,size,indexi,ghosts,p2);
        
        /* flow variables in the interior */
        _NavierStokes3DGetFlowVar_((phi+nvars*p2),_NavierStokes3D_stride_,rho,uvel,vvel,wvel,energy,pressure,gamma);

        /* set the ghost point values */
        double rho_gpt, uvel_gpt, vvel_gpt, wvel_gpt, energy_gpt, pressure_gpt;
        if (vel_normal > 0) {
          /* inflow */
          rho_gpt = boundary->FlowDensity;
          pressure_gpt = pressure;
          uvel_gpt = boundary->FlowVelocity[0];
          vvel_gpt = boundary->FlowVelocity[1];
          wvel_gpt = boundary->FlowVelocity[2];
        } else {
          /* outflow */
          rho_gpt = rho;
          pressure_gpt = boundary->FlowPressure;
          uvel_gpt = uvel;
          vvel_gpt = vvel;
          wvel_gpt = wvel;
        }
        energy_gpt = inv_gamma_m1*pressure_gpt
                    + 0.5 * rho_gpt 
                    * (uvel_gpt*uvel_gpt + vvel_gpt*vvel_gpt + wvel_gpt*wvel_gpt);

        phi[nvars*p1+0] = rho_gpt;
        phi[nvars*p1+1] = rho_gpt * uvel_gpt;
        phi[nvars*p1+2] = rho_gpt * vvel_gpt;
        phi[nvars*p1+3] = rho_gpt * wvel_gpt;
        phi[nvars*p1+4] = energy_gpt;

        _ArrayIncrementIndex_(ndims,bounds,indexb,done);
      }
    }

  }
  return(0);
}

BCSupersonicInflowU()

int BCSupersonicInflowU	(	void *	b,
		void *	m,
		int	ndims,
		int	nvars,
		int *	size,
		int	ghosts,
		double *	phi,
		double	waqt
	)

Supersonic inflow boundary conditions for the solution vector U

Applies the supersonic (steady) inflow boundary condition: All the flow variables (density, pressure, velocity) are specified at the physical boundary ghost points, since it is supersonic inflow. This boundary condition is specific to two and three dimensional Euler/Navier-Stokes systems (Euler2D, NavierStokes2D, NavierStokes3D).

Note: the Dirichlet boundary condition (_DIRICHLET_) could be used as well for supersonic inflow; however the specified Dirichlet state should be in terms of the conserved variables, while the specified supersonic inflow state here is in terms of the flow variables.

Parameters

b	Boundary object of type DomainBoundary
m	MPI object of type MPIVariables
ndims	Number of spatial dimensions
nvars	Number of variables/DoFs per grid point
size	Integer array with the number of grid points in each spatial dimension
ghosts	Number of ghost points
phi	The solution array on which to apply the boundary condition
waqt	Current solution time

Definition at line 24 of file BCSupersonicInflow.c.

{
  DomainBoundary *boundary = (DomainBoundary*) b;

  if (ndims == 2) {

    double gamma;
    gamma = boundary->gamma;
    double inv_gamma_m1 = 1.0/(gamma-1.0);

    if (boundary->on_this_proc) {
      int bounds[ndims], indexb[ndims];
      _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
      _ArraySetValue_(indexb,ndims,0);
      int done = 0;
      while (!done) {
        int p1; _ArrayIndex1DWO_(ndims,size,indexb,boundary->is,ghosts,p1);
        
        /* set the ghost point values */
        double rho_gpt, uvel_gpt, vvel_gpt, energy_gpt, pressure_gpt;
        rho_gpt      = boundary->FlowDensity;
        pressure_gpt = boundary->FlowPressure;
        uvel_gpt     = boundary->FlowVelocity[0];
        vvel_gpt     = boundary->FlowVelocity[1];
        energy_gpt   = inv_gamma_m1*pressure_gpt
                       + 0.5 * rho_gpt * (uvel_gpt*uvel_gpt + vvel_gpt*vvel_gpt);

        phi[nvars*p1+0] = rho_gpt;
        phi[nvars*p1+1] = rho_gpt * uvel_gpt;
        phi[nvars*p1+2] = rho_gpt * vvel_gpt;
        phi[nvars*p1+3] = energy_gpt;

        _ArrayIncrementIndex_(ndims,bounds,indexb,done);
      }
    }

  } else if (ndims == 3) {

    double gamma;
    gamma = boundary->gamma;
    double inv_gamma_m1 = 1.0/(gamma-1.0);

    if (boundary->on_this_proc) {
      int bounds[ndims], indexb[ndims];
      _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
      _ArraySetValue_(indexb,ndims,0);
      int done = 0;
      while (!done) {
        int p1; _ArrayIndex1DWO_(ndims,size,indexb,boundary->is,ghosts,p1);
        
        /* set the ghost point values */
        double rho_gpt, uvel_gpt, vvel_gpt, wvel_gpt, energy_gpt, pressure_gpt;
        rho_gpt      = boundary->FlowDensity;
        pressure_gpt = boundary->FlowPressure;
        uvel_gpt     = boundary->FlowVelocity[0];
        vvel_gpt     = boundary->FlowVelocity[1];
        wvel_gpt     = boundary->FlowVelocity[2];
        energy_gpt   = inv_gamma_m1*pressure_gpt
                       + 0.5 * rho_gpt 
                       * (uvel_gpt*uvel_gpt + vvel_gpt*vvel_gpt + wvel_gpt*wvel_gpt);

        phi[nvars*p1+0] = rho_gpt;
        phi[nvars*p1+1] = rho_gpt * uvel_gpt;
        phi[nvars*p1+2] = rho_gpt * vvel_gpt;
        phi[nvars*p1+3] = rho_gpt * wvel_gpt;
        phi[nvars*p1+4] = energy_gpt;

        _ArrayIncrementIndex_(ndims,bounds,indexb,done);
      }
    }

  }
  return(0);
}

BCSupersonicOutflowU()

int BCSupersonicOutflowU	(	void *	b,
		void *	m,
		int	ndims,
		int	nvars,
		int *	size,
		int	ghosts,
		double *	phi,
		double	waqt
	)

Supersonic outflow boundary conditions for the solution vector U

Applies the supersonic outflow boundary condition: All flow variables (density, pressure, velocity) are extrapolated from the interior since the outflow is supersonic. This boundary condition is specific to two and three dimensional Euler/Navier-Stokes systems (Euler2D, NavierStokes2D, NavierStokes3D).

Note: The extrapolate boundary condition (_EXTRAPOLATE_) can be used as well for this boundary. I am not entirely sure why I wrote the code for this boundary in such a complicated fashion.

Parameters

b	Boundary object of type DomainBoundary
m	MPI object of type MPIVariables
ndims	Number of spatial dimensions
nvars	Number of variables/DoFs per grid point
size	Integer array with the number of grid points in each spatial dimension
ghosts	Number of ghost points
phi	The solution array on which to apply the boundary condition
waqt	Current solution time

Definition at line 24 of file BCSupersonicOutflow.c.

{
  DomainBoundary *boundary = (DomainBoundary*) b;

  int dim   = boundary->dim;
  int face  = boundary->face;

  if (ndims == 2) {

    /* create a fake physics object */
    Euler2D physics; 
    double gamma;
    gamma = physics.gamma = boundary->gamma;
    double inv_gamma_m1 = 1.0/(gamma-1.0);

    if (boundary->on_this_proc) {
      int bounds[ndims], indexb[ndims], indexi[ndims];
      _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
      _ArraySetValue_(indexb,ndims,0);
      int done = 0;
      while (!done) {
        int p1, p2;
        _ArrayCopy1D_(indexb,indexi,ndims);
        _ArrayAdd1D_(indexi,indexi,boundary->is,ndims);
        if      (face ==  1) indexi[dim] = ghosts-1-indexb[dim];
        else if (face == -1) indexi[dim] = size[dim]-indexb[dim]-1;
        else return(1);
        _ArrayIndex1DWO_(ndims,size,indexb,boundary->is,ghosts,p1);
        _ArrayIndex1D_(ndims,size,indexi,ghosts,p2);
        
        /* flow variables in the interior */
        double rho, uvel, vvel, energy, pressure;
        double rho_gpt, uvel_gpt, vvel_gpt, energy_gpt, pressure_gpt;
        _Euler2DGetFlowVar_((phi+nvars*p2),rho,uvel,vvel,energy,pressure,(&physics));
        /* set the ghost point values */
        rho_gpt       = rho;
        pressure_gpt  = pressure;
        uvel_gpt      = uvel;
        vvel_gpt      = vvel;
        energy_gpt    = inv_gamma_m1*pressure_gpt
                        + 0.5 * rho_gpt * (uvel_gpt*uvel_gpt + vvel_gpt*vvel_gpt);

        phi[nvars*p1+0] = rho_gpt;
        phi[nvars*p1+1] = rho_gpt * uvel_gpt;
        phi[nvars*p1+2] = rho_gpt * vvel_gpt;
        phi[nvars*p1+3] = energy_gpt;

        _ArrayIncrementIndex_(ndims,bounds,indexb,done);
      }
    }

  } else if (ndims == 3) {

    /* create a fake physics object */
    double gamma;
    gamma = boundary->gamma;
    double inv_gamma_m1 = 1.0/(gamma-1.0);

    if (boundary->on_this_proc) {
      int bounds[ndims], indexb[ndims], indexi[ndims];
      _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
      _ArraySetValue_(indexb,ndims,0);
      int done = 0;
      while (!done) {
        int p1, p2;
        _ArrayCopy1D_(indexb,indexi,ndims);
        _ArrayAdd1D_(indexi,indexi,boundary->is,ndims);
        if      (face ==  1) indexi[dim] = ghosts-1-indexb[dim];
        else if (face == -1) indexi[dim] = size[dim]-indexb[dim]-1;
        else return(1);
        _ArrayIndex1DWO_(ndims,size,indexb,boundary->is,ghosts,p1);
        _ArrayIndex1D_(ndims,size,indexi,ghosts,p2);
        
        /* flow variables in the interior */
        double rho, uvel, vvel, wvel, energy, pressure;
        double rho_gpt, uvel_gpt, vvel_gpt, wvel_gpt, energy_gpt, pressure_gpt;
        _NavierStokes3DGetFlowVar_((phi+nvars*p2),_NavierStokes3D_stride_,rho,uvel,vvel,wvel,energy,pressure,gamma);
        /* set the ghost point values */
        rho_gpt       = rho;
        pressure_gpt  = pressure;
        uvel_gpt      = uvel;
        vvel_gpt      = vvel;
        wvel_gpt      = wvel;
        energy_gpt    = inv_gamma_m1*pressure_gpt
                        + 0.5 * rho_gpt 
                        * (uvel_gpt*uvel_gpt + vvel_gpt*vvel_gpt + wvel_gpt*wvel_gpt);

        phi[nvars*p1+0] = rho_gpt;
        phi[nvars*p1+1] = rho_gpt * uvel_gpt;
        phi[nvars*p1+2] = rho_gpt * vvel_gpt;
        phi[nvars*p1+3] = rho_gpt * wvel_gpt;
        phi[nvars*p1+4] = energy_gpt;

        _ArrayIncrementIndex_(ndims,bounds,indexb,done);
      }
    }

  }
  return(0);
}

BCTurbulentSupersonicInflowU()

int BCTurbulentSupersonicInflowU	(	void *	b,
		void *	m,
		int	ndims,
		int	nvars,
		int *	size,
		int	ghosts,
		double *	phi,
		double	waqt
	)

Turbulent Supersonic inflow boundary conditions for the solution vector U

Applies the turbulent supersonic inflow boundary condition: The inflow consists of a mean supersonic inflow on which turbulent flow fluctuations are added. This boundary condition is specific to the 3D Navier-Stokes system (NavierStokes3D).

Note: Some parts of the code may be hardcoded for use with the shock-turbulence interaction problem (for which this boundary condition was written).

Parameters

b	Boundary object of type DomainBoundary
m	MPI object of type MPIVariables
ndims	Number of spatial dimensions
nvars	Number of variables/DoFs per grid point
size	Integer array with the number of grid points in each spatial dimension
ghosts	Number of ghost points
phi	The solution array on which to apply the boundary condition
waqt	Current solution time

Definition at line 21 of file BCTurbulentSupersonicInflow.c.

{
  DomainBoundary *boundary = (DomainBoundary*) b;

  int dim   = boundary->dim;

  double *inflow_data = boundary->UnsteadyDirichletData;
  int    *inflow_size = boundary->UnsteadyDirichletSize;

  if (ndims == 3) {

    /* create a fake physics object */
    double gamma;
    gamma = boundary->gamma;
    double inv_gamma_m1 = 1.0/(gamma-1.0);

    if (boundary->on_this_proc) {
      /* the following bit is hardcoded for the inflow data
       * representing fluctuations in a domain of length 2pi */
      double  xt = boundary->FlowVelocity[dim] * waqt;
      int     N  = inflow_size[dim];
      double  L  = 2.0 * (4.0*atan(1.0));
      int     it = ((int) ((xt/L) * ((double)N))) % N;

      int bounds[ndims], indexb[ndims];
      _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
      _ArraySetValue_(indexb,ndims,0);
      int done = 0;
      while (!done) {
        int p1; _ArrayIndex1DWO_(ndims,size,indexb,boundary->is,ghosts,p1);
        
        /* set the ghost point values - mean flow */
        double rho_gpt, uvel_gpt, vvel_gpt, wvel_gpt, energy_gpt, pressure_gpt;
        rho_gpt      = boundary->FlowDensity;
        pressure_gpt = boundary->FlowPressure;
        uvel_gpt     = boundary->FlowVelocity[0];
        vvel_gpt     = boundary->FlowVelocity[1];
        wvel_gpt     = boundary->FlowVelocity[2];

        /* calculate the turbulent fluctuations */
        double duvel , dvvel , dwvel ;
        int index1[ndims]; _ArrayCopy1D_(indexb,index1,ndims);
        index1[dim] = it;
        int q; _ArrayIndex1D_(ndims,inflow_size,index1,0,q);
        duvel = inflow_data[q*nvars+1];
        dvvel = inflow_data[q*nvars+2];
        dwvel = inflow_data[q*nvars+3];

        /* add the turbulent fluctuations to the velocity field */
        uvel_gpt      += duvel;
        vvel_gpt      += dvvel;
        wvel_gpt      += dwvel;

        /* set the ghost point values */
        energy_gpt   = inv_gamma_m1*pressure_gpt
                       + 0.5 * rho_gpt 
                       * (uvel_gpt*uvel_gpt + vvel_gpt*vvel_gpt + wvel_gpt*wvel_gpt);
        phi[nvars*p1+0] = rho_gpt;
        phi[nvars*p1+1] = rho_gpt * uvel_gpt;
        phi[nvars*p1+2] = rho_gpt * vvel_gpt;
        phi[nvars*p1+3] = rho_gpt * wvel_gpt;
        phi[nvars*p1+4] = energy_gpt;

        _ArrayIncrementIndex_(ndims,bounds,indexb,done);
      }
    }

  }
  return(0);
}

BCNoFluxU()

int BCNoFluxU	(	void *	b,
		void *	m,
		int	ndims,
		int	nvars,
		int *	size,
		int	ghosts,
		double *	phi,
		double	waqt
	)

No-Flux (inviscid wall) boundary conditions for the solution vector U

Applies the no-flux boundary conditions: This boundary condition is specific to the NUMA 2D/3D (Numa2D, Numa3D). Used for simulating inviscid walls or symmetry boundaries. It's equivalent to the slip-wall BC of the Euler/Navier- Stokes system.

The density, potential temperature, and tangential velocity are extrapolated, while the normal velocity at the ghost point is set to the negative of that in the interior (to enforce zero-normal velocity at the boundary face).

Parameters

b	Boundary object of type DomainBoundary
m	MPI object of type MPIVariables
ndims	Number of spatial dimensions
nvars	Number of variables/DoFs per grid point
size	Integer array with the number of grid points in each spatial dimension
ghosts	Number of ghost points
phi	The solution array on which to apply the boundary condition
waqt	Current solution time

Definition at line 22 of file BCNoFlux.c.

{
  DomainBoundary *boundary = (DomainBoundary*) b;

  int dim   = boundary->dim;
  int face  = boundary->face;

  if (boundary->on_this_proc) {
    int bounds[ndims], indexb[ndims], indexi[ndims];
    _ArraySubtract1D_ (bounds,boundary->ie,boundary->is,ndims);
    _ArraySetValue_   (indexb,ndims,0);
    int done = 0;
    while (!done) {
      int p1, p2;
      _ArrayCopy1D_ (indexb,indexi,ndims);
      _ArrayAdd1D_  (indexi,indexi,boundary->is,ndims);
      if      (face ==  1) indexi[dim] = ghosts-1-indexb[dim];
      else if (face == -1) indexi[dim] = size[dim]-indexb[dim]-1;
      else return(1);
      _ArrayIndex1DWO_  (ndims,size,indexb,boundary->is,ghosts,p1);
      _ArrayIndex1D_    (ndims,size,indexi,ghosts,p2);
      
      if (nvars == 4) {
        phi[nvars*p1+0] = phi[nvars*p2+0];
        phi[nvars*p1+1] = (dim == _XDIR_ ? -phi[nvars*p2+1] : phi[nvars*p2+1] );
        phi[nvars*p1+2] = (dim == _YDIR_ ? -phi[nvars*p2+2] : phi[nvars*p2+2] );
        phi[nvars*p1+3] = phi[nvars*p2+3];
      } else if (nvars == 5) {
        phi[nvars*p1+0] = phi[nvars*p2+0];
        phi[nvars*p1+1] = (dim == _XDIR_ ? -phi[nvars*p2+1] : phi[nvars*p2+1] );
        phi[nvars*p1+2] = (dim == _YDIR_ ? -phi[nvars*p2+2] : phi[nvars*p2+2] );
        phi[nvars*p1+3] = (dim == _ZDIR_ ? -phi[nvars*p2+3] : phi[nvars*p2+3] );
        phi[nvars*p1+4] = phi[nvars*p2+4];
      }

      _ArrayIncrementIndex_(ndims,bounds,indexb,done);
    }
  }
  return(0);
}

BCSWSlipWallU()

int BCSWSlipWallU	(	void *	b,
		void *	m,
		int	ndims,
		int	nvars,
		int *	size,
		int	ghosts,
		double *	phi,
		double	waqt
	)

Slip (inviscid) wall boundary conditions for the solution vector U

Applies the slip-wall boundary condition: This is specific to the one and two dimenstional shallow water equations (ShallowWater1D, ShallowWater2D). It is used for simulating inviscid walls or symmetric boundaries. The height, and tangential velocity at the ghost points are extrapolated from the interior, while the normal velocity at the ghost points is set such that the interpolated value at the boundary face is equal to the specified wall velocity.

Parameters

b	Boundary object of type DomainBoundary
m	MPI object of type MPIVariables
ndims	Number of spatial dimensions
nvars	Number of variables/DoFs per grid point
size	Integer array with the number of grid points in each spatial dimension
ghosts	Number of ghost points
phi	The solution array on which to apply the boundary condition
waqt	Current solution time

Definition at line 21 of file BCSWSlipWall.c.

{
  DomainBoundary *boundary = (DomainBoundary*) b;

  int dim   = boundary->dim;
  int face  = boundary->face;

  if (ndims == 1) {

    if (boundary->on_this_proc) {
      int bounds[ndims], indexb[ndims], indexi[ndims];
      _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
      _ArraySetValue_(indexb,ndims,0);
      int done = 0;
      while (!done) {
        int p1, p2;
        _ArrayCopy1D_(indexb,indexi,ndims);
        _ArrayAdd1D_(indexi,indexi,boundary->is,ndims);
        if      (face ==  1) indexi[dim] = ghosts-1-indexb[dim];
        else if (face == -1) indexi[dim] = size[dim]-indexb[dim]-1;
        else return(1);
        _ArrayIndex1DWO_(ndims,size,indexb,boundary->is,ghosts,p1);
        _ArrayIndex1D_(ndims,size,indexi,ghosts,p2);
        
        /* flow variables in the interior */
        double h, uvel;
        double h_gpt, uvel_gpt;
        _ShallowWater1DGetFlowVar_((phi+nvars*p2),h,uvel);
        /* set the ghost point values */
        h_gpt = h;
        uvel_gpt = 2.0*boundary->FlowVelocity[_XDIR_] - uvel;

        phi[nvars*p1+0] = h_gpt;
        phi[nvars*p1+1] = h_gpt * uvel_gpt;

        _ArrayIncrementIndex_(ndims,bounds,indexb,done);
      }
    }

  } else if (ndims == 2) {

    if (boundary->on_this_proc) {
      int bounds[ndims], indexb[ndims], indexi[ndims];
      _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
      _ArraySetValue_(indexb,ndims,0);
      int done = 0;
      while (!done) {
        int p1, p2;
        _ArrayCopy1D_(indexb,indexi,ndims);
        _ArrayAdd1D_(indexi,indexi,boundary->is,ndims);
        if      (face ==  1) indexi[dim] = ghosts-1-indexb[dim];
        else if (face == -1) indexi[dim] = size[dim]-indexb[dim]-1;
        else return(1);
        _ArrayIndex1DWO_(ndims,size,indexb,boundary->is,ghosts,p1);
        _ArrayIndex1D_(ndims,size,indexi,ghosts,p2);
        
        /* flow variables in the interior */
        double h, uvel, vvel;
        double h_gpt, uvel_gpt, vvel_gpt;
        _ShallowWater2DGetFlowVar_((phi+nvars*p2),h,uvel,vvel);
        /* set the ghost point values */
        h_gpt = h;
        if (dim == _XDIR_) {
          uvel_gpt = 2.0*boundary->FlowVelocity[_XDIR_] - uvel;
          vvel_gpt = vvel;
        } else if (dim == _YDIR_) {
          uvel_gpt = uvel;
          vvel_gpt = 2.0*boundary->FlowVelocity[_YDIR_] - vvel;
        } else {
          uvel_gpt = 0.0;
          vvel_gpt = 0.0;
        }

        phi[nvars*p1+0] = h_gpt;
        phi[nvars*p1+1] = h_gpt * uvel_gpt;
        phi[nvars*p1+2] = h_gpt * vvel_gpt;

        _ArrayIncrementIndex_(ndims,bounds,indexb,done);
      }
    }

  }
  return(0);
}

BCSpongeSource()

int BCSpongeSource	(	void *	b,
		int	ndims,
		int	nvars,
		int	ghosts,
		int *	size,
		double *	grid,
		double *	u,
		double *	source
	)

a special BC enforcement - an absorbent sponge - enforced through a source term

Applies the sponge boundary condition: This function computes the source term required to apply a sponge boundary condition that gradually relaxes the solution to a specified state. This boundary condition is different from other boundary conditions in the sense that it is applied at interior grid points (but within the defined sponge zone).

The source term for the sponge is computed as:

\begin{align} {\bf S}_i &= \sigma_i \left( {\bf U}_i - {\bf U}_{\rm ref} \right),\\ \sigma_i &= \frac {x_i - x_{\rm start}} {x_{\rm end} - x_{\rm start}} \end{align}

where \(i\) is the grid index along the spatial dimension of the sponge, \({\bf U}_{\rm ref}\) is the specified state to which the solution is relaxed, and \(x_i\), \(x_{\rm start}\), and \(x_{\rm end}\) are the spatial coordinates of the grid point, sponge start, and sponge end, respectively, along the spatial dimension of the sponge.

Parameters

b	Boundary object of type DomainBoundary
ndims	Number of spatial dimensions
nvars	Number of variables/DoFs per grid point
ghosts	Number of ghost points
size	Integer array with the number of grid points in each spatial dimension
grid	1D array with the spatial coordinates of the grid points, one dimension after the other
u	Solution
source	Source term to which the sponge term is added

Definition at line 26 of file BCSponge.c.

{
  DomainBoundary *boundary = (DomainBoundary*) b;
  int            dim       = boundary->dim;
  int            face      = boundary->face;
  double         *uref     = boundary->SpongeValue;
  double         *xmin     = boundary->xmin;
  double         *xmax     = boundary->xmax;
  int            v;

  if (boundary->on_this_proc) {
    int bounds[ndims], indexb[ndims];
    _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
    _ArraySetValue_(indexb,ndims,0); 
    int done = 0;
    while (!done) {
      int i = indexb[dim] + boundary->is[dim];
      double x, xstart, xend;
      _GetCoordinate_(dim,i,size,ghosts,grid,x);
      xstart = xmin[dim];
      xend   = xmax[dim];
      /* calculate sigma */
      double sigma;
      if (face > 0) sigma = (x - xstart) / (xend - xstart);
      else          sigma = (x - xend  ) / (xstart - xend);
      /* add to the source term */
      int p; _ArrayIndex1DWO_(ndims,size,indexb,boundary->is,ghosts,p);
      for (v=0; v<nvars; v++) source[nvars*p+v] -= (sigma * (u[nvars*p+v]-uref[v]));
      _ArrayIncrementIndex_(ndims,bounds,indexb,done);
    }
  }
  return(0);
}

BCSpongeUDummy()

int BCSpongeUDummy	(	void *	b,
		void *	m,
		int	ndims,
		int	nvars,
		int *	size,
		int	ghosts,
		double *	phi,
		double	waqt
	)

dummy function that get called during applying BCs - they don't do anything

Dummy function to ensure consistency with the overall boundary condition implementation. The actual sponge boundary condition is implemented by BCSpongeSource()

Parameters

b	Boundary object of type DomainBoundary
m	MPI object of type MPIVariables
ndims	Number of spatial dimensions
nvars	Number of variables/DoFs per grid point
size	Integer array with the number of grid points in each spatial dimension
ghosts	Number of ghost points
phi	The solution array on which to apply the boundary condition
waqt	Current solution time

Definition at line 73 of file BCSponge.c.

{
  return(0);
}

BCReadTemperatureData()

int BCReadTemperatureData	(	void *	b,
		void *	m,
		int	ndims,
		int	nvars,
		int *	DomainSize
	)

Function to read in unsteady temperature data for thermal slip wall BC (BCThermalSlipWallU())

Read in the temperature data: The temperature data needs to be provided as a binary file. For parallel runs, only rank 0 reads the file, and then distributes the data to the other processors.

This function needs to be better documented.

Definition at line 147 of file BCIO.c.

{
  DomainBoundary *boundary = (DomainBoundary*) b;
  MPIVariables   *mpi      = (MPIVariables*)   m;

  char    *filename = boundary->UnsteadyTemperatureFilename;
  int     *temperature_field_size  = NULL;
  double  *time_level_data         = NULL;
  double  *temperature_field_data  = NULL;
  double  *time_buffer = NULL;
  double  *data_buffer = NULL;
  
  int dim = boundary->dim;
  int face= boundary->face;
  int d;

  if (!mpi->rank) {

    printf("Reading boundary temperature data from %s.\n",filename);

    FILE *in;
    int  ferr;

    /* calculate the number of processors that sit on this temperature boundary */
    int nproc = 1;
    for (d=0; d<ndims; d++) nproc *= mpi->iproc[d]; nproc /= mpi->iproc[dim];

    in = fopen(filename,"rb");
    if (!in) {
      fprintf(stderr,"Error in BCReadTemperatureData(): cannot open boundary temperature data file %s.\n",filename);
      return(1);
    }

    int count = 0;
    while ((!feof(in)) && (count < nproc)) {

      int rank[ndims], size[ndims];
      ferr = fread(rank,sizeof(int),ndims,in);
      if (ferr != ndims) {
        fprintf(stderr,"Error in BCReadTemperatureData(): Error (1) in file reading, count %d.\n",count);
        return(1);
      }
      if (rank[dim] != (face > 0 ? 0 : mpi->iproc[dim]-1) ) {
        fprintf(stderr,"Error in BCReadTemperatureData(): Error (2) in file reading, count %d.\n",count);
        return(1);
      }
      ferr = fread(size,sizeof(int),ndims,in);
      if (ferr != ndims) {
        fprintf(stderr,"Error in BCReadTemperatureData(): Error (3) in file reading, count %d.\n",count);
        return(1);
      }

      int n_data = size[dim]; /* number of time levels for which temperature data is provided */
      time_buffer = (double*) calloc (n_data,sizeof(double));
      ferr = fread(time_buffer,sizeof(double),n_data,in);
      if (ferr != n_data) {
        fprintf(stderr,"Error in BCReadTemperatureData(): Error (5) in file reading, count %d.\n",count);
        return(1);
      }

      int data_size = 1;
      for (d=0; d<ndims; d++) data_size *= size[d];
      data_buffer = (double*) calloc (data_size,sizeof(double));
      ferr = fread(data_buffer,sizeof(double),data_size,in);
      if (ferr != data_size) {
        fprintf(stderr,"Error in BCReadTemperatureData(): Error (6) in file reading, count %d.\n",count);
        return(1);
      }

      int rank1D = MPIRank1D(ndims,mpi->iproc,rank);

      if (!rank1D) {

        int index[ndims];

        temperature_field_size = (int*) calloc (ndims, sizeof(int));
        _ArrayCopy1D_(size,temperature_field_size,ndims);

        time_level_data = (double*) calloc (size[dim], sizeof(double));
        _ArrayCopy1D_(time_buffer,time_level_data,size[dim]);

        temperature_field_data = (double*) calloc (data_size, sizeof(double));
        ArrayCopynD(ndims,data_buffer,temperature_field_data,size,0,0,index,1);

      } else {

#ifndef serial
        MPI_Request req[3] = {MPI_REQUEST_NULL,MPI_REQUEST_NULL,MPI_REQUEST_NULL};
        MPI_Isend(size,ndims,MPI_INT,rank1D,2152,mpi->world,&req[0]);
        MPI_Isend(time_buffer,size[dim],MPI_DOUBLE,rank1D,2154,mpi->world,&req[2]);
        MPI_Isend(data_buffer,data_size,MPI_DOUBLE,rank1D,2153,mpi->world,&req[1]);
        MPI_Status status_arr[3];
        MPI_Waitall(3,&req[0],status_arr);
#else
        fprintf(stderr,"Error in BCReadTemperatureData(): This is a serial run. Invalid (non-zero) rank read.\n");
#endif

      }

      free(time_buffer);
      free(data_buffer);
      count++;
    }

    if (count < nproc) {
      fprintf(stderr,"Error in BCReadTemperatureData(): missing data in unsteady boundary data file %s.\n",filename);
      fprintf(stderr,"Error in BCReadTemperatureData(): should contain data for %d processors, ", nproc);
      fprintf(stderr,"Error in BCReadTemperatureData(): but contains data for %d processors!\n", count);
      return(1);
    }

    fclose(in);

  } else {

#ifndef serial
    if (mpi->ip[dim] == (face > 0 ? 0 : mpi->iproc[dim]-1) ) {

      MPI_Request req = MPI_REQUEST_NULL;

      temperature_field_size = (int*) calloc (ndims,sizeof(int));
      MPI_Irecv(temperature_field_size,ndims,MPI_INT,0,2152,mpi->world,&req);
      MPI_Wait(&req,MPI_STATUS_IGNORE);

      int flag = 1;
      for (d=0; d<ndims; d++) if ((d != dim) && (temperature_field_size[d] != DomainSize[d])) flag = 0;
      if (!flag) {
        fprintf(stderr,"Error in BCReadTemperatureData(): Error (4) (dimension mismatch) in file reading, rank %d.\n",mpi->rank);
        return(1);
      }

      time_level_data = (double*) calloc (temperature_field_size[dim],sizeof(double));
      MPI_Irecv(time_level_data, temperature_field_size[dim], MPI_DOUBLE,0,2154,mpi->world,&req);
      MPI_Wait(&req,MPI_STATUS_IGNORE);

      int data_size = 1;
      for (d=0; d<ndims; d++) data_size *= temperature_field_size[d];
      temperature_field_data = (double*) calloc (data_size,sizeof(double));
      MPI_Irecv(temperature_field_data,data_size,MPI_DOUBLE,0,2153,mpi->world,&req);
      MPI_Wait(&req,MPI_STATUS_IGNORE);

    }
#else
    fprintf(stderr,"Error in BCReadTemperatureData(): Serial code should not be here!.\n");
#endif

  }
  
  boundary->UnsteadyTemperatureSize = temperature_field_size;
  boundary->UnsteadyTimeLevels      = time_level_data;
  boundary->UnsteadyTemperatureData = temperature_field_data;

  return(0);
}

BCReadTurbulentInflowData()

int BCReadTurbulentInflowData	(	void *	b,
		void *	m,
		int	ndims,
		int	nvars,
		int *	DomainSize
	)

Function to read in unsteady boundary data for turbulent inflow

Read in the turbulent inflow data: The turbulent inflow data needs to be provided as a binary file. For parallel runs, only rank 0 reads the file, and then distributes the data to the other processors.

This function needs to be better documented.

Definition at line 19 of file BCIO.c.

{
  DomainBoundary *boundary = (DomainBoundary*) b;
  MPIVariables   *mpi      = (MPIVariables*)   m;

  char    *filename     = boundary->UnsteadyDirichletFilename;
  int     *inflow_size  = NULL;
  double  *inflow_data  = NULL;
  double  *buffer       = NULL;
  
  int dim = boundary->dim;
  int face= boundary->face;
  int d;

  if (!mpi->rank) {

    printf("Reading turbulent inflow boundary data from %s.\n",filename);

    FILE *in;
    int  ferr;

    /* calculate the number of processors that sit on unsteady boundary */
    int nproc = 1;
    for (d=0; d<ndims; d++) nproc *= mpi->iproc[d]; nproc /= mpi->iproc[dim];

    in = fopen(filename,"rb");
    if (!in) {
      fprintf(stderr,"Error in BCReadTurbulentInflowData(): cannot open unsteady boundary data file %s.\n",filename);
      return(1);
    }
    int count = 0;
    while ((!feof(in)) && (count < nproc)) {
      int rank[ndims], size[ndims];
      ferr = fread(rank,sizeof(int),ndims,in);
      if (ferr != ndims) {
        fprintf(stderr,"Error in BCReadTurbulentInflowData(): Error (1) in file reading, count %d.\n",count);
        return(1);
      }
      if (rank[dim] != (face > 0 ? 0 : mpi->iproc[dim]-1) ) {
        fprintf(stderr,"Error in BCReadTurbulentInflowData(): Error (2) in file reading, count %d.\n",count);
        return(1);
      }
      ferr = fread(size,sizeof(int),ndims,in);
      if (ferr != ndims) {
        fprintf(stderr,"Error in BCReadTurbulentInflowData(): Error (3) in file reading, count %d.\n",count);
        return(1);
      }
      int flag = 1;
      for (d=0; d<ndims; d++) if ((d != dim) && (size[d] != DomainSize[d])) flag = 0;
      if (!flag) {
        fprintf(stderr,"Error in BCReadTurbulentInflowData(): Error (4) (dimension mismatch) in file reading, count %d.\n",count);
        return(1);
      }

      int data_size = nvars;
      for (d=0; d<ndims; d++) data_size *= size[d];
      buffer = (double*) calloc (data_size,sizeof(double));
      ferr = fread(buffer,sizeof(double),data_size,in);
      if (ferr != data_size) {
        fprintf(stderr,"Error in BCReadTurbulentInflowData(): Error (6) in file reading, count %d.\n",count);
        return(1);
      }

      int rank1D = MPIRank1D(ndims,mpi->iproc,rank);

      if (!rank1D) {

        int index[ndims];
        inflow_size = (int*) calloc (ndims, sizeof(int));
        _ArrayCopy1D_(size,inflow_size,ndims);
        inflow_data = (double*) calloc (data_size, sizeof(double));
        ArrayCopynD(ndims,buffer,inflow_data,size,0,0,index,nvars);

      } else {
#ifndef serial
        MPI_Request req[2] = {MPI_REQUEST_NULL,MPI_REQUEST_NULL};
        MPI_Isend(size,ndims,MPI_INT,rank1D,2152,mpi->world,&req[0]);
        MPI_Isend(buffer,data_size,MPI_DOUBLE,rank1D,2153,mpi->world,&req[1]);
        MPI_Status status_arr[3];
        MPI_Waitall(2,&req[0],status_arr);
#else
        fprintf(stderr,"Error in BCReadTurbulentInflowData(): This is a serial run. Invalid (non-zero) rank read.\n");
#endif
      }

      free(buffer);
      count++;
    }

    if (count < nproc) {
      fprintf(stderr,"Error in BCReadTurbulentInflowData(): missing data in unsteady boundary data file %s.\n",filename);
      fprintf(stderr,"Error in BCReadTurbulentInflowData(): should contain data for %d processors, ", nproc);
      fprintf(stderr,"Error in BCReadTurbulentInflowData(): but contains data for %d processors!\n", count);
      return(1);
    }

    fclose(in);

  } else {
#ifndef serial
    if (mpi->ip[dim] == (face > 0 ? 0 : mpi->iproc[dim]-1) ) {
      MPI_Request req = MPI_REQUEST_NULL;
      inflow_size = (int*) calloc (ndims,sizeof(int));
      MPI_Irecv(inflow_size,ndims,MPI_INT,0,2152,mpi->world,&req);
      MPI_Wait(&req,MPI_STATUS_IGNORE);
      int data_size = nvars;
      for (d=0; d<ndims; d++) data_size *= inflow_size[d];
      inflow_data = (double*) calloc (data_size,sizeof(double));
      MPI_Irecv(inflow_data,data_size,MPI_DOUBLE,0,2153,mpi->world,&req);
      MPI_Wait(&req,MPI_STATUS_IGNORE);
    }
#else
    fprintf(stderr,"Error in BCReadTurbulentInflowData(): Serial code should not be here!.\n");
#endif
  }
  
  boundary->UnsteadyDirichletSize = inflow_size;
  boundary->UnsteadyDirichletData = inflow_data;

  return(0);
}

gpuBCPeriodicU()

int gpuBCPeriodicU	(	void *	,
		void *	,
		int	,
		int	,
		int *	,
		int	,
		double *	,
		double
	)

Periodic boundary conditions for the solution vector U

gpuBCSlipWallU()

int gpuBCSlipWallU	(	void *	,
		void *	,
		int	,
		int	,
		int *	,
		int	,
		double *	,
		double
	)

Slip (inviscid) wall boundary conditions for the solution vector U