# Boundary Conditions

PeleC manages boundary conditions in a form consistent with many AMReX codes. Ghost cell data are updated over an AMR level during a `FillPatch`

operation and fluxes are then computed over the entire box without specifically recognizing boundary cells. A generic boundary filler function fills standard boundary condition types that do not require user input, including:

*Interior*- Copy-in-intersect in index space (same as periodic boundary conditions). Periodic boundaries are set in the PeleC inputs file*Symmetry*- All conserved quantities and the tangential momentum component are reflected from interior cells without sign change (REFLECT_EVEN) while the normal component is reflected with a sign change (REFLECT_ODD)*NoSlipWall*- REFLECT_EVEN is applied to all conserved quantities except for both tangential and normal momentum components which are updated using REFLECT_ODD*SlipWall*- SlipWall is identical to Symmetry*FOExtrap*- First-order extrapolation: the value in the ghost-cells are a copy of the last interior cell.

More complex boundary conditions require user input that is prescribed explicitly. Boundaries identified as `UserBC`

or `Hard`

in the inputs will be tagged as `EXT_DIR`

in `pc_hypfill`

. Users will then fill the Dirichlet boundary values, typically by calling the helper function, `bcnormal`

.

Special care should be taken when prescribing subsonic `Inflow`

or an `Outflow`

boundary conditions. It might be tempting to directly impose target values in the boundary filler function (for `Inflow`

), or to perform a simple extrapolation (for `Outflow`

). However, this approach would fail to correctly respect the flow of information along solution characteristics - the system would be ill-posed and would lead to unphysical behavior. In particular, at a subsonic inflow boundary, at a subsonic inlet there is one outgoing characteristic, so one flow variable must be specified using information from inside the domain. Similarly, there is one incoming characteristic at outflow boundaries. The NSCBC method, described below, is the preffered method to account for this, but has not been ported to the all C++ version of PeleC. In the meantime, the recommended strategy for subsonic inflow and outflow boundaries for confined geometries such as nozzles and combustors is as follows:

Subsonic Inflows: Specify the desired temperature, velocity, and composition (if relevant) in the ghost cells. Take the pressure from the domain interior. Based on these values, compute the density, internal energy, and total energy for the ghost cells.

Subsonic Outflows: Specify the desired outlet pressure and extrapolate the other flow quantities. In particular, we recommend following the simple characteristic-based extrapolation proposed by Whitfield and Janus (Three-Dimensional Unsteady Euler Equations Solution Using Flux Vector Splitting. AIAA Paper 84-1552, 1984.) and described in Ch. 8 of Blazek’s textbook (Computational Fluid Dunamics - Principles and Applications). Implementations of this method can be found in the

`bcnormal`

function in`prob.H`

for various test cases, including EB-C10 and EB-ConvergingNozzle.

A detailed analysis comparing various boundary condition strategies and demonstrating their implementation is available for the Converging Nozzle Case case.

Warning

This following is currently deprecated as the GS-NSCBC boundary condition has not been ported from Fortran to C++.

A well-known approach to the subsonic problem is the Navier-Stokes Characteristic Boundary Conditions (NSCBC) strategy, and is described in the paper Poinsot and Lele (1992) JCP. In the method, the hyperbolic structure is decomposed to identify incoming and outgoing waves, given a statement of the “external” state outside the domain, and then to construct a model that gives “desired” behavior at the interface. One issue with direct application of the NSCBC treatment in PeleC is that it is formulated to impose boundary fluxes directly. In PeleC however, the Godunov approach that is implemented makes use of boundary data specified via grow cell values, and reconstructs fluxes at faces when required. Thus, the NSCBC strategy has been reformulated to provide the grow cell data required in PeleC. The strategy, the Ghost-Cells Navier-Stokes Boundary Conditions (GC-NSCBC) method, is described in Motheau et al. (2017) AIAA Journal.

For the characteristics-based boundary condition implementation, the solution is rewritten in terms of one-dimensional hyperbolic wave propagation. The waves leaving the domain are computed numerically, while the waves entering into the domain are provided by a model that is based on a “target state”. With the help of numerical relaxation parameters, the contribution of entering waves can be controlled to afford some freedom at the boundary to “push” toward a target state while also allowing acoustic waves to leave the domain. The approach allows some control to minimize the effects of reflected waves, which would otherwise result from the “hard” imposition of the external conditions. Description of the relevant theory is beyond the scope of this documentation (see the paper Motheau et al. (2017) AIAA Journal, which also contains examples demonstrating why imposing directly target values in ghost-cells does not work as expected, and why the NSCBC theory helps to get a more “desirable” solution).

In order to understand the impact of the GC-NSCBC treatment, we give an example that imposes “hard” values in the ghost-cells to represent external conditions, and uses first-order extrapolation at the outflow boundary. A precomputed 1D flame profile is interpolated onto a uniform PeleC grid. Because the solution has to adapt to the new grid and to the PeleC numerical discretization, it creates an unphysical acoustic bump that moves through the domain as an acoustic disturbance. With “hard” boundary conditions, this disturbance is reflected from the outflow boundary back into the domain, and interacts with the flame upstream. A steady solution to this system would require the propagation of this wave back and forth until numerical diffusion eventually reduces its magnitude below some threshold. With the GC-NSCBC boundary treatment, the acoustic wave simply leaves the computational domain. Often times, the latter is the desired behavior of the code.

No GC-NSCBC treatment, hard values set at the left boundary for the inflow, and first order extrapolation in the right boundary to mimic an outflow. The unphysical reflections of the acoustic wave at boundary can be clearly seen.

With the GC-NSCBC, the spurious acoustic wave simply leaves the domain with no unphysical reflection.

In PeleC, the subroutine `bcnormal`

is used to provide the target state for the GC-NSCBC treatment as well as the numerical parameters used by the GC-NSCBC method to efficiently “damp” the reflected waves. Note the signature and the content of the `bcnormal`

routine:

```
subroutine bcnormal(x,u_int,u_ext,dir,sgn,time,bc_type,bc_params,bc_target)
...
integer, optional, intent(out) :: bc_type
double precision, optional, intent(out) :: bc_params(6)
double precision, optional, intent(out) :: bc_target(5)
...
double precision :: relax_U, relax_V, relax_W, relax_T, beta, sigma_out
integer :: flag_nscbc, which_bc_type
flag_nscbc = 0
! When optional arguments are present, GC-NSCBC is activated
! Generic values are auto-filled for numerical parameters,
! but should be set by the user for each BC
! Note that in the impose_NSCBC_xD.f90 routine, not all parameters are used in same time
if (present(bc_type).and.present(bc_params).and.present(bc_target)) then
flag_nscbc = 1
relax_U = 0.5d0 ! For inflow only, relax parameter for x_velocity
relax_V = 0.5d0 ! For inflow only, relax parameter for y_velocity
relax_W = 0.5d0 ! For inflow only, relax parameter for z_velocity
relax_T = -0.2d0 ! For inflow only, relax parameter for temperature
beta = 1.0d0 ! Control the contribution of transverse terms, here they will be discarded
sigma_out = -0.6d0 ! For outflow only, relax parameter. A negative value means that the local Mach number will be used
which_bc_type = Interior ! This is to ensure that nothing will be done if the user don't set anything
endif
```

When `bc_type`

, `bc_params`

and `bc_target`

parameters are present, the routine is likely being called from `impose_NSCBC_(dir)d.F90`

. In this case the flag `flag_nscbc`

is activated to fill optional arrays with the requisite data. Note however that the `FillPatch`

operation called in the AMReX framework also calls `pc_hypfill`

, which then also calls `bcnormal`

. In this case, the GC-NSCBC parameters are not directly relevant. In order to make `bc_normal`

sufficiently generic for both purposes, only the target state is returned to `pc_hypfill`

and the parameters associated to the GC-NSCBC method are ignored. By default, the GC-NSCBC method is activated for all subsonic flow boundaries. It can be turned off by setting the flags `nscbc_adv`

and `nscbc_diff`

to zero. In that case, the ghost-cells will be filled directly with the target state (although, as mentioned, this will likely lead to undesired behavior in the solution!).

The use of `bc_type`

, `bc_params`

and `bc_target`

will be described in detail in other sections of this documentation, but let us focus here on the parameter, `bc_type`

. The `bc_type`

(an integer) is a coded form of the physical boundary condition that we want to impose, and this is done point-wise. This means that along a face of the domain, different physical boundary conditions
can be combined. For example, one may wish to impose an inflow in the middle of a wall in order to represent a localized inlet jet or an open boundary. Four physical boundary conditions are implemented in the GC-NSCBC framework: `Inflow`

, `Outflow`

, `SlipWall`

, `NoSlipWall`

.

`Inflow`

and `Outflow`

conditions rely on different models for the waves entering into the domain, and are computed in the routine `compute_waves`

in `impose_NSCBC_(dir)d.F90`

.
For example in 2D, `Inflow`

requires models for three incoming waves. Thus, three relaxation parameters are needed: `relax_U`

, `relax_V`

and `relax_T`

. Also, three state target
values are needed: `TARGET_VX`

, `TARGET_VY`

and `TARGET_TEMPERATURE`

. For an `Outflow`

, only one wave is leaving the domain, so only `TARGET_PRESSURE`

is needed, and
the relaxation parameter is controlled with `sigma_out`

. Note that transverse terms can be included in the computation of the waves, and the amount of contribution is controlled
by the parameter `beta`

, with values between 0 (full contribution) and 1 (no contribution). A negative input value of `beta`

indicates that `beta`

will be adjusted dynamically with the Mach number of the local flow (see Motheau et al. (2017) AIAA Journal and other references therein for details).

The `impose_NSCBC_(dir)d.F90`

routine is organized as follows:

First, data in ghost-cells along the direction at corners are treated. This is because we have to use a one-sided derivative to compute transverse terms at corners.

For each cell, we compute derivatives in the normal and tangential directions of the face.

We call bcnormal to get: the physical boundary (

`bc_type`

), the target state values (`bc_target`

), and the associated numerical parameters (`bc_params`

).Then we compute the NSCBC waves.

The last step is GC-NSCBC procedure to recompute the values in ghost-cells according to the characteristic waves that have been computed in the previous step.

This procedure is done for each face of the domain.

Below is an example to achieve an inflow/outflow along the x-axis of a channel, periodic in y. Note how the `bc_params`

and `bc_target`

arrays are constructed at the end of the routine.

```
subroutine bcnormal(x,u_int,u_ext,dir,sgn,time,bc_type,bc_params,bc_target)
use probdata_module
use eos_type_module
use eos_module
use meth_params_module, only : URHO, UMX, UMY, UMZ, UTEMP, UEDEN, UEINT, UFS
use network, only: nspecies, naux, molec_wt
use prob_params_module, only : Interior, Inflow, Outflow, SlipWall, NoSlipWall, &
problo, probhi
use bl_constants_module, only: M_PI
implicit none
double precision :: x(3), time
double precision :: u_int(*),u_ext(*)
integer :: dir,sgn
integer, optional, intent(out) :: bc_type
double precision, optional, intent(out) :: bc_params(6)
double precision, optional, intent(out) :: bc_target(5)
type (eos_t) :: eos_state
double precision :: u(3)
double precision :: y
double precision :: relax_U, relax_V, relax_W, relax_T, beta, sigma_out
integer :: flag_nscbc, which_bc_type
flag_nscbc = 0
! When optional arguments are present, GC-NSCBC is activated
! Generic values are auto-filled for numerical parameters,
! but should be set by the user for each BC
! Note that in the impose_NSCBC_xD.f90 routine, not all parameters are used in same time
if (present(bc_type).and.present(bc_params).and.present(bc_target)) then
flag_nscbc = 1
relax_U = 0.5d0 ! For inflow only, relax parameter for x_velocity
relax_V = 0.5d0 ! For inflow only, relax parameter for y_velocity
relax_W = 0.5d0 ! For inflow only, relax parameter for z_velocity
relax_T = 0.2d0 ! For inflow only, relax parameter for temperature
beta = 0.2d0 ! Control the contribution of transverse terms
sigma_out = 0.25d0 ! For outflow only, relax parameter
which_bc_type = Interior ! This is to ensure that nothing will be done if the user don't set anything
endif
call build(eos_state)
! at low X
if (dir == 1) then
if (sgn == 1) then
relax_U = 10.0d0
relax_V = 2.0d0
relax_T = - relax_V
beta = 0.6d0
which_bc_type = Inflow
u(1) = u_ref
u(2) = 0.0d0
u(3) = 0.0d0
eos_state % massfrac(1) = 1.d0
eos_state % p = p_ref
eos_state % T = T_ref
call eos_tp(eos_state)
end if
! at hi X
if (sgn == -1) then
! Set outflow pressure
which_bc_type = Outflow
sigma_out = 0.28d0
beta = -0.60d0
u(1:3) = 0.d0
eos_state % massfrac(1) = 1.d0
eos_state % p = p_ref
eos_state % T = T_ref
call eos_tp(eos_state)
end if
end if
! at low Y
if (dir == 2) then
if (sgn == 1) then
! Do nothing, this is periodic
end if
! at hi Y
if (sgn == -1) then
! Do nothing, this is periodic
end if
end if
u_ext(UFS:UFS+nspecies-1) = eos_state % massfrac * eos_state % rho
u_ext(URHO) = eos_state % rho
u_ext(UMX) = eos_state % rho * u(1)
u_ext(UMY) = eos_state % rho * u(2)
u_ext(UMZ) = eos_state % rho * u(3)
u_ext(UTEMP) = eos_state % T
u_ext(UEINT) = eos_state % rho * eos_state % e
u_ext(UEDEN) = eos_state % rho * (eos_state % e + 0.5d0 * (u(1)**2 + u(2)**2) + u(3)**2)
! Here the optional parameters are filled by the local variables if they were present
if (flag_nscbc == 1) then
bc_type = which_bc_type
bc_params(1) = relax_T! For inflow only, relax parameter for temperature
bc_params(2) = relax_U ! For inflow only, relax parameter for x_velocity
bc_params(3) = relax_V ! For inflow only, relax parameter for y_velocity
bc_params(4) = relax_W ! For inflow only, relax parameter for z_velocity
bc_params(5) = beta ! Control the contribution of transverse terms.
bc_params(6) = sigma_out ! For outflow only, relax parameter
bc_target(1) = U_ext(UMX)/U_ext(URHO) ! Target for Inflow
bc_target(2) = U_ext(UMY)/U_ext(URHO) ! Target for Inflow
bc_target(3) = U_ext(UMZ)/U_ext(URHO) ! Target for Inflow
bc_target(4) = U_ext(UTEMP) ! Target for Inflow
bc_target(5) = eos_state%p ! Target for Outflow
end if
call destroy(eos_state)
end subroutine bcnormal
```

The choice of the relaxation parameters in `bc_params`

is case-dependent, unfortunately. Some trial-and-error is often required to find the best values. However, we suggest the the following based on literature and practical experience:

`relax_U`

,`relax_V`

and`relax_W`

should have values near 0.2. Higher values will impose the velocity more “strongly”, but will likely lead to more unphysical waves reflection.`relax_T`

must be a negative value, typically near -0.2.For outflow boundaries,

`sigma_out`

= 0.25 is often reported to be a good choice.The

`beta`

must be between 0 and 1; it controls the contribution of transverse terms. The choice for this parameter is more complicated. For outflows, it should be close to the Mach number. For some cases, a spatially averaged Mach number will provide good results, while for other cases, the point-wise local Mach number is better.`beta`

will be set to the local Mach number if it is set to a negative value in the inputs. For inflows, it has been found that a value of 0.5 provides good results, but it may lead to instabilities, and for some case turning off the transverse terms (beta=1) will be better.