Boundary Condition Datafile Guide

Overview

The boundary condition datafile consists of two parts. First, the spatial segmentation of the boundaries is specified. Then, the actual boundary conditions are stated for all of the specified segments. They include information on the boundary condition type and the associated value, as well as the time from which on the stated condition is applied (time stamp).

Note

The semantics of the datafile changes if boundary.interpolateBCvalues = true is invoked in the parameter file.

Boundary Condition Types

DORiE currently supports four types of boundary conditions:

Dirichlet

dirichlet <head>

Set a fixed matric head [m] at the boundary.

Neumann

neumann <flux>

Set a fixed volumetric flux [m/s] at the boundary. Positive values imply fluxes out of the domain, negative values imply fluxes into the domain.

Limited Influx

limited_influx <flux>

Set a volumetric flux [m/s] at the boundary (Neumann BC). Fluxes into the domain are stopped as soon as the matric head at the respective boundary segment reaches \(h_m \geq 0.0 \, \text{m}\).

Evaporation

evaporation <flux> <head>

Set a volumetric flux [m/s] (Neumann BC) and a fallback matric head [m] (Dirichlet BC) at the boundary. The Neumann BC is invoked as long as the water content inside the soil can sustain the flux. Before every time step, the numeric flux caused by the given Dirichlet BC is estimated. If it is lower than the given Neumann flux, the Dirichlet BC is applied.

Datafile Structure

Boundary Segmentation

spatial_resolution_<position> <value_count> <values>...

In 2D, the spatial segmentation is fairly simple. The boundary is split into four parts, North, South, West and East. For each of these parts, we specify how many segmentation values will be stated, followed by the values themselves. The segmentation is applied along the respective spatial base vector (e.g. from West to East for the Northern and Southern boundaries).

Note

value_count = 0 means that the respective boundary is one segment for which an explicit BC has to be specified.

To split a northern boundary with a length of 3 m into three equidistant segments, we need to define

spatial_resolution_north 2 1.0 2.0
No Flow Boundary

To apply a Neumann BC with flux \(j=0 \, \text{ms}^{-1}\) to a boundary without any segmentation, set value_count = -1. This boundary will not be considered in the Boundary Condition Specification.

In 3D, every boundary is defined by two adjacent borders. By segmenting the borders, the boundary is split up into rectangles, for each of which a boundary condition has to be specified. To split a northern boundary of a cubic domain with an edge length of 2 m into four squares, we need to define

spatial_resolution_north_fb 1 1.0
spatial_resolution_north_we 1 1.0

No Flow Boundaries are possible here as well. Notice that this definition has to be consistent for every boundary. DORiE will throw an error if you try to segment a border of one boundary if the other border is defined as No Flow Boundary.

Number of Time Stamps

number_BC_change_times <int>

Specify the number of time stamps for which BCs will be defined.

Note

Every line is counted as one time stamp, even if two subsequent lines define BCs for the same time, as in the case of BC value interpolation.

Boundary Condition Specification

The boundary conditions have to be defined for every segment and in every time stamp. Boundary conditions and their respective values are separated by single whitespaces. All BC specifications for one time stamp have to be stated in the same line. BCs between adjacent borders are not interpolated and thus not steady. These lines follow a simple grammar:

bc_line  ::=  time { group }*
time     ::=  float
group    ::=  ( bc_type1 float ) | ( bc_type2 float float )
bc_type1 ::=  "neumann" | "dirichlet" | "limited_influx"
bc_type2 ::=  "evaporation"

The boundary conditions defined here are parsed in the same order as the boundary segments have been specified. In 3D, the rectangular boundary segments are parsed in a tabular fashion, where columns run faster than rows. Columns are defined along the first direction specified in the Boundary Segmentation, and rows are defined along the second direction.

_images/quickstart-bcfile-segm1.png

Spatial segmentation of the northern boundary for the example stated in Boundary Segmentation (top view). x points towards the east-west direction and y towards the front-back direction. Red numbers indicate the order in which the boundary conditions for the respective segment have to be stated.

If boundary.interpolateBCvalues = true is invoked in the parameter file, subsequent BCs of the same type will be linearly interpolated. In order to change the BC type for a certain boundary segment, the respective time stamp has to be stated twice:

0 neumann -1e-6 [...]
2000 neumann -5e-6 [...]
2000 dirichlet 0.0 [...]

This will linearly interpolate the Neumann BC values for \(0 \leq t < 2000\) between \(j(t=0) = -1 \cdot 10^{-6} \, \text{ms}^{-1}\) and \(j(t=2000) = -5 \cdot 10^{-6} \, \text{ms}^{-1}\), and switch the BC type to Dirichlet at \(t=2000\). DORiE will throw an error if you choose subsequent time stamps which cannot be interpolated. Notice that the second line is syntactically valid for disabled interpolation as well, although it has no effect then.

Example Files

Simple Boundary Condition datafiles can be created with the command dorie create.

Exemplary BC File for Infiltration (2D)

spatial_resolution_north 0
spatial_resolution_south 0
spatial_resolution_west -1
spatial_resolution_east -1
number_BC_change_times 1
0 neumann -5.55e-6 dirichlet 0

Exemplary BC File for Infiltration (3D)

spatial_resolution_north_we 0
spatial_resolution_north_fb 0
spatial_resolution_south_we 0
spatial_resolution_south_fb 0
spatial_resolution_west_sn -1
spatial_resolution_west_fb -1
spatial_resolution_east_sn -1
spatial_resolution_east_fb -1
spatial_resolution_front_sn -1
spatial_resolution_front_we -1
spatial_resolution_back_sn -1
spatial_resolution_back_we -1
number_BC_change_times 1
0 neumann -5.55e-6 dirichlet 0

Code Documentation

The readout of the boundary condition datafile and the handling of the associated data are managed by a single class.

template <typename Traits>
class RectangularGrid : public Dune::Dorie::BCReadoutInterface<Traits>

Boundary Condition Input Adapter for Rectangular Grid input.

This class inherits the BCReadoutInterface class and provides the input interface for a rectangular grid input: The boundary is segmented into rectangles. Different BCs and condition values can be specified for every rectangle. The BCs are not steady over the borders of the rectangles.

Public Functions

RectangularGrid(const ParameterTree &config)

Read necessary parameters. Then read out BC data file with rectangular boundary decomposition.

Exceptions
  • Dune::IOError: BC File path not specified

RF getNextTimeStamp(const RF &time) const

Returns the time value of the next BC change. Returns -1 if there is none.

BoundaryCondition::Type getBCtype(const Domain &xGlobal, const RF &time) const

Return Boundary Condition Type at given global position and time.

  • Determine boundary side
  • Check for noFlowBoundary
  • Determine boundary part index
  • Determine time index
  • Read out BC data
    See
    Dune::Dorie::RectangularGrid::bcSectionData
    Return
    BoundaryCondition::Type
    Parameters
    • xGlobal: Global coordinates
    • time: Current time

RF getHeadValue(const Domain &xGlobal, const RF &time) const

Return Matric Head Value at given global position and time.

  • Determine boundary side
  • Determine boundary part index
  • Determine time index
  • Interpolate (if enabled)
  • Read out BC data
    See
    Dune::Dorie::RectangularGrid::bcSectionData
    Return
    Matric Head Value
    Parameters
    • xGlobal: Global coordinates
    • time: Current time
    Exceptions
    • Dune::Exception: Boundary side not determined successfully
    • Dune::Exception: Dirichlet BC is queried at no-flow boundary

RF getFluxValue(const Domain &xGlobal, const RF &time) const

Return Flux Value at given global position and time.

  • Determine boundary side
  • Determine boundary part index
  • Determine time index
  • Interpolate (if enabled)
  • Read out BC data
    See
    Dune::Dorie::RectangularGrid::bcSectionData
    Return
    Flux Value
    Parameters
    • xGlobal: Global coordinates
    • time: Current time
    Exceptions
    • Dune::Exception: Boundary side not determined successfully

Private Types

enum [anonymous]

Values:

dim = Traits::dim
typedef Traits::RangeField RF
typedef Traits::Domain Domain
typedef Traits::DomainField DF
typedef Dune::Dorie::BCReadoutInterface<Traits>::BC BC

Private Functions

RF interpolateValues(const RF &y0, const RF &y1, const RF &t0, const RF &t1, const RF &t) const

Linear interpolation between the points (t0,y0) and (t1,y1), given t.

unsigned int getBCtime(const RF &time) const

Returns time stamp index according to real time value.

BC::Side getBCside(const Domain &xGlobal) const

Returns boundary index as enum (can be interpreted as int), given the global Coordinate.

int getBCindex(const Domain &xGlobal, const int &bcSide) const

Returns the BC index, given the boundary side and the coordinate.

For two dimensions we simply iterate through one spatial resolution array until a resolution value is larger than the coordinate. For three dimensions, we need to do this for two spatial resolution arrays. Afterwards, we have to calculate the index. The first vector of spatial resolution always runs fastest.

Exceptions
  • Dune::Exception: Position query exceeds last segment of spatial resolution (This should not happen if specified extensions are >= true grid extensions)

void readBCfile()

Handle readout of the BC file specified by the user.

  • Read spatial segmentation
  • Read number of time stamps
  • Read time stamps and BCs
    Exceptions
    • Dune::IOError: BC data file cannot be opened
    • Dune::IOError: BC data file has more lines than expected

void readSpatialResolution(std::istringstream &instream, const int &counter)

Read the spatial segmentation scheme.

Exceptions
  • Dune::IOError: Previously defined No Flow Boundary receives segmentation
  • Dune::IOError: Previously segmented boundary receives No Flow Boundary
  • Dune::IOError: Segmentation does not match number of segments
  • Dune::IOError: Spatial resolution values are not stated in increasing order
  • Dune::IOError: Number of segments is integer < -1
  • Dune::IOError: Spatial resolution keywords do not match the mandatory order

unsigned int getNumberTimeStamps(std::istringstream &instream, const int &counter) const

Read out the number of time stamps for BC changes (needs to be called only once)

Exceptions
  • Dune::IOError: Number of time stamps is stated incorrectly

void readBCresolution(std::istringstream &instream, const int &counter)

Read out the time resolution.

Exceptions
  • Dune::IOError: First time stamp does not coincide with time.start key
  • Dune::IOError: Time stamps are not stated in increasing order
  • Dune::IOError: Number of BC declarations does not match spatial segmentation
  • Dune::IOError: BC type is unknown
  • Dune::IOError: Activated interpolation only: Time stamps are not the same if any BC changes in type

void finalizeSpatialResolution()

Add the extension of the grid at the end of every segmentation entry.

void initBCdata()

Initialize the bcData matrix.

void userInputCheck(const unsigned int &changeTimes) const

Make sure saved BC data makes sense.

Exceptions
  • Dune::IOError: Number of saved timestamps does not match number_BC_change_times key
  • Dune::IOError: Last time stamp is larger than time.end key
  • Dune::IOError: BC arrays have incorrect size (i.e. size /= number of time stamps)
  • Dune::IOError: Segmentation does not match saved number of BC indices

void plotBCdata() const

Write out saved BC data matrix for Debugging.

Private Members

std::vector<RF> bcTimeStamps

Vector containig all times at which BCs change.

const Domain extensions

spatial dimensions of the domain

const DF width = extensions[0]

Domain width.

const DF height = extensions[dim-1]

Domain height. This must be set to the maximum height of the domain!

const DF depth = extensions[1]

Domain depth. depth=height for dim=2.

const std::string bcFilePath

Path to file containing BC data.

const RF tStart

Starting time of simulation.

const RF tEnd

Finish time of simulation.

const bool interpolate

Boolean whether consecutive BC values of the same type will be interpolated.

const unsigned int bcSides = (dim==2 ? 4 : 6)

Number of boundaries. Only necessary for setup.

const unsigned int headerLines = (dim==2 ? 4 : 12)

Number of lines in the BC data file containing the spatial resolution. Only necessary for setup.

std::array<std::vector<DF>, (dim==2 ? 4 : 12)> spatialResolution

Array containing the vectors with entries of spatial resolution. Row(2D:4 - 3D:12): BC border. Entries (col): Segmentation.

std::array<unsigned int, (dim==2 ? 4 : 6)> bcIndices

Array with the number of BC indices for every side.

std::array<bool, (dim==2 ? 4 : 6)> noFlowBoundary

Array stating if a boundary has no segmentation and no user specified BC values. Boundary will be set to j=0 (no flow)

std::array<std::vector<bcSectionData>, (dim==2 ? 4 : 6)> bcData

Matrix containing all BC information. Row: BC side. Col: BC index on this side.

const RF eps

Error margin for spatial query.

struct bcSectionData

Container for BC data of a single segment (2D) or rectangle (3D)

Every vector should have as many entries as there are changes in the BCs, i.e. time stamps. If the BC type does not account for a respective value, it is set to 0.0. This means that, for example, a flux value could be called for a Dirichlet BC. This has to be avoided at the value query itself.

For the special case of a No-Flow-Boundary, the arrays have the size 1. Which boundaries are declared as such is saved in the noFlowBoundary variable.

See
initBCdata()
See
readBCresolution()

Public Members

template<>
std::vector<BoundaryCondition::Type> type

Type of Boundary Condition.

See
Dune::Dorie::BoundaryCondition

template<>
std::vector<RF> fluxValue

Value of flux.

template<>
std::vector<RF> headValue

Value of matric head.