HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
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
bBoundary object of type DomainBoundary
flag_gpuFlag to indicate if GPU is being used

Definition at line 12 of file BCInitialize.c.

14 {
15  DomainBoundary *boundary = (DomainBoundary*) b;
16 
17 #if defined(HAVE_CUDA)
18  if (flag_gpu) {
19 
20  if (!strcmp(boundary->bctype,_PERIODIC_)) {
21  boundary->BCFunctionU = gpuBCPeriodicU;
22  } else if (!strcmp(boundary->bctype,_SLIP_WALL_)) {
23  boundary->BCFunctionU = gpuBCSlipWallU;
24  } else {
25  fprintf(stderr,"[GPU] Error in BCInitialize(): \"%s\" is not a supported boundary condition.\n",
26  boundary->bctype);
27  return 1;
28  }
29 
30  } else {
31 #endif
32 
33  if (!strcmp(boundary->bctype,_PERIODIC_ )) boundary->BCFunctionU = BCPeriodicU;
34  else if (!strcmp(boundary->bctype,_EXTRAPOLATE_ )) boundary->BCFunctionU = BCExtrapolateU;
35  else if (!strcmp(boundary->bctype,_DIRICHLET_ )) boundary->BCFunctionU = BCDirichletU;
36  else if (!strcmp(boundary->bctype,_REFLECT_ )) boundary->BCFunctionU = BCReflectU;
37  else if (!strcmp(boundary->bctype,_SPONGE_ )) boundary->BCFunctionU = BCSpongeUDummy;
38  else if (!strcmp(boundary->bctype,_NOSLIP_WALL_ )) boundary->BCFunctionU = BCNoslipWallU;
39  else if (!strcmp(boundary->bctype,_SLIP_WALL_ )) boundary->BCFunctionU = BCSlipWallU;
40  else if (!strcmp(boundary->bctype,_THERMAL_SLIP_WALL_ )) boundary->BCFunctionU = BCThermalSlipWallU;
41  else if (!strcmp(boundary->bctype,_THERMAL_NOSLIP_WALL_ )) boundary->BCFunctionU = BCThermalNoslipWallU;
42  else if (!strcmp(boundary->bctype,_SW_SLIP_WALL_ )) boundary->BCFunctionU = BCSWSlipWallU;
43  else if (!strcmp(boundary->bctype,_SUBSONIC_OUTFLOW_ )) boundary->BCFunctionU = BCSubsonicOutflowU;
44  else if (!strcmp(boundary->bctype,_SUBSONIC_INFLOW_ )) boundary->BCFunctionU = BCSubsonicInflowU;
45  else if (!strcmp(boundary->bctype,_SUBSONIC_AMBIVALENT_ )) boundary->BCFunctionU = BCSubsonicAmbivalentU;
46  else if (!strcmp(boundary->bctype,_SUPERSONIC_OUTFLOW_ )) boundary->BCFunctionU = BCSupersonicOutflowU;
47  else if (!strcmp(boundary->bctype,_SUPERSONIC_INFLOW_ )) boundary->BCFunctionU = BCSupersonicInflowU;
48  else if (!strcmp(boundary->bctype,_TURBULENT_SUPERSONIC_INFLOW_ )) boundary->BCFunctionU = BCTurbulentSupersonicInflowU;
49  else if (!strcmp(boundary->bctype,_NO_FLUX_BC_ )) boundary->BCFunctionU = BCNoFluxU;
50  else {
51  fprintf(stderr,"Error in BCInitialize(): \"%s\" is not a supported boundary condition.\n",
52  boundary->bctype);
53  return(1);
54  }
55 
56 #if defined(HAVE_CUDA)
57  }
58 #endif
59 
60  return 0;
61 }
#define _SUPERSONIC_INFLOW_
int BCSupersonicInflowU(void *, void *, int, int, int *, int, double *, double)
int(* BCFunctionU)(void *, void *, int, int, int *, int, double *, double)
#define _SUBSONIC_AMBIVALENT_
int BCNoFluxU(void *, void *, int, int, int *, int, double *, double)
Definition: BCNoFlux.c:22
int BCSpongeUDummy(void *, void *, int, int, int *, int, double *, double)
Definition: BCSponge.c:73
int gpuBCPeriodicU(void *, void *, int, int, int *, int, double *, double)
#define _DIRICHLET_
#define _THERMAL_SLIP_WALL_
#define _THERMAL_NOSLIP_WALL_
#define _SPONGE_
Structure containing the variables and function pointers defining a boundary.
#define _NO_FLUX_BC_
int BCPeriodicU(void *, void *, int, int, int *, int, double *, double)
Definition: BCPeriodic.c:19
int BCThermalNoslipWallU(void *, void *, int, int, int *, int, double *, double)
int BCSlipWallU(void *, void *, int, int, int *, int, double *, double)
Definition: BCSlipWall.c:22
int BCSubsonicAmbivalentU(void *, void *, int, int, int *, int, double *, double)
int BCNoslipWallU(void *, void *, int, int, int *, int, double *, double)
Definition: BCNoslipWall.c:20
int BCReflectU(void *, void *, int, int, int *, int, double *, double)
Definition: BCReflect.c:14
#define _SW_SLIP_WALL_
int BCSubsonicInflowU(void *, void *, int, int, int *, int, double *, double)
char bctype[_MAX_STRING_SIZE_]
int BCSWSlipWallU(void *, void *, int, int, int *, int, double *, double)
Definition: BCSWSlipWall.c:21
#define _SUBSONIC_INFLOW_
int BCSubsonicOutflowU(void *, void *, int, int, int *, int, double *, double)
#define _EXTRAPOLATE_
#define _REFLECT_
#define _SUBSONIC_OUTFLOW_
int BCDirichletU(void *, void *, int, int, int *, int, double *, double)
Definition: BCDirichlet.c:13
int gpuBCSlipWallU(void *, void *, int, int, int *, int, double *, double)
int BCSupersonicOutflowU(void *, void *, int, int, int *, int, double *, double)
#define _NOSLIP_WALL_
#define _SUPERSONIC_OUTFLOW_
int BCExtrapolateU(void *, void *, int, int, int *, int, double *, double)
Definition: BCExtrapolate.c:13
int BCThermalSlipWallU(void *, void *, int, int, int *, int, double *, double)
int BCTurbulentSupersonicInflowU(void *, void *, int, int, int *, int, double *, double)
#define _SLIP_WALL_
#define _PERIODIC_
#define _TURBULENT_SUPERSONIC_INFLOW_

◆ 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
bBoundary object of type DomainBoundary
flag_gpuFlag indicating if GPU is being used

Definition at line 12 of file BCCleanup.c.

14 {
15  DomainBoundary *boundary = (DomainBoundary*) b;
16  free(boundary->xmin);
17  free(boundary->xmax);
18  free(boundary->is);
19  free(boundary->ie);
20  if (boundary->DirichletValue) free(boundary->DirichletValue);
21  if (boundary->SpongeValue ) free(boundary->SpongeValue );
22  if (boundary->FlowVelocity ) free(boundary->FlowVelocity );
23 
24  if (boundary->UnsteadyDirichletSize) free(boundary->UnsteadyDirichletSize);
25  if (boundary->UnsteadyDirichletData) free(boundary->UnsteadyDirichletData);
26 
27  if (boundary->UnsteadyTemperatureSize) free(boundary->UnsteadyTemperatureSize);
28  if (boundary->UnsteadyTimeLevels) free(boundary->UnsteadyTimeLevels);
29  if (boundary->UnsteadyTemperatureData) free(boundary->UnsteadyTemperatureData);
30 
31 #if defined(HAVE_CUDA)
32  if (flag_gpu) {
33  gpuFree(boundary->gpu_bounds);
34  gpuFree(boundary->gpu_is);
35  gpuFree(boundary->gpu_ie);
36  if (boundary->FlowVelocity) gpuFree(boundary->gpu_FlowVelocity );
37  }
38 #endif
39 
40  return 0;
41 }
Structure containing the variables and function pointers defining a boundary.
double * UnsteadyDirichletData
double * UnsteadyTimeLevels
double * UnsteadyTemperatureData
void gpuFree(void *)

◆ 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
bBoundary object of type DomainBoundary
mMPI object of type MPIVariables
ndimsNumber of spatial dimensions
nvarsNumber of variables/DoFs per grid point
sizeInteger array with the number of grid points in each spatial dimension
ghostsNumber of ghost points
phiThe solution array on which to apply the boundary condition
waqtCurrent solution time

Definition at line 19 of file BCPeriodic.c.

29 {
30  DomainBoundary *boundary = (DomainBoundary*) b;
31  MPIVariables *mpi = (MPIVariables*) m;
32 
33  int dim = boundary->dim;
34  int face = boundary->face;
35 
36  if ((boundary->on_this_proc) && (mpi->iproc[dim] == 1)) {
37  int bounds[ndims], index1[ndims], index2[ndims];
38  _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
39  _ArraySetValue_(index1,ndims,0);
40  _ArraySetValue_(index2,ndims,0);
41  int done = 0;
42  while (!done) {
43  int p1 = 0, p2 = 0;
44  _ArrayCopy1D_(index1,index2,ndims);
45  if (face == 1) {
46  index2[dim] = index1[dim] + size[dim]-ghosts;
47  _ArrayIndex1DWO_(ndims,size,index1,boundary->is,ghosts,p1);
48  _ArrayIndex1D_(ndims,size,index2,ghosts,p2);
49  } else if (face == -1) {
50  _ArrayIndex1DWO_(ndims,size,index1,boundary->is,ghosts,p1);
51  _ArrayIndex1D_(ndims,size,index1,ghosts,p2);
52  }
53  _ArrayCopy1D_((phi+nvars*p2),(phi+nvars*p1),nvars);
54  _ArrayIncrementIndex_(ndims,bounds,index1,done);
55  }
56  }
57 
58  return(0);
59 }
Structure containing the variables and function pointers defining a boundary.
#define _ArraySubtract1D_(x, a, b, size)
#define _ArrayIndex1D_(N, imax, i, ghost, index)
#define _ArrayIndex1DWO_(N, imax, i, offset, ghost, index)
#define _ArraySetValue_(x, size, value)
#define _ArrayIncrementIndex_(N, imax, i, done)
Structure of MPI-related variables.
#define _ArrayCopy1D_(x, y, size)

◆ 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
bBoundary object of type DomainBoundary
mMPI object of type MPIVariables
ndimsNumber of spatial dimensions
nvarsNumber of variables/DoFs per grid point
sizeInteger array with the number of grid points in each spatial dimension
ghostsNumber of ghost points
phiThe solution array on which to apply the boundary condition
waqtCurrent solution time

Definition at line 13 of file BCExtrapolate.c.

23 {
24  DomainBoundary *boundary = (DomainBoundary*) b;
25 
26  int dim = boundary->dim;
27  int face = boundary->face;
28 
29  if (boundary->on_this_proc) {
30  int bounds[ndims], indexb[ndims], indexi[ndims];
31  _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
32  _ArraySetValue_(indexb,ndims,0);
33  int done = 0;
34  while (!done) {
35  _ArrayCopy1D_(indexb,indexi,ndims);
36  _ArrayAdd1D_(indexi,indexi,boundary->is,ndims);
37  if (face == 1) indexi[dim] = ghosts-1-indexb[dim];
38  else if (face == -1) indexi[dim] = size[dim]-indexb[dim]-1;
39  else return(1);
40  int p1,p2;
41  _ArrayIndex1DWO_(ndims,size,indexb,boundary->is,ghosts,p1);
42  _ArrayIndex1D_(ndims,size,indexi,ghosts,p2);
43  _ArrayCopy1D_((phi+nvars*p2),(phi+nvars*p1),nvars);
44  _ArrayIncrementIndex_(ndims,bounds,indexb,done);
45  }
46  }
47  return(0);
48 }
Structure containing the variables and function pointers defining a boundary.
#define _ArraySubtract1D_(x, a, b, size)
#define _ArrayAdd1D_(x, a, b, size)
#define _ArrayIndex1D_(N, imax, i, ghost, index)
#define _ArrayIndex1DWO_(N, imax, i, offset, ghost, index)
#define _ArraySetValue_(x, size, value)
#define _ArrayIncrementIndex_(N, imax, i, done)
#define _ArrayCopy1D_(x, y, size)

◆ 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
bBoundary object of type DomainBoundary
mMPI object of type MPIVariables
ndimsNumber of spatial dimensions
nvarsNumber of variables/DoFs per grid point
sizeInteger array with the number of grid points in each spatial dimension
ghostsNumber of ghost points
phiThe solution array on which to apply the boundary condition
waqtCurrent solution time

Definition at line 13 of file BCDirichlet.c.

23 {
24  DomainBoundary *boundary = (DomainBoundary*) b;
25 
26  if (boundary->on_this_proc) {
27  int bounds[ndims], indexb[ndims];
28  _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
29  _ArraySetValue_(indexb,ndims,0);
30  int done = 0;
31  while (!done) {
32  int p; _ArrayIndex1DWO_(ndims,size ,indexb,boundary->is,ghosts,p);
33  _ArrayCopy1D_((boundary->DirichletValue),(phi+nvars*p),nvars);
34  _ArrayIncrementIndex_(ndims,bounds,indexb,done);
35  }
36  }
37  return(0);
38 }
Structure containing the variables and function pointers defining a boundary.
#define _ArraySubtract1D_(x, a, b, size)
#define _ArrayIndex1DWO_(N, imax, i, offset, ghost, index)
#define _ArraySetValue_(x, size, value)
#define _ArrayIncrementIndex_(N, imax, i, done)
#define _ArrayCopy1D_(x, y, size)

◆ 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
bBoundary object of type DomainBoundary
mMPI object of type MPIVariables
ndimsNumber of spatial dimensions
nvarsNumber of variables/DoFs per grid point
sizeInteger array with the number of grid points in each spatial dimension
ghostsNumber of ghost points
phiThe solution array on which to apply the boundary condition
waqtCurrent solution time

Definition at line 14 of file BCReflect.c.

24 {
25  DomainBoundary *boundary = (DomainBoundary*) b;
26 
27  int dim = boundary->dim;
28  int face = boundary->face;
29 
30  if (boundary->on_this_proc) {
31  int bounds[ndims], indexb[ndims], indexi[ndims];
32  _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
33  _ArraySetValue_(indexb,ndims,0);
34  int done = 0;
35  while (!done) {
36  int p1, p2;
37  _ArrayCopy1D_(indexb,indexi,ndims);
38  _ArrayAdd1D_(indexi,indexi,boundary->is,ndims);
39  if (face == 1) indexi[dim] = ghosts-1-indexb[dim];
40  else if (face == -1) indexi[dim] = size[dim]-indexb[dim]-1;
41  else return(1);
42  _ArrayIndex1DWO_(ndims,size,indexb,boundary->is,ghosts,p1);
43  _ArrayIndex1D_(ndims,size,indexi,ghosts,p2);
44  _ArrayScaleCopy1D_((phi+nvars*p2),(-1.0),(phi+nvars*p1),nvars);
45  _ArrayIncrementIndex_(ndims,bounds,indexb,done);
46  }
47  }
48  return(0);
49 }
Structure containing the variables and function pointers defining a boundary.
#define _ArraySubtract1D_(x, a, b, size)
#define _ArrayScaleCopy1D_(x, a, y, size)
#define _ArrayAdd1D_(x, a, b, size)
#define _ArrayIndex1D_(N, imax, i, ghost, index)
#define _ArrayIndex1DWO_(N, imax, i, offset, ghost, index)
#define _ArraySetValue_(x, size, value)
#define _ArrayIncrementIndex_(N, imax, i, done)
#define _ArrayCopy1D_(x, y, size)

◆ 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
bBoundary object of type DomainBoundary
mMPI object of type MPIVariables
ndimsNumber of spatial dimensions
nvarsNumber of variables/DoFs per grid point
sizeInteger array with the number of grid points in each spatial dimension
ghostsNumber of ghost points
phiThe solution array on which to apply the boundary condition
waqtCurrent solution time

Definition at line 20 of file BCNoslipWall.c.

30 {
31  DomainBoundary *boundary = (DomainBoundary*) b;
32 
33  int dim = boundary->dim;
34  int face = boundary->face;
35 
36  if (ndims == 2) {
37 
38  /* create a fake physics object */
39  double gamma;
40  gamma = boundary->gamma;
41  double inv_gamma_m1 = 1.0/(gamma-1.0);
42 
43  if (boundary->on_this_proc) {
44  int bounds[ndims], indexb[ndims], indexi[ndims];
45  _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
46  _ArraySetValue_(indexb,ndims,0);
47  int done = 0;
48  while (!done) {
49  int p1, p2;
50  _ArrayCopy1D_(indexb,indexi,ndims);
51  _ArrayAdd1D_(indexi,indexi,boundary->is,ndims);
52  if (face == 1) indexi[dim] = ghosts-1-indexb[dim];
53  else if (face == -1) indexi[dim] = size[dim]-indexb[dim]-1;
54  else return(1);
55  _ArrayIndex1DWO_(ndims,size,indexb,boundary->is,ghosts,p1);
56  _ArrayIndex1D_(ndims,size,indexi,ghosts,p2);
57 
58  /* flow variables in the interior */
59  double rho, uvel, vvel, energy, pressure;
60  double rho_gpt, uvel_gpt, vvel_gpt, energy_gpt, pressure_gpt;
61  _NavierStokes2DGetFlowVar_((phi+nvars*p2),rho,uvel,vvel,energy,pressure,gamma);
62  /* set the ghost point values */
63  rho_gpt = rho;
64  pressure_gpt = pressure;
65  uvel_gpt = 2.0*boundary->FlowVelocity[0] - uvel;
66  vvel_gpt = 2.0*boundary->FlowVelocity[1] - vvel;
67  energy_gpt = inv_gamma_m1*pressure_gpt
68  + 0.5 * rho_gpt * (uvel_gpt*uvel_gpt + vvel_gpt*vvel_gpt);
69 
70  phi[nvars*p1+0] = rho_gpt;
71  phi[nvars*p1+1] = rho_gpt * uvel_gpt;
72  phi[nvars*p1+2] = rho_gpt * vvel_gpt;
73  phi[nvars*p1+3] = energy_gpt;
74 
75  _ArrayIncrementIndex_(ndims,bounds,indexb,done);
76  }
77  }
78 
79  } else if (ndims == 3) {
80 
81  double gamma;
82  gamma = boundary->gamma;
83  double inv_gamma_m1 = 1.0/(gamma-1.0);
84 
85  if (boundary->on_this_proc) {
86  int bounds[ndims], indexb[ndims], indexi[ndims];
87  _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
88  _ArraySetValue_(indexb,ndims,0);
89  int done = 0;
90  while (!done) {
91  int p1, p2;
92  _ArrayCopy1D_(indexb,indexi,ndims);
93  _ArrayAdd1D_(indexi,indexi,boundary->is,ndims);
94  if (face == 1) indexi[dim] = ghosts-1-indexb[dim];
95  else if (face == -1) indexi[dim] = size[dim]-indexb[dim]-1;
96  else return(1);
97  _ArrayIndex1DWO_(ndims,size,indexb,boundary->is,ghosts,p1);
98  _ArrayIndex1D_(ndims,size,indexi,ghosts,p2);
99 
100  /* flow variables in the interior */
101  double rho, uvel, vvel, wvel, energy, pressure;
102  double rho_gpt, uvel_gpt, vvel_gpt, wvel_gpt, energy_gpt, pressure_gpt;
103  _NavierStokes3DGetFlowVar_((phi+nvars*p2),_NavierStokes3D_stride_,rho,uvel,vvel,wvel,energy,pressure,gamma);
104  /* set the ghost point values */
105  rho_gpt = rho;
106  pressure_gpt = pressure;
107  uvel_gpt = 2.0*boundary->FlowVelocity[0] - uvel;
108  vvel_gpt = 2.0*boundary->FlowVelocity[1] - vvel;
109  wvel_gpt = 2.0*boundary->FlowVelocity[2] - wvel;
110  energy_gpt = inv_gamma_m1*pressure_gpt
111  + 0.5 * rho_gpt
112  * (uvel_gpt*uvel_gpt + vvel_gpt*vvel_gpt + wvel_gpt*wvel_gpt);
113 
114  phi[nvars*p1+0] = rho_gpt;
115  phi[nvars*p1+1] = rho_gpt * uvel_gpt;
116  phi[nvars*p1+2] = rho_gpt * vvel_gpt;
117  phi[nvars*p1+3] = rho_gpt * wvel_gpt;
118  phi[nvars*p1+4] = energy_gpt;
119 
120  _ArrayIncrementIndex_(ndims,bounds,indexb,done);
121  }
122  }
123 
124  }
125  return(0);
126 }
#define _NavierStokes2DGetFlowVar_(u, rho, vx, vy, e, P, gamma)
Structure containing the variables and function pointers defining a boundary.
#define _ArraySubtract1D_(x, a, b, size)
#define _ArrayAdd1D_(x, a, b, size)
#define _ArrayIndex1D_(N, imax, i, ghost, index)
#define _ArrayIndex1DWO_(N, imax, i, offset, ghost, index)
static const int _NavierStokes3D_stride_
#define _ArraySetValue_(x, size, value)
#define _ArrayIncrementIndex_(N, imax, i, done)
#define _ArrayCopy1D_(x, y, size)
#define _NavierStokes3DGetFlowVar_(u, stride, rho, vx, vy, vz, e, P, gamma)

◆ 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
bBoundary object of type DomainBoundary
mMPI object of type MPIVariables
ndimsNumber of spatial dimensions
nvarsNumber of variables/DoFs per grid point
sizeInteger array with the number of grid points in each spatial dimension
ghostsNumber of ghost points
phiThe solution array on which to apply the boundary condition
waqtCurrent solution time

Definition at line 27 of file BCThermalNoslipWall.c.

37 {
38  DomainBoundary *boundary = (DomainBoundary*) b;
39 
40  int dim = boundary->dim;
41  int face = boundary->face;
42 
43  if (ndims == 3) {
44 
45  /* create a fake physics object */
46  double gamma;
47  gamma = boundary->gamma;
48  double inv_gamma_m1 = 1.0/(gamma-1.0);
49 
50  if (boundary->on_this_proc) {
51 
52  int *temperature_field_size = boundary->UnsteadyTemperatureSize;
53  int n_time_levels = temperature_field_size[dim];
54  double *time_levels = boundary->UnsteadyTimeLevels;
55  double *temperature_data = boundary->UnsteadyTemperatureData;
56 
57  int it = n_time_levels - 1;
58  while ((time_levels[it] > waqt) && (it > 0)) it--;
59 
60  int bounds[ndims], indexb[ndims], indexi[ndims];
61  _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
62  _ArraySetValue_(indexb,ndims,0);
63 
64  int done = 0;
65  while (!done) {
66 
67  int p1, p2;
68  _ArrayCopy1D_(indexb,indexi,ndims);
69  _ArrayAdd1D_(indexi,indexi,boundary->is,ndims);
70  if (face == 1) indexi[dim] = ghosts-1-indexb[dim];
71  else if (face == -1) indexi[dim] = size[dim]-indexb[dim]-1;
72  else return(1);
73  _ArrayIndex1DWO_(ndims,size,indexb,boundary->is,ghosts,p1);
74  _ArrayIndex1D_(ndims,size,indexi,ghosts,p2);
75 
76  /* get the specified temperature */
77  int index1[ndims]; _ArrayCopy1D_(indexb,index1,ndims);
78  index1[dim] = it;
79  int q; _ArrayIndex1D_(ndims,temperature_field_size,index1,0,q);
80  double temperature_b = temperature_data[q];
81 
82  /* flow variables in the interior */
83  double rho, uvel, vvel, wvel, energy, pressure;
84  _NavierStokes3DGetFlowVar_((phi+nvars*p2),_NavierStokes3D_stride_,rho,uvel,vvel,wvel,energy,pressure,gamma);
85  /* set the ghost point values */
86  double rho_gpt, uvel_gpt, vvel_gpt, wvel_gpt, energy_gpt, pressure_gpt;
87  rho_gpt = rho;
88  uvel_gpt = 2.0*boundary->FlowVelocity[0] - uvel;
89  vvel_gpt = 2.0*boundary->FlowVelocity[1] - vvel;
90  wvel_gpt = 2.0*boundary->FlowVelocity[2] - wvel;
91  pressure_gpt = rho_gpt * temperature_b;
92  energy_gpt = inv_gamma_m1*pressure_gpt
93  + 0.5 * rho_gpt
94  * (uvel_gpt*uvel_gpt + vvel_gpt*vvel_gpt + wvel_gpt*wvel_gpt);
95 
96  phi[nvars*p1+0] = rho_gpt;
97  phi[nvars*p1+1] = rho_gpt * uvel_gpt;
98  phi[nvars*p1+2] = rho_gpt * vvel_gpt;
99  phi[nvars*p1+3] = rho_gpt * wvel_gpt;
100  phi[nvars*p1+4] = energy_gpt;
101 
102  _ArrayIncrementIndex_(ndims,bounds,indexb,done);
103  }
104  }
105 
106  }
107  return(0);
108 }
Structure containing the variables and function pointers defining a boundary.
#define _ArraySubtract1D_(x, a, b, size)
#define _ArrayAdd1D_(x, a, b, size)
double * UnsteadyTimeLevels
#define _ArrayIndex1D_(N, imax, i, ghost, index)
#define _ArrayIndex1DWO_(N, imax, i, offset, ghost, index)
static const int _NavierStokes3D_stride_
#define _ArraySetValue_(x, size, value)
#define _ArrayIncrementIndex_(N, imax, i, done)
double * UnsteadyTemperatureData
#define _ArrayCopy1D_(x, y, size)
#define _NavierStokes3DGetFlowVar_(u, stride, rho, vx, vy, vz, e, P, gamma)

◆ 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
bBoundary object of type DomainBoundary
mMPI object of type MPIVariables
ndimsNumber of spatial dimensions
nvarsNumber of variables/DoFs per grid point
sizeInteger array with the number of grid points in each spatial dimension
ghostsNumber of ghost points
phiThe solution array on which to apply the boundary condition
waqtCurrent solution time

Definition at line 22 of file BCSlipWall.c.

32 {
33  DomainBoundary *boundary = (DomainBoundary*) b;
34 
35  int dim = boundary->dim;
36  int face = boundary->face;
37 
38  if (ndims == 1) {
39 
40  /* create a fake physics object */
41  Euler1D physics;
42  double gamma;
43  gamma = physics.gamma = boundary->gamma;
44  double inv_gamma_m1 = 1.0/(gamma-1.0);
45 
46  if (boundary->on_this_proc) {
47  int bounds[ndims], indexb[ndims], indexi[ndims];
48  _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
49  _ArraySetValue_(indexb,ndims,0);
50  int done = 0;
51  while (!done) {
52  int p1, p2;
53  _ArrayCopy1D_(indexb,indexi,ndims);
54  _ArrayAdd1D_(indexi,indexi,boundary->is,ndims);
55  if (face == 1) indexi[dim] = ghosts-1-indexb[dim];
56  else if (face == -1) indexi[dim] = size[dim]-indexb[dim]-1;
57  else return(1);
58  _ArrayIndex1DWO_(ndims,size,indexb,boundary->is,ghosts,p1);
59  _ArrayIndex1D_(ndims,size,indexi,ghosts,p2);
60 
61  /* flow variables in the interior */
62  double rho, uvel, energy, pressure;
63  double rho_gpt, uvel_gpt, energy_gpt, pressure_gpt;
64  _Euler1DGetFlowVar_((phi+nvars*p2),rho,uvel,energy,pressure,(&physics));
65  /* set the ghost point values */
66  rho_gpt = rho;
67  pressure_gpt = pressure;
68  uvel_gpt = 2.0*boundary->FlowVelocity[_XDIR_] - uvel;
69  energy_gpt = inv_gamma_m1*pressure_gpt
70  + 0.5 * rho_gpt * uvel_gpt*uvel_gpt;
71 
72  phi[nvars*p1+0] = rho_gpt;
73  phi[nvars*p1+1] = rho_gpt * uvel_gpt;
74  phi[nvars*p1+2] = energy_gpt;
75 
76  _ArrayIncrementIndex_(ndims,bounds,indexb,done);
77  }
78  }
79 
80  } else if (ndims == 2) {
81 
82  /* create a fake physics object */
83  Euler2D physics;
84  double gamma;
85  gamma = physics.gamma = boundary->gamma;
86  double inv_gamma_m1 = 1.0/(gamma-1.0);
87 
88  if (boundary->on_this_proc) {
89  int bounds[ndims], indexb[ndims], indexi[ndims];
90  _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
91  _ArraySetValue_(indexb,ndims,0);
92  int done = 0;
93  while (!done) {
94  int p1, p2;
95  _ArrayCopy1D_(indexb,indexi,ndims);
96  _ArrayAdd1D_(indexi,indexi,boundary->is,ndims);
97  if (face == 1) indexi[dim] = ghosts-1-indexb[dim];
98  else if (face == -1) indexi[dim] = size[dim]-indexb[dim]-1;
99  else return(1);
100  _ArrayIndex1DWO_(ndims,size,indexb,boundary->is,ghosts,p1);
101  _ArrayIndex1D_(ndims,size,indexi,ghosts,p2);
102 
103  /* flow variables in the interior */
104  double rho, uvel, vvel, energy, pressure;
105  double rho_gpt, uvel_gpt, vvel_gpt, energy_gpt, pressure_gpt;
106  _Euler2DGetFlowVar_((phi+nvars*p2),rho,uvel,vvel,energy,pressure,(&physics));
107  /* set the ghost point values */
108  rho_gpt = rho;
109  pressure_gpt = pressure;
110  if (dim == _XDIR_) {
111  uvel_gpt = 2.0*boundary->FlowVelocity[_XDIR_] - uvel;
112  vvel_gpt = vvel;
113  } else if (dim == _YDIR_) {
114  uvel_gpt = uvel;
115  vvel_gpt = 2.0*boundary->FlowVelocity[_YDIR_] - vvel;
116  } else {
117  uvel_gpt = 0.0;
118  vvel_gpt = 0.0;
119  }
120  energy_gpt = inv_gamma_m1*pressure_gpt
121  + 0.5 * rho_gpt * (uvel_gpt*uvel_gpt + vvel_gpt*vvel_gpt);
122 
123  phi[nvars*p1+0] = rho_gpt;
124  phi[nvars*p1+1] = rho_gpt * uvel_gpt;
125  phi[nvars*p1+2] = rho_gpt * vvel_gpt;
126  phi[nvars*p1+3] = energy_gpt;
127 
128  _ArrayIncrementIndex_(ndims,bounds,indexb,done);
129  }
130  }
131 
132  } else if (ndims == 3) {
133 
134  /* create a fake physics object */
135  double gamma;
136  gamma = boundary->gamma;
137  double inv_gamma_m1 = 1.0/(gamma-1.0);
138 
139  if (boundary->on_this_proc) {
140  int bounds[ndims], indexb[ndims], indexi[ndims];
141  _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
142  _ArraySetValue_(indexb,ndims,0);
143  int done = 0;
144  while (!done) {
145  int p1, p2;
146  _ArrayCopy1D_(indexb,indexi,ndims);
147  _ArrayAdd1D_(indexi,indexi,boundary->is,ndims);
148  if (face == 1) indexi[dim] = ghosts-1-indexb[dim];
149  else if (face == -1) indexi[dim] = size[dim]-indexb[dim]-1;
150  else return(1);
151  _ArrayIndex1DWO_(ndims,size,indexb,boundary->is,ghosts,p1);
152  _ArrayIndex1D_(ndims,size,indexi,ghosts,p2);
153 
154  /* flow variables in the interior */
155  double rho, uvel, vvel, wvel, energy, pressure;
156  double rho_gpt, uvel_gpt, vvel_gpt, wvel_gpt, energy_gpt, pressure_gpt;
157  _NavierStokes3DGetFlowVar_((phi+nvars*p2),_NavierStokes3D_stride_,rho,uvel,vvel,wvel,energy,pressure,gamma);
158  /* set the ghost point values */
159  rho_gpt = rho;
160  pressure_gpt = pressure;
161  if (dim == _XDIR_) {
162  uvel_gpt = 2.0*boundary->FlowVelocity[_XDIR_] - uvel;
163  vvel_gpt = vvel;
164  wvel_gpt = wvel;
165  } else if (dim == _YDIR_) {
166  uvel_gpt = uvel;
167  vvel_gpt = 2.0*boundary->FlowVelocity[_YDIR_] - vvel;
168  wvel_gpt = wvel;
169  } else if (dim == _ZDIR_) {
170  uvel_gpt = uvel;
171  vvel_gpt = vvel;
172  wvel_gpt = 2.0*boundary->FlowVelocity[_ZDIR_] - wvel;
173  } else {
174  uvel_gpt = 0.0;
175  vvel_gpt = 0.0;
176  wvel_gpt = 0.0;
177  }
178  energy_gpt = inv_gamma_m1*pressure_gpt
179  + 0.5 * rho_gpt
180  * (uvel_gpt*uvel_gpt + vvel_gpt*vvel_gpt + wvel_gpt*wvel_gpt);
181 
182  phi[nvars*p1+0] = rho_gpt;
183  phi[nvars*p1+1] = rho_gpt * uvel_gpt;
184  phi[nvars*p1+2] = rho_gpt * vvel_gpt;
185  phi[nvars*p1+3] = rho_gpt * wvel_gpt;
186  phi[nvars*p1+4] = energy_gpt;
187 
188  _ArrayIncrementIndex_(ndims,bounds,indexb,done);
189  }
190  }
191 
192  }
193  return(0);
194 }
double gamma
Definition: euler1d.h:274
Structure containing variables and parameters specific to the 1D Euler equations. This structure cont...
Definition: euler1d.h:273
#define _ZDIR_
Structure containing the variables and function pointers defining a boundary.
#define _Euler2DGetFlowVar_(u, rho, vx, vy, e, P, p)
Definition: euler2d.h:44
#define _ArraySubtract1D_(x, a, b, size)
#define _ArrayAdd1D_(x, a, b, size)
#define _YDIR_
Definition: euler2d.h:41
double gamma
Definition: euler2d.h:245
#define _ArrayIndex1D_(N, imax, i, ghost, index)
#define _ArrayIndex1DWO_(N, imax, i, offset, ghost, index)
static const int _NavierStokes3D_stride_
#define _ArraySetValue_(x, size, value)
#define _ArrayIncrementIndex_(N, imax, i, done)
#define _XDIR_
Definition: euler1d.h:75
#define _Euler1DGetFlowVar_(u, rho, v, e, P, p)
Definition: euler1d.h:81
#define _ArrayCopy1D_(x, y, size)
#define _NavierStokes3DGetFlowVar_(u, stride, rho, vx, vy, vz, e, P, gamma)

◆ 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
bBoundary object of type DomainBoundary
mMPI object of type MPIVariables
ndimsNumber of spatial dimensions
nvarsNumber of variables/DoFs per grid point
sizeInteger array with the number of grid points in each spatial dimension
ghostsNumber of ghost points
phiThe solution array on which to apply the boundary condition
waqtCurrent solution time

Definition at line 27 of file BCThermalSlipWall.c.

37 {
38  DomainBoundary *boundary = (DomainBoundary*) b;
39 
40  int dim = boundary->dim;
41  int face = boundary->face;
42 
43  if (ndims == 3) {
44 
45  /* create a fake physics object */
46  double gamma;
47  gamma = boundary->gamma;
48  double inv_gamma_m1 = 1.0/(gamma-1.0);
49 
50  if (boundary->on_this_proc) {
51 
52  int *temperature_field_size = boundary->UnsteadyTemperatureSize;
53  int n_time_levels = temperature_field_size[dim];
54  double *time_levels = boundary->UnsteadyTimeLevels;
55  double *temperature_data = boundary->UnsteadyTemperatureData;
56 
57  int it = n_time_levels - 1;
58  while ((time_levels[it] > waqt) && (it > 0)) it--;
59 
60  int bounds[ndims], indexb[ndims], indexi[ndims];
61  _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
62  _ArraySetValue_(indexb,ndims,0);
63 
64  int done = 0;
65  while (!done) {
66 
67  int p1, p2;
68  _ArrayCopy1D_(indexb,indexi,ndims);
69  _ArrayAdd1D_(indexi,indexi,boundary->is,ndims);
70  if (face == 1) indexi[dim] = ghosts-1-indexb[dim];
71  else if (face == -1) indexi[dim] = size[dim]-indexb[dim]-1;
72  else return(1);
73  _ArrayIndex1DWO_(ndims,size,indexb,boundary->is,ghosts,p1);
74  _ArrayIndex1D_(ndims,size,indexi,ghosts,p2);
75 
76  /* get the specified temperature */
77  int index1[ndims]; _ArrayCopy1D_(indexb,index1,ndims);
78  index1[dim] = it;
79  int q; _ArrayIndex1D_(ndims,temperature_field_size,index1,0,q);
80  double temperature_b = temperature_data[q];
81 
82  /* flow variables in the interior */
83  double rho, uvel, vvel, wvel, energy, pressure;
84  _NavierStokes3DGetFlowVar_((phi+nvars*p2),_NavierStokes3D_stride_,rho,uvel,vvel,wvel,energy,pressure,gamma);
85  /* set the ghost point values */
86  double rho_gpt, uvel_gpt, vvel_gpt, wvel_gpt, energy_gpt, pressure_gpt;
87  rho_gpt = rho;
88  if (dim == _XDIR_) {
89  uvel_gpt = 2.0*boundary->FlowVelocity[_XDIR_] - uvel;
90  vvel_gpt = vvel;
91  wvel_gpt = wvel;
92  } else if (dim == _YDIR_) {
93  uvel_gpt = uvel;
94  vvel_gpt = 2.0*boundary->FlowVelocity[_YDIR_] - vvel;
95  wvel_gpt = wvel;
96  } else if (dim == _ZDIR_) {
97  uvel_gpt = uvel;
98  vvel_gpt = vvel;
99  wvel_gpt = 2.0*boundary->FlowVelocity[_ZDIR_] - wvel;
100  } else {
101  uvel_gpt = 0.0;
102  vvel_gpt = 0.0;
103  wvel_gpt = 0.0;
104  }
105  pressure_gpt = rho_gpt * temperature_b;
106  energy_gpt = inv_gamma_m1*pressure_gpt
107  + 0.5 * rho_gpt
108  * (uvel_gpt*uvel_gpt + vvel_gpt*vvel_gpt + wvel_gpt*wvel_gpt);
109 
110  phi[nvars*p1+0] = rho_gpt;
111  phi[nvars*p1+1] = rho_gpt * uvel_gpt;
112  phi[nvars*p1+2] = rho_gpt * vvel_gpt;
113  phi[nvars*p1+3] = rho_gpt * wvel_gpt;
114  phi[nvars*p1+4] = energy_gpt;
115 
116  _ArrayIncrementIndex_(ndims,bounds,indexb,done);
117  }
118  }
119 
120  }
121  return(0);
122 }
#define _ZDIR_
Structure containing the variables and function pointers defining a boundary.
#define _ArraySubtract1D_(x, a, b, size)
#define _ArrayAdd1D_(x, a, b, size)
#define _YDIR_
Definition: euler2d.h:41
double * UnsteadyTimeLevels
#define _ArrayIndex1D_(N, imax, i, ghost, index)
#define _ArrayIndex1DWO_(N, imax, i, offset, ghost, index)
static const int _NavierStokes3D_stride_
#define _ArraySetValue_(x, size, value)
#define _ArrayIncrementIndex_(N, imax, i, done)
#define _XDIR_
Definition: euler1d.h:75
double * UnsteadyTemperatureData
#define _ArrayCopy1D_(x, y, size)
#define _NavierStokes3DGetFlowVar_(u, stride, rho, vx, vy, vz, e, P, gamma)

◆ 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
bBoundary object of type DomainBoundary
mMPI object of type MPIVariables
ndimsNumber of spatial dimensions
nvarsNumber of variables/DoFs per grid point
sizeInteger array with the number of grid points in each spatial dimension
ghostsNumber of ghost points
phiThe solution array on which to apply the boundary condition
waqtCurrent solution time

Definition at line 20 of file BCSubsonicInflow.c.

30 {
31  DomainBoundary *boundary = (DomainBoundary*) b;
32 
33  int dim = boundary->dim;
34  int face = boundary->face;
35 
36  if (ndims == 2) {
37 
38  /* create a fake physics object */
39  Euler2D physics;
40  double gamma;
41  gamma = physics.gamma = boundary->gamma;
42  double inv_gamma_m1 = 1.0/(gamma-1.0);
43 
44  if (boundary->on_this_proc) {
45  int bounds[ndims], indexb[ndims], indexi[ndims];
46  _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
47  _ArraySetValue_(indexb,ndims,0);
48  int done = 0;
49  while (!done) {
50  int p1, p2;
51  _ArrayCopy1D_(indexb,indexi,ndims);
52  _ArrayAdd1D_(indexi,indexi,boundary->is,ndims);
53  if (face == 1) indexi[dim] = ghosts-1-indexb[dim];
54  else if (face == -1) indexi[dim] = size[dim]-indexb[dim]-1;
55  else return(1);
56  _ArrayIndex1DWO_(ndims,size,indexb,boundary->is,ghosts,p1);
57  _ArrayIndex1D_(ndims,size,indexi,ghosts,p2);
58 
59  /* flow variables in the interior */
60  double rho, uvel, vvel, energy, pressure;
61  double rho_gpt, uvel_gpt, vvel_gpt, energy_gpt, pressure_gpt;
62  _Euler2DGetFlowVar_((phi+nvars*p2),rho,uvel,vvel,energy,pressure,(&physics));
63  /* set the ghost point values */
64  rho_gpt = boundary->FlowDensity;
65  pressure_gpt = pressure;
66  uvel_gpt = boundary->FlowVelocity[0];
67  vvel_gpt = boundary->FlowVelocity[1];
68  energy_gpt = inv_gamma_m1*pressure_gpt
69  + 0.5 * rho_gpt * (uvel_gpt*uvel_gpt + vvel_gpt*vvel_gpt);
70 
71  phi[nvars*p1+0] = rho_gpt;
72  phi[nvars*p1+1] = rho_gpt * uvel_gpt;
73  phi[nvars*p1+2] = rho_gpt * vvel_gpt;
74  phi[nvars*p1+3] = energy_gpt;
75 
76  _ArrayIncrementIndex_(ndims,bounds,indexb,done);
77  }
78  }
79 
80  } else if (ndims == 3) {
81 
82  /* create a fake physics object */
83  double gamma;
84  gamma = boundary->gamma;
85  double inv_gamma_m1 = 1.0/(gamma-1.0);
86 
87  if (boundary->on_this_proc) {
88  int bounds[ndims], indexb[ndims], indexi[ndims];
89  _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
90  _ArraySetValue_(indexb,ndims,0);
91  int done = 0;
92  while (!done) {
93  int p1, p2;
94  _ArrayCopy1D_(indexb,indexi,ndims);
95  _ArrayAdd1D_(indexi,indexi,boundary->is,ndims);
96  if (face == 1) indexi[dim] = ghosts-1-indexb[dim];
97  else if (face == -1) indexi[dim] = size[dim]-indexb[dim]-1;
98  else return(1);
99  _ArrayIndex1DWO_(ndims,size,indexb,boundary->is,ghosts,p1);
100  _ArrayIndex1D_(ndims,size,indexi,ghosts,p2);
101 
102  /* flow variables in the interior */
103  double rho, uvel, vvel, wvel, energy, pressure;
104  double rho_gpt, uvel_gpt, vvel_gpt, wvel_gpt, energy_gpt, pressure_gpt;
105  _NavierStokes3DGetFlowVar_((phi+nvars*p2),_NavierStokes3D_stride_,rho,uvel,vvel,wvel,energy,pressure,gamma);
106  /* set the ghost point values */
107  rho_gpt = boundary->FlowDensity;
108  pressure_gpt = pressure;
109  uvel_gpt = boundary->FlowVelocity[0];
110  vvel_gpt = boundary->FlowVelocity[1];
111  wvel_gpt = boundary->FlowVelocity[2];
112  energy_gpt = inv_gamma_m1*pressure_gpt
113  + 0.5 * rho_gpt
114  * (uvel_gpt*uvel_gpt + vvel_gpt*vvel_gpt + wvel_gpt*wvel_gpt);
115 
116  phi[nvars*p1+0] = rho_gpt;
117  phi[nvars*p1+1] = rho_gpt * uvel_gpt;
118  phi[nvars*p1+2] = rho_gpt * vvel_gpt;
119  phi[nvars*p1+3] = rho_gpt * wvel_gpt;
120  phi[nvars*p1+4] = energy_gpt;
121 
122  _ArrayIncrementIndex_(ndims,bounds,indexb,done);
123  }
124  }
125 
126  }
127  return(0);
128 }
Structure containing the variables and function pointers defining a boundary.
#define _Euler2DGetFlowVar_(u, rho, vx, vy, e, P, p)
Definition: euler2d.h:44
#define _ArraySubtract1D_(x, a, b, size)
#define _ArrayAdd1D_(x, a, b, size)
double gamma
Definition: euler2d.h:245
#define _ArrayIndex1D_(N, imax, i, ghost, index)
#define _ArrayIndex1DWO_(N, imax, i, offset, ghost, index)
static const int _NavierStokes3D_stride_
#define _ArraySetValue_(x, size, value)
#define _ArrayIncrementIndex_(N, imax, i, done)
#define _ArrayCopy1D_(x, y, size)
#define _NavierStokes3DGetFlowVar_(u, stride, rho, vx, vy, vz, e, P, gamma)

◆ 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
bBoundary object of type DomainBoundary
mMPI object of type MPIVariables
ndimsNumber of spatial dimensions
nvarsNumber of variables/DoFs per grid point
sizeInteger array with the number of grid points in each spatial dimension
ghostsNumber of ghost points
phiThe solution array on which to apply the boundary condition
waqtCurrent solution time

Definition at line 19 of file BCSubsonicOutflow.c.

29 {
30  DomainBoundary *boundary = (DomainBoundary*) b;
31 
32  int dim = boundary->dim;
33  int face = boundary->face;
34 
35  if (ndims == 2) {
36 
37  /* create a fake physics object */
38  Euler2D physics;
39  double gamma;
40  gamma = physics.gamma = boundary->gamma;
41  double inv_gamma_m1 = 1.0/(gamma-1.0);
42 
43  if (boundary->on_this_proc) {
44  int bounds[ndims], indexb[ndims], indexi[ndims];
45  _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
46  _ArraySetValue_(indexb,ndims,0);
47  int done = 0;
48  while (!done) {
49  int p1, p2;
50  _ArrayCopy1D_(indexb,indexi,ndims);
51  _ArrayAdd1D_(indexi,indexi,boundary->is,ndims);
52  if (face == 1) indexi[dim] = ghosts-1-indexb[dim];
53  else if (face == -1) indexi[dim] = size[dim]-indexb[dim]-1;
54  else return(1);
55  _ArrayIndex1DWO_(ndims,size,indexb,boundary->is,ghosts,p1);
56  _ArrayIndex1D_(ndims,size,indexi,ghosts,p2);
57 
58  /* flow variables in the interior */
59  double rho, uvel, vvel, energy, pressure;
60  double rho_gpt, uvel_gpt, vvel_gpt, energy_gpt, pressure_gpt;
61  _Euler2DGetFlowVar_((phi+nvars*p2),rho,uvel,vvel,energy,pressure,(&physics));
62  /* set the ghost point values */
63  rho_gpt = rho;
64  pressure_gpt = pressure; /* useless statement to avoid compiler warning */
65  pressure_gpt = boundary->FlowPressure;
66  uvel_gpt = uvel;
67  vvel_gpt = vvel;
68  energy_gpt = inv_gamma_m1*pressure_gpt
69  + 0.5 * rho_gpt * (uvel_gpt*uvel_gpt + vvel_gpt*vvel_gpt);
70 
71  phi[nvars*p1+0] = rho_gpt;
72  phi[nvars*p1+1] = rho_gpt * uvel_gpt;
73  phi[nvars*p1+2] = rho_gpt * vvel_gpt;
74  phi[nvars*p1+3] = energy_gpt;
75 
76  _ArrayIncrementIndex_(ndims,bounds,indexb,done);
77  }
78  }
79 
80  } else if (ndims == 3) {
81 
82  /* create a fake physics object */
83  double gamma;
84  gamma = boundary->gamma;
85  double inv_gamma_m1 = 1.0/(gamma-1.0);
86 
87  if (boundary->on_this_proc) {
88  int bounds[ndims], indexb[ndims], indexi[ndims];
89  _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
90  _ArraySetValue_(indexb,ndims,0);
91  int done = 0;
92  while (!done) {
93  int p1, p2;
94  _ArrayCopy1D_(indexb,indexi,ndims);
95  _ArrayAdd1D_(indexi,indexi,boundary->is,ndims);
96  if (face == 1) indexi[dim] = ghosts-1-indexb[dim];
97  else if (face == -1) indexi[dim] = size[dim]-indexb[dim]-1;
98  else return(1);
99  _ArrayIndex1DWO_(ndims,size,indexb,boundary->is,ghosts,p1);
100  _ArrayIndex1D_(ndims,size,indexi,ghosts,p2);
101 
102  /* flow variables in the interior */
103  double rho, uvel, vvel, wvel, energy, pressure;
104  double rho_gpt, uvel_gpt, vvel_gpt, wvel_gpt, energy_gpt, pressure_gpt;
105  _NavierStokes3DGetFlowVar_((phi+nvars*p2),_NavierStokes3D_stride_,rho,uvel,vvel,wvel,energy,pressure,gamma);
106  /* set the ghost point values */
107  rho_gpt = rho;
108  pressure_gpt = pressure; /* useless statement to avoid compiler warning */
109  pressure_gpt = boundary->FlowPressure;
110  uvel_gpt = uvel;
111  vvel_gpt = vvel;
112  wvel_gpt = wvel;
113  energy_gpt = inv_gamma_m1*pressure_gpt
114  + 0.5 * rho_gpt
115  * (uvel_gpt*uvel_gpt + vvel_gpt*vvel_gpt + wvel_gpt*wvel_gpt);
116 
117  phi[nvars*p1+0] = rho_gpt;
118  phi[nvars*p1+1] = rho_gpt * uvel_gpt;
119  phi[nvars*p1+2] = rho_gpt * vvel_gpt;
120  phi[nvars*p1+3] = rho_gpt * wvel_gpt;
121  phi[nvars*p1+4] = energy_gpt;
122 
123  _ArrayIncrementIndex_(ndims,bounds,indexb,done);
124  }
125  }
126 
127  }
128  return(0);
129 }
Structure containing the variables and function pointers defining a boundary.
#define _Euler2DGetFlowVar_(u, rho, vx, vy, e, P, p)
Definition: euler2d.h:44
#define _ArraySubtract1D_(x, a, b, size)
#define _ArrayAdd1D_(x, a, b, size)
double gamma
Definition: euler2d.h:245
#define _ArrayIndex1D_(N, imax, i, ghost, index)
#define _ArrayIndex1DWO_(N, imax, i, offset, ghost, index)
static const int _NavierStokes3D_stride_
#define _ArraySetValue_(x, size, value)
#define _ArrayIncrementIndex_(N, imax, i, done)
#define _ArrayCopy1D_(x, y, size)
#define _NavierStokes3DGetFlowVar_(u, stride, rho, vx, vy, vz, e, P, gamma)

◆ 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
bBoundary object of type DomainBoundary
mMPI object of type MPIVariables
ndimsNumber of spatial dimensions
nvarsNumber of variables/DoFs per grid point
sizeInteger array with the number of grid points in each spatial dimension
ghostsNumber of ghost points
phiThe solution array on which to apply the boundary condition
waqtCurrent solution time

Definition at line 28 of file BCSubsonicAmbivalent.c.

38 {
39  DomainBoundary *boundary = (DomainBoundary*) b;
40 
41  int dim = boundary->dim;
42  int face = boundary->face;
43 
44  if (ndims == 2) {
45 
46  /* create a fake physics object */
47  Euler2D physics;
48  double gamma;
49  gamma = physics.gamma = boundary->gamma;
50  double inv_gamma_m1 = 1.0/(gamma-1.0);
51 
52  /* boundary normal (pointing into the domain) */
53  double nx, ny;
54  if (dim == 0) {
55  nx = 1.0;
56  ny = 0.0;
57  } else if (dim == 1) {
58  nx = 0.0;
59  ny = 1.0;
60  }
61  nx *= (double) face;
62  ny *= (double) face;
63 
64  if (boundary->on_this_proc) {
65  int bounds[ndims], indexb[ndims], indexi[ndims], indexj[ndims];
66  _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
67  _ArraySetValue_(indexb,ndims,0);
68  int done = 0;
69  while (!done) {
70  int p1, p2;
71  double rho, uvel, vvel, energy, pressure;
72 
73  /* compute boundary face velocity - 2nd order */
74  _ArrayCopy1D_(indexb,indexi,ndims);
75  _ArrayAdd1D_(indexi,indexi,boundary->is,ndims);
76  _ArrayCopy1D_(indexi,indexj,ndims);
77  if (face == 1) {
78  indexi[dim] = 0;
79  indexj[dim] = indexi[dim] + 1;
80  } else if (face == -1) {
81  indexi[dim] = size[dim]-1;
82  indexj[dim] = indexi[dim] - 1;
83  }
84  _ArrayIndex1D_(ndims,size,indexi,ghosts,p1);
85  _ArrayIndex1D_(ndims,size,indexj,ghosts,p2);
86  double uvel1, uvel2, uvelb,
87  vvel1, vvel2, vvelb;
88  _Euler2DGetFlowVar_((phi+nvars*p1),rho,uvel1,vvel1,energy,pressure,(&physics));
89  _Euler2DGetFlowVar_((phi+nvars*p2),rho,uvel2,vvel2,energy,pressure,(&physics));
90  uvelb = 1.5*uvel1 - 0.5*uvel2;
91  vvelb = 1.5*vvel1 - 0.5*vvel2;
92  double vel_normal = uvelb*nx + vvelb*ny;
93 
94  _ArrayCopy1D_(indexb,indexi,ndims);
95  _ArrayAdd1D_(indexi,indexi,boundary->is,ndims);
96  if (face == 1) indexi[dim] = ghosts-1-indexb[dim];
97  else if (face == -1) indexi[dim] = size[dim]-indexb[dim]-1;
98  else return(1);
99  _ArrayIndex1DWO_(ndims,size,indexb,boundary->is,ghosts,p1);
100  _ArrayIndex1D_(ndims,size,indexi,ghosts,p2);
101 
102  /* flow variables in the interior */
103  _Euler2DGetFlowVar_((phi+nvars*p2),rho,uvel,vvel,energy,pressure,(&physics));
104 
105  /* set the ghost point values */
106  double rho_gpt, uvel_gpt, vvel_gpt, energy_gpt, pressure_gpt;
107  if (vel_normal > 0) {
108  /* inflow */
109  rho_gpt = boundary->FlowDensity;
110  pressure_gpt = pressure;
111  uvel_gpt = boundary->FlowVelocity[0];
112  vvel_gpt = boundary->FlowVelocity[1];
113  } else {
114  /* outflow */
115  rho_gpt = rho;
116  pressure_gpt = boundary->FlowPressure;
117  uvel_gpt = uvel;
118  vvel_gpt = vvel;
119  }
120  energy_gpt = inv_gamma_m1*pressure_gpt
121  + 0.5 * rho_gpt * (uvel_gpt*uvel_gpt + vvel_gpt*vvel_gpt);
122 
123  phi[nvars*p1+0] = rho_gpt;
124  phi[nvars*p1+1] = rho_gpt * uvel_gpt;
125  phi[nvars*p1+2] = rho_gpt * vvel_gpt;
126  phi[nvars*p1+3] = energy_gpt;
127 
128  _ArrayIncrementIndex_(ndims,bounds,indexb,done);
129  }
130  }
131 
132  } else if (ndims == 3) {
133 
134  /* create a fake physics object */
135  double gamma;
136  gamma = boundary->gamma;
137  double inv_gamma_m1 = 1.0/(gamma-1.0);
138 
139  /* boundary normal (pointing into the domain) */
140  double nx, ny, nz;
141  if (dim == 0) {
142  nx = 1.0;
143  ny = 0.0;
144  nz = 0.0;
145  } else if (dim == 1) {
146  nx = 0.0;
147  ny = 1.0;
148  nz = 0.0;
149  } else if (dim == 2) {
150  nx = 0.0;
151  ny = 0.0;
152  nz = 1.0;
153  }
154  nx *= (double) face;
155  ny *= (double) face;
156  nz *= (double) face;
157 
158  if (boundary->on_this_proc) {
159  int bounds[ndims], indexb[ndims], indexi[ndims], indexj[ndims];
160  _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
161  _ArraySetValue_(indexb,ndims,0);
162  int done = 0;
163  while (!done) {
164  int p1, p2;
165  double rho, uvel, vvel, wvel, energy, pressure;
166 
167  /* compute boundary face velocity - 2nd order */
168  _ArrayCopy1D_(indexb,indexi,ndims);
169  _ArrayAdd1D_(indexi,indexi,boundary->is,ndims);
170  _ArrayCopy1D_(indexi,indexj,ndims);
171  if (face == 1) {
172  indexi[dim] = 0;
173  indexj[dim] = indexi[dim] + 1;
174  } else if (face == -1) {
175  indexi[dim] = size[dim]-1;
176  indexj[dim] = indexi[dim] - 1;
177  }
178  _ArrayIndex1D_(ndims,size,indexi,ghosts,p1);
179  _ArrayIndex1D_(ndims,size,indexj,ghosts,p2);
180  double uvel1, uvel2, uvelb,
181  vvel1, vvel2, vvelb,
182  wvel1, wvel2, wvelb;
183  _NavierStokes3DGetFlowVar_((phi+nvars*p1),_NavierStokes3D_stride_,rho,uvel1,vvel1,wvel1,energy,pressure,gamma);
184  _NavierStokes3DGetFlowVar_((phi+nvars*p2),_NavierStokes3D_stride_,rho,uvel2,vvel2,wvel2,energy,pressure,gamma);
185  uvelb = 1.5*uvel1 - 0.5*uvel2;
186  vvelb = 1.5*vvel1 - 0.5*vvel2;
187  wvelb = 1.5*wvel1 - 0.5*wvel2;
188  double vel_normal = uvelb*nx + vvelb*ny + wvelb*nz;
189 
190  _ArrayCopy1D_(indexb,indexi,ndims);
191  _ArrayAdd1D_(indexi,indexi,boundary->is,ndims);
192  if (face == 1) indexi[dim] = ghosts-1-indexb[dim];
193  else if (face == -1) indexi[dim] = size[dim]-indexb[dim]-1;
194  else return(1);
195  _ArrayIndex1DWO_(ndims,size,indexb,boundary->is,ghosts,p1);
196  _ArrayIndex1D_(ndims,size,indexi,ghosts,p2);
197 
198  /* flow variables in the interior */
199  _NavierStokes3DGetFlowVar_((phi+nvars*p2),_NavierStokes3D_stride_,rho,uvel,vvel,wvel,energy,pressure,gamma);
200 
201  /* set the ghost point values */
202  double rho_gpt, uvel_gpt, vvel_gpt, wvel_gpt, energy_gpt, pressure_gpt;
203  if (vel_normal > 0) {
204  /* inflow */
205  rho_gpt = boundary->FlowDensity;
206  pressure_gpt = pressure;
207  uvel_gpt = boundary->FlowVelocity[0];
208  vvel_gpt = boundary->FlowVelocity[1];
209  wvel_gpt = boundary->FlowVelocity[2];
210  } else {
211  /* outflow */
212  rho_gpt = rho;
213  pressure_gpt = boundary->FlowPressure;
214  uvel_gpt = uvel;
215  vvel_gpt = vvel;
216  wvel_gpt = wvel;
217  }
218  energy_gpt = inv_gamma_m1*pressure_gpt
219  + 0.5 * rho_gpt
220  * (uvel_gpt*uvel_gpt + vvel_gpt*vvel_gpt + wvel_gpt*wvel_gpt);
221 
222  phi[nvars*p1+0] = rho_gpt;
223  phi[nvars*p1+1] = rho_gpt * uvel_gpt;
224  phi[nvars*p1+2] = rho_gpt * vvel_gpt;
225  phi[nvars*p1+3] = rho_gpt * wvel_gpt;
226  phi[nvars*p1+4] = energy_gpt;
227 
228  _ArrayIncrementIndex_(ndims,bounds,indexb,done);
229  }
230  }
231 
232  }
233  return(0);
234 }
Structure containing the variables and function pointers defining a boundary.
#define _Euler2DGetFlowVar_(u, rho, vx, vy, e, P, p)
Definition: euler2d.h:44
#define _ArraySubtract1D_(x, a, b, size)
#define _ArrayAdd1D_(x, a, b, size)
double gamma
Definition: euler2d.h:245
#define _ArrayIndex1D_(N, imax, i, ghost, index)
#define _ArrayIndex1DWO_(N, imax, i, offset, ghost, index)
static const int _NavierStokes3D_stride_
#define _ArraySetValue_(x, size, value)
#define _ArrayIncrementIndex_(N, imax, i, done)
#define _ArrayCopy1D_(x, y, size)
#define _NavierStokes3DGetFlowVar_(u, stride, rho, vx, vy, vz, e, P, gamma)

◆ 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
bBoundary object of type DomainBoundary
mMPI object of type MPIVariables
ndimsNumber of spatial dimensions
nvarsNumber of variables/DoFs per grid point
sizeInteger array with the number of grid points in each spatial dimension
ghostsNumber of ghost points
phiThe solution array on which to apply the boundary condition
waqtCurrent solution time

Definition at line 24 of file BCSupersonicInflow.c.

34 {
35  DomainBoundary *boundary = (DomainBoundary*) b;
36 
37  if (ndims == 2) {
38 
39  double gamma;
40  gamma = boundary->gamma;
41  double inv_gamma_m1 = 1.0/(gamma-1.0);
42 
43  if (boundary->on_this_proc) {
44  int bounds[ndims], indexb[ndims];
45  _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
46  _ArraySetValue_(indexb,ndims,0);
47  int done = 0;
48  while (!done) {
49  int p1; _ArrayIndex1DWO_(ndims,size,indexb,boundary->is,ghosts,p1);
50 
51  /* set the ghost point values */
52  double rho_gpt, uvel_gpt, vvel_gpt, energy_gpt, pressure_gpt;
53  rho_gpt = boundary->FlowDensity;
54  pressure_gpt = boundary->FlowPressure;
55  uvel_gpt = boundary->FlowVelocity[0];
56  vvel_gpt = boundary->FlowVelocity[1];
57  energy_gpt = inv_gamma_m1*pressure_gpt
58  + 0.5 * rho_gpt * (uvel_gpt*uvel_gpt + vvel_gpt*vvel_gpt);
59 
60  phi[nvars*p1+0] = rho_gpt;
61  phi[nvars*p1+1] = rho_gpt * uvel_gpt;
62  phi[nvars*p1+2] = rho_gpt * vvel_gpt;
63  phi[nvars*p1+3] = energy_gpt;
64 
65  _ArrayIncrementIndex_(ndims,bounds,indexb,done);
66  }
67  }
68 
69  } else if (ndims == 3) {
70 
71  double gamma;
72  gamma = boundary->gamma;
73  double inv_gamma_m1 = 1.0/(gamma-1.0);
74 
75  if (boundary->on_this_proc) {
76  int bounds[ndims], indexb[ndims];
77  _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
78  _ArraySetValue_(indexb,ndims,0);
79  int done = 0;
80  while (!done) {
81  int p1; _ArrayIndex1DWO_(ndims,size,indexb,boundary->is,ghosts,p1);
82 
83  /* set the ghost point values */
84  double rho_gpt, uvel_gpt, vvel_gpt, wvel_gpt, energy_gpt, pressure_gpt;
85  rho_gpt = boundary->FlowDensity;
86  pressure_gpt = boundary->FlowPressure;
87  uvel_gpt = boundary->FlowVelocity[0];
88  vvel_gpt = boundary->FlowVelocity[1];
89  wvel_gpt = boundary->FlowVelocity[2];
90  energy_gpt = inv_gamma_m1*pressure_gpt
91  + 0.5 * rho_gpt
92  * (uvel_gpt*uvel_gpt + vvel_gpt*vvel_gpt + wvel_gpt*wvel_gpt);
93 
94  phi[nvars*p1+0] = rho_gpt;
95  phi[nvars*p1+1] = rho_gpt * uvel_gpt;
96  phi[nvars*p1+2] = rho_gpt * vvel_gpt;
97  phi[nvars*p1+3] = rho_gpt * wvel_gpt;
98  phi[nvars*p1+4] = energy_gpt;
99 
100  _ArrayIncrementIndex_(ndims,bounds,indexb,done);
101  }
102  }
103 
104  }
105  return(0);
106 }
Structure containing the variables and function pointers defining a boundary.
#define _ArraySubtract1D_(x, a, b, size)
#define _ArrayIndex1DWO_(N, imax, i, offset, ghost, index)
#define _ArraySetValue_(x, size, value)
#define _ArrayIncrementIndex_(N, imax, i, done)

◆ 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
bBoundary object of type DomainBoundary
mMPI object of type MPIVariables
ndimsNumber of spatial dimensions
nvarsNumber of variables/DoFs per grid point
sizeInteger array with the number of grid points in each spatial dimension
ghostsNumber of ghost points
phiThe solution array on which to apply the boundary condition
waqtCurrent solution time

Definition at line 24 of file BCSupersonicOutflow.c.

34 {
35  DomainBoundary *boundary = (DomainBoundary*) b;
36 
37  int dim = boundary->dim;
38  int face = boundary->face;
39 
40  if (ndims == 2) {
41 
42  /* create a fake physics object */
43  Euler2D physics;
44  double gamma;
45  gamma = physics.gamma = boundary->gamma;
46  double inv_gamma_m1 = 1.0/(gamma-1.0);
47 
48  if (boundary->on_this_proc) {
49  int bounds[ndims], indexb[ndims], indexi[ndims];
50  _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
51  _ArraySetValue_(indexb,ndims,0);
52  int done = 0;
53  while (!done) {
54  int p1, p2;
55  _ArrayCopy1D_(indexb,indexi,ndims);
56  _ArrayAdd1D_(indexi,indexi,boundary->is,ndims);
57  if (face == 1) indexi[dim] = ghosts-1-indexb[dim];
58  else if (face == -1) indexi[dim] = size[dim]-indexb[dim]-1;
59  else return(1);
60  _ArrayIndex1DWO_(ndims,size,indexb,boundary->is,ghosts,p1);
61  _ArrayIndex1D_(ndims,size,indexi,ghosts,p2);
62 
63  /* flow variables in the interior */
64  double rho, uvel, vvel, energy, pressure;
65  double rho_gpt, uvel_gpt, vvel_gpt, energy_gpt, pressure_gpt;
66  _Euler2DGetFlowVar_((phi+nvars*p2),rho,uvel,vvel,energy,pressure,(&physics));
67  /* set the ghost point values */
68  rho_gpt = rho;
69  pressure_gpt = pressure;
70  uvel_gpt = uvel;
71  vvel_gpt = vvel;
72  energy_gpt = inv_gamma_m1*pressure_gpt
73  + 0.5 * rho_gpt * (uvel_gpt*uvel_gpt + vvel_gpt*vvel_gpt);
74 
75  phi[nvars*p1+0] = rho_gpt;
76  phi[nvars*p1+1] = rho_gpt * uvel_gpt;
77  phi[nvars*p1+2] = rho_gpt * vvel_gpt;
78  phi[nvars*p1+3] = energy_gpt;
79 
80  _ArrayIncrementIndex_(ndims,bounds,indexb,done);
81  }
82  }
83 
84  } else if (ndims == 3) {
85 
86  /* create a fake physics object */
87  double gamma;
88  gamma = boundary->gamma;
89  double inv_gamma_m1 = 1.0/(gamma-1.0);
90 
91  if (boundary->on_this_proc) {
92  int bounds[ndims], indexb[ndims], indexi[ndims];
93  _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
94  _ArraySetValue_(indexb,ndims,0);
95  int done = 0;
96  while (!done) {
97  int p1, p2;
98  _ArrayCopy1D_(indexb,indexi,ndims);
99  _ArrayAdd1D_(indexi,indexi,boundary->is,ndims);
100  if (face == 1) indexi[dim] = ghosts-1-indexb[dim];
101  else if (face == -1) indexi[dim] = size[dim]-indexb[dim]-1;
102  else return(1);
103  _ArrayIndex1DWO_(ndims,size,indexb,boundary->is,ghosts,p1);
104  _ArrayIndex1D_(ndims,size,indexi,ghosts,p2);
105 
106  /* flow variables in the interior */
107  double rho, uvel, vvel, wvel, energy, pressure;
108  double rho_gpt, uvel_gpt, vvel_gpt, wvel_gpt, energy_gpt, pressure_gpt;
109  _NavierStokes3DGetFlowVar_((phi+nvars*p2),_NavierStokes3D_stride_,rho,uvel,vvel,wvel,energy,pressure,gamma);
110  /* set the ghost point values */
111  rho_gpt = rho;
112  pressure_gpt = pressure;
113  uvel_gpt = uvel;
114  vvel_gpt = vvel;
115  wvel_gpt = wvel;
116  energy_gpt = inv_gamma_m1*pressure_gpt
117  + 0.5 * rho_gpt
118  * (uvel_gpt*uvel_gpt + vvel_gpt*vvel_gpt + wvel_gpt*wvel_gpt);
119 
120  phi[nvars*p1+0] = rho_gpt;
121  phi[nvars*p1+1] = rho_gpt * uvel_gpt;
122  phi[nvars*p1+2] = rho_gpt * vvel_gpt;
123  phi[nvars*p1+3] = rho_gpt * wvel_gpt;
124  phi[nvars*p1+4] = energy_gpt;
125 
126  _ArrayIncrementIndex_(ndims,bounds,indexb,done);
127  }
128  }
129 
130  }
131  return(0);
132 }
Structure containing the variables and function pointers defining a boundary.
#define _Euler2DGetFlowVar_(u, rho, vx, vy, e, P, p)
Definition: euler2d.h:44
#define _ArraySubtract1D_(x, a, b, size)
#define _ArrayAdd1D_(x, a, b, size)
double gamma
Definition: euler2d.h:245
#define _ArrayIndex1D_(N, imax, i, ghost, index)
#define _ArrayIndex1DWO_(N, imax, i, offset, ghost, index)
static const int _NavierStokes3D_stride_
#define _ArraySetValue_(x, size, value)
#define _ArrayIncrementIndex_(N, imax, i, done)
#define _ArrayCopy1D_(x, y, size)
#define _NavierStokes3DGetFlowVar_(u, stride, rho, vx, vy, vz, e, P, gamma)

◆ 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
bBoundary object of type DomainBoundary
mMPI object of type MPIVariables
ndimsNumber of spatial dimensions
nvarsNumber of variables/DoFs per grid point
sizeInteger array with the number of grid points in each spatial dimension
ghostsNumber of ghost points
phiThe solution array on which to apply the boundary condition
waqtCurrent solution time

Definition at line 21 of file BCTurbulentSupersonicInflow.c.

31 {
32  DomainBoundary *boundary = (DomainBoundary*) b;
33 
34  int dim = boundary->dim;
35 
36  double *inflow_data = boundary->UnsteadyDirichletData;
37  int *inflow_size = boundary->UnsteadyDirichletSize;
38 
39  if (ndims == 3) {
40 
41  /* create a fake physics object */
42  double gamma;
43  gamma = boundary->gamma;
44  double inv_gamma_m1 = 1.0/(gamma-1.0);
45 
46  if (boundary->on_this_proc) {
47  /* the following bit is hardcoded for the inflow data
48  * representing fluctuations in a domain of length 2pi */
49  double xt = boundary->FlowVelocity[dim] * waqt;
50  int N = inflow_size[dim];
51  double L = 2.0 * (4.0*atan(1.0));
52  int it = ((int) ((xt/L) * ((double)N))) % N;
53 
54  int bounds[ndims], indexb[ndims];
55  _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
56  _ArraySetValue_(indexb,ndims,0);
57  int done = 0;
58  while (!done) {
59  int p1; _ArrayIndex1DWO_(ndims,size,indexb,boundary->is,ghosts,p1);
60 
61  /* set the ghost point values - mean flow */
62  double rho_gpt, uvel_gpt, vvel_gpt, wvel_gpt, energy_gpt, pressure_gpt;
63  rho_gpt = boundary->FlowDensity;
64  pressure_gpt = boundary->FlowPressure;
65  uvel_gpt = boundary->FlowVelocity[0];
66  vvel_gpt = boundary->FlowVelocity[1];
67  wvel_gpt = boundary->FlowVelocity[2];
68 
69  /* calculate the turbulent fluctuations */
70  double duvel , dvvel , dwvel ;
71  int index1[ndims]; _ArrayCopy1D_(indexb,index1,ndims);
72  index1[dim] = it;
73  int q; _ArrayIndex1D_(ndims,inflow_size,index1,0,q);
74  duvel = inflow_data[q*nvars+1];
75  dvvel = inflow_data[q*nvars+2];
76  dwvel = inflow_data[q*nvars+3];
77 
78  /* add the turbulent fluctuations to the velocity field */
79  uvel_gpt += duvel;
80  vvel_gpt += dvvel;
81  wvel_gpt += dwvel;
82 
83  /* set the ghost point values */
84  energy_gpt = inv_gamma_m1*pressure_gpt
85  + 0.5 * rho_gpt
86  * (uvel_gpt*uvel_gpt + vvel_gpt*vvel_gpt + wvel_gpt*wvel_gpt);
87  phi[nvars*p1+0] = rho_gpt;
88  phi[nvars*p1+1] = rho_gpt * uvel_gpt;
89  phi[nvars*p1+2] = rho_gpt * vvel_gpt;
90  phi[nvars*p1+3] = rho_gpt * wvel_gpt;
91  phi[nvars*p1+4] = energy_gpt;
92 
93  _ArrayIncrementIndex_(ndims,bounds,indexb,done);
94  }
95  }
96 
97  }
98  return(0);
99 }
Structure containing the variables and function pointers defining a boundary.
#define _ArraySubtract1D_(x, a, b, size)
double * UnsteadyDirichletData
#define _ArrayIndex1D_(N, imax, i, ghost, index)
#define _ArrayIndex1DWO_(N, imax, i, offset, ghost, index)
#define _ArraySetValue_(x, size, value)
#define _ArrayIncrementIndex_(N, imax, i, done)
#define _ArrayCopy1D_(x, y, size)

◆ 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
bBoundary object of type DomainBoundary
mMPI object of type MPIVariables
ndimsNumber of spatial dimensions
nvarsNumber of variables/DoFs per grid point
sizeInteger array with the number of grid points in each spatial dimension
ghostsNumber of ghost points
phiThe solution array on which to apply the boundary condition
waqtCurrent solution time

Definition at line 22 of file BCNoFlux.c.

32 {
33  DomainBoundary *boundary = (DomainBoundary*) b;
34 
35  int dim = boundary->dim;
36  int face = boundary->face;
37 
38  if (boundary->on_this_proc) {
39  int bounds[ndims], indexb[ndims], indexi[ndims];
40  _ArraySubtract1D_ (bounds,boundary->ie,boundary->is,ndims);
41  _ArraySetValue_ (indexb,ndims,0);
42  int done = 0;
43  while (!done) {
44  int p1, p2;
45  _ArrayCopy1D_ (indexb,indexi,ndims);
46  _ArrayAdd1D_ (indexi,indexi,boundary->is,ndims);
47  if (face == 1) indexi[dim] = ghosts-1-indexb[dim];
48  else if (face == -1) indexi[dim] = size[dim]-indexb[dim]-1;
49  else return(1);
50  _ArrayIndex1DWO_ (ndims,size,indexb,boundary->is,ghosts,p1);
51  _ArrayIndex1D_ (ndims,size,indexi,ghosts,p2);
52 
53  if (nvars == 4) {
54  phi[nvars*p1+0] = phi[nvars*p2+0];
55  phi[nvars*p1+1] = (dim == _XDIR_ ? -phi[nvars*p2+1] : phi[nvars*p2+1] );
56  phi[nvars*p1+2] = (dim == _YDIR_ ? -phi[nvars*p2+2] : phi[nvars*p2+2] );
57  phi[nvars*p1+3] = phi[nvars*p2+3];
58  } else if (nvars == 5) {
59  phi[nvars*p1+0] = phi[nvars*p2+0];
60  phi[nvars*p1+1] = (dim == _XDIR_ ? -phi[nvars*p2+1] : phi[nvars*p2+1] );
61  phi[nvars*p1+2] = (dim == _YDIR_ ? -phi[nvars*p2+2] : phi[nvars*p2+2] );
62  phi[nvars*p1+3] = (dim == _ZDIR_ ? -phi[nvars*p2+3] : phi[nvars*p2+3] );
63  phi[nvars*p1+4] = phi[nvars*p2+4];
64  }
65 
66  _ArrayIncrementIndex_(ndims,bounds,indexb,done);
67  }
68  }
69  return(0);
70 }
#define _ZDIR_
Structure containing the variables and function pointers defining a boundary.
#define _ArraySubtract1D_(x, a, b, size)
#define _ArrayAdd1D_(x, a, b, size)
#define _YDIR_
Definition: euler2d.h:41
#define _ArrayIndex1D_(N, imax, i, ghost, index)
#define _ArrayIndex1DWO_(N, imax, i, offset, ghost, index)
#define _ArraySetValue_(x, size, value)
#define _ArrayIncrementIndex_(N, imax, i, done)
#define _XDIR_
Definition: euler1d.h:75
#define _ArrayCopy1D_(x, y, size)

◆ 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
bBoundary object of type DomainBoundary
mMPI object of type MPIVariables
ndimsNumber of spatial dimensions
nvarsNumber of variables/DoFs per grid point
sizeInteger array with the number of grid points in each spatial dimension
ghostsNumber of ghost points
phiThe solution array on which to apply the boundary condition
waqtCurrent solution time

Definition at line 21 of file BCSWSlipWall.c.

31 {
32  DomainBoundary *boundary = (DomainBoundary*) b;
33 
34  int dim = boundary->dim;
35  int face = boundary->face;
36 
37  if (ndims == 1) {
38 
39  if (boundary->on_this_proc) {
40  int bounds[ndims], indexb[ndims], indexi[ndims];
41  _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
42  _ArraySetValue_(indexb,ndims,0);
43  int done = 0;
44  while (!done) {
45  int p1, p2;
46  _ArrayCopy1D_(indexb,indexi,ndims);
47  _ArrayAdd1D_(indexi,indexi,boundary->is,ndims);
48  if (face == 1) indexi[dim] = ghosts-1-indexb[dim];
49  else if (face == -1) indexi[dim] = size[dim]-indexb[dim]-1;
50  else return(1);
51  _ArrayIndex1DWO_(ndims,size,indexb,boundary->is,ghosts,p1);
52  _ArrayIndex1D_(ndims,size,indexi,ghosts,p2);
53 
54  /* flow variables in the interior */
55  double h, uvel;
56  double h_gpt, uvel_gpt;
57  _ShallowWater1DGetFlowVar_((phi+nvars*p2),h,uvel);
58  /* set the ghost point values */
59  h_gpt = h;
60  uvel_gpt = 2.0*boundary->FlowVelocity[_XDIR_] - uvel;
61 
62  phi[nvars*p1+0] = h_gpt;
63  phi[nvars*p1+1] = h_gpt * uvel_gpt;
64 
65  _ArrayIncrementIndex_(ndims,bounds,indexb,done);
66  }
67  }
68 
69  } else if (ndims == 2) {
70 
71  if (boundary->on_this_proc) {
72  int bounds[ndims], indexb[ndims], indexi[ndims];
73  _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
74  _ArraySetValue_(indexb,ndims,0);
75  int done = 0;
76  while (!done) {
77  int p1, p2;
78  _ArrayCopy1D_(indexb,indexi,ndims);
79  _ArrayAdd1D_(indexi,indexi,boundary->is,ndims);
80  if (face == 1) indexi[dim] = ghosts-1-indexb[dim];
81  else if (face == -1) indexi[dim] = size[dim]-indexb[dim]-1;
82  else return(1);
83  _ArrayIndex1DWO_(ndims,size,indexb,boundary->is,ghosts,p1);
84  _ArrayIndex1D_(ndims,size,indexi,ghosts,p2);
85 
86  /* flow variables in the interior */
87  double h, uvel, vvel;
88  double h_gpt, uvel_gpt, vvel_gpt;
89  _ShallowWater2DGetFlowVar_((phi+nvars*p2),h,uvel,vvel);
90  /* set the ghost point values */
91  h_gpt = h;
92  if (dim == _XDIR_) {
93  uvel_gpt = 2.0*boundary->FlowVelocity[_XDIR_] - uvel;
94  vvel_gpt = vvel;
95  } else if (dim == _YDIR_) {
96  uvel_gpt = uvel;
97  vvel_gpt = 2.0*boundary->FlowVelocity[_YDIR_] - vvel;
98  } else {
99  uvel_gpt = 0.0;
100  vvel_gpt = 0.0;
101  }
102 
103  phi[nvars*p1+0] = h_gpt;
104  phi[nvars*p1+1] = h_gpt * uvel_gpt;
105  phi[nvars*p1+2] = h_gpt * vvel_gpt;
106 
107  _ArrayIncrementIndex_(ndims,bounds,indexb,done);
108  }
109  }
110 
111  }
112  return(0);
113 }
Structure containing the variables and function pointers defining a boundary.
#define _ArraySubtract1D_(x, a, b, size)
#define _ArrayAdd1D_(x, a, b, size)
#define _YDIR_
Definition: euler2d.h:41
#define _ArrayIndex1D_(N, imax, i, ghost, index)
#define _ArrayIndex1DWO_(N, imax, i, offset, ghost, index)
#define _ArraySetValue_(x, size, value)
#define _ArrayIncrementIndex_(N, imax, i, done)
#define _ShallowWater1DGetFlowVar_(u, h, v)
#define _XDIR_
Definition: euler1d.h:75
#define _ArrayCopy1D_(x, y, size)
#define _ShallowWater2DGetFlowVar_(u, h, uvel, vvel)

◆ 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
bBoundary object of type DomainBoundary
ndimsNumber of spatial dimensions
nvarsNumber of variables/DoFs per grid point
ghostsNumber of ghost points
sizeInteger array with the number of grid points in each spatial dimension
grid1D array with the spatial coordinates of the grid points, one dimension after the other
uSolution
sourceSource term to which the sponge term is added

Definition at line 26 of file BCSponge.c.

36 {
37  DomainBoundary *boundary = (DomainBoundary*) b;
38  int dim = boundary->dim;
39  int face = boundary->face;
40  double *uref = boundary->SpongeValue;
41  double *xmin = boundary->xmin;
42  double *xmax = boundary->xmax;
43  int v;
44 
45  if (boundary->on_this_proc) {
46  int bounds[ndims], indexb[ndims];
47  _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
48  _ArraySetValue_(indexb,ndims,0);
49  int done = 0;
50  while (!done) {
51  int i = indexb[dim] + boundary->is[dim];
52  double x, xstart, xend;
53  _GetCoordinate_(dim,i,size,ghosts,grid,x);
54  xstart = xmin[dim];
55  xend = xmax[dim];
56  /* calculate sigma */
57  double sigma;
58  if (face > 0) sigma = (x - xstart) / (xend - xstart);
59  else sigma = (x - xend ) / (xstart - xend);
60  /* add to the source term */
61  int p; _ArrayIndex1DWO_(ndims,size,indexb,boundary->is,ghosts,p);
62  for (v=0; v<nvars; v++) source[nvars*p+v] -= (sigma * (u[nvars*p+v]-uref[v]));
63  _ArrayIncrementIndex_(ndims,bounds,indexb,done);
64  }
65  }
66  return(0);
67 }
Structure containing the variables and function pointers defining a boundary.
#define _ArraySubtract1D_(x, a, b, size)
#define _ArrayIndex1DWO_(N, imax, i, offset, ghost, index)
#define _ArraySetValue_(x, size, value)
#define _ArrayIncrementIndex_(N, imax, i, done)
#define _GetCoordinate_(dir, i, dim, ghosts, x, coord)
Definition: basic.h:31

◆ 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
bBoundary object of type DomainBoundary
mMPI object of type MPIVariables
ndimsNumber of spatial dimensions
nvarsNumber of variables/DoFs per grid point
sizeInteger array with the number of grid points in each spatial dimension
ghostsNumber of ghost points
phiThe solution array on which to apply the boundary condition
waqtCurrent solution time

Definition at line 73 of file BCSponge.c.

83 {
84  return(0);
85 }

◆ 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.

148 {
149  DomainBoundary *boundary = (DomainBoundary*) b;
150  MPIVariables *mpi = (MPIVariables*) m;
151 
152  char *filename = boundary->UnsteadyTemperatureFilename;
153  int *temperature_field_size = NULL;
154  double *time_level_data = NULL;
155  double *temperature_field_data = NULL;
156  double *time_buffer = NULL;
157  double *data_buffer = NULL;
158 
159  int dim = boundary->dim;
160  int face= boundary->face;
161  int d;
162 
163  if (!mpi->rank) {
164 
165  printf("Reading boundary temperature data from %s.\n",filename);
166 
167  FILE *in;
168  int ferr;
169 
170  /* calculate the number of processors that sit on this temperature boundary */
171  int nproc = 1;
172  for (d=0; d<ndims; d++) nproc *= mpi->iproc[d]; nproc /= mpi->iproc[dim];
173 
174  in = fopen(filename,"rb");
175  if (!in) {
176  fprintf(stderr,"Error in BCReadTemperatureData(): cannot open boundary temperature data file %s.\n",filename);
177  return(1);
178  }
179 
180  int count = 0;
181  while ((!feof(in)) && (count < nproc)) {
182 
183  int rank[ndims], size[ndims];
184  ferr = fread(rank,sizeof(int),ndims,in);
185  if (ferr != ndims) {
186  fprintf(stderr,"Error in BCReadTemperatureData(): Error (1) in file reading, count %d.\n",count);
187  return(1);
188  }
189  if (rank[dim] != (face > 0 ? 0 : mpi->iproc[dim]-1) ) {
190  fprintf(stderr,"Error in BCReadTemperatureData(): Error (2) in file reading, count %d.\n",count);
191  return(1);
192  }
193  ferr = fread(size,sizeof(int),ndims,in);
194  if (ferr != ndims) {
195  fprintf(stderr,"Error in BCReadTemperatureData(): Error (3) in file reading, count %d.\n",count);
196  return(1);
197  }
198 
199  int n_data = size[dim]; /* number of time levels for which temperature data is provided */
200  time_buffer = (double*) calloc (n_data,sizeof(double));
201  ferr = fread(time_buffer,sizeof(double),n_data,in);
202  if (ferr != n_data) {
203  fprintf(stderr,"Error in BCReadTemperatureData(): Error (5) in file reading, count %d.\n",count);
204  return(1);
205  }
206 
207  int data_size = 1;
208  for (d=0; d<ndims; d++) data_size *= size[d];
209  data_buffer = (double*) calloc (data_size,sizeof(double));
210  ferr = fread(data_buffer,sizeof(double),data_size,in);
211  if (ferr != data_size) {
212  fprintf(stderr,"Error in BCReadTemperatureData(): Error (6) in file reading, count %d.\n",count);
213  return(1);
214  }
215 
216  int rank1D = MPIRank1D(ndims,mpi->iproc,rank);
217 
218  if (!rank1D) {
219 
220  int index[ndims];
221 
222  temperature_field_size = (int*) calloc (ndims, sizeof(int));
223  _ArrayCopy1D_(size,temperature_field_size,ndims);
224 
225  time_level_data = (double*) calloc (size[dim], sizeof(double));
226  _ArrayCopy1D_(time_buffer,time_level_data,size[dim]);
227 
228  temperature_field_data = (double*) calloc (data_size, sizeof(double));
229  ArrayCopynD(ndims,data_buffer,temperature_field_data,size,0,0,index,1);
230 
231  } else {
232 
233 #ifndef serial
234  MPI_Request req[3] = {MPI_REQUEST_NULL,MPI_REQUEST_NULL,MPI_REQUEST_NULL};
235  MPI_Isend(size,ndims,MPI_INT,rank1D,2152,mpi->world,&req[0]);
236  MPI_Isend(time_buffer,size[dim],MPI_DOUBLE,rank1D,2154,mpi->world,&req[2]);
237  MPI_Isend(data_buffer,data_size,MPI_DOUBLE,rank1D,2153,mpi->world,&req[1]);
238  MPI_Status status_arr[3];
239  MPI_Waitall(3,&req[0],status_arr);
240 #else
241  fprintf(stderr,"Error in BCReadTemperatureData(): This is a serial run. Invalid (non-zero) rank read.\n");
242 #endif
243 
244  }
245 
246  free(time_buffer);
247  free(data_buffer);
248  count++;
249  }
250 
251  if (count < nproc) {
252  fprintf(stderr,"Error in BCReadTemperatureData(): missing data in unsteady boundary data file %s.\n",filename);
253  fprintf(stderr,"Error in BCReadTemperatureData(): should contain data for %d processors, ", nproc);
254  fprintf(stderr,"Error in BCReadTemperatureData(): but contains data for %d processors!\n", count);
255  return(1);
256  }
257 
258  fclose(in);
259 
260  } else {
261 
262 #ifndef serial
263  if (mpi->ip[dim] == (face > 0 ? 0 : mpi->iproc[dim]-1) ) {
264 
265  MPI_Request req = MPI_REQUEST_NULL;
266 
267  temperature_field_size = (int*) calloc (ndims,sizeof(int));
268  MPI_Irecv(temperature_field_size,ndims,MPI_INT,0,2152,mpi->world,&req);
269  MPI_Wait(&req,MPI_STATUS_IGNORE);
270 
271  int flag = 1;
272  for (d=0; d<ndims; d++) if ((d != dim) && (temperature_field_size[d] != DomainSize[d])) flag = 0;
273  if (!flag) {
274  fprintf(stderr,"Error in BCReadTemperatureData(): Error (4) (dimension mismatch) in file reading, rank %d.\n",mpi->rank);
275  return(1);
276  }
277 
278  time_level_data = (double*) calloc (temperature_field_size[dim],sizeof(double));
279  MPI_Irecv(time_level_data, temperature_field_size[dim], MPI_DOUBLE,0,2154,mpi->world,&req);
280  MPI_Wait(&req,MPI_STATUS_IGNORE);
281 
282  int data_size = 1;
283  for (d=0; d<ndims; d++) data_size *= temperature_field_size[d];
284  temperature_field_data = (double*) calloc (data_size,sizeof(double));
285  MPI_Irecv(temperature_field_data,data_size,MPI_DOUBLE,0,2153,mpi->world,&req);
286  MPI_Wait(&req,MPI_STATUS_IGNORE);
287 
288  }
289 #else
290  fprintf(stderr,"Error in BCReadTemperatureData(): Serial code should not be here!.\n");
291 #endif
292 
293  }
294 
295  boundary->UnsteadyTemperatureSize = temperature_field_size;
296  boundary->UnsteadyTimeLevels = time_level_data;
297  boundary->UnsteadyTemperatureData = temperature_field_data;
298 
299  return(0);
300 }
char UnsteadyTemperatureFilename[_MAX_STRING_SIZE_]
Structure containing the variables and function pointers defining a boundary.
int MPIRank1D(int, int *, int *)
Definition: MPIRank1D.c:26
double * UnsteadyTimeLevels
MPI_Comm world
INLINE int ArrayCopynD(int, const double *, double *, int *, int, int, int *, int)
Structure of MPI-related variables.
double * UnsteadyTemperatureData
#define _ArrayCopy1D_(x, y, size)

◆ 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.

20 {
21  DomainBoundary *boundary = (DomainBoundary*) b;
22  MPIVariables *mpi = (MPIVariables*) m;
23 
24  char *filename = boundary->UnsteadyDirichletFilename;
25  int *inflow_size = NULL;
26  double *inflow_data = NULL;
27  double *buffer = NULL;
28 
29  int dim = boundary->dim;
30  int face= boundary->face;
31  int d;
32 
33  if (!mpi->rank) {
34 
35  printf("Reading turbulent inflow boundary data from %s.\n",filename);
36 
37  FILE *in;
38  int ferr;
39 
40  /* calculate the number of processors that sit on unsteady boundary */
41  int nproc = 1;
42  for (d=0; d<ndims; d++) nproc *= mpi->iproc[d]; nproc /= mpi->iproc[dim];
43 
44  in = fopen(filename,"rb");
45  if (!in) {
46  fprintf(stderr,"Error in BCReadTurbulentInflowData(): cannot open unsteady boundary data file %s.\n",filename);
47  return(1);
48  }
49  int count = 0;
50  while ((!feof(in)) && (count < nproc)) {
51  int rank[ndims], size[ndims];
52  ferr = fread(rank,sizeof(int),ndims,in);
53  if (ferr != ndims) {
54  fprintf(stderr,"Error in BCReadTurbulentInflowData(): Error (1) in file reading, count %d.\n",count);
55  return(1);
56  }
57  if (rank[dim] != (face > 0 ? 0 : mpi->iproc[dim]-1) ) {
58  fprintf(stderr,"Error in BCReadTurbulentInflowData(): Error (2) in file reading, count %d.\n",count);
59  return(1);
60  }
61  ferr = fread(size,sizeof(int),ndims,in);
62  if (ferr != ndims) {
63  fprintf(stderr,"Error in BCReadTurbulentInflowData(): Error (3) in file reading, count %d.\n",count);
64  return(1);
65  }
66  int flag = 1;
67  for (d=0; d<ndims; d++) if ((d != dim) && (size[d] != DomainSize[d])) flag = 0;
68  if (!flag) {
69  fprintf(stderr,"Error in BCReadTurbulentInflowData(): Error (4) (dimension mismatch) in file reading, count %d.\n",count);
70  return(1);
71  }
72 
73  int data_size = nvars;
74  for (d=0; d<ndims; d++) data_size *= size[d];
75  buffer = (double*) calloc (data_size,sizeof(double));
76  ferr = fread(buffer,sizeof(double),data_size,in);
77  if (ferr != data_size) {
78  fprintf(stderr,"Error in BCReadTurbulentInflowData(): Error (6) in file reading, count %d.\n",count);
79  return(1);
80  }
81 
82  int rank1D = MPIRank1D(ndims,mpi->iproc,rank);
83 
84  if (!rank1D) {
85 
86  int index[ndims];
87  inflow_size = (int*) calloc (ndims, sizeof(int));
88  _ArrayCopy1D_(size,inflow_size,ndims);
89  inflow_data = (double*) calloc (data_size, sizeof(double));
90  ArrayCopynD(ndims,buffer,inflow_data,size,0,0,index,nvars);
91 
92  } else {
93 #ifndef serial
94  MPI_Request req[2] = {MPI_REQUEST_NULL,MPI_REQUEST_NULL};
95  MPI_Isend(size,ndims,MPI_INT,rank1D,2152,mpi->world,&req[0]);
96  MPI_Isend(buffer,data_size,MPI_DOUBLE,rank1D,2153,mpi->world,&req[1]);
97  MPI_Status status_arr[3];
98  MPI_Waitall(2,&req[0],status_arr);
99 #else
100  fprintf(stderr,"Error in BCReadTurbulentInflowData(): This is a serial run. Invalid (non-zero) rank read.\n");
101 #endif
102  }
103 
104  free(buffer);
105  count++;
106  }
107 
108  if (count < nproc) {
109  fprintf(stderr,"Error in BCReadTurbulentInflowData(): missing data in unsteady boundary data file %s.\n",filename);
110  fprintf(stderr,"Error in BCReadTurbulentInflowData(): should contain data for %d processors, ", nproc);
111  fprintf(stderr,"Error in BCReadTurbulentInflowData(): but contains data for %d processors!\n", count);
112  return(1);
113  }
114 
115  fclose(in);
116 
117  } else {
118 #ifndef serial
119  if (mpi->ip[dim] == (face > 0 ? 0 : mpi->iproc[dim]-1) ) {
120  MPI_Request req = MPI_REQUEST_NULL;
121  inflow_size = (int*) calloc (ndims,sizeof(int));
122  MPI_Irecv(inflow_size,ndims,MPI_INT,0,2152,mpi->world,&req);
123  MPI_Wait(&req,MPI_STATUS_IGNORE);
124  int data_size = nvars;
125  for (d=0; d<ndims; d++) data_size *= inflow_size[d];
126  inflow_data = (double*) calloc (data_size,sizeof(double));
127  MPI_Irecv(inflow_data,data_size,MPI_DOUBLE,0,2153,mpi->world,&req);
128  MPI_Wait(&req,MPI_STATUS_IGNORE);
129  }
130 #else
131  fprintf(stderr,"Error in BCReadTurbulentInflowData(): Serial code should not be here!.\n");
132 #endif
133  }
134 
135  boundary->UnsteadyDirichletSize = inflow_size;
136  boundary->UnsteadyDirichletData = inflow_data;
137 
138  return(0);
139 }
Structure containing the variables and function pointers defining a boundary.
int MPIRank1D(int, int *, int *)
Definition: MPIRank1D.c:26
double * UnsteadyDirichletData
char UnsteadyDirichletFilename[_MAX_STRING_SIZE_]
MPI_Comm world
INLINE int ArrayCopynD(int, const double *, double *, int *, int, int, int *, int)
Structure of MPI-related variables.
#define _ArrayCopy1D_(x, y, size)

◆ 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