Native XML Input Specification for Amanzi-U, Version 1.6-dev

OVERVIEW

This is a continuously evolving specification format used by the code developers. Its main purpose is to develop and test new capabilities without disruption of end-users.

PARAMETER LIST

The Amanzi input file is an ASCII text XML-formatted file that must be framed at the beginning and end by the following statements:

<ParameterList name="transport">
  various parameters and sublists
</ParameterList>

The value of name can be anything (transport in this example). A ParameterList consists of just two types of entries: Parameter and ParameterList. ParameterLists are labeled with name [string], while Parameters have a separate fields called name [string], type [string] and value [TYPE], where TYPE can be any of the following: double, int, bool, string, Array(double), Array(int), and Array(string). The value of the parameter is given in quotes (e.g. value=”2.7e3”). Array data is specified as a single comma-delimited string bounded by brackets (e.g. value=”{linear, constant, linear}”).

<ParameterList name="transport">
  <Parameter name="cfl" type="double" value="0.9"/>
  <Parameter name="ratio" type="Array(int)" value="{2, 1, 4}"/>
</ParameterList>

In this example, the list transport has parameter cfl that is the double with value 0.9, and parameter ratio that is the integer array such that ratio[0] = 2, ratio[1]=1, and ratio[2]=4.

Syntax of the specification

Input specification for each ParameterList entry consists of two parts. First, a bulleted list defines the usage syntax and available options. This is followed by example snippets of XML code to demonstrate usage.

In many cases, the input specifies data for a particular parameterized model, and Amanzi supports a number of parameterizations. For example, initial data might be uniform (the value is required), or linear in y (the value and its gradient are required). Where Amanzi supports a number of parameterized models for parameter model, the available models will be listed by name, and then will be described in the subsequent section. In the manufactured example below, the specification looks as follows:

• SOIL [list] accepts parameters that describes properties of this soil.

  • “region” [string] defines a subdomain of the computational domain.

  • “model” [list] specifies a model for the soil. Available options are “van Genuchten” and “Brooks-Corey”.

Here SOIL is defined by a region and a model. The region is a string parameter but the model is given by a sublist with its own set of parameters. The parameter for model can be described in the same section or in a separate section of this document. For instance, the local description may look like:

• “model” [list] specifies a model for the soil. Available options are “van Genuchten” and “Brooks-Corey”. The option “van Genuchten” requires “m” [double]. The option “Brooks-Corey” requires “lambda” [double] and “alpha” [double].

Each part of the spec is illustrated by an example followed by optional comments:

<ParameterList name="water retention models">
  <ParameterList name="SOIL">
    <Parameter name="region" type="string" value="TOP_DOMAIN"/>
    <ParameterList name="Brooks-Corey">
      <Parameter name="lambda" type="double" value="0.7"/>
      <Parameter name="alpha" type="double" value="1e-3"/>
    </ParameterList>
  </ParameterList>
</ParameterList>

This defines soil properties in region TOP_DOMAIN using the Brooks-Corey model with parameters lambda=0.7 and alpha=1e-3.

Additional conventions:

• Reserved keywords and labels are italicized in discussions and “quoted and italicized” in the spec. These are usually labels or values of parameters in the input file and must match (using XML matching rules) the specified or allowable values.

• User-defined labels are marked with _ALL_CAPS in this document. In practice, no rules are imposed on these names.

• Lists with too many parameters are described using multiple sections and multiple examples. For most examples we show name of the parent sublist.

Naming convention rule

It is hard to overestimate importance of a reasonable naming convention rule for efficient code development and its daily usage in research.

• Camel-case names should not be used as names for fixed keywords (parameters and parameter lists). The following is a short list of allowed exceptions.

  • The names created by the user are not fixed/reserved keywords and are exempt from the above rule. In this documents, we always prefix user-defined names with the underscore symbol.

  • Proper names such as an individual person, place, or organization, including their derivatives should be spelled using capital letters. Examples: van Genuchten m, Brooks-Corey lambda, Jacobian matrix, and Newton correction.

  • Names of chemical species (inside fixed keywords) should be capitalized. Examples: CO2, H+.

  • A few well-established abbreviations. Their complete list is here: PK, MPC, BDF1, EOS, IEM, PFloTran, pH, TP. Note that names of linear and nonlinear solvers and preconditioners are not included in this list. Thus, we have to use pcg, gmres, nka, amg, ml, and ilu.

  • Units such as energy [J] and temperature [K].

  • The Hilbert spaces L2 and H1. Note that L2 and l2 are different spaces and should be used appropriately.

  • Trilinos parameters. There are a few camel-case parameters that go directly to Trilinos functions and therefore outside of our control, e.g. ML output.

Verbose output

Output of all components of Amanzi is controlled by a standard verbose object list. This list can be inserted in almost any significant component of this spec to produce a verbose output, see the embedded examples. If this list is not specified, the default verbosity value is used.

• “verbosity level” [string] Available options are none, low, medium, high, and extreme. Option extreme is used by the developers only. For communication between users and developers, the recommended option is *high.

• “hide line prefix” [bool] defines prefix for output messages. Default value is true.

• “name” [string] is the name of the prefix.

• “write on rank” [int] is processor rank on which the output is performed. Default is 0.

<ParameterList name="verbose object">
  <Parameter name="verbosity level" type="string" value="medium"/>
  <Parameter name="name" type="string" value="my header"/>
  <Parameter name="hide line prefix" type="bool" value="false"/>
  <Parameter name="write on rank" type="int" value="0"/>
</ParameterList>

Residual debugger

Some components (currently just nonlinear solver, this may change) leverage a residual debugger object for writing, to file, residuals, corrections, and internal iterates of a solve process for solver debugging/work. Control of when these iterates are written is controlled by a few parameters. This should be written sparingly – each attempt at a timestep and each cycle is its own file, and writes its own mesh file, so this should be considered i/o and runtime expensive.

• “cycles start period stop” [Array(int)] the first entry is the start cycle, the second is the cycle period, and the third is the stop cycle or -1 in which case there is no stop cycle. All iterations shall be written at such cycles that satisfy formula cycle = start + n*period, for n=0,1,2,… and cycle < stop if stop != -1.0.

• “cycles start period stop n” [Array(int)] if multiple cycles start-period-stop parameters are needed, then use these parameters with n=0,1,2,…, and not the single “cycles start period stop” parameter.

• “cycles” [Array(int)] an array of single cycles at which all iterations shall be written.

Note: cycle here means the current time integration step and not the global cycle.

<ParameterList name="BDF1">  <!-- parent list -->
<ParameterList name="residual debugger">
  <Parameter name="cycles start period stop" type="Array(int)" value="{0,100,-1}"/>
  <Parameter name="cycles" type="Array(int)" value="{999, 1001}"/>
</ParameterList>
</ParameterList>

Units

Amanzi’s internal default units are SI units except for the concentration.

Curently format of units is very rigid, we need to use multiplication and division operators, e.g. Pa*s and kg/s. However, this allows to mix and match, e.g. m/s^2 or ms^-2 or s^-2*m.

• “concentration” [string] defines units for concentration. Available options are “molar” (default) which is “mol/L” and “SI” which is “mol/m^3”.

<ParameterList>  <!-- parent list -->
<ParameterList name="units">
  <Parameter name="length" type="string" value="m"/>
  <Parameter name="time" type="string" value="s"/>
  <Parameter name="mass" type="string" value="kg"/>
  <Parameter name="temperature" type="string" value="K"/>
  <Parameter name="concentration" type="string" value="molar"/>
</ParameterList>
</ParameterList>

Symbol index

Symbol	Description
\(C\)	concentration \([mol/m^3]\)
\(\boldsymbol{D}\)	dispersion tensor \([m^2/s]\)
\(H\)	specific molar enthalpy \([J/mol]\)
\(\boldsymbol{K}\)	absolute permeability tensor \([m^2]\)
\(\boldsymbol{M}\)	molecular diffusion coefficient \([m^2/s]\)
\(Q\)	source, units depends on equation
\(T\)	temperature \([K]\)
\(U_\alpha\)	specific internal energy of phase \(\alpha\) \([J/mol]\)
\(\varepsilon\)	specific energy \([J/m^3]\)
\(g\)	gravitational acceleration vector \([m/s^2]\)
\(h\)	ponded depth, or the water head over the surface \([m]\)
\(k_r\)	relative permeability \([-]\)
\(p_\alpha\)	pressure of phase \(\alpha\) \([Pa]\)
\(\mathbf{q}\)	Darcy flux vector \([mol\ / m^2 / s]\)
\(s_\alpha\)	saturation of phase \(\alpha\) \([-]\)
\(x\), \(y\), \(z\)	spatial coordinates \([m]\)
\(z\)	elevation \([m]\)
\(\alpha_l\), \(\alpha_t\)	longitudinal and transverse dispersivities \([m]\)
\(\alpha\)	Biot-Willis coefficient \([-]\)
\(\beta_*\)	thermal dilation coefficients \([W/K]\)
\(\kappa\)	thermal conductivity \([W/m/K]\)
\(\mu_\alpha\)	dynamic viscosity of phase \(alpha\) \([Pa\cdot s]\)
\(\eta_\alpha\)	molar density of phase \(\alpha\) \([mol/m^3]\)
\(\rho_\alpha\)	mass density of phase \(\alpha\) \([kg / m^3]\)
\(\phi\)	porosity of the soil \([-]\)
\(\tau\)	tortuosity \([-]\)
\(\theta\)	water content \([mol]\)

CYCLE DRIVER

CycleDriver is a class to hold the cycle driver, which runs the overall, top level timestep loop. It instantiates states, ensures they are initialized, and runs the timestep loop including Vis and restart/checkpoint dumps. It contains one and only one PK – most likely this PK is an MPC of some type – to do the actual work.

The new multi-processor cycle driver provides more flexibility to handle multiphysics process kernels (PKs) and multiple time periods.

cycle_driver-spec

• “component names” [Array(string)] provides the list of species names. It is required for reactive transport.

• “component molar masses” [Array(string)] provides the list of molar masses of species. It is required for proper conversion to and from dimensionless units. Default is 1.

• “number of liquid components” [int] is the number of liquid components.

• “time periods” [list] contains the list of time periods involved in the simulation. The number of time periods is not limited.

  • “TP #” [list] defines a particular time period. The numbering should be sequential starting with 0.

    • “PK tree” [list] describes a hierarchical structure of the process kernels that reflect their weak and strong coupling.

      • “PKNAME” [list] name of PK which is used in the simulation. Name can be arbitrary but the sublist with the same name should exist in the list of PKs (see below).

      • “PK type” [string] specifies the type of PK supported by Amanzi. At the moment available options are (“darcy”, “richards”, “transport”, “one-phase energy”, “two-phase energy”, “reactive transport”, “flow reactive transport”, “thermal richards”, “chemistry”, “transport implicit”, “transport matrix fracture”, “transport matrix fracture implicit”, “flow”, and “darcy matrix fracture”).

      • “start period time” [double] is the start time of the current time period.

      • “end period time” [double] is the end time of the current time period.

      • “maximum cycle number” [int] is the maximum allowed number of cycles in the current time period. Special value -1 means unlimited number of cycles.

      • “initial timestep” [double] is the initial timestep for the current time period.

• “io frequency” [int] defines frequency of printing field statistics.

<ParameterList>  <!-- parent list -->
<ParameterList name="cycle driver">
  <Parameter name="component names" type="Array(string)" value="{H+, Na+, NO3-, Zn++}"/>
  <Parameter name="component molar masses" type="Array(double)" value="{1.0e-3, 23.0e-3, 62.0e-3, 65.4e-3}"/>
  <Parameter name="number of liquid components" type="int" value="4"/>
  <ParameterList name="time periods">
    <ParameterList name="TP 0">
      <ParameterList name="PK tree">
        <ParameterList name="_FLOW and REACTIVE TRANSPORT">
          <Parameter name="PK type" type="string" value="flow reactive transport"/>
          <ParameterList name="_REACTIVE TRANSPORT">
            <Parameter name="PK type" type="string" value="reactive transport"/>
            <ParameterList name="_TRANSPORT">
              <Parameter name="PK type" type="string" value="transport"/>
            </ParameterList>
            <ParameterList name="_CHEMISTRY">
              <Parameter name="PK type" type="string" value="chemistry"/>
            </ParameterList>
          </ParameterList>
          <ParameterList name="_FLOW">
            <Parameter name="PK type" type="string" value="darcy"/>
          </ParameterList>
        </ParameterList>
      </ParameterList>
      <Parameter name="start period time" type="double" value="0.0"/>
      <Parameter name="end period time" type="double" value="1.5778463e+09"/>
      <Parameter name="maximum cycle number" type="int" value="-1"/>
      <Parameter name="initial timestep" type="double" value="1.57680e+05"/>
    </ParameterList>

    <ParameterList name="TP 1">
    ...
    </ParameterList>
  </ParameterList>
</ParameterList>
</ParameterList>

In this simulation, we use the PK labeled as flow reactive transport. It is defined internally as sequential application of two PKs, flow and reactive transport. The latter is defined as sequential application of two PKs, transport and chemistry. Process kernel reactive transport can susbcycle with respect to flow. Process kernel chemistry can susbcycle with respect to transport.

Time period control

A set of times that simulation hits exactly can be used to avoid problems with sudden change of boundary conditions or source/sink terms. This list must NOT include start times for time periods TP #.

time_period_controls-spec

• “start times” [Array(double)] is the list of particular times that we want to hit exactly.

• “initial time step” [Array(double)] is the size of the first time step after we hit a special time specified above.

• “maximum time step” [Array(double)] allows the user to limit the time step between two particular times.

<ParameterList name="cycle driver">  <!-- parent list -->
<ParameterList name="time period control">
  <Parameter name="start times" type="Array(double)" value="{3.16e+10, 6.32e+10}"/>
  <Parameter name="initial time step" type="Array(double)" value="{100.0, 100.0}"/>
  <Parameter name="maximum time step" type="Array(double)" value="{3.2e+8, 4e+17}"/>
 </ParameterList>
 </ParameterList>

Between approximately 1000 and 2000 years, we limit the maximum time step to 10 years.

Restart from checkpoint data file

A user may request to restart a simulation from a checkpoint data file by creating list restart. In this scenario, the code will overwrite data initialized using the input XML file. The purpose of restart is to continue the simulation that has been terminated before for some reasons, e.g. because its allocation of time ran out. The value for the current time and current cycle is read from the checkpoint file.

restart_from_file-spec

• “restart” [list]

  • “file name” [string] provides name of the existing checkpoint data file to restart from.

<ParameterList name="cycle driver">  <!-- parent list -->
<ParameterList name="restart">
  <Parameter name="file name" type="string" value="_CHECK00123.h5"/>
</ParameterList>
</ParameterList>

In this example, Amanzi is restarted with all state data initialized from file CHECK00123.h5.

STATE

State is a simple data-manager, allowing PKs to require, read, and write various fields.

• Acts as a factory for data through the various require methods.

• Provides some data protection by providing both const and non-const data pointers to PKs.

• Provides some initialization capability – this is where all independent variables can be initialized (as independent variables are owned by state, not by any PK).

state-spec

• “evaluators” [evaluator-typedinline-spec-list] A list of evaluators.

• “initial conditions” [list] A list of constant-in-time data. Note that “initial conditions” is not a particularly descriptive name here – PDE initial conditions are generally not here. This list consists of

• “model parameters” [list] A list of shared model parameters that can be used across all evaluators.

evaluator-typedinline-spec

• “evaluator type” [string] Type of the evaluator Included for convenience in defining data that is not in the dependency graph, constants are things (like gravity, or atmospheric pressure) which are stored in state but never change. Typically they’re limited to scalars and dense, local vectors.

Example:

<ParameterList name="state">
  <Parameter name="initialization filename" type="string" value="_CHECK00123.h5"/>
  <ParameterList name="evaluators">
    <ParameterList name="pressure">
      <Parameter name="evaluator type" type="string" value="primary variable" />
    </ParameterList>
  </ParameterList>

  <ParameterList name="initial conditions">
    <Parameter name="time" type="double" value="0.0">
    <ParameterList name="atmospheric pressure">
      <Parameter name="value" type="double" value="101325.0" />
    </ParameterList>
    <ParameterList name="gravity">
      <Parameter name="value" type="Array(double)" value="{0.0,0.0,-9.80665}" />
    </ParameterList>
  </ParameterList>
</ParameterList>

Warning

• “initialization filename” [string] (optional) provides name of the existing checkpoint data file. The initialization sequence is as follows. First, we try to initialize a field using the provided checkpoint file. Second, regardless of the outcome of the previous step, we try to initialize the field using the sublist “initial conditions”. By design, the second step allows us to overwrite only part for the field.

• “initial conditions” [list] There are several options available to initialize fields using this list: “restart file” - read field from existing HDF5 file, “exodus file initialization” - read field from existing Exodus file, “cells from file” - read cell components from HDF5 file, “constant” - set field values to constant, “initialize from 1D column” - initialize 1D column from file and “function” - field is initialized by function. Initialization time is specified with “time” that defaults to 0.

Field evaluators

There are five different types of field evaluators.

Independent variable field evaluator

An independent variable field evaluator has no dependencies and is specified by a function. Typically, it is evaluated once per simulation. The evaluator has the following control parameters.

• “field evaluator type” [string] The value of this parameter is used by the factory of evaluators. The available option are “independent variable”, “primary variable”, and “secondary variable”.

• “constant in time” [bool] specifies time-dependence nature of the field.

• “function” [list] defines a piecewise continuous function for calculating the independent variable. In may contain multiple sublists “_DOMAIN” with identical structure.

  • “_DOMAIN” [list]

    • “region” [string] specifies domain of the function, a single region.

    • “regions” [Array(string)] is the alternative to option “region”, domain on the function consists of many regions.

    • “component” [string] specifies geometric object associated with the mesh function. Available options are “cell”, “face”, and “node”.

    • “function” [list] defines an analytic function for calculation. Its structure is described in the Functions section below.

    • “initialize faces from cells” [bool] instructs state to initialize face-component and boundary face-component (if any) of a composite vector from a cell-component using simple averaging. Default is false.

<ParameterList name="field_evaluators">  <!-- parent list -->
<ParameterList name="saturation_liquid">
  <Parameter name="field evaluator type" type="string" value="independent variable"/>
  <Parameter name="constant in time" type="bool" value="true"/>
  <ParameterList name="function">
    <ParameterList name="_DOMAIN">
      <Parameter name="region" type="string" value="_ALL DOMAIN"/>
      <Parameter name="component" type="string" value="cell"/>
      <ParameterList name="function">
        <ParameterList name="function-constant">
          <Parameter name="value" type="double" value="0.8"/>
        </ParameterList>
      </ParameterList>
    </ParameterList>
  </ParameterList>
  <ParameterList name="verbose object">
    <Parameter name="verbosity level" type="string" value="extreme"/>
  </ParameterList>
</ParameterList>
</ParameterList>

In this example field saturation_liquid is defined as a cell-based variable with constant value 0.8. Note that the user-defined name for this field cannot have spaces.

Independent variable field evaluator from file

This evaluator is typically used for providing data that are functions of space and time. Data is provided, discretely (e.g. with one data point per cell/face/node), at a series of time slices. The time slices are interpolated linearly in time to provide the value.

Within the file, data is expected to meet the following (HDF5) layout:

/time : a 1D array of length NTIMES, providing the time in seconds.
/variable_name.ENTITY.DOF  (group)

   /0 : a 1D array of length NENTITIES, providing the values for each entity
        at time /time[0]
   /1 : ...
   /NTIMES-1 : 1D array at time /time[NTIMES-1]

This evaluator is used by providing the option:

“evaluator type” == “independent variable”

independent-variable-from-file-evaluator-spec

• “filename” [string] Path to the file.

• “variable name” [string] Name of the dataset to read from the file.

• “domain name” [string] domain Name of the domain on which the field is defined.

• “component name” [string] cell Name of the component in the field to populate.

• “mesh entity” [string] cell Name of the entity on which the component is defined.

• “number of dofs” [int] 1 Number of degrees of freedom to read.

• “time function” [function-spec] optional If provided, time is first manipulated by this function before interpolation. This is useful for things like cyclic data, which can use a modulo time function to repeat the same data.

<ParameterList name="field_evaluators">  <!-- parent list -->
<ParameterList name="porosity">
  <Parameter name="field evaluator type" type="string" value="independent variable from file"/>
  <Parameter name="filename" type="string" value="_DATA_FILE.h5"/>
  <Parameter name="domain name" type="string" value="domain"/>
  <Parameter name="variable name" type="string" value="porosity"/>
  <Parameter name="component name" type="string" value="cell"/>
  <Parameter name="mesh entity" type="string" value="cell"/>
  <Parameter name="number of dofs" type="int" value="1"/>

  <ParameterList name="time function">
    <Parameter name="times" type="Array(double)" value="{1.0, 2.0, 3.0}"/>
  </ParameterList>
</ParameterList>
</ParameterList>

The field porosity is defined as a cell-based variable and interpolated between three time intervals.

Constant variable field evaluator

Constant variable field evaluator as a simplified version of independent field evaluator from file which allows one to define constant in time field. Initialization of the field has to be done in the initial conditions sublist of state.

<ParameterList name="initial conditions">  <!-- parent list -->
<ParameterList name="porosity">
  <ParameterList name="function">
    <ParameterList name="_ANY NAME">
      <Parameter name="regions" type="Array(string)" value="{_ALL DOMAIN}"/>
      <Parameter name="component" type="string" value="cell"/>
      <ParameterList name="function">
        <ParameterList name="function-constant">
          <Parameter name="value" type="double" value="90000.0"/>
        </ParameterList>
      </ParameterList>
    </ParameterList>
  </ParameterList>
</ParameterList>
</ParameterList>

<ParameterList name="field_evaluators">  <!-- parent list -->
<ParameterList name="porosity">
  <Parameter name="field evaluator type" type="string" value="constant variable"/>
</ParameterList>
</ParameterList>

Primary variable field evaluator

An evaluator with no dependencies that serves as the primary variable to be solved for by a PK. Note that users almost never are required to write an input spec for these – they are controlled by the PK and therefore the input spec for this evaluator is written by that PK.

evaluator-primary-spec

• “tag” [string] “” Time tag at which this primary variable is used.

Secondary variable field evaluators

Secondary fields are derived either from primary fields or other secondary fields. There are two types of secondary fields evaluators. The first type is used to evaluate a single field. The second type is used to evaluate efficiently (in one call of an evaluator) multiple fields.

Typically, secondary fields are created by high-level PKs during the setup phase and inserted automatically in the list of evaluators. The related XML syntax can provide various parameters needed for evaluation as explained in two examples below. The developer can create a secondary field evaluator using common parameters as well as custom parameters (see the examples).

• “evaluator dependencies” [Array(string)] provides a list of fields on which this evaluator depends.

• “check derivatives” [bool] allows the develop to check derivatives with finite differences. This is the expensive option involving finite difference approximations and is recommended for code debugging only. Default is false.

• “finite difference epsilon” [double] defines the finite difference epsilon. Default is 1e-10.

<ParameterList name="field_evaluators">  <!-- parent list -->
<ParameterList name="molar_density_liquid">
  <Parameter name="field evaluator type" type="string" value="eos"/>
  <Parameter name="eos basis" type="string" value="both"/>
  <Parameter name="molar density key" type="string" value="molar_density_liquid"/>
  <Parameter name="mass density key" type="string" value="mass_density_liquid"/>
  <ParameterList name="EOS parameters">
    <Parameter name="eos type" type="string" value="liquid water"/>
  </ParameterList>
  <ParameterList name="verbose object">
    <Parameter name="verbosity level" type="string" value="extreme"/>
  </ParameterList>
</ParameterList>
</ParameterList>

In this example, the secondary field internal_energy_rock is evaluated using one of the internal iem evaluators. A particular evaluator is selected dynamically using parameter iem type.

Initial conditions

Constant scalar field

A constant field is the global (with respect to the mesh) constant. At the moment, the set of such fields includes fluid_density and fluid_viscosity. The initialization requires to provide a named sublist with a single parameter value.

<ParameterList name="initial conditions">  <!-- parent list -->
<ParameterList name="fluid_density">
  <Parameter name="value" type="double" value="998.0"/>
</ParameterList>
</ParameterList>

Constant vector field

A constant vector field is the global (with respect to the mesh) vector constant. At the moment, the set of such vector constants includes gravity. The initialization requires to provide a named sublist with a single parameter value. In three dimensions, it looks like

<ParameterList name="initial conditions">  <!-- parent list -->
<ParameterList name="gravity">
  <Parameter name="value" type="Array(double)" value="{0.0, 0.0, -9.81}"/>
</ParameterList>
</ParameterList>

A scalar field

A variable scalar field is defined by a few functions (labeled with _MESH BLOCK # in our example) with non-overlapping domains. The required parameters for each function are region, component, and function.

• “regions” [Array(string)] is list of mesh regions where the function should be applied, the domain of the function.

• “component” [string] specifies a mesh object on which the discrete field is defined.

Optional parameters are write checkpoint and write vis. These parameters define whether the field has to be written into checkpoints of visualization files. Default values are true.

<ParameterList name="initial conditions">  <!-- parent list -->
<ParameterList name="pressure">
  <Parameter name="write checkpoint" type="bool" value ="false">
  <Parameter name="write vis" type="bool" value ="true">
  <ParameterList name="function">
    <ParameterList name="_MESH BLOCK 1">
      <Parameter name="regions" type="Array(string)" value="_DOMAIN 1"/>
      <Parameter name="component" type="string" value="cell"/>
      <ParameterList name="function">
        <ParameterList name="function-constant">
          <Parameter name="value" type="double" value="90000.0"/>
        </ParameterList>
      </ParameterList>
    </ParameterList>
    <ParameterList name="_MESH BLOCK 2">
      <Parameter name="regions" type="Array(string)" value="{_DOMAIN 2, _DOMAIN 3}"/>
      ...
    </ParameterList>
  </ParameterList>
</ParameterList>
</ParameterList>

In this example, the field pressure has constant value 90000 [Pa] in each mesh cell of the first region. The second mesh block will define the pressure in the second mesh region and so on. Note that the names of functions may coincide with the names of regions.

A vector or tensor field

A variable tensor (or vector) field is defined similarly to a variable scalar field. The difference lies in the definition of the function which is now a multi-valued function.

tensor_field-spec

• “number of dofs” [int] is the number of components in the vector or tensor.

• “function type” [string] defines the function type. The only available option is “composite function”.

• “dot with normal” [bool] triggers the special initialization of a vector field such as the “volumetric_flow_rate”. This field is defined by projection of the velocity (a vector field) on face normals. Changing value to false will produce the vector field.

Optional parameters are write checkpoint, and write vis. These parameters define whether the field has to be written into checkpoints of vis files. Default values are true.

<ParameterList name="initial conditions">  <!-- parent list -->
<ParameterList name="volumetric_flow_rate">
  <Parameter name="write checkpoint" type="bool" value="true"/>
  <Parameter name="write vis" type="bool" value="false"/>
  <Parameter name="dot with normal" type="bool" value="true"/>

  <ParameterList name="function">
    <ParameterList name="_MESH BLOCK 1">
      <Parameter name="regions" type="Array(string)" value="{_ALL DOMAIN}"/>
      <Parameter name="component" type="string" value="face"/>
      <ParameterList name="function">
        <Parameter name="number of dofs" type="int" value="2"/>
        <Parameter name="function type" type="string" value="composite function"/>
        <ParameterList name="dof 1 function">
          <ParameterList name="function-constant">
            <Parameter name="value" type="double" value="0.002"/>
          </ParameterList>
        </ParameterList>
        <ParameterList name="dof 2 function">
          <ParameterList name="function-constant">
            <Parameter name="value" type="double" value="0.001"/>
          </ParameterList>
        </ParameterList>
      </ParameterList>
    </ParameterList>
  </ParameterList>
</ParameterList>
</ParameterList>

In this example the constant Darcy velocity (0.002, 0.001) [m/s] is dotted with the face normal producing one number per mesh face.

Geochemical constraint

We can define geochemical constraint as follows:

<ParameterList name="initial conditions">  <!-- parent list -->
<ParameterList name="geochemical conditions">
  <ParameterList name="initial">
    <Parameter name="regions" type="Array(string)" value="{_ENTIRE DOMAIN}"/>
  </ParameterList>
</ParameterList>
</ParameterList>

Initialization from Exodus II file

Some fields can be initialized from Exodus II files. For each field, an additional sublist has to be added to the named sublist of State list with the file name and the name of attributes. For a serial run, the file extension must be .exo. For a parallel run, it must be .par.

• “attributes” [Array(string)] defines names of attributes. The number of names must be equal to the number of components in the field. The names can be repeated. Scalar fields (e.g. porosity) require one name, tensorial fields (e.g. permeability) require two or three names.

<ParameterList name="initial conditions">  <!-- parent list -->
<ParameterList name="permeability">
  <ParameterList name="exodus file initialization">
    <Parameter name="file" type="string" value="_MESH_AND_DATA.exo"/>
    <Parameter name="attributes" type="Array(string)" value="{permx, permx, permz}"/>
  </ParameterList>
</ParameterList>
</ParameterList>

Initialization from HDF5 file

Some field can be initialized from HDF5 file. The field has to written to HDF5 file as 2D array (number_elements, number_of_components) and has to name as field_name.entity.component, e.g transport_porosity.cell.0. Parameter “cell from file” initializes only cell part of the field.

<ParameterList name="initial conditions">  <!-- parent list -->
<ParameterList name="transport_porosity">
  <Parameter name="restart file" type="string" value="_TEST1.h5"/>
  </ParameterList>
<ParameterList name="porosity">
  <Parameter name="cells from file" type="string" value="_TEST3.h5"/>
</ParameterList>
</ParameterList>

Initialization of 1D column

It is possible to initialize only 1D column portion of a particular field.

<ParameterList name="initial conditions">  <!-- parent list -->
<ParameterList name="temperature">
  <ParameterList name="initialize from 1D column">
    <Parameter name="file" type="string" value="_COLUMN_DATA.h5"/>
    <Parameter name="z header" type="string" value="/z"/>
    <Parameter name="f header" type="string" value="/temperature"/>
    <Parameter name="coordinate orientation" type="string" value="depth"/>
    <Parameter name="surface sideset" type="string" value="surface"/>
  </ParameterList>
</ParameterList>
</ParameterList>

Mesh partitioning

Amanzi’s state has a number of tools to verify completeness of initial data. This is done using list mesh partitions. Each sublist must have parameter region list specifying regions that define unique partition of the mesh.

<ParameterList name="state">  <!-- parent list -->
<ParameterList name="mesh partitions">
  <ParameterList name="_MATERIALS">
    <Parameter name="region list" type="Array(string)" value="{_SAND1,_CLAY,_SAND2}"/>
  </ParameterList>
</ParameterList>
</ParameterList>

In this example, we verify that three mesh regions representing sand and clay cover completely the mesh without overlaps. If so, all material fields, such as porosity and permeability, will be initialized properly.

Data IO control

Two parameters below allow us to control fields that will go into a visualization file. First, we remove all fields matching the patterns specified by blacklist. Second, we add all fields matching the patterns specified by whitelist. Both parameters are optional.

• “blacklist” [Array(string)] list of fields that should not be written to the visualization file. Standard regular expression rules can be used, e.g. “(secondary_)(.*)” skips all fields those names start with “secondary_”.

• “whitelist” [Array(string)] list of fields that should be written to the visualization file. Standard regular expression rules can be used, e.g. “(primary_)(.*)” adds all fields those names start with “primary_”.

Example

The complete example of a state initialization is below. Note that _MATERIAL1 and _MATERIAL2 must be labels of the existing regions that cover the computational domain. The fields porosity and pressure are constant over the whole domain. The field permeability is the piecewise constant diagonal tensor.

<ParameterList name="state">
<ParameterList name="field evaluators">
  <ParameterList name="porosity">
    <ParameterList name="function">
      <ParameterList name="_ANY NAME ">
        <Parameter name="regions" type="Array(string)" value="{_ALL DOMAIN}"/>
        <Parameter name="component" type="string" value="cell"/>
        <ParameterList name="function">
          <ParameterList name="function-constant">
            <Parameter name="value" type="double" value="0.408"/>
          </ParameterList>
        </ParameterList>
      </ParameterList>
    </ParameterList>
  </ParameterList>
</ParameterList>

<ParameterList name="initial conditions">
  <ParameterList name="fluid_density">
    <Parameter name="value" type="double" value="998.0"/>
  </ParameterList>

  <ParameterList name="gravity">
    <Parameter name="value" type="Array(double)" value="{0.0, -9.81}"/>
  </ParameterList>

  <ParameterList name="pressure">
    <ParameterList name="function">
      <ParameterList name="_ANY NAME">
        <Parameter name="regions" type="Array(string)" value="{_ALL DOMAIN}"/>
        <Parameter name="component" type="string" value="cell"/>
        <ParameterList name="function">
          <ParameterList name="function-constant">
            <Parameter name="value" type="double" value="90000.0"/>
          </ParameterList>
        </ParameterList>
      </ParameterList>
    </ParameterList>
  </ParameterList>

  <ParameterList name="permeability">
    <ParameterList name="function">
      <ParameterList name="_ANY NAME">
        <Parameter name="regions" type="Array(string)" value="_MATERIAL1"/>
        <Parameter name="component" type="string" value="cell"/>
        <ParameterList name="function">
          <Parameter name="function type" type="string" value="composite function"/>
          <Parameter name="number of dofs" type="int" value="2"/>
          <ParameterList name="dof 1 function">
            <ParameterList name="function-constant">
              <Parameter name="value" type="double" value="1e-12"/>
            </ParameterList>
          </ParameterList>
          <ParameterList name="dof 2 function">
            <ParameterList name="function-constant">
              <Parameter name="value" type="double" value="1e-13"/>
            </ParameterList>
          </ParameterList>
        </ParameterList>
      </ParameterList>
      <ParameterList name="_ANY_NAME">
        <Parameter name="regions" type="Array(string)" value="_MATERIAL2"/>
        <Parameter name="component" type="string" value="cell"/>
        <ParameterList name="function">
          <Parameter name="function type" type="string" value="composite function"/>
          <Parameter name="number of dofs" type="int" value="2"/>
          <ParameterList name="dof 1 function">
            <ParameterList name="function-constant">
              <Parameter name="value" type="double" value="2e-13"/>
            </ParameterList>
          </ParameterList>
          <ParameterList name="dof 2 function">
            <ParameterList name="function-constant">
              <Parameter name="value" type="double" value="2e-14"/>
            </ParameterList>
          </ParameterList>
        </ParameterList>
      </ParameterList>
    </ParameterList>
  </ParameterList>
</ParameterList>
</ParameterList>

PROCESS KERNELS (PKs)

A process kernel represents a single or system of partial/ordinary differential equation(s) or conservation law(s), and is used as the fundamental unit for coupling strategies.

Implementations of this interface typically are either an MPC (multi-process coupler) whose job is to heirarchically couple several other PKs and represent the system of equations, or a Physical PK, which represents a single equation.

Note there are two PK specs – the first is the “typed” spec, which appears in the “cycle driver” list in the PK tree and has no other parameters other than its type and its children. The second is the spec for the base class PK, which is inherited and included by each actual PK, lives in the “PKs” sublist of “main”, and has all needed parameters.

pk-spec

• “PK type” [string] One of the registered PK types. Note this must match the corresponding entry in the [pk-typed-spec]

• “verbose object” [verbose-object-spec] optional See Verbose Object

Example:

<ParameterList name="PK tree">
  <ParameterList name="my cool MPC PK">
    <Parameter name="PK type" type="string" value="my cool MPC PK"/>
    <ParameterList name="sub PK 1">
      ...
    </ParameterList>
    <ParameterList name="sub PK 2">
      ...
    </ParameterList>
    ...
  </ParameterList>
</ParameterList>

BASE PKs

There are several types of PKs, and each PK has its own valid input spec. However, there are three main types of PKs, from which nearly all PKs derive. Note that none of these are true PKs and cannot stand alone.

PK_Physical

Physical PKs represent the physics capability implemented within Amanzi.

pks-typed-spec

• “darcy” See PK: Flow

• “richards” See PK: Flow

• “transport” See PK: Transport

• “transport implicit” See PK: Transport

• “one-phase energy” See PK: Energy

• “chemistry alquimia” See PK: Chemistry Alquimia

• “chemistry amanzi” See PK: Chemistry Amanzi

• “shallow water” See PK: Shallow Water

• “navier stokes” See PK: Navier-Stokes

• “multiphase” See PK: Multiphase

• “mechanics” See PK: Mechanics

PK_BDF

This is a base class from which PKs that want to use the BDF series of implicit time integrators must derive. It specifies both the BDFFnBase interface and implements some basic functionality for BDF PKs.

PK_PhysicalBDF

A base class for all PKs that are both physical, in the sense that they implement an equation and are not couplers, and BDF, in the sense that they support the implicit integration interface. This largely just supplies a default error norm based on error in conservation relative to the extent of the conserved quantity.

By default, the error norm used by solvers is given by:

\[ENORM(u, du) = |du| / ( a_{tol} + r_{tol} * |u| )\]

The defaults here are typically good, or else good defaults are set in the code, so usually are not supplied by the user.

PK_DomainFunction

Prototype class for computing boundary and source terms. PK’s classes for source and boundary terms may inherit from this class or mimic its functionality.

Simple function

To be written.

Coupling function

This provides coupling of fields located on conforming manifold and 3D meshes. For the 3D mesh, the coupling generates a list of boundary conditions. For the manifold mesh, the coupling creates a list of sources. The manifold mesh has a child to parent map, e.g. a manifold cell to a 3D mesh face. For the 3D mesh, we need the reverse map. There are a few coupling submodels:

• “rate” computes data using the following formulas:

  \[\begin{split}\begin{array}{l} value[i][c] = value[i][c] + flux[f] * external\_field[i][cc] / V_c \\ value[N][c] = value[N][c] - \Delta t * flux[f] \end{array}\end{split}\]

  where cc is a space cell incident to face f, and N is the auxiliary value added to the result. Note that an internal face f (resp., boundary face f) is shared by two (resp. one) space cells.

• “field” computes data using the following formula:

  \[value[i][f] = external\_field[i][c]\]

• “conserved quantity” is not used in Amanzi.

Weight function

The total source Q is given for each domain. A weighted source distribution model is employed. The cell-based source density is calculated as (Q / W_D) * weight, where W_D is the weighted domain volume. The weight is defined globally, for the whole computational domain.

Field function

This leverages the Arcos DAG to evaluate a user-provided field and use that in the domain function.

domain-function-field-spec

• “field key” [string] Field used in the domain function.

• “component” [string] Component of the field. Default is “cell”.

• “scaling factor” [double] Constant multiplication factor. Default is 1.

PHYSICS PKs

Each PK represents a single physical process implemented within Amanzi such as mass transport, solute transport and dispersion, energy transport, rock deformation, and chemical reactions.

Flow PK

Mathematical models

A few PDE models can be instantiated using the parameters described below.

Fully saturated flow

The conceptual PDE model for the fully saturated flow is

\[\left(\frac{S_s}{g} + \frac{S_y}{Lg}\right)\frac{\partial p_l}{\partial t} = -\boldsymbol{\nabla} \cdot (\rho_l \boldsymbol{q}_l) + Q, \quad \boldsymbol{q}_l = -\frac{\boldsymbol{K}}{\mu} (\boldsymbol{\nabla} p - \rho_l \boldsymbol{g}),\]

where \(S_s\) and \(S_y\) are specific storage [1/m] and specific yield [-], respectively, \(L\) is characteristic length [m], \(\rho_l\) is fluid density [\(kg / m^3\)], \(Q\) is source or sink term [\(kg / m^3 / s\)], \(\boldsymbol{q}_l\) is the Darcy velocity [\(m/s\)], and \(\boldsymbol{g}\) is gravity [\(m/s^2\)]. The specific storgae can be defined using

\[S_s = \left(\phi\, \beta_f + \beta_m\right)\rho_l\,g\]

where \(\beta_f\) [1/Pa] and \(\beta_m\) [1/Pa] are fluid and matrix compressibilities, respectively.

Partially saturated flow with water vapor

The conceptual PDE model for the partially saturated flow with water vapor includes liquid phase (liquid water) and gas phase (water vapor):

\[\frac{\partial \theta}{\partial t} = - \boldsymbol{\nabla} \cdot (\eta_l \boldsymbol{q}_l) - \boldsymbol{\nabla} \cdot (\boldsymbol{K}_g \boldsymbol{\nabla} \big(\frac{p_v}{p_g}\big)) + Q, \quad \boldsymbol{q}_l = -\frac{\boldsymbol{K} k_r}{\mu} (\boldsymbol{\nabla} p - \rho_l \boldsymbol{g})\]

where \(\theta\) is total water storage [\(mol/m^3\)], \(\eta_l\) is molar density of liquid (water) [\(mol/m^3\)], \(\rho_l\) is fluid density [\(kg/m^3\)], \(Q\) is source or sink term [\(mol/m^3/s\)], \(\boldsymbol{q}_l\) is the Darcy velocity [\(m/s\)], \(k_r\) is relative permeability [-], \(\boldsymbol{g}\) is gravity [\(m/s^2\)], \(p_v\) is the vapor pressure [Pa], \(p_g\) is the gas pressure [Pa], and \(\boldsymbol{K}_g\) is the effective diffusion coefficient of the water vapor. We define

\[\theta = \phi \eta_l s_l + \phi \eta_g (1 - s_l) X_g\]

where \(s_l\) is liquid saturation [-], \(\phi\) is porosity [-], \(\eta_g\) is molar density of water vapor [\(mol/m^3\)], and \(X_g\) is molar fraction of water vapor. The effective diffusion coefficient of the water vapor is given by

\[\boldsymbol{K}_g = \phi s_g \tau_g \eta_g \boldsymbol{D}_g\]

where \(s_g\) is gas saturation [-], \(\tau_g\) is the tortuosity of the gas phase [-], \(\eta_g\) is the molar density of gas [\(kg/m^3\)], and \(\boldsymbol{D}_g\) is the diffusion coefficient of the gas phase [\(m^2/s\)], The gas pressure \(p_g\) is set to the atmosperic pressure and the vapor pressure model assumes thermal equlibrium of liquid and gas phases:

\[p_v = P_{sat}(T) \exp\left(\frac{P_{cgl}}{\eta_l R T}\right)\]

where \(R\) is the ideal gas constant [\(kg m^2/K/mol/s^2\)], \(P_{cgl}\) is the liquid-gas capillary pressure [Pa], \(P_{sat}\) is the saturated vapor pressure [Pa], and \(T\) is the temperature [K]. The diffusion coefficient is based of TOUGHT2 model

\[D_g = D_0 \frac{P_{ref}}{p} \left(\frac{T}{273.15}\right)^a\]

where \(D_0 = 2.14 \cdot 10^{-5}\), \(P_{ref}\) is atmospheric pressure, and \(a = 1.8\). finally we need a model for the gas tortuosity. We use the Millington and Quirk model:

\[\tau_g = \phi^\beta s_g^\gamma\]

where \(\beta = 1/3\) and \(\gamma = 7/3\).

Isothermal partially saturated flow with dual porosity model

The conceptual model for the partially saturated flow with dual porosity model assumes that water flow is restricted to the fractures and the water in the matrix does not move. The rock matrix represents immobile pockets that can exchange, retain and store water but do not permit convective flow. This leads to dual-porosity type flow and transport models that partition the liquid phase into mobile and immobile regions. The Richards equation in the mobile region is augmented by the water exchange term \(\Sigma_w\):

\[\frac{\partial \theta_{lf}}{\partial t} = -\boldsymbol{\nabla} \cdot (\eta_l \boldsymbol{q}_l) -\frac{K_m\,k_{rm}\,\eta_l}{\mu\,L_m}\, \nabla p_m + Q_f, \qquad \boldsymbol{q}_l = -\frac{\boldsymbol{K}_f\, k_{rf}}{\mu} (\boldsymbol{\nabla} p_f - \rho_l \boldsymbol{g})\]

where \(p_f\) is fracture pressure [Pa], \(p_m\) is matrix pressure [Pa], \(L_m\) is the characteristic matrix depth defined typically as the ratio of a matrix block [m], and \(Q_f\) is source or sink term [\(kg \cdot m^{-3} \cdot s^{-1}\)]. The equation for water balance in the matrix is

\[\frac{\partial \theta_{lm}}{\partial t} = Q_m +\nabla\cdot \left(\frac{K_m\, k_{rm}\,\eta_l}{\mu}\, \nabla p_{m}\right),\]

where \(Q_m\) is source or sink term [\(kg / m^3 / s\)]. The water storages are defined as

\[\theta_f = \phi_f\, \eta_l\, s_{lf},\quad \theta_m = \phi_m\, \eta_l\, s_{lm},\]

where saturations \(s_{lf}\) and \(s_{lm}\) may use different capillary pressure - saturation models. In the simplified model, the rate of water exchange between the fracture and matrix regions is proportional to the difference in hydraulic heads:

\[\frac{K_m\,k_{rm}\,\eta_l}{\mu\,L_m}\, \nabla p_m \approx \alpha_w (h_f - h_m),\]

where \(\alpha_w\) is the mass transfer coefficient. Since hydraulic heads are needed for both regions, this equation requires retention curves for both regions and therefore is nonlinear.

Physical models and assumptions

This list is used to summarize physical models and assumptions, such as coupling with other PKs. This list is often generated or extended by a high-level MPC PK. In the code development, this list plays a two-fold role. First, it provides necessary information for coupling different PKs such as flags for adding a vapor diffusion to Richards’ equation. Second, the developers may use it instead of a factory of evaluators such as creation of primary and secondary evaluators for rock porosity models. Combination of both approaches may lead to a more efficient code.

flow_assumptions-spec

• “vapor diffusion” [bool] is set up automatically by a high-level PK, e.g. by EnergyFlow PK. The default value is “false”.

• “flow and transport in fractures” [bool] indicates that Darcy flow is calculated in fractures. This option is ignored is mesh dimentionaly equals to manifold dimensionality.

• “multiscale model” [string] specifies a multiscale model. Available options are “single porosity” (default) and “dual continuum discontinum matrix”.

• “viscosity model” [string] changes the evaluator for liquid viscosity. Available options are “generic” and “constant viscosity” (default).

• “porosity model” [string] specifies an isothermal porosity model. Available options are “compressible” and “constant” (default).

• “coupled matrix fracture flow” [string] specifies PK’s role in the strong coupling of two flow PKs. The value is either “matrix” or “fracture”.

• “eos lookup table” [string] provides the name for optional EOS lookup table.

<ParameterList name="flow">  <!-- parent list -->
<ParameterList name="physical models and assumptions">
  <Parameter name="vapor diffusion" type="bool" value="false"/>
  <Parameter name="viscosity model" type="string" value="constant viscosity"/>
  <Parameter name="porosity model" type="string" value="compressible"/>
  <Parameter name="multiscale model" type="string" value="single porosity"/>
  <Parameter name="coupled matrix fracture flow" type="string" value="matrix"/>
  <Parameter name="eos lookup table" type="string" value="h2o.eos"/>
</ParameterList>
</ParameterList>

Global parameters

• “domain name” [string] specifies mesh name that defined domain of this PK. Default is “domain”.

Water retention models

User defines water retention models in sublist water retention models. It contains as many sublists, e.g. SOIL_1, SOIL_2, etc, as there are different soils. This list is required for the Richards problem only.

The water retention models are associated with non-overlapping regions. Each of the sublists (e.g. Soil 1) includes a few mandatory parameters: region name, model name, and parameters for the selected model.

wrm-spec

• “water retention model” [string] specifies a model for the soil. The available models are “van Genuchten”, “Brooks Corey”, and “fake”. The later is used only to set up a simple analytic solution for convergence study.

  • The model “van Genuchten” requires “van Genuchten alpha” [double], “van Genuchten m” [double], “van Genuchten l” [double], “residual saturation liquid” [double], and “relative permeability model” [string].

  • The model “Brooks-Corey” requires “Brooks Corey lambda” [double], “Brooks Corey alpha” [double], “Brooks Corey l” [double], “residual saturation liquid” [double], and “relative permeability model” [string].

  • The model “saturated” does not require any parameters. It is equaivalent to the Darcy model if all other parameters (e.g. porosity) are constant

  • The model “linear” requires “alpha” [double] and “residual saturation liquid” [double].

• “relative permeability model” [string] The available options are “Mualem” (default) and “Burdine”.

• “regularization interval” [double] removes the kink in the water retention curve at the saturation point using a cubic spline. The parameter specifies the regularization region [Pa]. Default value is 0.

Amanzi performs rudimentary checks of validity of the provided parameters. The relative permeability curves can be calculated and saved in an ASCI file if the list output is provided. This list has two mandatory parameters:

• “file” [string] is the user defined file name. It should be different for each soil.

• “number of points” [int] is the number of data points. Each file will contain a table with three columns: saturation, relative permeability, and capillary pressure. The data points are equidistributed between the residual saturation and 1.

<ParameterList name="flow">  <!-- parent list -->
<ParameterList name="water retention models">
  <ParameterList name="_SOIL_1">
    <Parameter name="regions" type="Array(string)" value="{_TOP}"/>
    <Parameter name="water retention model" type="string" value="van Genuchten"/>
    <Parameter name="van Genuchten alpha" type="double" value="0.000194"/>
    <Parameter name="van Genuchten m" type="double" value="0.28571"/>
    <Parameter name="van Genuchten l" type="double" value="0.5"/>
    <Parameter name="residual saturation liquid" type="double" value="0.103"/>
    <Parameter name="regularization interval" type="double" value="100.0"/>
    <Parameter name="relative permeability model" type="string" value="Mualem"/>
    <ParameterList name="output">
      <Parameter name="file" type="string" value="soil1.txt"/>
      <Parameter name="number of points" type="int" value="1000"/>
    </ParameterList>
  </ParameterList>

  <ParameterList name="_SOIL_2">
    <Parameter name="regions" type="Array(string)" value="{_BOTTOM}"/>
    <Parameter name="water retention model" type="string" value="Brooks Corey"/>
    <Parameter name="Brooks Corey lambda" type="double" value="0.0014"/>
    <Parameter name="Brooks Corey alpha" type="double" value="0.000194"/>
    <Parameter name="Brooks Corey l" type="double" value="0.51"/>
    <Parameter name="residual saturation liquid" type="double" value="0.103"/>
    <Parameter name="regularization interval" type="double" value="0.0"/>
    <Parameter name="relative permeability model" type="string" value="Burdine"/>
  </ParameterList>

  <ParameterList name="_SOIL_3">
    <Parameter name="regions" type="Array(string)" value="{_MIDDLE}"/>
    <Parameter name="alpha" type="double" value="5e-5"/>  <!-- Pa^-1 -->
    <Parameter name="residual saturation liquid" type="double" value="0.02"/>
  </ParameterList>
</ParameterList>
</ParameterList>

In this example, we define three different water retention models in three soils.

Porosity models

User defines porosity models in sublist porosity models. It contains as many sublists, e.g. _SOIL1 and _SOIL2, as there are different soils. The porosity models are associated with non-overlapping regions. Each of the sublists (e.g. _SOIL1) includes a few mandatory parameters: regions names, model name, and parameters for the selected model.

porosity_model-spec

• “porosity model” [string] specifies a model for the soil. The available models are “compressible” and “constant”.

  • The model “compressible” requires mandatory parameters:

    • “underformed soil porosity”” [double]

    • “reference pressure” [double]. Default value is 101325.0 Pa.

    • “pore compressibility” [double] [Pa^-1]

    Additional physics introduces more parameters:

    • “Biot coefficient” [double] [-]

    • “volumetric thermal dilation” [double] [K^-1]

    • “reference temperature” [double] [K]. Default value is 20 C.

  • The model “constant” requires “value” [double].

<ParameterList name="flow">  <!-- parent list -->
<ParameterList name="porosity models">
  <ParameterList name="_SOIL1">
    <Parameter name="regions" type="Array(string)" value="{_TOP HALF}"/>
    <Parameter name="porosity model" type="string" value="constant"/>
    <Parameter name="value" type="double" value="0.2"/>
  </ParameterList>

  <ParameterList name="_SOIL2">
    <Parameter name="regions" type="Array(string)" value="{_BOTTOM HALF}"/>
    <Parameter name="porosity model" type="string" value="compressible"/>
    <Parameter name="underformed soil porosity" type="double" value="0.2"/>
    <Parameter name="reference pressure" type="double" value="101325.0"/>
    <Parameter name="pore compressibility" type="double" value="1e-8"/>
  </ParameterList>
</ParameterList>
</ParameterList>

In this example, we define two different porosity models in two soils.

Multiscale continuum models

The list multiscale models is the place for various multiscale models. The list is extension of the list water retention models. Its ordered by soil regions and includes parameters for the multiscale, capillary pressure, and relative permebility models. This list is optional.

multiscale_porosity-spec

• “multiscale model” [string] is the model name. Available options are “dual porosity” and “generalized dual porosity”.

• “xxx parameters” [sublist] provides parameters for the model specified by variable “multiscale model”.

• “water retention model” [string] specifies a model for the soil. The available models are “van Genuchten” and “Brooks Corey”. Parameters for each model are described above.

• “relative permeability model” [string] The available options are “Mualem” (default) and “Burdine”.

Dual porosity model

A two-scale porosity model (fracture + matrix) aka dual porosity model. Current naming convention is that the fields used in the single-porosity model correspond now to the fracture continuum. Example: pressure = pressure in the fracture continuum; pressure_msp = pressure in the matrix continuum.

dual_porosity-spec

• “mass transfer coefficient” [double] is the mass transfer coefficient.

• “tolerance” [double] defines tolerance for iterative methods used to solve secondary equations. Default is 1e-8.

Generalized dual porosity model

A two-scale porosity model (fracture + matrix) aka dual porosity model. Current naming convention is that the fields used in the single-porosity model correspond now to the fracture continuum. Example: pressure = pressure in the fracture continuum; pressure_msp = pressure in the matrix continuum.

generalized_dual_porosity-spec

• “number of matrix nodes” [int] defines number of matrix layers.

• “matrix depth” [double] is the characteristic length for matrix continuum.

• “matrix volume fraction” [double] defines relative volume of matrix continuum.

<ParameterList name="flow">  <!-- parent list -->
<ParameterList name="multiscale models">
  <ParameterList name="_SOIL1">
    <Parameter name="regions" type="Array(string)" value="{_TOP HALF}"/>
    <Parameter name="multiscale model" type="string" value="dual porosity"/>
    <ParameterList name="dual porosity parameters">
      <Paramater name="mass transfer coefficient" type="double" value="4e-5"/>
      <Paramater name="tolerance" type="double" value="1e-8"/>
    </ParameterList>

    <Parameter name="water retention model" type="string" value="van Genuchten"/>
    <Parameter name="van Genuchten alpha" type="double" value="0.000194"/>
    <Parameter name="van Genuchten m" type="double" value="0.28571"/>
    <Parameter name="van Genuchten l" type="double" value="0.5"/>
    <Parameter name="residual saturation" type="double" value="0.103"/>
    <Parameter name="relative permeability model" type="string" value="Mualem"/>
  </ParameterList>
</ParameterList>
</ParameterList>

Permeability

The available options are

permeability-spec

• “coordinate system” [string] defines coordinate system for calculating absolute permeability. The available options are “cartesian” and “layer”.

• “off-diagonal components” [int] defines additional (typically off-diagonal) components of the absolute permeability. Default is 0.

<ParameterList name="flow">  <!-- parent list -->
<ParameterList name="absolute permeability">
  <Parameter name="coordinate system" type="string" value="cartesian"/>
  <Parameter name="off-diagonal components" type="int" value="0"/>
</ParameterList>
</ParameterList>

** DOC GENERATION ERROR: file not found ‘ Permeability ‘ **

Kozeny-Carman

** DOC GENERATION ERROR: file not found ‘ Permeability_KozenyCarman ‘ **

Power law

** DOC GENERATION ERROR: file not found ‘ Permeability_PowerLaw ‘ **

Relative permeability

This section discusses interface treatment of cell-centered fields such as relative permeability, density and viscosity.

rel_perm-spec

• “relative permeability” [list] collects information required for treatment of relative permeability, density and viscosity on mesh faces.

  • “upwind method” [string] defines a method for calculating the upwinded relative permeability. The available options are: “upwind: gravity”, “upwind: darcy velocity” (default), “upwind: second-order”, “upwind: amanzi” (experimental), “upwind: amanzi new” (experiemental), “other: harmonic average”, and “other: arithmetic average”.

  • “upwind frequency” [string] defines frequency of recalculating Darcy flux inside nonlinear solver. The available options are “every timestep” and “every nonlinear iteration”. The first option freezes the Darcy flux for the whole time step. The second option updates it on each iteration of a nonlinear solver. The second option is recommended for the Newton solver. It may impact significantly upwinding of the relative permeability and convergence rate of this solver.

  • “upwind parameters” [list] defines parameters for upwind method specified by “relative permeability”.

    • “tolerance” [double] specifies relative tolerance for almost zero local flux. In such a case the flow is assumed to be parallel to a mesh face. Default value is 1e-12.

    • “method” [string] specifies a reconstruction method. Available option is “cell-based” (default).

    • “polynomial order” [int] defines the polynomial order of a reconstructed function. Default is 1.

    • “limiter” [string] specifies limiting method for a high-order reconstruction. Available options are “Barth-Jespersen” (default), “Michalak-Gooch”, “tensorial”, and “Kuzmin”.

<ParameterList name="flow">  <!-- parent list -->
<ParameterList name="relative permeability">
  <Parameter name="upwind method" type="string" value="upwind: darcy velocity"/>
  <Parameter name="upwind frequency" type="string" value="every timestep"/>
  <ParameterList name="upwind parameters">
     <Parameter name="tolerance" type="double" value="1e-12"/>
     <Parameter name="method" type="string" value="cell-based"/>
     <Parameter name="polynomial order" type="int" value="1"/>
     <Parameter name="limiter" type="string" value="Barth-Jespersen"/>
  </ParameterList>
</ParameterList>
</ParameterList>

Fracture permeability models

A list of permeability models in a fracture network contains similar sublists that must cover all network. Each sublist has two paremeters.

• “region” [string] defines region where model applies.

• “model” [string] specifies the model name. Currently only one parameter is available, “cubic law”.

<ParameterList name="flow">  <!-- parent list -->
<ParameterList name="fracture permeability models">
  <ParameterList name="_ONE FRACTURE LEAVE">
    <Parameter name="region" type="string" value="fracture"/>
    <Parameter name="model" type="string" value="cubic law"/>
  </ParameterList>
</ParameterList>

Diffusion operators

List operators describes the PDE structure of the flow, specifies a discretization scheme, and selects assembling schemas for matrices and preconditioners.

diffusion_operator-spec

• “operators” [list]

  • “diffusion operator” [list] defines parameters for generating and assembling diffusion matrix.

    • “matrix” [list] defines parameters for generating and assembling diffusion matrix. See section describing operators. When the Richards problem is set up, Flow PK sets up proper value for parameter “upwind method” of this sublist.

    • “preconditioner” [list] defines parameters for generating and assembling diffusion matrix that is used to create preconditioner. This sublist is ignored for the saturated problem. Since update of preconditioner can be lagged, we need two objects called “matrix” and “preconditioner”. When the Richards problem is set up, Flow PK sets up proper value for parameter “upwind method” of this sublist.

<ParameterList name="flow">  <!-- parent list -->
<ParameterList name="operators">
  <ParameterList name="diffusion operator">
    <ParameterList name="matrix">
      <Parameter name="discretization primary" type="string" value="mfd: optimized for monotonicity"/>
      <Parameter name="discretization secondary" type="string" value="mfd: optimized for sparsity"/>
      <Parameter name="schema" type="Array(string)" value="{face, cell}"/>
      <Parameter name="preconditioner schema" type="Array(string)" value="{face}"/>
      <Parameter name="gravity" type="bool" value="true"/>
      <Parameter name="gravity term discretization" type="string" value="hydraulic head"/>
    </ParameterList>
    <ParameterList name="preconditioner">
      <Parameter name="discretization primary" type="string" value="mfd: optimized for monotonicity"/>
      <Parameter name="discretization secondary" type="string" value="mfd: optimized for sparsity"/>
      <Parameter name="schema" type="Array(string)" value="{face, cell}"/>
      <Parameter name="preconditioner schema" type="Array(string)" value="{face, cell}"/>
      <Parameter name="Newton correction" type="string" value="approximate Jacobian"/>
    </ParameterList>
  </ParameterList>
</ParameterList>
</ParameterList>

This example creates a p-lambda system, i.e. the pressure is discretized in mesh cells and on mesh faces. The preconditioner is defined on faces only, i.e. cell-based unknowns are eliminated explicitly and the preconditioner is applied to the Schur complement.

Boundary conditions

Boundary conditions are defined in sublist boundary conditions. Four types of boundary conditions are supported. Each type has a similar structure: a list of identical elements that contain information about a part of the boundary where it is prescribed, a function to calculate it, and optional parameters to modify it slightly. This modification is referred to as a submodel and requires additional parameters as described below.

boundary_conditions-spec

• “pressure” [list] is the Dirichlet boundary condition where the pressure is prescribed on a part of the boundary surface. No submodels is available.

• “mass flux” [list] is the Neumann boundary condition where an outward mass flux is prescribed on a part of the boundary surface. This is the default boundary condition. If no condition is specified on a mesh face, the zero flux boundary condition is used.

  • “rainfall” [bool] indicates the submodel where the mass flux is defined with respect to the gravity vector and the actual flux depends on the boundary slope. Default is “false”.

• “static head” [list] is the Dirichlet boundary condition where the hydrostatic pressure is prescribed on a part of the boundary.

  • “relative to top” [bool] indicates the submodel where the static head is defined with respect to the top boundary (a curve in 3D) of the specified regions. Support of 2D is turned off. Default value is “false”.

  • “relative to bottom” [bool] indicates the submodel where the static head is defined with respect to the bottom boundary (a curve in 3D) of the specified regions. Support of 2D is turned off. Default value is “false”.

  • “no flow above water table” [bool] indicates the submodel where the no-flow boundary condition has to be used above the water table. This switch uses the pressure value at a face centroid. Default is “false”.

• “seepage face” [list] is the seepage face boundary condition, a dynamic combination of the “pressure” and “mass flux” boundary conditions over the specified region. The atmospheric pressure is prescribed if internal pressure is higher it. Otherwise, the outward mass flux is prescribed.

  • “reference pressure” [double] defaults to the atmospheric pressure.

  • “rainfall” [bool] indicates the submodel where the mass flux is defined with respect to the gravity vector and the actual influx depends on the boundary slope. Default is “false”.

  • “submodel” [string] indicates different models for the seepage face boundary condition. It can take values “PFloTran” and “FACT”. The first option leads to a discontinuous change of the boundary condition type from the infiltration to pressure. The second option is described in the document on mathematical models. It employs a smooth transition from the infiltration to mixed boundary condition. The recommended value is “PFloTran”.

  • “seepage flux threshold” [double] sets up the threshold for switching from the pressure to influx boundary condition in submodel “PFloTran”. The pressure condition remains for a small influx value until it exceeds the certain fraction of the “mass flux” specified by this parameter. The admissible range is from 0 to 0.1. Default value is 0.

Each boundary condition accepts three parameters: “regions”, “use area fractions”, and “spatial distribution method”. Parameter “regions” specifies the list of regions where the boundary condition is defined. The boolean parameter “use area fractions” instructs the code to use all available volume fractions. Default value is false, it corresponds to \(f=1\) in the formulas below. Parameter “spatial distribution method” defines the method for distributing data (e.g. the total mass flux) over the specified regions. The available options are “volume”, “permeability”, “domain coupling”, “subgrid”, “simple well”, or “none”. For instance, for a given boundary function \(g(x)\), these options correspond to different boundary conditions for the Darcy velocity in the original PDE:

\[{\boldsymbol q} \cdot {\boldsymbol n} = g(x)\, f\, \frac{1}{|B|},\quad\mbox{and}\quad {\boldsymbol q} \cdot {\boldsymbol n} = g(x)\, f,\]

where \(f\) is the volume fraction function, and \(|B|\) is the area of the specified regions calculated using the volume fraction function.

<ParameterList name="flow">  <!-- parent list -->
<ParameterList name="boundary conditions">
  <ParameterList name="pressure">
    <ParameterList name="_BC 0">
      <Parameter name="regions" type="Array(string)" value="{_WEST_SIDE}"/>
      <Parameter name="spatial distribution method" type="string" value="none"/>
      <ParameterList name="boundary pressure">
        <ParameterList name="function-constant">
          <Parameter name="value" type="double" value="101325.0"/>
        </ParameterList>
      </ParameterList>
    </ParameterList>
  </ParameterList>

  <ParameterList name="mass flux">
    <ParameterList name="_BC 1">
      <Parameter name="regions" type="Array(string)" value="{_NORTH_SIDE, _SOUTH_SIDE}"/>
      <Parameter name="spatial distribution method" type="string" value="volume"/>
      <Parameter name="rainfall" type="bool" value="false"/>
      <ParameterList name="outward mass flux">
        <ParameterList name="function-constant">
          <Parameter name="value" type="double" value="0.0"/>
        </ParameterList>
      </ParameterList>
    </ParameterList>
  </ParameterList>

  <ParameterList name="static head">
    <ParameterList name="_BC 2">
      <Parameter name="regions" type="Array(string)" value="{_EAST_SIDE_TOP}"/>
      <Parameter name="spatial distribution method" type="string" value="none"/>
      <Parameter name="relative to top" type="bool" value="true"/>
      <Parameter name="relative to bottom" type="bool" value="true"/>
      <ParameterList name="static head">
        <ParameterList name="function-static-head">
          <Parameter name="p0" type="double" value="101325.0"/>
          <Parameter name="density" type="double" value="1000.0"/>
          <Parameter name="gravity" type="double" value="9.8"/>
          <Parameter name="space dimension" type="int" value="3"/>
          <ParameterList name="water table elevation">
            <ParameterList name="function-constant">
              <Parameter name="value" type="double" value="10.0"/>
            </ParameterList>
          </ParameterList>
        </ParameterList>
      </ParameterList>
    </ParameterList>
  </ParameterList>

  <ParameterList name="seepage face">
    <ParameterList name="_BC 3">
      <Parameter name="regions" type="Array(string)" value="{_EAST_SIDE_BOTTOM}"/>
      <Parameter name="spatial distribution method" type="string" value="none"/>
      <Parameter name="rainfall" type="bool" value="true"/>
      <Parameter name="submodel" type="string" value="PFloTran"/>
      <Parameter name="reference pressure" type="double" value="101325.0"/>
      <ParameterList name="outward mass flux">
        <ParameterList name="function-constant">
          <Parameter name="value" type="double" value="1.0"/>
        </ParameterList>
      </ParameterList>
    </ParameterList>
  </ParameterList>
</ParameterList>
</ParameterList>

This example includes all four types of boundary conditions. The boundary of a square domain is split into six pieces. Constant function is used for simplicity and can be replaced by any of the other available functions.

Sources and sinks

The sources and sinks are typically associated with wells. Negative source means a producing well. Positive source means an injecting well. The structure of list source terms mimics that of list boundary conditions. Again, constant functions can be replaced by any of the available functions.

flow_source-spec

• “regions” [Array(string)] is the list of regions where the source is defined.

• “spatial distribution method” [string] is the method for distributing source Q over the specified regions. The available options are “volume”, “none”, “permeability” and “simple well”. For option “none”, the source term function Q is measured in [kg/m^3/s]. For the other options, it is measured in [kg/s]. When the source function is defined over a few regions, Q is distributed over their union. Option “volume fraction” can be used when the region geometric model support volume fractions. Option “simple well” implements the Peaceman model. The well flux is defined as \(q_w = W (p - p_w)\) [kg/s], where \(W\) is the well index and \(p_w\) is the well pressure. The pressure in a well is assumed to be hydrostatic.

• “use volume fractions” instructs the code to use all available volume fractions. Note that the region geometric model supports volume fractions only for a few regions.

• “submodel” [string] refines definition of the source. Available options are

  • “rate” defines the source in a natural way as the rate of change, Q. It requires a rate function.

  • “integrated source” defines the indefinite integral, I, of the rate of change, i.e. the source term is calculated as Q = dI/dt. It requires a function of I.

  • “bhp” stands for the bottom hole pressure. It requires “depth”, “well radius” and “bhp” function.

  Default submodel is “rate”.

  For most distributions methods, two submodels are supported: “rate” and “integrated source”. For distribution method “simple well”, the following two submodels are supported: “rate” and “bhp”.

<ParameterList name="flow">  <!-- parent list -->
<ParameterList name="source terms">
  <ParameterList name="_SRC 0">
    <Parameter name="regions" type="Array(string)" value="{{_WELL_EAST}}"/>
    <Parameter name="spatial distribution method" type="string" value="volume"/>
    <Parameter name="submodel" type="string" value="rate"/>
    <ParameterList name="well">
      <ParameterList name="function-constant">
        <Parameter name="value" type="double" value="-0.1"/>
      </ParameterList>
    </ParameterList>
  </ParameterList>

  <ParameterList name="_SRC 1">
    <Parameter name="regions" type="Array(string)" value="{{_WELL_WEST}}"/>
    <Parameter name="spatial distribution method" type="string" value="permeability"/>
    <ParameterList name="well">
      <ParameterList name="function-constant">
        <Parameter name="value" type="double" value="-0.2"/>
      </ParameterList>
    </ParameterList>
  </ParameterList>

  <ParameterList name="_SRC 2">
    <Parameter name="regions" type="Array(string)" value="{{_WELL_NORTH}}"/>
    <Parameter name="spatial distribution method" type="string" value="simple well"/>
    <ParameterList name="well">
      <Parameter name="submodel" type="string" value="bhp"/>
      <Parameter name="depth" type="double" value="-2.5"/>
      <Parameter name="well radius" type="double" value="0.1"/>
      <ParameterList name="bhp">
        <ParameterList name="function-constant">
          <Parameter name="value" type="double" value="10.0"/>
        </ParameterList>
      </ParameterList>
    </ParameterList>
  </ParameterList>

  <ParameterList name="_SRC 3">
    <Parameter name="regions" type="Array(string)" value="{{_WELL_SOUTH}}"/>
      <Parameter name="spatial distribution method" type="string" value="simple well"/>
      <ParameterList name="well">
        <Parameter name="submodel" type="string" value="rate"/>
        <ParameterList name="rate">
          <ParameterList name="function-constant">
            <Parameter name="value" type="double" value="100.0"/>
          </ParameterList>
        </ParameterList>
      </ParameterList>
    </ParameterList>
  </ParameterList>
</ParameterList>
</ParameterList>

Time integrator

The list time integrator defines a generic time integrator used by the cycle driver. This driver assumes that each PK has only one time integrator. The list time integrator defines parameters controlling linear and nonlinear solvers during a time integration period. We break this long sublist into smaller parts.

Initialization and constraints

• “error control options” [Array(string)] lists various error control options. A nonlinear solver is terminated when all listed options are passed. The available options are “pressure”, “saturation”, and “residual”. All errors are relative, i.e. dimensionless. The error in pressure is compared with capillary pressure plus atmospheric pressure. The other two errors are compared with 1. The option “pressure” is always active during steady-state time integration. The option “saturation” is always active during transient time integration.

• “linear solver” [string] refers to a generic linear solver from list “solvers”. It is used in all cases except for “initialization” and “dae constraint”. Currently, it is used by the Darcy PK only.

• “preconditioner” [string] specifies preconditioner for linear and nonlinear solvers.

• “preconditioner enhancement” [string] specifies a linear solver that binds the above preconditioner to improve spectral properties. Default is “none”.

• “initialization” [list] defines parameters for calculating initial pressure guess. It can be used to obtain pressure field which is consistent with the boundary conditions. Default is empty list.

  • “method” [string] specifies an optional initialization methods. The available options are “picard” and “saturated solver”. The latter option leads to solving a Darcy problem. The former option uses sublist “picard parameters”. Picard works better if a bounded initial pressure guess is provided.

  • “active wells” [bool] specifies if wells are active or turned off. Default is false.

  • “picard parameters” [list] defines control parameters for the Picard solver.

    • “convergence tolerance” [double] specifies nonlinear convergence tolerance. Default is 1e-8.

    • “maximum number of iterations” [int] limits the number of iterations. Default is 400.

  • “linear solver” [string] refers to a solver sublist of the list “solvers”.

  • “clipping saturation value” [double] is an experimental option. It is used after pressure initialization to cut-off small values of pressure. The new pressure is calculated based of the provided saturation value. Default is 0.6.

  • “clipping pressure value” [double] is an experimental option. It is used after pressure initialization to cut-off small values of pressure below the provided value.

• “dae constraint” [list] each time the time integrator is restarted, we need to re-enforce the pressure-lambda relationship for new boundary conditions. Default is empty list.

  • “method” [string] is a placeholder for different algorithms. Now, the only available option is “projection” which is default.

  • “linear solver” [string] refers to a solver sublist of the list “solvers”.

  • “inflow krel correction” [bool] estimates relative permeability on inflow mesh faces. This estimate is more reliable than the upwinded relative permeability value, especially in steady-state calculations.

<ParameterList name="flow">  <!-- parent list -->
<ParameterList name="time integrator">
  <Parameter name="error control options" type="Array(string)" value="{pressure, saturation}"/>
  <Parameter name="linear solver" type="string" value="_GMRES_WITH_AMG"/>
   <Parameter name="preconditioner" type="string" value="_HYPRE_AMG"/>
   <Parameter name="preconditioner enhancement" type="string" value="none"/>

   <ParameterList name="initialization">  <!-- first method -->
     <Parameter name="method" type="string" value="saturated solver"/>
     <Parameter name="linear solver" type="string" value="_PCG_WITH_AMG"/>
     <Parameter name="clipping pressure value" type="double" value="50000.0"/>
   </ParameterList>

   <ParameterList name="initialization">  <!-- alternative method -->
     <Parameter name="method" type="string" value="picard"/>
     <Parameter name="linear solver" type="string" value="_PCG_WITH_AMG"/>
     <ParameterList name="picard parameters">
       <Parameter name="convergence tolerance" type="double" value="1e-8"/>
       <Parameter name="maximum number of iterations" type="int" value="20"/>
     </ParameterList>
   </ParameterList>

   <ParameterList name="dae constraint">
     <Parameter name="method" type="string" value="projection"/>
     <Parameter name="inflow krel correction" type="bool" value="false"/>
     <Parameter name="linear solver" type="string" value="_PCG_WITH_AMG"/>
   </ParameterList>
 </ParameterList>
 </ParameterList>

Time step controller and nonlinear solver

The time step is controlled by parameter time step controller type and the related list of options, see section TimeStepController for the list of supported parameter. Nonlinear solver is controlled by parameter solver type and related list of options. Amanzi supports a few nonlinear solvers described in details in a separate section.

<ParameterList name="flow">  <!-- parent list -->
<ParameterList name="time integrator">
  <Parameter name="time integration method" type="string" value="BDF1"/>
  <ParameterList name="BDF1">
    <Parameter name="max preconditioner lag iterations" type="int" value="5"/>
    <Parameter name="extrapolate initial guess" type="bool" value="true"/>
    <Parameter name="restart tolerance relaxation factor" type="double" value="1000.0"/>
    <Parameter name="restart tolerance relaxation factor damping" type="double" value="0.9"/>

    <Parameter name="timestep controller type" type="string" value="standard"/>
    <ParameterList name="timestep controller standard parameters">
      <Parameter name="min iterations" type="int" value="10"/>
      <Parameter name="max iterations" type="int" value="15"/>
      <Parameter name="time step increase factor" type="double" value="1.2"/>
      <Parameter name="time step reduction factor" type="double" value="0.5"/>
      <Parameter name="max time step" type="double" value="1e+9"/>
      <Parameter name="min time step" type="double" value="0.0"/>
    </ParameterList>

    <Parameter name="solver type" type="string" value="nka"/>
    <ParameterList name="nka parameters">
      <Parameter name="nonlinear tolerance" type="double" value="1e-5"/>
      <Parameter name="limit iterations" type="int" value="30"/>
      <Parameter name="diverged tolerance" type="double" value="1e+10"/>
      <Parameter name="diverged l2 tolerance" type="double" value="1e+10"/>
      <Parameter name="diverged pc tolerance" type="double" value="1e+10"/>
      <Parameter name="max du growth factor" type="double" value="1e+5"/>
      <Parameter name="max divergent iterations" type="int" value="3"/>
      <Parameter name="max nka vectors" type="int" value="10"/>
      <Parameter name="modify correction" type="bool" value="false"/>
      <ParameterList name="verbose object">
        <Parameter name="verbosity level" type="string" value="high"/>
      </ParameterList>
    </ParameterList>

    <!-- alternative solver
    <Parameter name="solver type" type="string" value="aa"/>
    <ParameterList name="aa parameters">
      <Parameter name="nonlinear tolerance" type="double" value="1e-5"/>
      <Parameter name="limit iterations" type="int" value="30"/>
      <Parameter name="diverged tolerance" type="double" value="1e+10"/>
      <Parameter name="diverged l2 tolerance" type="double" value="1e+10"/>
      <Parameter name="diverged pc tolerance" type="double" value="1e+10"/>
      <Parameter name="max du growth factor" type="double" value="1e+5"/>
      <Parameter name="max divergent iterations" type="int" value="3"/>
      <Parameter name="max aa vectors" type="int" value="10"/>
      <Parameter name="modify correction" type="bool" value="false"/>
      <Parameter name="relaxation parameter" type="double" value="0.75"/>
    </ParameterList-->
  </ParameterList>
</ParameterList>
</ParameterList>

In this example, the time step is increased by factor 1.2 when the nonlinear solver converges in 10 or less iterations. The time step is not changed when the number of nonlinear iterations is between 11 and 15. The time step will be cut twice if the number of nonlinear iterations exceeds 15.

Other parameters

The remaining flow parameters are

• “clipping parameters” [list] defines how solution increment calculated by a nonlinear solver is modified e.g., clipped.

  • “maximum saturation change” [double] Default is 0.25.

  • “pressure damping factor” [double] Default is 0.5.

• “plot time history” [bool] produces an ASCII file with the time history. Default is “false”.

• “algebraic water storage balance” [bool] uses algebraic correction to enforce consistency of water storage and Darcy fluxes. It leads to a monotone transport. Default is false.

<ParameterList name="flow">  <!-- parent list -->
<ParameterList name="clipping parameters">
   <Parameter name="maximum saturation change" type="double" value="0.25"/>
   <Parameter name="pressure damping factor" type="double" value="0.5"/>
</ParameterList>

<Parameter name="plot time history" type="bool" value="false"/>
<Parameter name="algebraic water storage balance" type="bool" value="false"/>
</ParameterList>

Explanation of verbose output

When verbosity is set to high, this PK reports information about current status of the simulation. Here after keyword global refers to the whole simulation including all time periods, keyword local refers to the current time period. The incomplete list is

• [global] cycle number, time before the step, and time step dt (in years)

• [local] step number, time T, and dT inside the time integrator (in seconds)

• [local] frequency of preconditioner updates

• [local] number of performed nonlinear steps and value of the nonlinear residual

• [local] total number of successful time steps (TS), failed time steps (FS), preconditioner updates (PC/1) and preconditioner applies (PC/2), linear solves insides preconditioner (LS)

• [local] amount of liquid (water) in the reservoir and amount of water entering and living domain through its boundary (based on darcy flux).

• [global] current simulation time (in years)

CycleDriver      |   Cycle 40: time(y) = 0.953452, dt(y) = 0.238395
TI::BDF1         |    step 40 T = 3.00887e+07 [sec]  dT = 7.52316e+06
TI::BDF1         |    preconditioner lag is 20 out of 20
TI::BDF1         |    success: 4 nonlinear itrs error=7.87642e-08
TI::BDF1         |    TS:40 FS:0 NS:64 PC:42 64 LS:0 dt:1.0000e+03 7.5232e+06
FlowPK::Richards |    reservoir water mass=1.36211e+06 [kg], total influx=897.175 [kg]
CycleDriver      |   New time(y) = 1.19185

Transport PK

Mathematical models

A few PDE models can be instantiated using the parameters described below.

Single-phase transport

The conceptual PDE model for the transport in partially saturated media is

\[\frac{\partial (\phi s_l C)}{\partial t} = - \boldsymbol{\nabla} \cdot (\boldsymbol{q}_l C) + \boldsymbol{\nabla} \cdot (\phi_e s_l\, (\boldsymbol{D}_l + \tau \boldsymbol{M}_l) \boldsymbol{\nabla} C) + Q,\]

where \(\phi\) is total porosity [-], \(\phi_e\) is effective transport porosity [-], \(s_l\) is liquid saturation [-], \(Q\) is source or sink term, \(\boldsymbol{q}_l\) is the Darcy velocity [m/s], \(\boldsymbol{D}_l\) is dispersion tensor, \(\boldsymbol{M}_l\) is diffusion coefficient, and \(\tau\) is tortuosity [-]. For an isotropic medium with no preferred axis of symmetry the dispersion tensor has the following form:

\[\boldsymbol{D}_l = \alpha_t \|\boldsymbol{v}\| \boldsymbol{I} + \left(\alpha_l-\alpha_t \right) \frac{\boldsymbol{v} \boldsymbol{v}}{\|\boldsymbol{v}\|}, \qquad \boldsymbol{v} = \frac{\boldsymbol{q}}{\phi_e}\]

where \(\alpha_l\) is longitudinal dispersivity [m], \(\alpha_t\) is transverse dispersivity [m], and \(\boldsymbol{v}\) is average pore velocity [m/s]. Amanzi supports two additional models for dispersivity with 3 and 4 parameters.

Single-phase transport with dual porosity model

The dual porosity formulation of the solute transport consists of two equations for the fracture and matrix regions. In the fracture region, we have citep{simunek-vangenuchten_2008}

\[\begin{split}\begin{array}{rcl} \frac{\partial (\phi_f\, s_{lf}\, C_f)}{\partial t} &=& - \boldsymbol{\nabla} \cdot (\boldsymbol{q}_l C_f) + \boldsymbol{\nabla} \cdot (\phi_f\, s_{lf}\, (\boldsymbol{D}_l + \tau_f M) \boldsymbol{\nabla} C_f)\\ && - \displaystyle\frac{\phi_m\,\tau_m}{L_m}\, M \nabla C_m - \Sigma_w C^* + Q_f, \end{array}\end{split}\]

where \(\phi_f\) is fracture porosity [-], \(\phi_m\) is matrix porosity [-], \(s_{lf}\) is liquid saturation in fracture [-], \(\boldsymbol{q}_l\) is the Darcy velocity [m/s], \(\boldsymbol{D}_l\) is dispersion tensor, \(\tau_f\) is fracture tortuosity [-], \(\tau_m\) is matrix tortuosity [-], \(M\) is molecular diffusion coefficient [\(m^2/s\)], and \(L_m\) is the characteristic matrix depth defined typically as the ratio of a matrix block [m], \(\Sigma_w\) is transfer rate due to flow from the matrix to the fracture, \(C^*\) is equal to \(C_f\) if \(\Sigma_w > 0\) and \(C_m\) is \(\Sigma_w < 0\), and \(Q_f\) is source or sink term. In the matrix region, we have

\[\frac{\partial (\phi_m\, s_{lm}\, C_m)}{\partial t} = \nabla\cdot (\phi_m\, \tau_m\, M_m \nabla C_m) + \Sigma_w C^* + Q_m,\]

where \(\phi_m\) is matrix porosity [-], \(s_{lm}\) is liquid saturation in matrix [-], \(Q_m\) is source or sink term. The simplified one-node dual porosity model uses a finite difference approximation of the solute gradient:

\[\nabla C_m \approx WR \, \frac{C_f - C_m}{L_m},\]

where \(WR\) is the Warren-Root coefficient that estimates the poro-space geometry, [-]

Physical models and assumptions

This list is used to summarize physical models and assumptions, such as coupling with other PKs. This list is often generated or extended by a high-level MPC PK.

transport-spec

• “gas diffusion” [bool] indicates that air-water partitioning coefficients are used to distribute components between liquid and as phases. Default is false.

• “permeability field is required” [bool] indicates if some transport features require absolute permeability. Default is false.

• “multiscale model” [string] specifies a multiscale model. Available options are “single porosity” (default) and “dual porosity”.

• “effective transport porosity” [bool] If true, effective transport porosity will be used by dispersive-diffusive fluxes instead of total porosity. Default is false.

• “eos lookup table” [string] provides the name for optional EOS lookup table.

• “use dispersion solver” [bool] instructs PK to instantiate a solver but do not call it. It is used now by MPC to form a global solver. Default is false.

<ParameterList name="_TRANSPORT">  <!-- parent list -->
<ParameterList name="physical models and assumptions">
  <Parameter name="gas diffusion" type="bool" value="false"/>
  <Parameter name="permeability field is required" type="bool" value="false"/>
  <Parameter name="multiscale model" type="string" value="single porosity"/>
  <Parameter name="effective transport porosity" type="bool" value="false"/>
  <Parameter name="eos lookup table" type="string" value="h2o.eos"/>
  <Parameter name="use dispersion solver" type="bool" value="false"/>
</ParameterList>
</ParameterList>

Global parameters

The transport component of Amanzi performs advection of aqueous and gaseous components and their dispersion and diffusion. The main parameters control temporal stability, spatial and temporal accuracy, and verbosity:

transport_global_params-spec

• “domain name” [string] specifies mesh name that defined domain of this PK. Default is “domain”.

• “cfl” [double] Time step limiter, a number less than 1. Default value is 1.

• “method” [string] defines flux method. Available options are “muscl” (default) and “fct”.

• “spatial discretization order” [int] defines accuracy of spatial discretization. It permits values 1 or 2. Default value is 1.

• “temporal discretization order” [int] defines accuracy of temporal discretization. It permits values 1 or 2 and values 3 or 4. Note that RK3 is not monotone. Default value is 1.

• “reconstruction” [list] collects reconstruction parameters. The available options are describe in the separate section below.

• “solver” [string] Specifies the dispersion/diffusion solver.

• “preconditioner” [string] specifies preconditioner for dispersion solver.

• “number of aqueous components” [int] The total number of aqueous components. Default value is the total number of components.

• “number of gaseous components” [int] The total number of gaseous components. Default value is 0.

<ParameterList>  <!-- parent list -->
<ParameterList name="_TRANSPORT">
  <Parameter name="domain name" type="string" value="domain"/>
  <Parameter name="cfl" type="double" value="1.0"/>
  <Parameter name="method" type="string" value="muscl"/>
  <Parameter name="spatial discretization order" type="int" value="1"/>
  <Parameter name="temporal discretization order" type="int" value="1"/>
  <Parameter name="solver" type="string" value="_PCG_SOLVER"/>

  <ParameterList name="reconstruction">
    <Parameter name="method" type="string" value="cell-based"/>
    <Parameter name="polynomial order" type="int" value="1"/>
    <Parameter name="limiter" type="string" value="tensorial"/>
    <Parameter name="limiter extension for transport" type="bool" value="true"/>
  </ParameterList>

  <ParameterList name="verbose object">
    <Parameter name="verbosity level" type="string" value="high"/>
  </ParameterList>
</ParameterList>
</ParameterList>

Material properties

The material properties include dispersivity model and diffusion parameters for aqueous and gaseous phases. The dispersivity is defined as a soil property. The diffusivity is defined independently for each solute.

mdm-spec

• _SOIL [list] Defines material properties.

  • “region” [Array(string)] Defines geometric regions for material SOIL.

  • “model” [string] Defines dispersivity model, choose exactly one of the following: “scalar”, “Bear”, “Burnett-Frind”, or “Lichtner-Kelkar-Robinson”.

  • “parameters for MODEL” [list] where “MODEL” is the model name. For model “scalar”, only one of the following options must be specified:

    • “alpha” [double] defines dispersivity in all directions, [m].

    • “dispersion coefficient” [double] defines dispersion coefficient [m^2/s].

    For model “Bear”, the following options must be specified:

    • “alpha_l” [double] defines dispersion in the direction of Darcy velocity, [m].

    • “alpha_t” [double] defines dispersion in the orthogonal direction, [m].

    For model “Burnett-Frind”, the following options must be specified:

    • “alphaL” [double] defines the longitudinal dispersion in the direction of Darcy velocity, [m].

    • “alpha_th” [double] Defines the transverse dispersion in the horizonla direction orthogonal directions, [m].

    • “alpha_tv” [double] Defines dispersion in the orthogonal directions, [m]. When “alpha_th” equals to “alpha_tv”, we obtain dispersion in the direction of the Darcy velocity. This and the above parameters must be defined for “Burnett-Frind” and “Lichtner-Kelkar-Robinson” models.

    For model “Lichtner-Kelker-Robinson”, the following options must be specified:

    • “alpha_lh” [double] defines the longitudinal dispersion in the horizontal direction, [m].

    • “alpha_lv” [double] Defines the longitudinal dispersion in the vertical direction, [m]. When “alpha_lh” equals to “alpha_lv”, we obtain dispersion in the direction of the Darcy velocity. This and the above parameters must be defined for “Burnett-Frind” and “Lichtner-Kelker-Robinson” models.

    • “alpha_th” [double] Defines the transverse dispersion in the horizontal direction orthogonal directions, [m].

    • “alpha_tv” ``[double]` Defines dispersion in the orthogonal directions. When “alpha_th” equals to “alpha_tv”, we obtain dispersion in the direction of the Darcy velocity. This and the above parameters must be defined for “Burnett-Frind” and “Lichtner-Kelker-Robinson” models.

  • “aqueous tortuosity” [double] Defines tortuosity for calculating diffusivity of liquid solutes, [-].

  • “gaseous tortuosity” [double] Defines tortuosity for calculating diffusivity of gas solutes, [-].

Three examples are below:

<ParameterList name="_TRANSPORT">  <!-- parent list -->
<ParameterList name="material properties">
  <ParameterList name="_WHITE SOIL">
    <Parameter name="regions" type="Array(string)" value="{_TOP_REGION, _BOTTOM_REGION}"/>
    <Parameter name="model" type="string" value="Bear"/>
    <ParameterList name="parameters for Bear">
      <Parameter name="alpha_l" type="double" value="1e-2"/>
      <Parameter name="alpha_t" type="double" value="1e-5"/>
    <ParameterList>
    <Parameter name="aqueous tortuosity" type="double" value="1.0"/>
    <Parameter name="gaseous tortuosity" type="double" value="1.0"/>
  </ParameterList>

  <ParameterList name="_GREY SOIL">
    <Parameter name="regions" type="Array(string)" value="{_MIDDLE_REGION}"/>
    <Parameter name="model" type="string" value="Burnett-Frind"/>
    <ParameterList name="parameters for Burnett-Frind">
      <Parameter name="alpha_l" type="double" value="1e-2"/>
      <Parameter name="alpha_th" type="double" value="1e-3"/>
      <Parameter name="alpha_tv" type="double" value="2e-3"/>
    <ParameterList>
    <Parameter name="aqueous tortuosity" type="double" value="0.5"/>
    <Parameter name="gaseous tortuosity" type="double" value="1.0"/>
  </ParameterList>
</ParameterList>
</ParameterList>

material_properties-spec

• “molecular diffusion” [list] defines names of solutes in aqueous and gaseous phases and related diffusivity values.

<ParameterList name="_TRANSPORT">  <!-- parent list -->
<ParameterList name="molecular diffusion">
  <Parameter name="aqueous names" type=Array(string)" value="{CO2(l),Tc-99}"/>
  <Parameter name="aqueous values" type=Array(double)" value="{1e-8,1e-9}"/>

  <Parameter name="gaseous names" type=Array(string)" value="{CO2(g)}"/>
  <Parameter name="gaseous values" type=Array(double)" value="{1e-8}"/>
  <Parameter name="air-water partitioning coefficient" type=Array(double)" value="{0.03}"/>
</ParameterList>
</ParameterList>

Dispersion operator

List operators describes the PDE structure of the flow, specifies a discretization scheme, and selects assembling schemas for matrices and preconditioners.

dispersion-spec

• “operators” [list]

  • “diffusion operator” [list]

    • “matrix” [list] defines parameters for generating and assembling dispersion matrix. See section describing operators.

<ParameterList name="_TRANSPORT">  <!-- parent list -->
<ParameterList name="operators">
  <ParameterList name="diffusion operator">
    <ParameterList name="matrix">
      <Parameter name="discretization primary" type="string" value="mfd: optimized for monotonicity"/>
      <Parameter name="discretization secondary" type="string" value="mfd: two-point flux approximation"/>
      <Parameter name="schema" type="Array(string)" value="{face, cell}"/>
      <Parameter name="preconditioner schema" type="Array(string)" value="{face}"/>
    </ParameterList>
  </ParameterList>
</ParameterList>
</ParameterList>

This example creates a p-lambda system, i.e. the concentration is discretized in mesh cells and on mesh faces. The later unknowns are auxiliary unknowns.

Multiscale continuum models

The list of multiscale models is the place for various subscale models that coul be mixed and matched. Its ordered by materials and includes parameters for the assigned multiscale model This list is optional.

transport_multiscale-spec

• “multiscale model” [string] is the model name. Available option is “dual porosity” and “generalized dual porosity”.

• “regions” [Array(string)] is the list of regions where this model should be applied.

• “xxx parameters” [list] provides parameters for the model specified by variable “multiscale model”.

Dual porosity model

A two-scale porosity model (fracture + matrix) aka dual porosity model. Current naming convention is that the fields used in the single-porosity model correspond now to the fracture continuum. Example: tcc = total component concentration in the fracture continuum; tcc_matrix = total component concentration in the matrix continuum.

transport_dual_porosity-spec

• “Warren Root parameter” [list] scales diffusive solute transport due to concentration gradient.

• “tortousity” [double] defines tortuosity to correct diffusivity of a liquid solute.

Generalized dual porosity model

A generalized dual porosity model, fracture + multiple matrix nodes. Current naming convention is that the fields used in the single-porosity model correspond now to the fracture continuum. Example: tcc = total component concentration in the fracture continuum; tcc_matrix = total component concentration in the matrix continuum.

transport_generalized_dual_porosity-spec

• “number of matrix nodes” [int] defines number of matrix layers.

• “matrix depth” [double] is the characteristic length for matrix continuum.

• “tortousity” [double] defines tortuosity to correct diffusivity of a liquid solute.

• “matrix volume fraction” [double] defines relative volume of matrix continuum.

<ParameterList name="_TRANSPORT">  <!-- parent list -->
<ParameterList name="multiscale models">
  <ParameterList name="_WHITE SOIL">
    <Parameter name="multiscale model" type="string" value="dual porosity"/>
    <Parameter name="regions" type="Array(string)" value="{_TOP_REGION, _BOTTOM_REGION}"/>
    <ParameterList name="dual porosity parameters">
      <Paramater name="Warren Root parameter" type="double" value="4.0e-5"/>
      <Paramater name="matrix tortuosity" type="double" value="0.95"/>
      <Paramater name="matrix volume fraction" type="double" value="0.9999"/>
    </ParameterList>
  </ParameterList>

  <ParameterList name="_GREY SOIL">
    <Parameter name="multiscale model" type="string" value="generalized dual porosity"/>
    <Parameter name="regions" type="Array(string)" value="{_MIDDLE_REGION}"/>
    <ParameterList name="generalized dual porosity parameters">
      <Paramater name="number of matrix nodes" type="int" value="2"/>
      <Paramater name="matrix depth" type="double" value="0.01"/>
      <Paramater name="matrix tortuosity" type="double" value="1.0"/>
    </ParameterList>
  </ParameterList>
</ParameterList>
</ParameterList>

Boundary conditions

For the advective transport, the boundary conditions must be specified on inflow parts of the boundary. If no value is prescribed through the XML input, the zero influx boundary condition is used. Note that the boundary condition is set up separately for each component. The structure of boundary conditions is aligned with that used for flow and allows us to define spatially variable boundary conditions.

transport_boundary_conditions-spec

• “boundary conditions” [list]

  • “concentration” [list] This is a reserved keyword.

    • “_SOLUTE” [list] contains a few sublists (e.g. _EAST_CRIB) for boundary conditions. The name _SOLUTE must be the name in the list of solutes.

      • “_BC1” [list] defines boundary conditions using arrays of boundary regions and attached functions.

      • “regions” [Array(string)] defines a list of boundary regions where a boundary condition must be applied.

      • “spatial distribution method” [string] defines the method for distributing data over the specified regions. The available options are “area” or “none”.

      • “boundary concentration” [list] defines a function for calculating boundary conditions. The function specification is described in subsection Functions.

The example below sets constant boundary condition 1e-5 for the duration of transient simulation.

<ParameterList name="_TRANSPORT">  <!-- parent list -->
<ParameterList name="boundary conditions">
  <ParameterList name="concentration">
    <ParameterList name="NO3-">
      <ParameterList name="_EAST_CRIB">   <!-- user defined name -->
        <Parameter name="regions" type="Array(string)" value="{_TOP, _LEFT}"/>
        <Parameter name="spatial distribution method" type="string" value="none"/>
        <ParameterList name="boundary concentration">
          <ParameterList name="function-constant">  <!-- any time function -->
            <Parameter name="value" type="double" value="1e-5"/>
          </ParameterList>
        </ParameterList>
      </ParameterList>
      <ParameterList name="_WEST CRIB">   <!-- user defined name -->
        ...
      </ParameterList>
    </ParameterList>

    <ParameterList name="CO2"> <!-- Next component -->
      ...
    </ParameterList>
  </ParameterList>
</ParameterList>
</ParameterList>

Geochemical boundary condition is the Dirichlet boundary condition which requires calculation of a geochemical balance. Note that the number of time functions below is one less than the number of times and geochemical conditions.

<ParameterList name="boundary conditions">  <!-- parent list -->
<ParameterList name="geochemical">
  <ParameterList name="_EAST_CRIB">
    <Parameter name="solutes" type="Array(string)" value={H+,HCO3-,Ca++}"/>
    <Parameter name="times" type="Array(double)" value="{0.0, 100.0}"/>
    <Parameter name="geochemical conditions" type="Array(string)" value="{cond1, cond2}"/>
    <Parameter name="time functions" type="Array(string)" value="{constant}"/>
    <Parameter name="regions" type="Array(string)" value="{EAST_CRIB}"/>
  </ParameterList>
</ParameterList>
</ParameterList>

Sources and sinks

The sources and sinks are typically located at wells. Stability condition requires to distinguish between injecting and producing wells. A source function used for an injecting well specifies concentration of solute. A source function used for a producing well specifies volumetric flow rate [m^3/s] of water.

The structure of list source terms includes only sublists named after components. Again, constant functions can be replaced by any available time-function. Note that the source values are set up separately for each component.

transport_sources-spec

• “concentration” [list] This is a reserved keyword.

• “_SOLUTE” [list] contains a few sublists (e.g. _SRC1 and _SRC2) for various sources and sinks. The name _SOLUTE must exist in the list of solutes.

• “_SRC1” [list] defines a source using arrays of domain regions, a function, and a distribution method.

• “regions” [Array(string)] defines a list of domain regions where a source term must be applied.

• “spatial distribution method” [string] identifies a method for distributing source Q over the specified regions. The available options are “volume”, “none”, and “permeability”. For option “none” the source term Q is measured in [mol/L/s] (if units for concentration is mol/L) or [mol/m^3/s] (otherwise). For the other options, it is measured in [mol/s]. When the source function is defined over a few regions, Q will be distributed over their union.

• “submodel” [string] refines definition of source. Available options are “rate” and “integrand”. The first option defines rate of change q, the second one defines integrand Q of a rate Q = dq/dt. Default is “rate”.

• “sink” [list] is a function for calculating a source term. The function specification is described in subsection Functions.

This example defines one well and one sink.

<ParameterList name="source terms"> <!-- parent list -->
<ParameterList name="concentration">
  <ParameterList name="H+">
    <ParameterList name="_SOURCE: EAST WELL">   <!-- user defined name -->
      <Parameter name="regions" type="Array(string)" value="{_EAST_WELL}"/>
      <Parameter name="spatial distribution method" type="string" value="volume"/>
      <Parameter name="submodel" type="string" value="rate"/>
      <ParameterList name="injector">   <!-- reserved keyword -->
        <ParameterList name="function-constant">
          <Parameter name="value" type="double" value="0.01"/>
        </ParameterList>
      </ParameterList>
    </ParameterList>
    <ParameterList name="_SOURCE: WEST_WELL">
       ...
    </ParameterList>
  </ParameterList>

  <ParameterList name="CO2(g)">   <!-- next component, a gas -->
    <ParameterList name="_SOURCE: WEST WELL">   <!-- user defined name -->
      <Parameter name="regions" type="Array(string)" value="{_WEST_WELL}"/>
      <Parameter name="spatial distribution method" type="string" value="permeability"/>
      <ParameterList name="injector">
        <ParameterList name="function-constant">
          <Parameter name="value" type="double" value="0.02"/>
        </ParameterList>
      </ParameterList>
    </ParameterList>
  </ParameterList>
</ParameterList>
</ParameterList>

Developer parameters

The remaining parameters that can be used by a developer include

• “enable internal tests” [bool] turns on various internal tests during run time. Default value is “false”.

• “internal tests tolerance” [double] tolerance for internal tests such as the divergence-free condition. The default value is 1e-6.

• “runtime diagnostics: solute names” [Array(string)] defines solutes that will be tracked closely each time step if verbosity “high”. Default value is the first solute in the global list of “aqueous names” and the first gas in the global list of “gaseous names”.

• “runtime diagnostics: regions” [Array(string)] defines a boundary region for tracking solutes. Default value is a seepage face boundary, see Flow PK.

Explanation of verbose output

When verbosity is set to high, this PK reports information about current status of the simulation. Here after keyword global refers to the whole simulation including all time periods, keyword local refers to the current time period. The incomplete list is

• [global] cycle number, time before step, and time step dt (in years)

• [local] cell id and position with the smallest time step

• [local] convergence of a linear solver for dispersion, PCG here

• [local] number of subcycles, stable time step, and global time step (in seconds)

• [local] species’s name, concentration extrema, and total amount of it in the reservoir

• [global] current simulation time (in years)

CycleDriver      |   Cycle 10: time(y) = 0.803511, dt(y) = 0.089279
TransportPK      |    cell 0 has smallest dt, (-270, -270)
TransportPK      |    dispersion solver (pcg) ||r||=8.33085e-39 itrs=2
TransportPK      |    1 sub-cycles, dt_stable=2.81743e+06 [sec]  dt_MPC=2.81743e+06 [sec]
TransportPK      |    Tc-99: min=8.08339e-06 mol/L max=0.0952948 mol/L, total=9.07795 mol
CycleDriver      |   New time(y) = 0.89279

Chemistry PK

The chemistry header includes three parameters:

chemistry_pk-spec

• “chemistry model” [string] defines chemical model. The available options are “Alquimia” and “Amanzi” (default).

<ParameterList name="_CHEMISTRY">
  <Parameter name="chemistry model" type="string" value="Amanzi"/>
</ParameterList>

Geochemical engines

Here we specify either the default or the third-party geochemical engine.

Common parameters

The following parameters are common for all supported engines.

chemistry_params-spec

• “timestep control method” [string] specifies timestep control method for chemistry subcycling. Choose either “fixed” (default) or “simple”. For option “fixed”, timestep is fixed. For option “simple”, the timestep is adjusted in response to stiffness of system of equations based on a simple scheme. This option require the following parameters: “timestep cut threshold”, “timestep cut factor”, “timestep increase threshold”, and “timestep increase factor”.

• “timestep cut threshold” [int] is the number of Newton iterations that if exceeded will trigger a timestep cut. Default is 8.

• “max timestep (s)” [double] is the maximum timestep that chemistry will allow the MPC to take.

• “initial timestep (s)” [double] is the initial timestep that chemistry will ask the MPC to take.

• “timestep cut factor” [double] is the factor by which the timestep is cut. Default is 2.0

• “timestep increase threshold” [int] is the number of consecutive successful timesteps that will trigger a timestep increase. Default is 4.

• “timestep increase factor” [double] is the factor by which the timestep is increased. Default is 1.2

• “free ion initial guess” [double] provides an estimate of the free ion concentration for solutes. It used to help convergence of the initial solution of the chemistry. If this parameter is absent, a fraction (10%) of the total component concentration is used.

• “initial conditions time” [double] specifies time for applying initial conditions. This parameter is useful for simulation restart. Default value is the state time when chemistry PK is instantiated.

Alquimia

The Alquimia chemistry process kernel only requires the Engine and Engine Input File entries, but will also accept and respect the value given for max timestep (s). Most details are provided in the trimmed PFloTran file 1d-tritium-trim.in.

alquimia-spec

• “minerals” [Array(string)] is the list of mineral names.

• “sorption sites” [Array(string)] is the list of sorption sites.

• “auxiliary data” [Array(string)] defines additional chemistry related data that the user can request be saved to vis files.

• “min timestep (s)” [double] is the minimum timestep that chemistry will allow the MPC to take.

<ParameterList>  <!-- parent list -->
<ParameterList name="_CHEMISTRY">
  <Parameter name="engine" type="string" value="PFloTran"/>
  <Parameter name="engine input file" type="string" value="_TRITIUM.in"/>
  <Parameter name="minerals" type="Array(string)" value="{quartz, kaolinite, goethite, opal}"/>
  <Parameter name="min timestep (s)" type="double" value="1.5778463e-07"/>
  <Parameter name="max timestep (s)" type="double" value="1.5778463e+07"/>
  <Parameter name="initial timestep (s)" type="double" value="1.0e-02"/>
  <Parameter name="timestep control method" type="string" value="simple"/>
  <Parameter name="timestep cut threshold" type="int" value="8"/>
  <Parameter name="timestep cut factor" type="double" value="2.0"/>
  <Parameter name="timestep increase threshold" type="int" value="4"/>
  <Parameter name="timestep increase factor" type="double" value="1.2"/>
</ParameterList>
</ParameterList>

Amanzi

The Amanzi chemistry process kernel uses the following parameters.

chemistry_amanzi-spec

• “thermodynamic database” [list]

  • “file” [string] is the name of the chemistry database file, relative to the execution directory.

  • “format” [string] is the format of the database file. Actual database format is not XML and is the same as described for the 2010 demo with additions for the new chemical processes. Valid values: “simple”.

• “minerals” [Array(string)] is the list of mineral names.

• “sorption sites” [Array(string)] is the list of sorption sites.

• “activity model” [string] is the type of model used for activity corrections. Valid options are “unit”, “debye-huckel”, and “pitzer-hwm”,

• “tolerance” [double] defines tolerance in Newton solves inside the chemistry library.

• “maximum Newton iterations” [int] is the maximum number of iteration the chemistry library can take.

• “auxiliary data” [Array(string)] defines additional chemistry related data that the user can request be saved to vis files. Currently “pH” is the only variable supported.

<ParameterList>  <!-- parent list -->
<ParameterList name="_CHEMISTRY">
  <ParameterList name="thermodynamic database">
    <Parameter name="file" type="string" value="_TRITIUM.bgd"/>
    <Parameter name="format" type="string" value="simple"/>
  </ParameterList>
  <Parameter name="activity model" type="string" value="unit"/>
  <Parameter name="tolerance" type="double" value="1.5e-12"/>
  <Parameter name="maximum Newton iterations" type="int" value="25"/>
  <Parameter name="max timestep (s)" type="double" value="1.5e+07"/>
  <Parameter name="auxiliary data" type="Array(string)" value="{pH}"/>
  <Parameter name="number of component concentrations" type="int" value="1"/>
  <Parameter name="timestep control method" type="string" value="simple"/>
</ParameterList>
</ParameterList>

Energy PK

The conceptual PDE model for the energy equation is

\[\frac{\partial \varepsilon}{\partial t} = \boldsymbol{\nabla} \cdot (\kappa \nabla T) - \boldsymbol{\nabla} \cdot (\eta_l H_l \boldsymbol{q}_l) + Q\]

where \(\varepsilon\) is the energy density [\(J/m^3\)], \(\eta_l\) is molar density of liquid [\(mol/m^3\)], \(Q\) is heat source term, \(\boldsymbol{q}_l\) is the Darcy velocity [m/s], \(\kappa\) is thermal conductivity, and \(H_l\) is molar enthalpy of liquid [J/mol]. We define

\[\varepsilon = \phi (\eta_l s_l U_l + \eta_g s_g U_g) + (1 - \phi) \rho_r c_r T\]

where \(s_l\) is liquid saturation [-], \(s_g\) is gas saturation (water vapor), \(\eta_l\) is molar density of liquid [\(mol/m^3\)], \(\eta_g\) is molar density of gas, \(U_l\) is molar internal energy of liquid [J/mol], \(U_g\) is molar internal energy of gas (water vapor) [J/mol], \(\phi\) is porosity [-], \(\rho_r\) is rock density [\(kg/m^3\)], \(c_r\) is specific heat of rock [J/kg/K], and \(T\) is temperature [K].

Physical models and assumptions

This list is used to summarize physical models and assumptions, such as coupling with other PKs. This list is often generated on a fly by a high-level MPC PK.

energy_pk-spec

• “vapor diffusion” [bool] is set up automatically by a high-level PK, e.g. by EnergyFlow PK. The default value is “false”.

• “eos lookup table” [string] provides the name for optional EOS lookup table.

<ParameterList>  <!-- parent list -->
<ParameterList name="_ENERGY">
  <ParameterList name="physical models and assumptions">
    <Parameter name="vapor diffusion" type="bool" value="false"/>
    <Parameter name="eos lookup table" type="string" value="h2o.eos"/>
  </ParameterList>
</ParameterList>
</ParameterList>

Internal energy

Internal energy model is function of temperature only. Units are J/{mol/kg}. Internal energy list has a few parameters that allows us to run this PK in a variety of regimes, e.g. with or without gas phase.

iem-spec

• “energy key” [string] specifies name for the internal energy field. The default value is “energy”.

• “evaluator type” [string] changes the evaluator for internal energy. Available options are “generic” and “constant liquid density” (default).

• “vapor diffusion” [bool] specifies presence of a gas phase. The default value is “true”.

<ParameterList name="_ENERGY">  <!-- parent list -->
<ParameterList name="energy evaluator">
  <Parameter name="energy key" type="string" value="energy"/>
  <Parameter name="evaluator type" type="string" value="constant liquid density"/>
  <Parameter name="vapor diffusion" type="bool" value="true"/>
  <ParameterList name="verbose object">
    <Parameter name="verbosity level" type="string" value="high"/>
  </ParameterList>
</ParameterList>
</ParameterList>

Molar enthalpy

Field evaluator for specific molar enthalpy,

\[h = u + \displaystyle\frac{p}{\rho}.\]

<ParameterList name="_ENERGY">  <!-- parent list -->
<ParameterList name="enthalpy evaluator">
  <Parameter name="enthalpy key" type="string" value="enthalpy_liquid"/>
  <Parameter name="internal energy key" type="string" value="internal_energy_liquid"/>

  <Parameter name="include work term" type="bool" value="true"/>
  <Parameter name="pressure key" type="string" value="pressure"/>
  <Parameter name="molar density key" type="string" value="molar_density_liquid"/>
</ParameterList>
</ParameterList>

Thermal conductivity

Evaluator for thermal conductivity allows us to select a proper model. The variety of available models allows to run the energy PK by itself or in coupling with flow PK. The single-phase model accepts the following parameters:

tcm_one_phase-spec

• “TCM0” [list] defines a reginal model and its parameters.

• “thermal conductivity type” [string] is the name of conductivity model for a mixture of phases. Available two-phase models are “two-phase Peters-Lidard”, and “two-phase wet/dry”. Available one-phase model is “one-phase polynomial”.

• “regions” [Array(string)] specifies regions where model is applied.

• “solid phase” [list] define a thermal conductivity model for solid phase.

  • “eos type” [stirng] defines EOS model. Available options are “constant” and “salt”.

  • “reference conductivity” [double] defines constant or reference conductivity.

  • “reference temperature” [double] defines temperature at which reference conductivity of liquid is calculated. Default value is 298.15 [K].

  • “Sutherland constant” [double] defines parameter in temperature dependence.

• “gas phase” [list] define a thermal conductivity model for gas phase.

  • “eos type” [stirng] defines EOS model. Available options are “constant” and “ideal gas”.

  • “reference conductivity” [double] defines constant or reference conductivity.

  • “Sutherland constant” [double] defines parameter in temperature dependence

• “liquid phase” [list] define a thermal conductivity model for liquid phase.

  • “eos type” [stirng] defines EOS model. Available options are “constant” and “liquid water”.

  • “reference conductivity” [double] defines constant or reference conductivity. Default value is 0.6065 [W/m/K].

  • “reference temperature” [double] defines temperature at which reference conductivity of liquid is calculated. Default value is 298.15 [K].

  • “polynomial expansion” [Array(double)] collect coefficients in the quadratic representation of the thermal conductivity of liquid with respect to the dimensionless parameter T/Tref.

<ParameterList name="_ENERGY">  <!-- parent list -->
<ParameterList name="thermal conductivity evaluator">
  <ParameterList name="TCM0">
    <Parameter name="regions" type="Array(string)" value="{EntireDomain}"/>
    <Parameter name="thermal conductivity type" type="string" value="one-phase polynomial"/>
    <ParameterList name="solid phase">
      <Parameter name="eos type" type="string" value="constant"/>
      <Parameter name="reference conductivity" type="double" value="0.2"/>
      <Parameter name="reference temperature" type="double" value="298.15"/>
    </ParameterList>
    <ParameterList name="liquid phase">
      <Parameter name="eos type" type="string" value="polynomial"/>
      <Parameter name="reference conductivity" type="double" value="0.6"/>
      <Parameter name="polynomial expansion" type="Array(double)" value="{-1.48445, 4.12292, -1.63866}"/>
    </ParameterList>
    <ParameterList name="gas phase">
    </ParameterList>
  </ParameterList>
</ParameterList>
</ParameterList>

Operators

This section contains sublist for diffusion and advection operators. It also has one global parameter.

operators-spec

• “operators” [list]

  • “include enthalpy in preconditioner” [bool] allows us to study impact (usually positive) of including enthalpy term in the preconditioner. Default value is true.

  • “diagonal shift” [double] allows for a constant shift to be applied to the diagonal of the assembled operator, which can b useful for dealing with singular or near-singular matrices. Default is 0.0.

Diffusion operator

Operators sublist describes the PDE structure of the flow, specifies a discretization scheme, and selects assembling schemas for matrices and preconditioners.

diffusion_op-spec

• “diffusion operator” [list] defines parameters for generating and assembling diffusion matrix.

  • “matrix” [list] defines parameters for generating and assembling diffusion matrix. See section describing operators.

  • “preconditioner” [list] defines parameters for generating and assembling diffusion matrix that is used to create preconditioner. Since update of preconditioner can be lagged, we need two objects called “matrix” and “preconditioner”.

<ParameterList name="operators">
  <Parameter name="include enthalpy in preconditioner" type="boll" value="true"/>
  <ParameterList name="diffusion operator">
    <ParameterList name="matrix">
      <Parameter name="discretization primary" type="string" value="mfd: optimized for monotonicity"/>
      <Parameter name="discretization secondary" type="string" value="mfd: optimized for sparsity"/>
      <Parameter name="schema" type="Array(string)" value="{face, cell}"/>
      <Parameter name="preconditioner schema" type="Array(string)" value="{face}"/>
      <Parameter name="gravity" type="bool" value="false"/>
      <Parameter name="upwind method" type="string" value="standard: cell"/>
    </ParameterList>
    <ParameterList name="preconditioner">
      <Parameter name="discretization primary" type="string" value="mfd: optimized for monotonicity"/>
      <Parameter name="discretization secondary" type="string" value="mfd: optimized for sparsity"/>
      <Parameter name="schema" type="Array(string)" value="{face, cell}"/>
      <Parameter name="preconditioner schema" type="Array(string)" value="{face}"/>
      <Parameter name="gravity" type="bool" value="true"/>
      <Parameter name="Newton correction" type="string" value="approximate Jacobian"/>
      <Parameter name="upwind method" type="string" value="standard: cell"/>
    </ParameterList>
  </ParameterList>
</ParameterList>
</ParameterList>

This example uses cell-centered discretization for

Advection operator

This section to be written.

<ParameterList name="operators">  <!-- parent list -->
<ParameterList name="advection operator">
  <Parameter name="method" type="string" value="upwind"/>
</ParameterList>
</ParameterList>

Sources and sinks

The sources and sinks for injecting and removing energy from the system. Negative source removes energy. Positive source inject energy. The structure of list source terms mimics that of list boundary conditions. Again, constant functions can be replaced by any of the available functions.

sources-spec

• “regions” [Array(string)] is the list of regions where the source is defined.

• “spatial distribution method” [string] is the method for distributing source Q over the specified regions. The available options are “volume” and “none”. For option “none”, the source term function Q is measured in [J/m^3/s]. For option “volume”, it is measured in [J/s]. When the source function is defined over a few regions, Q is distributed over their union. Option “volume fraction” can be used when the region geometric model support volume fractions.

• “use volume fractions” instructs the code to use all available volume fractions. Note that the region geometric model supports volume fractions only for a few regions.

• “submodel” [string] refines definition of the source. Available options are “rate”, “integrated source”. The first option defines the source in a natural way as the rate of change q. The second option defines the indefinite integral Q of the rate of change, i.e. the source term is calculated as q = dQ/dt. Default is “rate”.

<ParameterList name="energy">  <!-- parent list -->
  <ParameterList name="source terms">
    <ParameterList name="_SRC 0">
      <Parameter name="regions" type="Array(string)" value="{_WELL_EAST}"/>
      <Parameter name="spatial distribution method" type="string" value="volume"/>
      <Parameter name="submodel" type="string" value="rate"/>
      <ParameterList name="source">
        <ParameterList name="function-constant">
          <Parameter name="value" type="double" value="-0.1"/>
        </ParameterList>
      </ParameterList>
    </ParameterList>
  </ParameterList>
</ParameterList>

Navier Stokes PK

The conceptual PDE model for the incompressible Navier Stokes equations are

\[\frac{\partial (\rho \boldsymbol{u})}{\partial t} + \boldsymbol{\nabla} \cdot (\rho \boldsymbol{u} \otimes \boldsymbol{u}) = - \boldsymbol{\nabla} p + \boldsymbol{\nabla} \cdot \boldsymbol{\sigma} + \rho \boldsymbol{g}\]

where \(\rho\) is the fluid density [kg/m^3], \(p\) is the pressure [Pa], \(\boldsymbol{\sigma}\) is the deviatoric stress tensor, \(\boldsymbol{g}\) is the gravity vector [\(m/s^2\)], and \(u \otimes v = u \times v^T\). The Stokes stress contitutive law for incompressible viscous fluid is

\[\boldsymbol{\sigma} = \mu \left(\boldsymbol{\nabla} \boldsymbol{u} + \boldsymbol{\nabla} \boldsymbol{u}^{T}\right),\]

where \(\mu\) is the dynamic viscosity [\(Pa \cdot s\)]. It can depend on density and pressure.

Physical models and assumptions

This list is used to summarize physical models and assumptions, such as coupling with other PKs. This list is often generated on a fly by a high-level MPC PK.

• “gravity” [bool] is set up automatically by a high-level PK. The default value is “false”.

<ParameterList name="_NAVIER_STOKES">  <!-- parent list -->
<ParameterList name="physical models and assumptions">
  <Parameter name="gravity" type="bool" value="false"/>
</ParameterList>
</ParameterList>

Operators

This section contains sublist for diffusion and advection operators. It also has one global parameter.

ns_operators-spec

• “operators” [list]

Elasticity operator

<ParameterList name="operators">  <!-- parent list -->
<ParameterList name="elasticity operator">
  <Parameter name="method" type="string" value="BernardiRaugel"/>
  <ParameterList name="schema">
    <Parameter name="base" type="string" value="cell"/>
    <Parameter name="location" type="Array(string)" value="{node, face}"/>
    <Parameter name="type" type="Array(string)" value="{vector, normal component}"/>
    <Parameter name="number" type="Array(int)" value="{2, 1}"/>
  </ParameterList>
</ParameterList>
</ParameterList>

Convection operator

This section to be written.

<ParameterList name="operators">  <!-- parent list -->
<ParameterList name="convection operator">
  <Parameter name="flux formula" type="string" value="Rusanov"/>
</ParameterList>
</ParameterList>

Divergence operator

This section to be written.

<ParameterList name="operators">  <!-- parent list -->
<ParameterList name="divergence operator">
  <Parameter name="method" type="string" value="BernardiRaugel"/>
  <ParameterList name="schema domain">
    <Parameter name="base" type="string" value="cell"/>
    <Parameter name="location" type="Array(string)" value="{node, face}"/>
    <Parameter name="type" type="Array(string)" value="{scalar, normal component}"/>
    <Parameter name="number" type="Array(int)" value="{2, 1}"/>
  </ParameterList>
  <ParameterList name="schema range">
    <Parameter name="base" type="string" value="cell"/>
    <Parameter name="location" type="Array(string)" value="{cell}"/>
    <Parameter name="type" type="Array(string)" value="{scalar}"/>
    <Parameter name="number" type="Array(int)" value="{1}"/>
  </ParameterList>
</ParameterList>
</ParameterList>

Multiphase PK

Mathematical models

The conceptual PDE model for the isothermal multiphase flow include transport equations for components and nonlinear algebraic constraints for the phase presence. At the moment we consider two phases (liquid and gas), multiple components, and one constraint. Each transport equation has the following form:

\[\frac{\partial \Theta}{\partial t} + \nabla \cdot \boldsymbol{\Psi} = Q,\]

where \(\Theta\) is the storage and \(\boldsymbol{\Psi}\) is the flux. The storage term sums up component amount across two phases, \(\alpha=l\) for liquid phase and \(\alpha=g\) for gas phase:

\[\Theta = \phi \sum_\alpha \eta_\alpha x_\alpha s_\alpha\]

where \(\phi\) is porosity [-], \(\eta\) is molar density [mol/m^3], \(x\) is molar fraction of component [-], and \(s\) is phase saturation [-].

The flux includes advective (wrt volumetric flux) and diffusive terms:

\[\boldsymbol{\Psi} = -\sum_\alpha \eta_\alpha \left(\boldsymbol{q}_\alpha + D_\alpha \nabla x \right)\]

where \(\boldsymbol{q}\) is Darcy phase velocity, \(D\) is molecular duffusion coefficient.

The nonlinear algebraic constraint may have different forms. One of the available forms is

\[min (s_g, 1 - \sum\limits_i x_g^i) = 0.\]

It implies that if gas compounent is present, then the sum of component mole fractions on the gas phase must be equal to 1.

The test examples illustrate three choices of primary variables for various models. The first one includes pressure liquid, mole gas fraction, and saturation liquid. The second one includes pressure liquid, molar gas density, and saturation liquid. The third one is based on the model in Jaffre’s paper. This model describes the two-phase two-component system with water and hydrogen.

The structure of a system of multiphase equations is described by sublist “system” which contains as many blocks as there are equations. The names of blocks are currently reserved and cannot be changed. Each block except for the “constraint eqn” has the following parameters:

multiphase_pk-spec

• “primary unknown” [string] defines a name of the primary variable which goes into a solution vector.

• “accumulation” [string] defines an evaluator for the accumulattion (or storage) term in this equation.

• “terms” [list] specifies details of underluing PDE operators. The list constains as many sublists as there are operators. Some of the operators can be considered either as diffusion operators with respect to particular fields or as advection operators associate with a Darcy velocities. We input parameters follow the diffusion viewpoint for all operators.

  • “coefficient” [string] defined the diffusion coefficient.

  • “argument” [string] defines a field for which the operator has the diffusion structure.

  • “scaling factor” [double] defines a scalar multiplier for the operator.

  • “phase” [int] specifies a phase accosiated with the operator. It is used to upwind the diffusion coefficient w.r.t. to the corresponding Darcy velocity.

<ParameterList name="system">
  <ParameterList name="pressure eqn">
    <Parameter name="primary unknown" type="string" value="pressure_liquid"/>
    <Parameter name="accumulation" type="string" value="total_water_storage"/>

    <ParameterList name="terms">
      <ParameterList name="advection">
        <Parameter name="coefficient" type="string" value="advection_water"/>
        <Parameter name="argument" type="string" value="pressure_liquid"/>
        <Parameter name="scaling factor" type="double" value="1.0"/>
        <Parameter name="phase" type="int" value="0"/>
      </ParameterList>
      <ParameterList name="diffusion">
        <Parameter name="coefficient" type="string" value="diffusion_liquid"/>
        <Parameter name="argument" type="string" value="molar_density_liquid"/>
        <Parameter name="scaling factor" type="double" value="-0.2"/>
        <Parameter name="phase" type="int" value="0"/>
      </ParameterList>
    </ParameterList>
  </ParameterList>

  <ParameterList name="constraint eqn">
    <Parameter name="primary unknown" type="string" value="saturation_liquid"/>
    <Parameter name="ncp evaluators" type="Array(string)" value="{ ncp_f, ncp_g }"/>
  </ParameterList>
</ParameterList>

The available boundary conditions include the prescibed pressure, total influx, and saturation.

Shallow water PK

Reconstruction and limiters

The control of the second-order numerical scheme is done via “reconstruction” sublist, described in Reconstruction. Here is the example:

<ParameterList name="shallow water">  <!-- parent list -->
<ParameterList name="reconstruction">
  <Parameter name="method" type="string" value="cell-based"/>
  <Parameter name="polynomial order" type="int" value="1"/>
  <Parameter name="weight" type="string" value="constant"/>
  <Parameter name="limiter" type="string" value="Barth-Jespersen"/>
  <Parameter name="limiter stencil" type="string" value="cell to closest cells"/>
  <Parameter name="limiter location" type="string" value="node"/>
  <Parameter name="limiter points" type="int" value="0"/>
  <Parameter name="limiter cfl" type="double" value="0.1"/>
</ParameterList>
</ParameterList>

Mechanics PK

The conceptual PDE model for the quasi-static elastic deformation is

\[\boldsymbol{\nabla} \cdot (C \varepsilon(\boldsymbol{d})) = \rho \boldsymbol{g}\]

where \(\boldsymbol{d}\) is the displacement [m], \(\rho\) is the rock density [kg/m^3], \(\boldsymbol{\varepsilon}\) is the strain tensor, and \(\boldsymbol{g}\) is the gravity vector [\(m/s^2\)].

For a linear elasticity problem, the stress tensor \(C\) is a linear operator acting on the strain tensor. In general, the stress tensor is a nonlinear operator.

Global parameters

Global parameters are placed in the sublist “mechanics”. The list of global parameters include:

mechanics_pk-spec

• “domain name” [string] specifies mesh name that defined domain of this PK. Default is “domain”.

Physical models and assumptions

This list is used to summarize physical models and assumptions, such as coupling with other PKs. This list is often generated or extended by a high-level MPC PK.

mechanics_assumptions-spec

• “use gravity” [bool] defines non-zero source term. Default is false.

• “biot scheme: undrained split” [bool] defines iterative coupling with a flow PK where mechanics is solved first.

• “biot scheme: fixed-stress split” [bool] defines iterative coupling with a flow PK where flow is solved first.

<ParameterList name="mechanics">  <!-- parent list -->
<ParameterList name="physical models and assumptions">
  <Parameter name="use gravity" type="bool" value="false"/>
  <Parameter name="biot scheme: undrained split" type="bool" value="false"/>
  <Parameter name="biot scheme: fixed stress split" type="bool" value="false"/>
</ParameterList>

Main sublists

The following sublists are needed to create and control this PK:

<ParameterList name="mechanics">  <!-- parent list -->
  <ParameterList name="time integrator">
  <ParameterList name="operators">
  <ParameterList name="physical models and assumptions">
  <ParameterList name="boundary conditions">

Nonlinear material models

Currently, the only available option is the Hardin-Drnevich hyperbolic model. Recall that

\[\boldsymbol{\sigma} = 2 G {\rm dev}(\varepsilon) + K {\rm trace}(\varepsilon) \boldsymbol{I}\]

where \(K\) is the bulk modulus and \(G\) is the shear modulus. The shear modulus is a nonlinear function of shear strain \(\gamma\):

\[G = G_{max} \left(1 + \displaystyle\frac{\gamma}{\gamma_{ref}}\right)^{-1}\]

where \(G_{max}\) is the maximum shear modulus and \(\gamma_{ref}\) is the reference shear strain. The model parameters are placed in the following sublist:

<ParameterList name="small strain models">
  <ParameterList name="All">
    <Parameter name="regions" type="Array(string)" value="{All}"/>
    <Parameter name="model" type="string" value="Hardin Drnevich"/>
    <Parameter name="reference shear strain" type="double" value="0.8"/>
    <Parameter name="maximum shear stress" type="double" value="70.0"/>
  </ParameterList>
</ParameterList>

Boundary conditions

Boundary conditions are defined in sublist boundary conditions. Two types of boundary conditions are supported. Each type has a similar structure: a list of identical elements that contain information about a part of the boundary where it is prescribed and a function to calculate it.

mechanics_boundary_conditions-spec

• “displacement” [list] is the Dirichlet boundary condition where the displacement is prescribed on a part of the boundary surface.

• “traction” [list] is the Neumann boundary condition where an outward traction is prescribed on a part of the boundary surface. This is the default boundary condition.

• “kinematic” [list] is the kinematic boundary condition, the essential condition for the normal velocity.

<ParameterList name="mechanics">  <!-- parent list -->
<ParameterList name="boundary conditions">
  <ParameterList name="displacement">
    <ParameterList name="BC_0">
      <Parameter name="regions" type="Array(string)" value="{_SURFACE_0}"/>
      <Parameter name="spatial distribution method" type="string" value="none"/>
      <ParameterList name="no slip">
        <Parameter name="number of dofs" type="int" value="2"/>
        <Parameter name="function type" type="string" value="composite function"/>
        <ParameterList name="dof 1 function">
          <ParameterList name="function-constant">
            <Parameter name="value" type="double" value="0.0"/>
          </ParameterList>
        </ParameterList>
        <ParameterList name="dof 2 function">
          <ParameterList name="function-constant">
            <Parameter name="value" type="double" value="0.0"/>
          </ParameterList>
        </ParameterList>
      </ParameterList>
    </ParameterList>
  </ParameterList>
  <ParameterList name="traction">
    <ParameterList name="BC_1">
      <Parameter name="regions" type="Array(string)" value="{_SURFACE_1}"/>
      <Parameter name="spatial distribution method" type="string" value="none"/>
      <ParameterList name="traction">
        <!-- similar to the displacement boundary condition -->
      </ParameterList>
    </ParameterList>
  </ParameterList>

COUPLED PKs (MPCs)

Coupling of process kernels requires additional parameters for PK described above.

mpcs-typed-spec

• “reactive transport” See MPC: Reactive Transport

• “flow reactive transport” See MPC: Flow and Reactive Transport

• “thermal flow” See MPC: Thermal Flow

• “darcy matrix fracture” See MPC: Flow in Fractured Rock

Reactive transport PK

Process kernel for coupling of Transport_PK and Chemistry_PK. Reactive transport can be setup using a steady-state flow. The two PKs are executed consequitively. The input spec requires new keyword reactive transport.

<ParameterList name="PK tree">  <!-- parent list -->
<ParameterList name="_REACTIVE TRANSPORT">
  <Parameter name="PK type" type="string" value="reactive transport"/>
  <ParameterList name="_TRANSPORT">
    <Parameter name="PK type" type="string" value="transport"/>
  </ParameterList>
  <ParameterList name="_CHEMISTRY">
    <Parameter name="PK type" type="string" value="chemistry amanzi"/>
  </ParameterList>
</ParameterList>
</ParameterList>

Flow and reactive transport PK

PK for coupling of Flow PK with Transport_PK and Chemistry_PK Amanzi uses operator splitting approach for coupled physical kernels. The coupling of PKs is described as a tree where flow and reactive transport are executed consequitively. The input spec requires new keyword flow reactive transport.

<ParameterList name="PK tree">  <!-- parent list -->
<ParameterList name="_FLOW and REACTIVE TRANSPORT">
  <Parameter name="PK type" type="string" value="flow reactive transport"/>
  <ParameterList name="_FLOW">
    <Parameter name="PK type" type="string" value="darcy"/>
  </ParameterList>
  <ParameterList name="_REACTIVE TRANSPORT">
    <Parameter name="PK type" type="string" value="reactive transport"/>
    <ParameterList name="_TRANSPORT">
    <Parameter name="PK type" type="string" value="transport"/>
    </ParameterList>
    <ParameterList name="_CHEMISTRY">
      <Parameter name="PK type" type="string" value="chemistry amanzi"/>
    </ParameterList>
  </ParameterList>
</ParameterList>
</ParameterList>

This example describe four PKs identified by keywords darcy, reactive transport, transport, and chemistry amanzi. The flow is fully saturated. The transport of reactive chemicals is based on the native chemistry package chemistry amanzi.

Details of PKs are organized as a plain list of ParameterLists. Note that reactive transport is MPC-PK and hence its description is short.

<ParameterList name="PKs">
<ParameterList name="_FLOW and REACTIVE TRANSPORT">
  <Parameter name="PK type" type="string" value="flow reactive transport"/>
  <Parameter name="PKs order" type="Array(string)" value="{_FLOW, _REACTIVE TRANSPORT}"/>
  <Parameter name="master PK index" type="int" value="0"/>
</ParameterList>

<ParameterList name="_REACTIVE TRANSPORT">
  <Parameter name="PK type" type="string" value="reactive transport"/>
  <Parameter name="PKs order" type="Array(string)" value="{_CHEMISTRY, _TRANSPORT}"/>
</ParameterList>

<ParameterList name="_FLOW">
  ...
</ParameterList>

<ParameterList name="_TRANSPORT">
  ...
</ParameterList>

<ParameterList name="_CHEMISTRY">
  ...
</ParameterList>
</ParameterList>

Thermal flow PK

Process kernel that strongly couples Flow PK with Energy PK. The conceptual PDE model of the coupled flow and energy equations is

\[\begin{split}\begin{array}{l} \frac{\partial \theta}{\partial t} = - \boldsymbol{\nabla} \cdot (\eta_l \boldsymbol{q}_l) - \boldsymbol{\nabla} \cdot (\phi s_g \tau_g D_g \boldsymbol{\nabla} X_g) + Q_1, \quad \boldsymbol{q}_l = -\frac{\boldsymbol{K} k_r}{\mu} (\boldsymbol{\nabla} p - \rho_l \boldsymbol{g}) \\ % \frac{\partial \varepsilon}{\partial t} = \boldsymbol{\nabla} \cdot (\kappa \nabla T) - \boldsymbol{\nabla} \cdot (\eta_l H_l \boldsymbol{q}_l) + Q_2 \end{array}\end{split}\]

In the first equation, \(\theta\) is total water storage (we use non-conventional definition) [\(mol/m^3\)], \(\eta_l\) is molar density of liquid [\(mol/m^3\)], \(\rho_l\) is fluid density [\(kg/m^3\)], \(Q_1\) is source or sink term, \(\boldsymbol{q}_l\) is the Darcy velocity [m/s], \(k_r\) is relative permeability [-], \(\boldsymbol{g}\) is gravity [\(m/s^2\)], \(\phi\) is porosity [-], \(s_g\) is gas saturation (water vapor) [-], \(\tau_g\) is tortuosity of gas [-], \(D_g\) is diffusion coefficient, and \(X_g\) is molar fraction of water in the gas phase [-]. We define

\[\theta = \phi (s_g \eta_g X_g + s_l \eta_l)\]

where \(s_l\) is liquid saturation [-], and \(\eta_g\) is molar density of gas.

In the second equation, \(\varepsilon\) is the energy density [\(J/mol^3\)], \(Q_2\) is source or sink term, \(\kappa\) is thermal conductivity [W/m/K], \(H_l\) is molar enthalphy of liquid [J/mol], and \(T\) is temperature [K]. We define

\[\varepsilon = \phi (\eta_l s_l U_l + \eta_g s_g U_g) + (1 - \phi) \rho_r c_r T\]

where \(U_l\) is molar internal energy of liquid [J/mol], \(U_g\) is molar internal energy of gas (water vapor) [J/mol], \(\rho_r\) is rock density [kg/m^3], and \(c_r\) is specific heat of rock [J/kg/K].

Diffusion operator

<ParameterList name="_NAVIER_STOKES">  <!-- parent lists -->
<ParameterList name="operator">
  <ParameterList name="diffusion operator">
   <ParameterList name="vapor matrix">
     <Parameter name="discretization primary" type="string" value="mfd: optimized for sparsity"/>
     <Parameter name="discretization secondary" type="string" value="mfd: optimized for sparsity"/>
     <Parameter name="schema" type="Array(string)" value="{face, cell}"/>
     <Parameter name="nonlinear coefficient" type="string" value="standard: cell"/>
     <Parameter name="exclude primary terms" type="bool" value="false"/>
     <Parameter name="scaled constraint equation" type="bool" value="false"/>
     <Parameter name="gravity" type="bool" value="false"/>
     <Parameter name="Newton correction" type="string" value="none"/>
   </ParameterList>
 </ParameterList>
 </ParameterList>

Coupled matrix-fracture Darcy flow PK

The process kernel that couples Darcy flow in matrix and fracture network.

Mathematical models

Let subscripts \(m\) and \(f\) correspond to matrix and fracture, respectively. The conceptual PDE model of the stationary coupled matrix-fracture flow is

\[\begin{split}\begin{array}{l} \phi_m \frac{S_{s,m}}{g} \frac{\partial p_m}{\partial t} + \boldsymbol{\nabla} \cdot (\rho_l \boldsymbol{q}_m) = Q_m, \quad \boldsymbol{q}_m = -\frac{\boldsymbol{K}_m}{\mu} (\boldsymbol{\nabla} p_m - \rho_l \boldsymbol{g}) \\ % \phi_f \frac{S_{s,f}}{g} \frac{\partial p_f}{\partial t} + \boldsymbol{\nabla} \cdot (\rho_l \boldsymbol{q}_f) = -\rho_l [[ \tilde{\boldsymbol{q}}_m \cdot \boldsymbol{n} ]], \quad \boldsymbol{q}_f = -\frac{\boldsymbol{K}_f}{\mu} (\boldsymbol{\nabla} p_f - \rho_l \boldsymbol{g}) \\ % \tilde{\boldsymbol{q}}_m \cdot \boldsymbol{n} = \frac{k}{g} (p_f - p_m) \end{array}\end{split}\]

subject to convential boundary conditions for both matrix and fracture domains expect for the matrix-fracture boundary where the boundary condition is

\[\boldsymbol{q}_m \cdot \boldsymbol{n} = \tilde{\boldsymbol{q}}_m \cdot \boldsymbol{n}\]

Here \(\rho_l\) is fluid density [kg/m^3], \(\phi\) is porosity [-], \(S_s\) is ispecific storage [m], \(p\) is aqueous pressure [Pa], \(\boldsymbol{K}\) is absolute permeability [m^2] for matrix domain and [m^3] for fracture domain, \(Q_m\) is source or sink term, \(\boldsymbol{q}\) is the Darcy velocity [m/s] for matrix domain and [m^2/s] for fracture domain, \(k\) is effective normal premeability [s^-1], and \(\boldsymbol{g}\) is gravity [\(m/s^2\)].

Main parameters and sublists

• “PKs order” [array(string)] defines user names for two flow PKs. The matrix PK should be defined first.

• “time integrator” [list] defines a generic time integrator used by the cycle driver.

<ParameterList name="PKs">  <!-- parent list -->
<ParameterList name="_COUPLED DARCY FLOW">
  <Parameter name="PKs order" type="Array(string)" value="{_FLOW MATRIX, _FLOW FRACTURE}"/>
  <Parameter name="master PK index" type="int" value="0"/>
  <ParameterList name="time integrator">
    ...
  </ParameterList>
</ParameterList>
</ParameterList>

Aperture models

The Barton-Bandis aperture-stress relationship is defined by initial aperture \(a_0\), overbuden pressure/normal stress \(p_{ov}\), fluid pressure \(p_f\), and two fitting parameters A and B:

\[a = a_0 - \frac{A (p_{ov} - p_f)}{1 + B (p_{ov} - p_f)}\]

• “undeformed aperture” [double] aperture at zero effective normal stress

• “overburden pressure” [double] overburden pressure/total normal stress

• “BartonBandis A” [double] fitting parameter, \([m / Pa]\)

• “BartonBandis B” [double] fitting parameter, \([Pa^{-1}]\)

<ParameterList name="fracture aperture models"> <!-- parent list -->
  <ParameterList name="FAM 0">
     <Parameter name="fracture aperture model" type="string" value="Barton Bandis"/>
     <Parameter name="regions" type="Array(string)" value="{RegionBottom}"/>
     <Parameter name="undeformed aperture" type="double" value="1e-5"/>
     <Parameter name="overburden pressure" type="double" value="1e+10"/>
     <Parameter name="BartonBandis A" type="double" value="2.22e-5"/>
     <Parameter name="BartonBandis B" type="double" value="3.47e-2"/>
  </ParameterList>

The exponential aperture-stress relationship is defined by initial aperture \(a_0\), overbuden pressure/normal stress \(p_{ov}\), fluid pressure \(p_f\), and fracture compressibility (fitting parameter) \(alpha\):

\[a = a_0 \exp{\alpha (p_{ov} - p_f)}\]

• “undeformed aperture” [double] aperture at zero effective normal stress

• “overburden pressure” [double] overburden pressure/total normal stress

• “compressibility” [double] fracture compressibility, \([Pa^{-1}]\)

<ParameterList name="fracture aperture models"> <!-- parent list -->
  <ParameterList name="FAM 0">
     <Parameter name="fracture aperture model" type="string" value="exponential law"/>
     <Parameter name="regions" type="Array(string)" value="{RegionBottom}"/>
     <Parameter name="undeformed aperture" type="double" value="1e-5"/>
     <Parameter name="overburden pressure" type="double" value="1e+10"/>
     <Parameter name="compressibility" type="double" value="1.0e-11"/>
  </ParameterList>

Thermoporoelasticty PK

Sequential coupling of fully flow and mechanics PKs via the fixed stress split algorithm. The conceptual PDE model of the coupled flow and mechanics equations is

\[\begin{split}\begin{array}{l} -\boldsymbol{\nabla} \cdot (C \varepsilon(\boldsymbol{d})) + \alpha \nabla p + \beta_r \nabla T = \rho \boldsymbol{g} \\[1ex] % \displaystyle\frac{\partial (\eta \phi)}{\partial t} = - \boldsymbol{\nabla} \cdot (\eta \boldsymbol{q}) + Q_1\\[1ex] % \displaystyle\frac{\partial (\phi\, \eta\, U_l + (1 - \phi) \rho_r U_r)}{\partial t} = \boldsymbol{\nabla} \cdot (\kappa \nabla T) - \boldsymbol{\nabla} \cdot (\eta\, H\, \boldsymbol{q}) + Q_2. \end{array}\end{split}\]

In the first equation, \(C\) is elasticity tensor [Pa], \(\varepsilon\) is stain tensor [-], \(\boldsymbol{d}\) is displacement [m], \(\rho\) is density of liquid [\(kg/m^3\)], \(\alpha\) is Biot coefficient [-], \(\beta_r\) is thermal stress coefficient [\(Pa/K\)], \(\boldsymbol{g}\) is gravity [\(m/s^2\)], \(p\) is pressure [Pa], and \(T\) is temperature [K].

In the second equation, \(\eta\) is molar density of liquid [\(mol/m^3\)], \(\phi\) is porosity [-], \(Q_1\) is source or sink term, and \(\boldsymbol{q}\) is Darcy velocity [m/s]. The porosity model is

\[\phi = \phi_0 + \alpha\, {\rm div}\, \varepsilon + c_0\,(p - p_0) - \beta_l (T - T_0).\]

where \(\beta_l\) is thermal dilation coefficient [\(Pa/K\)], \(c_0\) is pore compressibility [\(Pa^{-1}\)], \(p_0\) is reference pressure [\(Pa\)], and \(T_0\) is reference temperature [\(K\)]. Currently \(T_0 = 0^\circ C\).

In the third equation, \(\kappa\) is thermal conductivity, \(U_l\) is internal energy of liquid [J/mol], \(U_r\) is internal energy of rock [J/kg], \(H\) is molar enthalpy of liquid [J/mol], \(\rho_r\) is rock density [\(kg/m^3\)], and \(Q_2\) is heat source term. The internal energies are based on linear models:

\[U_i = c_i\,(T - T_0),\quad i=l,r,\]

where \(c_i\) are specific heat capacities.

GENERIC CAPABILITIES

Collection of generic tools used by PKs.

Operators

Operators are discrete forms of linearized PDEs operators. They form a layer between physical process kernels and solvers and include accumulation, diffusion, advection, elasticity, reaction, and source operators. The residual associated with an operator \(L_h\) helps to understand the employed sign convention:

\[r = f - L_h u.\]

Operator represents a map from linear space X to linear space Y. Typically, this map is a linear map, and encapsulates much of the discretization involved in moving from continuous to discrete equations. The spaces X and Y are described by CompositeVectors (CV). A few maps X->Y are supported.

An operator provides an interface for applying both the forward and inverse linear map (assuming the map is invertible).

Typically the Operator is never seen by the user; instead the user provides input information for helper classes based on the continuous mathematical operator and the desired discretization. These helpers build the needed Operator, which may include information from multiple helpers (i.e. in the case of Jacobian Operators for a PDE).

However, one option may be provided by the user, which is related to dealing with nearly singular operators:

operators-spec

• “diagonal shift” [double] 0.0 Adds a scalar shift to the diagonal of the Operator, which can be useful if the Operator is singular or near-singular.

A PK decides how to bundle operators in a collection of operators. For example, an advection-diffusion problem may benefit from using a single operator that combines two operators representing diffusion and advection process. Collection of operators must be used for implicit solvers and for building preconditioners. In such a case, the collections acts as a single operator.

Operators use a few tools that are generic in nature and can be used independently by PKs. The list includes reconstruction and limiting algorithms.

Schema

The operators use notion of schema to describe operator’s abstract structure. Old operators use a simple schema which is simply the list of geometric objects where scalar degrees of freedom are defined. New operators use a list to define location, type, and number of degrees of freedom. In addition, the base of local stencil is either face or cell. A rectangular operator needs two schemas do describe its domain (called “schema domain”) and its range (called “schema range”). A square operator may use either two identical schema lists or a single list called “schema”.

<ParameterList name="pks operator name">  <!-- parent list-->
<ParameterList name="schema domain">
  <Parameter name="base" type="string" value="cell"/>
  <Parameter name="location" type="Array(string)" value="{node, face}"/>
  <Parameter name="type" type="Array(string)" value="{scalar, normal component}"/>
  <Parameter name="number" type="Array(int)" value="{2, 1}"/>
</ParameterList>
<ParameterList name="schema domain">
  <Parameter name="base" type="string" value="cell"/>
  <Parameter name="location" type="Array(string)" value="{node, face}"/>
  <Parameter name="type" type="Array(string)" value="{scalar, normal component}"/>
  <Parameter name="number" type="Array(int)" value="{2, 1}"/>
</ParameterList>
</ParameterList>

This example describes a square operator with two degrees of freedom per mesh node and one degree of freedom per mesh face. The face-based degree of freedom is the normal component of a vector field. Such set of degrees of freedom is used in the Bernardi-Raugel element for discretizing Stokes equations. Parameter “base” indicates that local matrices are associated with mesh cells.

Diffusion operator

Diffusion is the most frequently used operator. It employs the old schema.

diffusion_op-spec

• “pks operator name” [list] a PK specific name for the diffusion operator.

  • “discretization primary” [string] specifies an advanced discretization method that has useful properties under some a priori conditions on the mesh and/or permeability tensor. The available options are “mfd: optimized for sparsity”, “mfd: optimized for monotonicity”, “mfd: default”, “mfd: support operator”, “mfd: two-point flux approximation”, “fv: default”, and “nlfv: default”. The first option is recommended for general meshes. The second option is recommended for orthogonal meshes and diagonal absolute permeability tensor.

  • “discretization secondary” [string] specifies the most robust discretization method that is used when the primary selection fails to satisfy all a priori conditions. Default value is equal to that for the primary discretization.

  • “diffusion tensor” [string] specifies additional properties of the diffusion tensor. It allows us to solve problems with non-symmetric but positive definite tensors. Available options are symmetric (default) and nonsymmetric.

  • “nonlinear coefficient” [string] specifies a method for treating nonlinear diffusion coefficient, if any. Available options are “none”, “upwind: face”, “divk: cell-face” (default), “divk: face”, “standard: cell”, and “divk: cell-face-twin”. Symmetry preserving methods are the divk-family of methods and the classical cell-centered method (“standard: cell”). The first part of the name indicates the base scheme. The second part (after the semi-column) indicates required components of the composite vector that must be provided by a physical PK. Default is “none”.

  • “schema” [Array(string)] defines the operator stencil. It is a collection of geometric objects. It equals to “{cell}” for finite volume schemes. It is typically “{face, cell}” for mimetic discretizations.

  • “preconditioner schema” [Array(string)] defines the preconditioner stencil. It is needed only when the default assembling procedure is not desirable. If skipped, the “schema” is used instead.

  • “gravity” [bool] specifies if flow is driven also by the gravity.

  • “gravity term discretization” [string] selects a model for discretizing the gravity term. Available options are “hydraulic head” [default] and “finite volume”. The first option starts with equation for the shifted solution, i.e. the hydraulic head, and derives gravity discretization by the reserve shifting. The second option is based on the divergence formula.

  • “gravity magnitude” [double] defined magnitude of the gravity vector.

  • “Newton correction” [string] specifies a model for correction (non-physical) terms that must be added to the preconditioner. These terms approximate some Jacobian terms. Available options are “true Jacobian” and “approximate Jacobian”. The FV scheme accepts only the first options. The othre schemes accept only the second option.

  • “scaled constraint equation” [bool] rescales flux continuity equations on mesh faces. These equations are divided by the nonlinear coefficient. This option allows us to treat the case of zero nonlinear coefficient. At moment this feature does not work with non-zero gravity term. Default is false.

  • “constraint equation scaling cutoff”” [double] specifies the cutoff value for applying rescaling strategy described above.

  • “consistent faces” [list] may contain a “preconditioner” and “linear operator” list (see sections Preconditioners and LinearSolvers respectively). If these lists are provided, and the “discretization primary” is of type “mfd: *”, then the diffusion method UpdateConsistentFaces() can be used. This method, given a set of cell values, determines the faces constraints that satisfy the constraint equation in MFD by assembling and inverting the face-only system. This is not currently used by any Amanzi PKs.

  • “fracture” [Array(string)] provides list of regions that defines a fracture network. This parameter is used only by the coupled flow PK.

Example:

<ParameterList name="pks operator name">
  <Parameter name="discretization primary" type="string" value="mfd: optimized for monotonicity"/>
  <Parameter name="discretization secondary" type="string" value="mfd: two-point flux approximation"/>
  <Parameter name="schema" type="Array(string)" value="{face, cell}"/>
  <Parameter name="preconditioner schema" type="Array(string)" value="{face}"/>
  <Parameter name="gravity" type="bool" value="true"/>
  <Parameter name="gravity term discretization" type="string" value="hydraulic head"/>
  <Parameter name="gravity magnitude" type="double" value="9.81"/>
  <Parameter name="nonlinear coefficient" type="string" value="upwind: face"/>
  <Parameter name="Newton correction" type="string" value="true Jacobian"/>

  <ParameterList name="consistent faces">
    <ParameterList name="linear solver">
      ...
    </ParameterList>
    <ParameterList name="preconditioner">
      ...
    </ParameterList>
  </ParameterList>
</ParameterList>

This example creates a p-lambda system, i.e. the pressure is discretized in mesh cells and on mesh faces. The preconditioner is defined on faces only, i.e. cell-based unknowns are eliminated explicitly and the preconditioner is applied to the Schur complement.

Additional options available only for the MFD family of discretizations include:

diffusion_op_addon-spec

• “nonlinear coefficient” [string] specifies a method for treating nonlinear diffusion coefficient, if any. Available options are “none”, “upwind: face”, “divk: cell-face” (default), “divk: face”, “standard: cell”, “divk: cell-face-twin” and “divk: cell-grad-face-twin”. Symmetry preserving methods are the divk-family of methods and the classical cell-centered method (“standard: cell”). The first part of the name indicates the base scheme. The second part (after the semi-column) indicates required components of the composite vector that must be provided by a physical PK.

• “discretization secondary” [string] specifies the most robust discretization method that is used when the primary selection fails to satisfy all a priori conditions. This is typically “mfd: default”, and is used only when an MFD “discretization primary” is used.

• “schema” [Array(string)] defines the operator stencil. It is a collection of geometric objects. Typically this is set by the implementation and is not provided.

• “preconditioner schema” [Array(string)] {face,cell} Defines the preconditioner stencil. It is needed only when the default assembling procedure is not desirable. If skipped, the “schema” is used instead. In addition to the default, {face} may be used, which forms the Schur complement.

• “consistent faces” [list] may contain a “preconditioner” and “linear operator” list (see sections Preconditioners and LinearSolvers respectively). If these lists are provided, and the “discretization primary” is of type “mfd: *”, then the diffusion method UpdateConsistentFaces() can be used. This method, given a set of cell values, determines the faces constraints that satisfy the constraint equation in MFD by assembling and inverting the face-only system. This is not currently used by any Amanzi PKs.

• “diffusion tensor” [string] allows us to solve problems with symmetric and non-symmetric (but positive definite) tensors. Available options are symmetric (default) and nonsymmetric.

• “use manifold flux” [bool] false Computes the flux using algorithms and data structures for manifolds or fracture networks.

Advection operator

A high-order advection operator may have different domain and range and therefore requires two schemas. The structure of the new schema is described in the previous section. A high-order advection operator has two terms in a weak formulation, corresponding to volume and surface integrals. These two terms are discretixed using two operators with matrix of types advection and flux, respectively.

advection_op-spec

• “pks operator name” [list] a PK specific name for the advection operator.

  • “method” [string] defines a discretization method. The available option is “dg modal”.

  • “method order” [int] defines method order. For example, the classical low-order finite volume scheme is equivalent to DG of order 0.

  • “matrix type” [string] defines matrix type. The supported options are “advection” and “flux”.

  • “dg basis” [string] defines bases for DG schemes. The available options are “regularized” (recommended), “normalized”, “orthonormalized”, and “natural” (not recommended).

  • “gradient operator on test function” [bool] defines place of the gradient operator. For integration by parts schemes, the gradient is transfered to a test function. This option is needed for discretizing volumetric integrals.

  • “jump operator on test function” [bool] defines place of the jump operator. For integration by parts schemes, the jump operator is applied to a test function. This option is needed for discretizing surface fluxes.

  • “flux formula” [string] defines type of the flux. The available options are “Rusanov” (default), “upwind”, “downwind”, and “NavierStokes”.

  • “schema domain” [list] defines a discretization schema for the operator domain.

  • “schema range” [list] defines a discretization schema for the operator range.

<ParameterList name="pks operator name">
  <Parameter name="method" type="string" value="dg modal"/>
  <Parameter name="method order" type="int" value="2"/>
  <Parameter name="flux formula" type="string" value="Rusanov"/>
  <Parameter name="matrix type" type="string" value="flux"/>
  <Parameter name="jump operator on test function" type="bool" value="true"/>

  <ParameterList name="schema domain">
    <Parameter name="base" type="string" value="cell"/>
    <Parameter name="location" type="Array(string)" value="{node, face}"/>
    <Parameter name="type" type="Array(string)" value="{scalar, normal component}"/>
    <Parameter name="number" type="Array(int)" value="{2, 1}"/>
  </ParameterList>
  <ParameterList name="schema range">
    <Parameter name="base" type="string" value="cell"/>
    <Parameter name="location" type="Array(string)" value="{cell}"/>
    <Parameter name="type" type="Array(string)" value="{scalar}"/>
    <Parameter name="number" type="Array(int)" value="{1}"/>
  </ParameterList>
</ParameterList>

In this example, we construct an operator for volumetric integrals in a weak formulation of advection problem.

The only low-order advection operator in Amanzi is the upwind operator. It employes the old schema.

<ParameterList name="pks operator name">
  <Parameter name="base" type="string" value="face"/>
  <Parameter name="schema" type="Array(string)" value="{cell}"/>
  <Parameter name="method order" type="int" value="0"/>
  <Parameter name="matrix type" type="string" value="advection"/>
</ParameterList>

Reaction operator

A reaction operator may represent either reaction of identity operator. It is symmetric so far and requires one schema. The structure of the schema is described in the previous section.

reaction_op-spec

• “pks operator name” [list] a PK specific name for the advection operator.

  • “method” [string] defines a discretization method. The only supported option is “dg nodal”.

  • “schema” [list] defines a discretization schema for the operator domain.

<ParameterList name="pks operator name">
  <Parameter name="method" type="string" value="dg modal"/>
  <Parameter name="method order" type="int" value="1"/>
  <Parameter name="matrix type" type="string" value="mass"/>
  <ParameterList name="schema">
    <Parameter name="base" type="string" value="cell"/>
    <Parameter name="location" type="Array(string)" value="{cell}"/>
    <Parameter name="type" type="Array(string)" value="{scalar}"/>
    <Parameter name="number" type="Array(int)" value="{3}"/>
  </ParameterList>
</ParameterList>

Elasticity operator

Elasticity operator is used for describing soil deformation or fluid flow (Stokes and Navier-Stokes).

elasticity_op-spec

• “method” [string] defines a discretization method. The available options are “BernardiRaugel”.

• “schema” [list] defines a discretization schema.

  • “location” [Array(string)] defines geometric location of degrees of freedom.

  • “type” [Array(string)] defines type of degrees of freedom. The available options are “scalar” and “normal component”.

  • “number” [Array(int)] indicates how many time this degree of freedom is repeated.

<ParameterList name="elasticity operator">
  <Parameter name="method" type="string" value="BernardiRaugel"/>
  <ParameterList name="schema">
    <Parameter name="base" type="string" value="cell"/>
    <Parameter name="location" type="Array(string)" value="{node, face}"/>
    <Parameter name="type" type="Array(string)" value="{scalar, normal component}"/>
    <Parameter name="number" type="Array(int)" value="{2, 1}"/>
  </ParameterList>
</ParameterList>

Abstract operator

An abstract operator is designed for testing new discretization methods. It uses the factory of discretization methods and a few control parameters required by this factory and/or particular method in it.

abstract_op-spec

• “method” [string] defines a discretization method. The available options are “diffusion”, “diffusion generalized”, “BernardiRaugel”, “CrouzeixRaviart”, “CrouzeixRaviart serendipity”, “Lagrange”, “Lagrange serendipity”, and “dg modal”.

• “method order” [int] defines disretization order. It is used by high-order discretization methods such as the discontinuous Galerkin.

• “matrix type” [string] defines type of local matrix. Available options are “mass”, “mass inverse”, “stiffness”, “divergence”, and “advection”.

<ParameterList name="_ABSTRACT OPERATOR">
  <Parameter name="method" type="string" value="dg modal"/>
  <Parameter name="method order" type="int" value="2"/>
  <Parameter name="dg basis" type="string" value="regularized"/>
  <Parameter name="matrix type" type="string" value="flux"/>

  <ParameterList name="schema domain">
    ...
  </ParameterList>
  <ParameterList name="schema range">
    ...
  </ParameterList>
</ParameterList>

Diffusion is the most frequently used operator. Examples of usage this operator are in test/operators_stokes.cc and test/operators_diffusion_curved.cc In the first example, we set up a discrete divergence operator that corersponds to a rectangular matrix. In the second example, we set up an elliptic operator when Hermite-type degrees of freedom are used on curved faces.

Reconstruction and limiters

A reconstruction of discrete fields is used to increase accuracy of discrete models. The reconstruction can be either unconstrained or limited. Amanzi supports a variety of state-of-the-art reconstruction and limiting algorithms and their extensions for various PKs.

reconstruction-spec

• “reconstruction” [list] describes parameters used by reconstruction algorithms.

• “method” [string] specifies a reconstruction method. Available option is “cell-based” (default).

• “polynomial order” [int] defines the polynomial order of the reconstructed function. Default is 1.

• “weight” [string] defines weight for reconstruction. Available options are “constant” (default) and “inverse distance”.

• “limiter” [string] specifies limiting method. Available options are “Barth-Jespersen” (default), “Michalak-Gooch”, “tensorial”, and “Kuzmin”.

• “limiter stencil” [string] specifies stencil for calculating local bounds. Available options are “face to cells”, “cell to closets cells”, “cell to all cells”, and “node to cells”. For a square mesh, the above options define stencils of size 2, 5, 9, and 4, respectively. Option “face to cells” is default for “Barth-Jespersen”, “Michalak-Gooch”, and “tensorial”. Option “node to cells” is default for “Kuzmin”.

• “limiter points” [int] specifies the number of integration points (Gauss points in 2D) on face where limiting occurs. Default is 1. Limited to 2D.

• “limiter location” [string] defines geometry entity where the limiter points are located. Available options are “node”, “face”, and “cell”. Option “node” is default for “node to cells” stencil. Option “face” is default for other stencils.

• “limiter cfl” [double] is a safety factor (less than 1) applied to the limiter. Default value is 1.

• “use external bounds” [bool] specifies if bounds for limiters are provided by the hosting application. Default is “false”.`

• “limiter extension for transport” [bool] adds additional corrections to limiters required by the transport PK. Default value is false.

<ParameterList name="reconstruction">
  <Parameter name="method" type="string" value="cell-based"/>
  <Parameter name="polynomial order" type="int" value="1"/>
  <Parameter name="weight" type="string" value="inverse distance"/>
  <Parameter name="limiter" type="string" value="tensorial"/>
  <Parameter name="limiter extension for transport" type="bool" value="false"/>
  <Parameter name="limiter stencil" type="string" value="face to cells"/>
  <Parameter name="limiter points" type="int" value="0"/>
</ParameterList>

Time integrator

Backward Euler is the simplest of the implicit methods. It solves time integration schemes by evaluating all time derivatives at the new time. This makes it unconditionally stable, though potentially not very accurate. This unconditional stability tends to make it the workhorse of the types of stiff, nonlinear parabolic equations such as Richards equation and the diffusion wave approximation.

In this method, we look to solve:

\[\frac{\partial \mathbf{u}}{\partial t} = f(\mathbf{u},\mathbf{x},t)\]

via the time discretization scheme:

\[\frac{\mathbf{u}^{t + \Delta t} - \mathbf{u}^{t}}{\Delta t} = f(\mathbf{u}^{t + \Delta t}, \mathbf{x}, t + \Delta t)\]

bdf1-ti-spec

• “verbose object” [verbose-object-spec] A Verbose Object

• “residual debugger” [residual-debugger-spec] A Residual Debugger object.

• “max preconditioner lag iterations” [int] 0 specifies frequency of preconditioner recalculation.

• “freeze preconditioner” [bool] false enforces preconditioner to be updated only once per non-linear solver. When set to true, the above parameter is ignored.

• “extrapolate initial guess” [bool] true identifies forward time extrapolation of the initial guess.

• “nonlinear iteration initial guess extrapolation order” [int] 1 defines extrapolation algorithm. Zero value implies no extrapolation.

• “restart tolerance relaxation factor” [double] 1 Changes the nonlinear tolerance on restart. The time integrator is usually restarted when a boundary condition changes drastically. It may be beneficial to loosen the nonlinear tolerance on the first several timesteps after the time integrator restart. The default value is 1, while a reasonable value may be as large as 1000.

• “restart tolerance relaxation factor damping” [double] 1 Controls how fast the loosened nonlinear tolerance will revert back to the one specified in “nonlinear tolerance”. If the nonlinear tolerance is “tol”, the relaxation factor is “factor”, and the damping is “d”, and the timestep count is “n” then the actual nonlinear tolerance is “tol * max(1.0, factor * d ** n)”. Reasonable values are between 0 and 1.

INCLUDES - [solver-typed-spec] Uses a Solver. - [timestep-controller-typed-spec] Uses a Timestep Controller

Note this also accepts an object that provides the BDF1 Solver Interface.

<ParameterList name="time integrator">
  <Parameter name="time integration method" type="string" value="BDF1"/>
  <ParameterList name="BDF1">
    <Parameter name="max preconditioner lag iterations" type="int" value="5"/>
    <Parameter name="freeze preconditioner" type="bool" value="false"/>
    <Parameter name="extrapolate initial guess" type="bool" value="true"/>
    <Parameter name="nonlinear iteration initial guess extrapolation order" type="int" value="1"/>
    <Parameter name="restart tolerance relaxation factor" type="double" value="1.0"/>
    <Parameter name="restart tolerance relaxation factor damping" type="double" value="1.0"/>

    <Parameter name="timestep controller type" type="string" value="standard"/>
    <ParameterList name="timestep controller standard parameters">
      ...
    </ParameterList>

    <Parameter name="solver type" type="string" value="nka"/>
    <ParameterList name="nka parameters">
      ...
    </ParameterList>
  </ParameterList>
</ParameterList>

Time step controller

A timestep controller is used to abstract across PKs algorithms for computing a physical timestep size that is valid for that PK. Several controller types are available which implement common choices.

The timestep controller fixed provides a single, uniform dt.

The timestep controller standard is a simple timestep control mechanism which sets the next timestep based upon the previous timestep and how many nonlinear iterations the previous timestep took to converge. The next timestep is given by the following rule:

• if \(N_k > N^{max}\) then \(\Delta t_{k+1} = f_{reduction} \Delta t_{k}\)

• if \(N_k < N^{min}\) then \(\Delta t_{k+1} = f_{increase} \Delta t_{k}\)

• otherwise \(\Delta t_{k+1} = \Delta t_{k}\)

where \(\Delta t_{k}\) is the previous timestep and \(N_k\) is the number of nonlinear iterations required to solve step \(k\).

The timestep controller smart is based on standard, but also tries to be a bit smarter to avoid repeated increase/decrease loops where the step size decreases, converges in few iterations, increases, but then fails again. It also tries to grow the timestep geometrically to more quickly recover from tricky nonlinearities.

The timestep controller from file loads a timestep history from a file, then advances the step size with those values. This is mostly used for testing purposes, where we need to force the same timestep history as previous runs to do regression testing. Otherwise roundoff errors can eventually alter number of iterations enough to alter the timestep history, resulting in solutions which are enough different to cause doubt over their correctness.

Note, the timestep controller is also responsible for throwing an error if the timestep becomes invalid.

Standard controller

This is a simple timestep control mechanism which sets the next timestep based upon the previous timestep and how many nonlinear iterations the previous timestep took to converge.

The timestep for step \(k+1\), \(\Delta t_{k+1}\), is given by:

• if \(N_k > N^{max}\) then \(\Delta t_{k+1} = f_{reduction} * \Delta t_{k}\)

• if \(N_k < N^{min}\) then \(\Delta t_{k+1} = f_{increase} * \Delta t_{k}\)

• otherwise \(\Delta t_{k+1} = \Delta t_{k}\)

where \(\Delta t_{k}\) is the previous timestep and \(N_k\) is the number of nonlinear iterations required to solve step \(k\):.

timestep-controller-standard-spec

• “max iterations” [int] \(N^{max}\), decrease the timestep if the previous step took more than this.

• “min iterations” [int] \(N^{min}\), increase the timestep if the previous step took less than this.

• “timestep reduction factor” [double] \(f_{reduction}\), reduce the previous timestep by this multiple.

• “timestep increase factor” [double] \(f_{increase}\), increase the previous timestep by this multiple.

INCLUDES:

• [timestep-controller-recoverable-spec]

<ParameterList name="BDF1"> <!-- parent list -->
  <Parameter name="timestep controller type" type="string" value="standard"/>
  <ParameterList name="timestep controller standard parameters">
    <Parameter name="min iterations" type="int" value="10"/>
    <Parameter name="max iterations" type="int" value="15"/>
    <Parameter name="timestep increase factor" type="double" value="1.2"/>
    <Parameter name="timestep reduction factor" type="double" value="0.5"/>
    <Parameter name="max timestep [s]" type="double" value="1e+9"/>
    <Parameter name="min timestep [s]" type="double" value="0.0"/>
    <Parameter name="initial timestep [s]" type="double" value="86400"/>
  </ParameterList>
</ParameterList>

In this example, the timestep is increased by factor 1.2 when the nonlinear solver converges in 10 or less iterations. The timestep is not changed when the number of nonlinear iterations is between 11 and 15. The timestep will be cut twice if the number of nonlinear iterations exceeds 15.

Smarter controller

This is based on Timestep Controller Standard, but also tries to be a bit smarter to avoid repeated increase/decrease loops where the step size decreases, converges in few iterations, increases, but then fails again. It also tries to grow the step geometrically to more quickly recover from tricky nonlinearities.

timestep-controller-smarter-spec

• “max iterations” [int] \(N^{max}\), decrease the timestep if the

  previous step took more than this.

• “min iterations” [int] \(N^{min}\), increase the timestep if the

  previous step took less than this.

• “timestep reduction factor” [double] \(f_{reduction}\), reduce the previous timestep by this multiple.

• “timestep increase factor” [double] \(f_{increase}\), increase the previous timestep by this multiple. Note that this can be modified geometrically in the case of repeated successful steps.

• “max timestep increase factor” [double] 10. The max \(f_{increase}\) will ever get.

• “growth wait after fail” [int] Wait at least this many timesteps before attempting to grow the timestep after a failed timestep.

• “count before increasing increase factor” [int] Require this many successive increasions before multiplying \(f_{increase}\) by itself.

Adaptive controller

This is under development and is based on a posteriori error estimates.

Fixed controller

A fixed timestep controller simply sets a constant timestep size.

“timestep controller type” = “fixed”

timestep-controller-fixed-spec

• “initial timestep [s]” [double] The fixed timestep size.

File-based controller

This loads a timestep history from a file, then advances the step size with those values. This is mostly used for testing purposes, where we need to force the same timestep history as previous runs to do regression testing. Otherwise even machine roundoff can eventually alter number of iterations enough to alter the timestep history, resulting in solutions which are enough different to cause doubt over their correctness.

timestep-controller-from-file-spec

• “file name” [string] Path to hdf5 file containing timestep information.

• “timestep header” [string] Name of the dataset containing the history of timestep sizes.

Functions

To set up non-trivial boundary conditions and/or initial fields, Amanzi supports a few mathematical functions. New function types can added easily. Each function is defined by a list:

<ParameterList name="function name">
  function-specification
</ParameterList>

The parameter list name string NAME is arbitrary and meaningful only to the parent parameter list. This list is given as input to the Amanzi::FunctionFactory::Create method which instantiates a new Amanzi::Function object. The function-specification is one of the following parameter lists.

Analytic, algabraic functions of space and time are used for a variety of purposes, including boundary conditions, initial conditions, and independent variables.

For initial conditions, functions are prescribed of space only, i.e.

\[u = f(x,y,z)\]

For boundary conditions and independent variables, functions are also a function of time:

\[u = f(t,x,y,z)\]

Additive function

An additive function simply adds two other function results together.

\[f(x) = f_1(x) + f_2(x)\]

where \(f_1\) is defined by the “function1” sublist, and \(f_2\) by the “function2” sublist.

function-additive-spec

• “function1” [function-typedinline-spec] \(f_1\) in \(f(x) = f_1(x) + f_2(x)\)

• “function2” [function-typedinline-spec] \(f_2\) in \(f(x) = f_1(x) + f_2(x)\)

Example:

<ParameterList name="function-additive">
  <ParameterList name="function1">
    function-specification
  </ParameterList>
  <ParameterList name="function2">
    function-specification
  </ParameterList>
</ParameterList>

Bilinear function

A piecewise bilinear function extends the linear form of the tabular function to two variables.

Define \(i(x) = i : x_i < x <= x_{{i+1}}\) and similarly \(j(y) = j : y_j < y <= y_{{j+1}}\) for monotonically increasing \(x_i\) and \(y_j\).

Given a two-dimensional array \(u_{i,j}\), \(f\) is then defined by bilinear interpolation on \(u_{i(x),j(y)}, u_{i(x)+1,j(y)}, u_{i(x),j(y)+1}, u_{i(x)+1,j(y)+1}\), if \((x,y)\) is in \([x_0,x_n] \times [y_0,y_m]\), linear interpolation if one of \(x,y\) are out of those bounds, and constant at the corner value if both are out of bounds.

function-bilinear-spec

• “file” [string] HDF5 or NetCDF filename of the data

• “row header” [string] name of the row dataset, the \(x_i\)

• “row coordinate” [string] one of “t”,`”x`”,`”y`”,`”z`”

• “column header” [string] name of the column dataset, the \(y_i\)

• “column coordinate” [string] one of “t”,`”x`”,`”y`”,`”z`”

• “value header” [string] name of the values dataset, the \(u_{{i,j}}\)

Example:

<ParameterList name="function-bilinear">
  <Parameter name="file" type="string" value="pressure.h5"/>
  <Parameter name="row header" type="string" value="/time"/>
  <Parameter name="row coordinate" type="string" value="t"/>
  <Parameter name="column header" type="string" value="/x"/>
  <Parameter name="column coordinate" type="string" value="x"/>
  <Parameter name="value header" type="string" value="/pressure"/>
</ParameterList>

Bilinear is space-time function

Define \(i(x) = i : x_i < x <= x_{{i+1}}\) and similarly \(j(y) = j : y_j < y <= y_{{j+1}}\) for monotonically increasing \(x_i\) and \(y_j\).

Given a two-dimensional array \(u_{{i,j}}\), \(f\) is then defined by bilinear interpolation on \(u_{i(x),j(y)}, u_{i(x)+1,j(y)}, u_{i(x),j(y)+1}, u_{i(x)+1,j(y)+1}\), if \((x,y)\) is in \([x_0,x_n] \times [y_0,y_m]\), linear interpolation if one of \(x,y\) are out of those bounds, and constant at the corner value if both are out of bounds.

function-bilinear-and-time-spec

• “file” [string] HDF5 or NetCDF filename of the data

• “time header” [string] time Name of the temporal dimension indices, the \(t_i\).

• “row header” [string] x name of the row dataset, the \(x_i\)

• “row coordinate” [string] x one of “x”,`”y`”,`”z`”

• “column header” [string] y name of the column dataset, the \(y_i\)

• “column coordinate” [string] y one of “x”,`”y`”,`”z`”

• “value header” [string] name of the values dataset, the \(u_{{i,j}}\)

• “forms” [string] linear Describes the temporal interpolant, one of “linear” or “constant”, where “linear” is therefore trilinear interpolation (2x space and time) and “constant” indicates that the value on an interval is provided by the left point’s (earlier in time) value.

Example:

<ParameterList name="function-bilinear-and-time">
  <Parameter name="file" type="string" value="head.h5"/>
  <Parameter name="time header" type="string" value="time"/>
  <Parameter name="row header" type="string" value="x"/>
  <Parameter name="column header" type="string" value="y"/>
  <Parameter name="value header" type="string" value="rain"/>
</ParameterList>

An example HDF5 file, called head.h5, might then look like:

time: array(4) = [0, 60, 120, 180]

x: array(3) = [0, 1, 2]

y: array(2) = [0, 1]

rain: (group)

| 0: array(3,2) = … # values at time 0

| 1: array(3,2) = … # values at time 60

| 2: array(3,2) = … # values at time 120

| 3: array(3,2) = … # values at time 180

An example NetCDF file, called head.nc, might look like:

time: array(4) = [0, 60, 120, 180]

x: array(3) = [0, 1, 2]

y: array(2) = [0, 1]

rain: array(NTIMES, NX, NY) = …

Composition function

Function composition simply applies one function to the result of another.

\[f(x) = f_1( f_2(x) )\]

where \(f_1\) is defined by the “function1” sublist, and \(f_2\) by the “function2” sublist.

function-composition-spec

• “function1” [function-typedinline-spec] \(f_1\) in \(f(x) = f_1(f_2(x))\)

• “function2” [function-typedinline-spec] \(f_2\) in \(f(x) = f_1(f_2(x))\)

<ParameterList name="function-composition">
  <ParameterList name="function1">
    function-specification
  </ParameterList>
  <ParameterList name="function2">
    function-specification
  </ParameterList>
</ParameterList>

Constant function

Constant function is defined as \(f(x) = a\), for all \(x\).

function_constant-spec

• “value” [double] The constant to be applied.

Example:

<ParameterList name="function-constant">
  <Parameter name="value" type="double" value="1.0"/>
</ParameterList>

Distance function

A distance function calculates distance from reference point \(x_0\) using by the following expression:

\[f(x) = \sqrt( \sum_{j=0}^{n} m_j (x_j - x_{0,j})^2 )\]

Note that the first parameter in \(x\) can be time.

function-distance-spec

• “x0” [Array(double)] Point from which distance is measured.

• “metric” [Array(double)] Linear scaling metric, typically all 1s.

Here is an example of a distance function using isotropic metric:

Example:

<ParameterList name="function-distance">
  <Parameter name="x0" type="Array(double)" value="{1.0, 3.0, 0.0}"/>
  <Parameter name="metric" type="Array(double)" value="{1.0, 1.0, 1.0}"/>
</ParameterList>

Expression function

This function parses a string expression. The function has min(N, D + 1) arguments t, x, y, and z. The argument t is required. D is the space dimension, and N is the user specified number of arguments which could be less than D + 1.

Example:

<ParameterList name="function-exprtk">
  <Parameter name="number of arguments" type="int" value="3"/>
  <Parameter name="formula" type="string" value="t + x + 2 * y"/>
</ParameterList>

Multiplicative function

A multiplicative function simply multiplies two other function results together.

\[f(x) = f_1(x) * f_2(x)\]

where \(f_1\) is defined by the “function1” sublist, and \(f_2\) by the “function2” sublist.

function-multiplicative-spec

• “function1” [function-typedinline-spec] \(f_1\) in \(f(x) = f_1(x) * f_2(x)\)

• “function2” [function-typedinline-spec] \(f_2\) in \(f(x) = f_1(x) * f_2(x)\)

Example:

<ParameterList name="function-multiplicative">
  <ParameterList name="function1">
    function-specification
  </ParameterList>
  <ParameterList name="function2">
    function-specification
  </ParameterList>
</ParameterList>

Multi-variable linear function

A multi-variable linear function is formally defined by

\[f(x) = y_0 + \sum_{{j=0}}^{{n-1}} g_j (x_j - x_{{0,j}})\]

with the constant term “math:y_0 and gradient \(g_0,\, g_1\,..., g_{{n-1}}\). If the reference point \(x_0\) is specified, it must have the same number of values as the gradient. Otherwise, it defaults to zero. Note that one of the parameters in a multi-valued linear function can be time.

function-linear-spec

• “y0” [double] y_0 in f = y0 + g * (x - x0)

• “gradient” [Array(double)] g in f = y0 + g * (x - x0)

• “x0” [Array(double)] x0 in f = y0 + g * (x - x0)

Conditions:

len(x0) == len(gradient)

Example:

<ParameterList name="function-linear">
  <Parameter name="y0" type="double" value="1.0"/>
  <Parameter name="gradient" type="Array(double)" value="{1.0, 2.0, 3.0}"/>
  <Parameter name="x0" type="Array(double)" value="{2.0, 3.0, 1.0}"/>
</ParameterList>

Multi-variable monomial function

A multi-variable monomial function is given by the following expression:

\[f(x) = c \prod_{j=0}^{n} (x_j - x_{0,j})^{p_j}\]

with the constant factor \(c\), the reference point \(x_0\), and integer exponents \(p_j\). Note that the first parameter in \(x\) can be time.

function-monomial-spec

• “c” [double] c in \(f = c \prod_{j=0}^{n} (x_j - x_{0,j})^{p_j}\)

• “x0” [Array(double)] x0 in \(f = c \prod_{j=0}^{n} (x_j - x_{0,j})^{p_j}\)

• “exponents” [Array(int)] p in \(f = c \prod_{j=0}^{n} (x_j - x_{0,j})^{p_j}\)

Conditions:

len(x0) == len(exponents)

Here is an example of monomial of degree 6 in three variables:

<ParameterList name="function-monomial">
  <Parameter name="c" type="double" value="1.0"/>
  <Parameter name="x0" type="Array(double)" value="{1.0, 3.0, 0.0}"/>
  <Parameter name="exponents" type="Array(int)" value="{2, 3, 1}"/>
</ParameterList>

Polynomial function

A generic polynomial function is given by the following expression:

\[f(x) = \sum_{{j=0}}^n c_j (x - x_0)^{{p_j}}\]

where \(c_j\) are coefficients of monomials, \(p_j\) are integer exponents, and \(x_0\) is the reference point.

function-polynomial-spec

• “coefficients” [Array(double)] c_j polynomial coefficients

• “exponents” [Array(int)] p_j polynomail exponents

• “reference point” [double] x0 to which polynomial argument is normalized.

Example:

<ParameterList name="function-polynomial">
  <Parameter name="coefficients" type="Array(double)" value="{1.0, 1.0}"/>
  <Parameter name="exponents" type="Array(int)" value="{2, 4}"/>
  <Parameter name="reference point" type="double" value="0.0"/>
</ParameterList>

Separable function

A separable function is defined as the product of other functions such as

\[f(x_0, x_1,...,x_{{n-1}}) = f_1(x_0)\, f_2(x_1,...,x_{{n-1}})\]

where \(f_1\) is defined by the “function1” sublist, and \(f_2\) by the “function2” sublist.

function-separable-spec

• “function1” [function-spec] \(f_1\) in \(f(x) = f_1(x0) * f_2(x1...)\)

• “function2” [function-spec] \(f_2\) in \(f(x) = f_1(x0) * f_2(x1...)\)

<ParameterList name="function-separable">
  <ParameterList name="function1">
    function-specification
  </ParameterList>
  <ParameterList name="function2">
    function-specification
  </ParameterList>
</ParameterList>

Smooth step function

A smooth \(C^2\) function f(x) on interval \([x_0,\,x_1]\) is defined such that f(x) = y_0 for x < x0, f(x) = y_1 for x > x_1, and monotonically increasing for \(x \in [x_0, x_1]\) through cubic interpolation.

function-smooth-step-spec

• “x0” [double] First fitting point

• “y0” [double] First fitting value

• “x1” [double] Second fitting point

• “y1” [double] Second fitting value

Example:

<ParameterList name="function-smooth-step">
  <Parameter name="x0" type="double" value="0.0"/>
  <Parameter name="y0" type="double" value="0.0"/>
  <Parameter name="x1" type="double" value="1.0"/>
  <Parameter name="y1" type="double" value="2.0"/>
</ParameterList>

Standard math functions

These functions allow to set up non-trivial time-dependent boundary conditions which increases a set of analytic solutions that can be used in convergence analysis tests.

\[f(x) = A * operator( p * (x - s) )\]

or

\[f(x) = A * operator(x-s, p)\]

Note that these operate only on the first coordinate, which is often time. Function composition can be used to apply these to other coordinates (or better yet a dimension could/should be added upon request).

function-standard-math-spec

• “operator” [string] specifies the name of a standard mathematical function. Available options are:

  • trigonometric operators: “cos”, “sin”, “tan”, “acos”, “asin”, “atan”

  • hyperbolic trig operators: “cosh”, “sinh”, “tanh”

  • power/log operators: “pow”, “exp”, “log”, “log10”, “sqrt”,

  • integral operators: “ceil”, “floor”, “mod”,

  • “abs”, “fabs”, “positive” (0 for negative values), “negative” (0 for positive values), “heaviside”, “sign”

• “amplitude” [double] specifies a multiplication factor a in formula a f(x). The multiplication factor is ignored by function mod. Default value is 1.

• “parameter” [double] 1.0 specifies additional parameter p for math functions with two arguments. These functions are “a pow(x[0], p)” and “a mod(x[0], p)”. Alternative, scales the argument before application, for use in changing the period of trig functions.

• “shift” [double] specifies a shift of the function argument. Default is 0.

Example:

<ParameterList name="function-standard-math">
  <Parameter name="operator" type="string" value="sqrt"/>
  <Parameter name="amplitude" type="double" value="1e-7"/>
  <Parameter name="shift" type="double" value="0.1"/>
</ParameterList>

This example defines function 1e-7 sqrt(t-0.1).

Static head function

\(f(z) = p0 + rho * g * (z0 - z)\)

Note that dimension 0 is usually time.

static_head-spec

• “p0” [double] Pressure at z0

• “density” [double] Density of water

• “gravity” [double] Gravity

• “space dimension” [int] Dimensionality, usually 3

• “water table elevation” [function-spec] Water table elevation function.

Example:

<ParameterList name="function-static-head">
  <Parameter name="p0" type="double" value="101325.0"/>
  <Parameter name="density" type="double" value="1000.0"/>
  <Parameter name="gravity" type="double" value="9.8"/>
  <Parameter name="space dimension" type="int" value="3"/>
  <ParameterList name="water table elevation">
    <ParameterList name="function-constant">
      <Parameter name="value" type="double" value="1.0"/>
    </ParameterList>
  </ParameterList>
</ParameterList>

Tabular function

A piecewise function of one variable.

A tabular function is tabulated on a series of intervals; given values \({{x_i}}, {{y_i}},, i=0, ... n-1\) and functional forms \({{f_j}},, j=0, ... n-2\) a tabular function \(f(x)\) is defined as:

\[\begin{split}\begin{matrix} f(x) &=& y_0, & x \le x_0,\\ f(x) &=& f_{{i-1}}(x) & x \in (x_{{i-1}}, x_i],\\ f(x) &=& y_{{n-1}}, & x > x_{{n-1}}. \end{matrix}\end{split}\]

The functional forms \({f_j}\) may be constant, which uses the left endpoint, i.e.

\(f_i(x) = y_i\),

linear, i.e.

\(f_i(x) = ( y_i * (x - x_i) + y_{{i+1}} * (x_{{i+1}} - x) ) / (x_{{i+1}} - x_i)\)

or arbitrary, in which the \(f_j\) must be provided.

The \(x_i\) and \(y_i\) may be provided in one of two ways – explicitly in the input spec or from an HDF5 file. The length of these must be equal, and the \(x_i\) must be monotonically increasing. Forms, as defined on intervals, must be of length equal to the length of the \(x_i\) less one.

Explicitly specifying the data:

function-tabular-spec

• “x values” [Array(double)] the \(x_i\)

• “y values” [Array(double)] the \(y_i\)

• “forms” [Array(string)] optional Form of the interpolant, either “constant”, “linear”, or “USER_DEFINED” Default is linear for each * interval. Note the length of this array should be one per interval, or one less than len of x and y values.

• “USER_DEFINED” [function-spec] optional user-provided functional forms on the interval

• “x coordinate” [string] t, “x”, “y”, “z” defines which coordinate direction the \(x_i\) are formed, defaulting to time.

The below example defines a function that is zero on interval \((-\infty,\,0]\), linear on interval \((0,\,1]\), constant (f(x)=1) on interval \((1,\,2]\), square root of t on interval \((2,\,3]\), and constant (f(x)=2) on interval \((3,\,\infty]\).

Example:

<ParameterList name="function-tabular">
  <Parameter name="x values" type="Array(double)" value="{0.0, 1.0, 2.0, 3.0}"/>
  <Parameter name="x coordinate" type="string" value="t"/>
  <Parameter name="y values" type="Array(double)" value="{0.0, 1.0, 2.0, 2.0}"/>
  <Parameter name="forms" type="Array(string)" value="{linear, constant, USER_FUNC}"/>

  <ParameterList name="USER_FUNC">
    <ParameterList name="function-standard-math">
      <Parameter name="operator" type="string" value="sqrt"/>
    </ParameterList>
  </ParameterList>
</ParameterList>

Loading table from file. (Note that “USER_DEFINED” is not an option here, but could be made so if requested).

function-tabular-fromfile-spec

• “file” [string] filename of either HDF5 or NetCDF file.

• “x header” [string] name of the dataset for the \(x_i\) in the file

• “y header” [string] name of the dataset for the \(y_i\) in the file

• “forms” [string] optional, Form of the interpolant, either “constant” or “linear”

The example below would perform linear-interpolation on the intervals provided by data within an hdf5 file “my_data.h5”.

Example:

<ParameterList name="function-tabular">
  <Parameter name="file" type="string" value="my_data.h5"/>
  <Parameter name="x coordinate" type="string" value="t"/>
  <Parameter name="x header" type="string" value="/time"/>
  <Parameter name="y header" type="string" value="/data"/>
</ParameterList>

Time functions

Boundary condition functions utilize a parameterized model for time variations that is either piecewise constant or piecewise linear. For example:

<Parameter name="times" type="Array(double)" value="{1, 2, 3}"/>
<Parameter name="time values" type="Array(double)" value="{10, 20, 30}"/>
<Parameter name="time functions" type="Array(string)" value="{constant, linear}"/>

This defines four time intervals: (-inf,1), (1,2), (2,3), (3,+inf). By assumption the function is constant over the first and last intervals. The remaining two intervals are specified by the parameter time functions. Thus, the value here is 10 anytime prior to t=2. The value increases linearly from 10 to 20 over the interval t=2 to t=3, and then is constant at 30 for t>3.

SOLVERS

Nonlinear solvers are used within implicit time integration schemes to drive the residual to zero and thereby solve for the primary variable at the new time.

solver-typed-spec

• “solver type” [string] Type of the solver. One of:

  • “Newton” See Solver: Newton and Inexact Newton

  • “JFNK” See Solver: Jacobian-Free Newton Krylov

  • “line search” See Solver: Newton with Line Search

  • “continuation” See Solver: Nonlinear Continuation

  • “nka” See Solver: Nonlinear Krylov Acceleration

  • “aa” See Solver: Anderson Acceleration

  • “nka line search” See Solver: NKA with Line Search”

  • “nka_ls_ats” See Solver: NKA with Line Search, ATS

  • “nka_bt_ats” See Solver: NKA with backtracking, ATS

  • “nox” See Solver: NOX

• “_solver_type_ parameters” [_solver_type_-spec] A sublist containing parameters specific to the type.

Warning

“JFNK”, “line search”, and “continuation” methods have not been beaten on as much as other methods. “nka_ls_ats” is somewhat deprecated and probably shouldn’t be used. Prefer “nka” for simple problems, “nka_bt_ats” for freeze-thaw problems or other problems with strong nonlinearities, and “Newton” when you have a good Jacobian. While “nox” hasn’t been used extensively, it may be quite useful.

Iterative solvers

This list contains sublists for various linear algebra solvers such as PCG, GMRES, and NKA. Note that only PK can provide a preconditioner for a linear solver; hence, we cannot specify it here.

• “iterative method” [string] defines a Krylov-based method. The available options include “pcg” and “gmres”.

• “direct method” [string] defines a direct method. The available option is “amesos”.

• “xxx parameters” [list] provides parameters for the iterative method specified by variable “iterative method”.

<ParameterList>  <!-- parent list -->
<ParameterList name="solvers">
  <ParameterList name="_GMRES WITH HYPRE AMG">
    <Parameter name="iterative method" type="string" value="gmres"/>

    <ParameterList name="gmres parameters">
      ...
    </ParameterList>
  </ParameterList>

  <ParameterList name="_PCG with HYPRE AMG">
    <Parameter name="iterative method" type="string" value="pcg"/>
    <ParameterList name="pcg parameters">
      ...
    </ParameterList>
  </ParameterList>
</ParameterList>
</ParameterList>

Generalized minimal residuals (GMRES)

Based on the methods of Yu. Kuznetsov, 1968; Y.Saad, 1986. Deflated version of GMRES is due to R.Morgan, GMRES with deflated restarting, 2002 SISC; S.Rollin, W.Fichtner, Improving accuracy of GMRES with deflated restarting, 2007 SISC.

iterative-method-gmres-spec

• “error tolerance” [double] 1.e-6 Tolerance on which to declare success.

• “maximum number of iterations” [int] 100 Maximum iterations before declaring failure.

• “overflow tolerance” [double] 3.e50 Error above this value results in failure.

• “convergence criterial” [Array(string)] {relative rhs} A list of criteria, any of which can be applied. Valid include:

  • “relative rhs” : measure error relative to the norm of the RHS vector

  • “relative residual” : measure error relative to the norm of the residual

  • “absolute residual” : measure error directly, norm of error

  • “make one iteration” : require at least one iteration to be performed before declaring success

• “size of Krylov space” [int] 10 Size of the Krylov space used to span the residual.

• “controller training start” [int] 0 Start iteration for determining convergence rates. (Add more please!)

• “controller training end” [int] 3 Start iteration for determining convergence rates. (Add more please!)

• “preconditioning strategy” [string] left Valid are “left” and “right”-type preconditioning (see Saad 1986)

• “maximum size of deflation space” [int] 0 Size of the deflation space, see Rollin et al.

<ParameterList name="_GMRES with HYPRE AMG">  <!-- parent list -->
<ParameterList name="gmres parameters">
  <Parameter name="error tolerance" type="double" value="1e-12"/>
  <Parameter name="maximum number of iterations" type="int" value="400"/>
  <Parameter name="convergence criteria" type="Array(string)" value="{relative residual}"/>
  <Parameter name="size of Krylov space" type="int" value="10"/>
  <Parameter name="overflow tolerance" type="double" value="3.0e+50"/>
  <Parameter name="maximum size of deflation space" type="int" value="0"/>
  <Parameter name="preconditioning strategy`" type="string" value="left"/>
  <Parameter name="release Krylov vectors" type="bool" value="false"/>

  <ParameterList name="verbose object">
    <Parameter name="verbosity level" type="string" value="high"/>
  </ParameterList>
</ParameterList>

<!-- Alternative implementation
<ParameterList name="belos gmres parameters">
  <Parameter name="error tolerance" type="double" value="1e-12"/>
  <Parameter name="maximum number of iterations" type="int" value="400"/>
  <Parameter name="convergence criteria" type="Array(string)" value="{relative residual}"/>
  <Parameter name="size of Krylov space" type="int" value="10"/>
  <Parameter name="overflow tolerance" type="double" value="3.0e+50"/>
</ParameterList-->
</ParameterList>

Preconditioner conjugate gradient (PCG)

iterative-method-pcg-spec

• “error tolerance” [double] 1.e-6 Tolerance on which to declare success.

• “maximum number of iterations” [int] 100 Maximum iterations before declaring failure.

• “overflow tolerance” [double] 3.e50 Error above this value results in failure.

• “convergence criterial” [Array(string)] {relative rhs} A list of criteria, any of which can be applied. Valid include:

  • “relative rhs” : measure error relative to the norm of the RHS vector

  • “relative residual” : measure error relative to the norm of the residual

  • “absolute residual” : measure error directly, norm of error

  • “make one iteration” : require at least one iteration to be performed before declaring success

<ParameterList name="_PCG with HYPRE AMG">  <!-- parent list -->
<ParameterList name="pcg parameters">
  <Parameter name="error tolerance" type="double" value="1e-12"/>
  <Parameter name="maximum number of iterations" type="int" value="400"/>
  <Parameter name="convergence criteria" type="Array(string)" value="{relative residual,make one iteration}"/>
  <Parameter name="overflow tolerance" type="double" value="3.0e+50"/>

  <ParameterList name="verbose object">
    <Parameter name="verbosity level" type="string" value="high"/>
  </ParameterList>
</ParameterList>
</ParameterList>

Newton-Krylov acceleration (NKA)

This is effectively equivalent to GMRES with a rolling restart, where vectors fall off the end of the space.

iterative-method-nka-spec

• “error tolerance” [double] 1.e-6 Tolerance on which to declare success.

• “maximum number of iterations” [int] 100 Maximum iterations before declaring failure.

• “overflow tolerance” [double] 3.e50 Error above this value results in failure.

• “convergence criterial” [Array(string)] {relative rhs} A list of criteria, any of which can be applied. Valid include:

  • “relative rhs” : measure error relative to the norm of the RHS vector

  • “relative residual” : measure error relative to the norm of the residual

  • “absolute residual” : measure error directly, norm of error

  • “make one iteration” : require at least one iteration to be performed before declaring success

• “max nka vectors” [int] 10 Size of the NKA space used to span the residual, conceptually equivalent to the size of the Krylov space.

• “nka vector tolerance” [double] 0.05 Vectors whose dot product are within this tolerance are considered parallel, and therefore the old vector is thrown out.

<ParameterList name="_NKA">  <!-- parent list -->
<ParameterList name="nka parameters">
  <Parameter name="error tolerance" type="double" value="1e-12"/>
  <Parameter name="maximum number of iterations" type="int" value="400"/>
  <Parameter name="convergence criteria" type="Array(string)" value="{relative residual}"/>
  <Parameter name="overflow tolerance" type="double" value="3.0e+50"/>
  <Parameter name="max nka vectors" type="int" value="10"/>
  <Parameter name="nka vector tolerance" type="double" value="0.05"/>

  <ParameterList name="verbose object">
    <Parameter name="verbosity level" type="string" value="high"/>
  </ParameterList>
</ParameterList>
</ParameterList>

Direct solvers from Amesos library

Amesos library of Trilinos package provides interfaces to a few direct solvers. List “amesos parameters” contains parameters that understood by this library. These parameters may violate the camel-case convention employed by this spec. Additional parameters are:

solvers_amesos-spec

• “solver name” [string] declares name of one of the supported direct solvers. Available options are “klu”, “superludist”, “basker”, etc, see Amesos and Amesos2 manuals for details. The default value is serial solver “klu”.

• “amesos version” [int] specifies version of Amesos. Available options are 1 and 2. The default value is 1.

<ParameterList name="_AMESOS KLU">  <!-- parent list -->
<Parameter name="direct method" type="string" value="amesos"/>
<ParameterList name="amesos parameters">
  <Parameter name="solver name" type="string" value="klu"/>
  <Parameter name="amesos version" type="int" value="1"/>
</ParameterList>
</ParameterList>

Nonlinear solvers

Amanzi supports a few nonlinear solvers. Typically, a process kernel uses a factory to select a nonlinear solver. This factory uses parameter solver type to find parameters for the selected solver.

Newton-Krylov acceleration (NKA)

Uses the Nonlinear Krylov acceleration method of Carlson and Miller to do effectively a multivariant secant method, accelerating the solution of a nonlinear solve. This method can be significantly faster than Newton, especially with an approximate Jacobian.

Calef et al. “Nonlinear Krylov acceleration applied to a discrete ordinates formulation of the k-eigenvalue problem.” JCP 238 (2013): 188-209.

N. N. Carlson, K. Miller, Design and application of a gradient-weighted moving finite element code II: In two dimensions, SIAM J. Sci. Comput. 19 (3) (1998) 766–798.

solver-nka-spec

• “nonlinear tolerance” [double] 1.e-6 Defines the required error tolerance. The error is calculated by a PK.

• “monitor” [string] monitor update Specifies control of the nonlinear residual. The available options are “monitor update”, “monitor residual”, “monitor preconditioned residual”, “monitor l2 residual”, and “monitor preconditioned l2 residual”.

• “limit iterations” [int] 20 Defines the maximum allowed number of iterations.

• “diverged tolerance” [double] 1.e10 Defines the error level indicating divergence of the solver. The error is calculated by a PK. Set to a negative value to ignore this check.

• “diverged l2 tolerance” [double] 1.e10 Defines another way to identify divergence of the solver. If the relative l2 (little l) norm of the solution increment is above this value, the solver is terminated. Set to a negative value to ignore this check.

• “diverged pc tolerance” [double] 1e10 Defines another way to identify divergence of the solver. If the relative maximum norm of the solution increment (with respect to the initial increment) is above this value, the solver is terminated. Set to a negative value to ignore this check.

• “diverged residual tolerance” [double] 1e10 Defines another way to identify divergence of the solver. If the relative l2 norm of the residual (with respect to the initial residual) is above this value, the solver is terminated. Set to a negative value to ignore this check.

• “max du growth factor” [double] 1e5 Allows the solver to identify divergence pattern on earlier iterations. If the maximum norm of the solution increment changes drastically on two consecutive iterations, the solver is terminated.

• “max error growth factor” [double] 1e5 Defines another way to identify divergence pattern on earlier iterations. If the PK-specific error changes drastically on two consecutive iterations, the solver is terminated.

• “max divergent iterations” [int] 3 Defines another way to identify divergence pattern on earlier iterations. If the maximum norm of the solution increment grows on too many consecutive iterations, the solver is terminated.

• “make one iteration” [bool] false require at least one iteration to be performed before declaring success. This options makes any effect only when “monitor residual” is choose.

• “modify correction” [bool] false Allows a PK to modify the solution increment. One example is a physics-based clipping of extreme solution values.

• “lag iterations” [int] 0 Delays the NKA acceleration, but updates the Krylov space.

• “max nka vectors” [int] 10 Defines the maximum number of consecutive vectors used for a local space.

• “nka vector tolerance” [double] 0.05 Defines the minimum allowed orthogonality between vectors in the local space. If a new vector does not satisfy this requirement, the space is modified.

Anderson acceleration (AA)

This is a variation of the GMRES solver for nonlinear problems.

solver-aa-spec

• “nonlinear tolerance” [double] 1.e-6 Defines the required error tolerance. The error is calculated by a PK.

• “limit iterations” [int] 20 Defines the maximum allowed number of iterations.

• “diverged tolerance” [double] 1.e10 Defines the error level indicating divergence of the solver. The error is calculated by a PK. Set to a negative value to ignore this check.

• “diverged l2 tolerance” [double] 1.e10 Defines another way to identify divergence of the solver. If the relative L2 norm of the solution increment is above this value, the solver is terminated. Set to a negative value to ignore this check.

• “max du growth factor” [double] 1e5 Allows the solver to identify divergence pattern on earlier iterations. If the maximum norm of the solution increment changes drastically on two consecutive iterations, the solver is terminated.

• “max divergent iterations” [int] 3 Defines another way to identify divergence pattern on earlier iterations. If the maximum norm of the solution increment grows on too many consecutive iterations, the solver is terminated.

• “max aa vectors” [int] 10 Defines the maximum number of consecutive vectors used for a local space.

• “modify correction” [bool] false Allows a PK to modify the solution increment. One example is a physics-based clipping of extreme solution values.

• “relaxation parameter” [double] 1 Damping factor for increment.

Newton

The classical Newton method works well for cases where Jacobian is available and corresponds to a stable (e.g. upwind) discretization. The inexact Newton methods work for cases where the discrete Jacobian is either not available, or not stable, or computationally expensive. The discrete Jacobian is replaced by a stable approximation of the continuum Jacobian. The choice between exact and inexact is not made by the Solver, but instead by the PK. Both use the ApplyPreconditioner() method – if this applies the true Jacobian, then the method is Newton. If it applies an appoximation, it is inexact Newton.

solver-newton-spec

• “nonlinear tolerance” [double] 1.e-6 defines the required error tolerance. The error is calculated by a PK.

• “monitor” [string] monitor update specifies control of the nonlinear residual. The available options are “monitor update” and “monitor residual”.

• “limit iterations” [int] 50 defines the maximum allowed number of iterations.

• “diverged tolerance” [double] 1.e10 defines the error level indicating divergence of the solver. The error is calculated by a PK.

• “max du growth factor” [double] 1.e5 allows the solver to identify divergence pattern on earlier iterations. If the maximum norm of the solution increment changes drastically on two consecutive iterations, the solver is terminated.

• “max error growth factor” [double] 1.e5 defines another way to identify divergence pattern on earlier iterations. If the PK-specific error changes drastically on two consecutive iterations, the solver is terminated.

• “max divergent iterations” [int] 3 defines another way to identify divergence pattern on earlier iterations. If the maximum norm of the solution increment grows on too many consecutive iterations, the solver is terminated.

• “make one iteration” [bool] false require at least one iteration to be performed before declaring success. This options makes any effect only when “monitor residual” is choose.

• “modify correction” [bool] true allows a PK to modify the solution increment. One example is a physics-based clipping of extreme solution values.

• “stagnation iteration check” [int] 8 determines the number of iterations before the stagnation check is turned on. The stagnation happens when the current l2-error exceeds the initial l2-error.

Inexact Newton

The inexact Newton methods work for cases where the discrete Jacobian is either not available, or not stable, or computationally expensive. The discrete Jacobian is replaced by a stable approximation of the continuum Jacobian. This solver has the same list of parameters as the Newton solver.

The difference between these solvers is in the preconditioner parameters. Here is the list of selected parameters for the Newton-Picard solver.

<ParameterList name="operators">
  <ParameterList name="diffusion operator">
    <ParameterList name="preconditioner">
      <Parameter name="discretization primary" type="string" value="mfd: optimized for monotonicity"/>
      <Parameter name="discretization secondary" type="string" value="mfd: optimized for sparsity"/>
      <Parameter name="schema" type="Array(string)" value="{face, cell}"/>
      <Parameter name="preconditioner schema" type="Array(string)" value="{face, cell}"/>
      <Parameter name="Newton correction" type="string" value="approximate Jacobian"/>
    </ParameterList>
  </ParameterList>
</ParameterList>

<Parameter name="solver type" type="string" value="Newton"/>
<ParameterList name="Newton parameters">
  <Parameter name="nonlinear tolerance" type="double" value="1.0e-05"/>
  <Parameter name="diverged tolerance" type="double" value="1.0e+10"/>
  <Parameter name="max du growth factor" type="double" value="1.0e+03"/>
  <Parameter name="max divergent iterations" type="int" value="3"/>
  <Parameter name="limit iterations" type="int" value="20"/>
  <Parameter name="modify correction" type="bool" value="false"/>
</ParameterList>

Jacobian-free Newton-Krylov (JFNK)

Jacobian-Free Newton Krylov uses a finite difference scheme to approximate the action of the Jacobian matrix, then uses a Krylov method (which only needs the action of the Jacobian and not the Jacobian itself) to calculate the action of the inverse of the Jacobian, thereby providing a Newton-like update. As the linear Krylov scheme converges to the inverse action, the nonlinear solution converges to the same solution as a true Newton method.

This implementation simply replaces a SolverFnBase’s ApplyPreconditioner() with a new ApplyPreconditioner() which uses the Krylov method with the action of the forward operator to (hopefully) improve, relative to the supplied approximate inverse, the estimate of the inverse.

solver-jfnk-spec

• “nonlinear solver” [solver-typed-spec] The outer nonlinear solver to use.

• “inverse” [inverse-typed-spec] The Krylov method to use.

• “JF matrix parameters” [jf-matrix-spec] See jf-matrix-spec

<Parameter name="solver type" type="string" value="JFNK"/>
<ParameterList name="JFNK parameters">
  <Parameter name="typical solution value" type="double" value="1.0"/>

  <ParameterList name="JF matrix parameters">
    <Parameter name="finite difference epsilon" type="double" value="1.0e-8"/>
    <Parameter name="method for epsilon" type="string" value="Knoll-Keyes L2"/>
  </ParameterList>

  <ParameterList name="nonlinear solver">
    <Parameter name="nonlinear tolerance" type="double" value="1.0e-05"/>
    <Parameter name="diverged tolerance" type="double" value="1.0e+10"/>
    <Parameter name="limit iterations" type="int" value="20"/>
    <Parameter name="max divergent iterations" type="int" value="3"/>
  </ParameterList>

  <ParameterList name="linear operator">
    <Parameter name="iterative method" type="string" value="gmres"/>
    <ParameterList name="gmres parameters">
      ...
    </ParameterList>
  </ParameterList>
</ParameterList>
</ParameterList>

Nonlinear Object-Oriented Solution (NOX)

The interface to Trilinos NOX solver is as follows:

<Parameter name="solver type" type="string" value="nox"/>
  <ParameterList name="nox parameters">
    <Parameter name="typical solution value" type="double" value="1.0"/>

    <ParameterList name="JF matrix parameters">
      <Parameter name="finite difference epsilon" type="double" value="1.0e-8"/>
      <Parameter name="method for epsilon" type="string" value="Knoll-Keyes L2"/>
    </ParameterList>

    <ParameterList name="nonlinear solver">
      <Parameter name="solver type" type="string" value="Newton"/>
      <ParameterList name="Newton parameters">
        ...
      </ParameterList>
    </ParameterList>

    <ParameterList name="linear operator">
      <Parameter name="iterative method" type="string" value="gmres"/>
      <ParameterList name="gmres parameters">
        ...
      </ParameterList>
    </ParameterList>
  </ParameterList>
</ParameterList>

NKA with line search

Does NKA, then checks if that correction has reduced the residual by at least a tolerance. If not, this uses a Brendt minimization algorithm to try and find an \(\alpha\) that minimizes the reduction in the residual.

Note, this always monitors the residual.

Nonlinear continuation

Continuation methods are useful when the nonlinearity can be controlled by a single simple parameter. In this method, the nonlinear problem is solved with a less-nonlinear value of the parameter, and the solution of that is used as the initial guess to solve a harder problem. As each successive problem is solved, the continuation parameter is changed closer and closer to the true value.

Few if any PKs support this method currently – it requires the PK to provide more interface about how to update the continuation parameter.

solver-continuation-spec

• “nonlinear tolerance” [double] 1.e-6 defines the required error tolerance. The error is calculated by a PK.

• “number of continuation steps” [int] 5 How many steps to take from initial parameter to final parameter.

• “inner solver” [solver-typed-spec] A Solver, used at each step.

Preconditioners

This sublist contains entries for various preconditioners required by a simulation. At the moment, we support Trilinos multilevel preconditioner, Hypre BoomerAMG preconditioner, ILU preconditioner, Hypre’s ILU preconditioner, and identity preconditioner.

• “preconditioning method” [string] defines preconditioner algorithm.

• “xxx parameters” [list] provides parameters for the preconditioner specified by parameter “preconditioning method”.

<ParameterList>  <!-- parent list -->
<ParameterList name="preconditioners">
  <ParameterList name="_TRILINOS ML">
    <Parameter name="preconditioning method" type="string" value="ml"/>
    <ParameterList name="ml parameters">
      ...
    </ParameterList>
  </ParameterList>

  <ParameterList name="_HYPRE AMG">
    <Parameter name="preconditioning method" type="string" value="boomer amg"/>
    <ParameterList name="boomer amg parameters">
      ...
    </ParameterList>
  </ParameterList>

  <ParameterList name="_BLOCK ILU">
    <Parameter name="preconditioning method" type="string" value="block ilu"/>
    <ParameterList name="block ilu parameters">
      ...
    </ParameterList>
  </ParameterList>

  <ParameterList name="_DIAGONAL">
    <Parameter name="preconditioning method" type="string" value="diagonal"/>
  </ParameterList>
</ParameterList>

Hypre’s algebraic multigrid (AMG) and Euclid (ILU)

Boomer AMG is a HYPRE product consisting of a variety of Algebraic Multigrid methods. It is accessed through Ifpack.

This is provided when using the “preconditioning method”=`”boomer amg`” or “preconditioning method” = “hypre: boomer amg” in the Preconditioner spec.

preconditioner-boomer-amg-spec:

• “tolerance” [double] 0. If is not zero, the preconditioner is dynamic and approximate the inverse matrix with the prescribed tolerance (in the energy norm ???).

• “smoother sweeps” [int] 3 defines the number of smoothing loops. Default is 3.

• “cycle applications” [int] 5 defines the number of V-cycles.

• “strong threshold” [double] 0.5 defines the number of V-cycles. Default is 5.

• “relaxation type” [int] 6 defines the smoother to be used. Default is 6 which specifies a symmetric hybrid Gauss-Seidel / Jacobi hybrid method. TODO: add others!

• “coarsen type” [int] 0 defines the coarsening strategy to be used. Default is 0 which specifies a Falgout method. TODO: add others!

• “max multigrid levels” [int] optionally defined the maximum number of multigrid levels.

• “use block indices” [bool] false If true, uses the “systems of PDEs” code with blocks given by the SuperMap, or one per DoF per entity type.

• “number of functions” [int] 1 Any value > 1 tells Boomer AMG to use the “systems of PDEs” code with strided block type. Note that, to use this approach, unknowns must be ordered with DoF fastest varying (i.e. not the native Epetra_MultiVector order). By default, it uses the “unknown” approach in which each equation is coarsened and interpolated independently.

• “nodal strength of connection norm” [int] tells AMG to coarsen such that each variable has the same coarse grid - sometimes this is more “physical” for a particular problem. The value chosen here for nodal determines how strength of connection is determined between the coupled system. I suggest setting nodal = 1, which uses a Frobenius norm. This does NOT tell AMG to use nodal relaxation. Default is 0.

• “verbosity” [int] 0 prints a summary of run time settings and timing information to stdout. “1” prints coarsening info, “2” prints smoothing info, and “3’” prints both.

Example:

<ParameterList name="boomer amg parameters">
  <Parameter name="tolerance" type="double" value="0.0"/>
  <Parameter name="smoother sweeps" type="int" value="3"/>
  <Parameter name="cycle applications" type="int" value="5"/>
  <Parameter name="strong threshold" type="double" value="0.5"/>
  <Parameter name="coarsen type" type="int" value="0"/>
  <Parameter name="relaxation type" type="int" value="3"/>
  <Parameter name="verbosity" type="int" value="0"/>
  <Parameter name="number of functions" type="int" value="1"/>
</ParameterList>

ILU is a Parallel Incomplete LU, provided as part of the HYPRE project through the Ifpack interface.

This is provided when using the “preconditioning method”=`”ILU`” or =`”hypre: ILU`” in the Preconditioner spec.

preconditioner-ILU-spec:

• “ilu(k) fill level” [int] 1 The factorization level.

• “ilut drop tolerance” [double] 0 Defines a drop tolerance relative to the largest absolute value of any entry in the row being factored.

• “verbosity” [int] 0 Prints a summary of runtime settings and timing information to stdout.

Hypre’s Block ILU

The internal parameters for block ILU are as follows (see https://docs.trilinos.org/dev/packages/ifpack/doc/html/index.html for more detail):

• “fact: level-of-fill” [int] sets the level of fill used by the level-based ilu(k) strategy. This is based on powers of the graph, so the value 0 means no-fill. Default is 0.

• “fact: absolute threshold” [double] defines the value to add to each diagonal element (times the sign of the actual diagonal element). Default is 0.

• “fact: relative threshold” [double] multiplies the diagonal by this value before checking the threshold. Default is 1.

• “fact: relax value” [double] if nonzero, dropped values are added to the diagonal (times this factor). Default is 0.

• “overlap” [int] defines overlap of the additive Schwarz. Default is 0.

• “schwarz: combine mode” [string] defines how values corresponding to overlapping nodes are handled. Users should set this value to “Add” if interested in a symmetric preconditioner. Otherwise, the default value of “Zero” usually results in better convergence. If the level of overlap is set to zero, the rows of the user matrix that are stored on a given processor are treated as a self-contained local matrix and all column entries that reach to off-processor entries are ignored. Setting the level of overlap to one tells Ifpack to increase the size of the local matrix by adding rows that are reached to by rows owned by this processor. Increasing levels of overlap are defined recursively in the same way. For sufficiently large levels of overlap, the entire matrix would be part of each processor’s local ILU factorization process. Default is “Add”.

<ParameterList name="_MY ILU">  <!-- parent list -->
<ParameterList name="block ilu parameters">
  <Parameter name="fact: relax value" type="double" value="1.0"/>
  <Parameter name="fact: absolute threshold" type="double" value="0.0"/>
  <Parameter name="fact: relative threshold" type="double" value="1.0"/>
  <Parameter name="fact: level-of-fill" type="int" value="10"/>
  <Parameter name="overlap" type="int" value="0"/>
  <Parameter name="schwarz: combine mode" type="string" value="Add"/>
</ParameterList>
</ParameterList>

Trilinos ML

This is provided when using the “preconditioning method”=`”ml`” in the Preconditioner spec.

Warning

no input spec defined

See also: https://trilinos.github.io/pdfs/mlguide5.pdf

Example:

<ParameterList name="ml parameters">
  <Parameter name="ML output" type="int" value="0"/>
  <Parameter name="aggregation: damping factor" type="double" value="1.33"/>
  <Parameter name="aggregation: nodes per aggregate" type="int" value="3"/>
  <Parameter name="aggregation: threshold" type="double" value="0.0"/>
  <Parameter name="aggregation: type" type="string" value="Uncoupled"/>
  <Parameter name="coarse: type" type="string" value="Amesos-KLU"/>
  <Parameter name="coarse: max size" type="int" value="128"/>
  <Parameter name="coarse: damping factor" type="double" value="1.0"/>
  <Parameter name="cycle applications" type="int" value="2"/>
  <Parameter name="eigen-analysis: iterations" type="int" value="10"/>
  <Parameter name="eigen-analysis: type" type="string" value="cg"/>
  <Parameter name="max levels" type="int" value="40"/>
  <Parameter name="prec type" type="string" value="MGW"/>
  <Parameter name="smoother: damping factor" type="double" value="1.0"/>
  <Parameter name="smoother: pre or post" type="string" value="both"/>
  <Parameter name="smoother: sweeps" type="int" value="2"/>
  <Parameter name="smoother: type" type="string" value="Gauss-Seidel"/>
</ParameterList>

Identity

Simply copies the input vector to the output – uses the Identity matrix as a preconditioner.

This is provided when using the “preconditioning method”=`”identity`” in the Preconditioner spec.

No parameters are required.

Diagonal

Simply applys the pointwise inverse of the diagonal of the matrix as an extremely cheap matrix.

This is provided when using the “preconditioning method”=`”diagonal`” in the Preconditioner spec.

No parameters are required.

MESH

Amanzi supports both structured and unstructured numerical solution approaches. This flexibility has a direct impact on the selection and design of the underlying numerical algorithms, the style of the software implementations, and, ultimately, the complexity of the user-interface. This specification format uses and describes the unstructured mesh only.

mesh-spec

• “mesh” [list] accepts “unstructured” to indicate the meshing option that Amanzi will use. This instructs Amanzi to use data structures provided in the Trilinos or MSTK software frameworks. To the extent possible, the discretization algorithms implemented under this option are largely independent of the shape and connectivity of the underlying cells. As a result, this option supports an arbitrarily complex computational mesh structure that enables users to work with numerical meshes that can be aligned with geometrically complex man-made or geostatigraphical features. Under this option, the user typically provides a mesh file that was generated with an external software package. The following mesh file formats are currently supported: “Exodus II” (see example), “MSTK” (see example), and “MOAB” (obsolete). Amanzi also provides a rudimentary capability to generate unstructured meshes automatically.

  • “unstructured” [list] accepts instructions to either (1) read or, (2) generate an unstructured mesh.

    • “read mesh file” [list] accepts name, format of pre-generated mesh file

      • “file” [string] name of pre-generated mesh file. Note that in the case of an Exodus II mesh file, the suffix of the serial mesh file must be .exo and the suffix of the parallel mesh file must be .par. When running in serial the code will read this the indicated file directly. When running in parallel and the suffix is .par, the code will instead read the partitioned files, that have been generated with a Nemesis tool and named as filename.par.N.r where N is the number of processors and r is the rank. When running in parallel and the suffix is .exo, the code will partition automatically the serial file.

      • “format” [string] format of pre-generated mesh file (“MSTK”, “MOAB”, or “Exodus II”)

    • “generate mesh” [list] accepts parameters of generated mesh

      • “domain low coordinate” [Array(double)] Location of low corner of domain

      • “domain high coordinate” [Array(double)] Location of high corner of domain

      • “number of cells” [Array(int)] the number of uniform cells in each coordinate direction

    • “expert” [list] accepts parameters that control which particular mesh framework is to be used.

      • “framework” [string] one of “stk::mesh”, “MSTK”, “MOAB” or “Simple”.

      • “verify mesh” [bool] true or false.

      • “partitioner” [string] defines the partitioning algorithm for parallel unstructured meshes. The available options are “metis” (default), “zoltan_graph” and “zoltan_rcb”. “metis” and “zoltan_graph” perform a graph partitioning of the mesh with no regard to the geometry of the mesh. “zoltan_rcb” partitions meshes using Recursive Coordinate Bisection which can lead to better partitioning in meshes that are thin in a particular direction. Additionally, the use of “zoltan_rcb” with the MSTK framework triggers an option to detect columns of elements in a mesh and adjust the partitioning such that no column is split over multiple partitions. If no partitioner is specified, the default one is used.

      • “request edges” [bool] builds support for mesh edges. Only in 3D.

      • “contiguous global ids” [bool] enforces contiguous global ids. Default is true.

    • “submesh” [list] parameters for extracted meshes

      • “domain name” [string] specifies name of the domain. Available options are “fracture” for the fracture network or “surface” for surface models.

      • “extraction method” [string] specifies the extraction method. The only available option is “manifold mesh”. If this parameter is missing, the parent mesh framework is used for submesh extraction..

      • “regions” [Array(string)] defines a list of regions for submesh. Parameter “extraction method” requires a single name in this list.

Example of Unstructured mesh generated internally:

<ParameterList>  <!-- parent list -->
<ParameterList name="mesh">
  <ParameterList name="unstructured"/>
    <ParameterList name="generate mesh">
      <Parameter name="number of cells" type="Array(int)" value="{100, 1, 100}"/>
      <Parameter name="domain low coordinate" type="Array(double)" value="{0.0, 0.0, 0.0}"/>
      <Parameter name="domain high coordinate" type="Array(double)" value="{103.2, 1.0, 103.2}"/>
    </ParameterList>

    <ParameterList name="expert">
      <Parameter name="framework" type="string" value="MSTK"/>
      <Parameter name="partitioner" type="string" value="metis"/>
    </ParameterList>
  </ParameterList>
</ParameterList>

Example of Unstructured mesh read from an external file:

<ParameterList name="mesh">  <!-- parent list -->
<ParameterList name="unstructured">
  <ParameterList name="read mesh file">
    <Parameter name="file" type="string" value="mesh_filename"/>
    <Parameter name="format" type="string" value="Exodus II"/>
  </ParameterList>
</ParameterList>
</ParameterList>

GEOMETRIC MODEL

Domain

It is not always possible to extract space dimension from provided data. Therefore, we require the user to provide simple list domain with only one parameter spatial dimension.

geometric_model-spec

• “spatial dimension” [int] defined space dimension. The available values are 2 or 3.

<ParameterList name="domain">
  <Parameter name="spatial dimension" type="int" value="2"/>
</ParameterList>

Regions

Regions are geometrical constructs used to define subsets of the computational domain in order to specify the problem to be solved, and the output desired. Regions may represents zero-, one-, two- or three-dimensional subsets of physical space. For a three-dimensional problem, the simulation domain will be a three-dimensional region bounded by a set of two-dimensional regions. If the simulation domain is N-dimensional, the boundary conditions must be specified over a set of regions are (N-1)-dimensional.

Warning

Surface files contain labeled triangulated face sets. The user is responsible for ensuring that the intersections with other surfaces in the problem, including the boundaries, are exact (i.e. that surface intersections are watertight where applicable), and that the surfaces are contained within the computational domain. If nodes in the surface fall outside the domain, the elements they define are ignored.

Examples of surface files are given in the Exodus II file format here.

Warning

Region names must NOT be repeated.

Example:

<ParameterList>  <!-- parent list -->
  <ParameterList name="regions">
    <ParameterList name="TOP SECTION">
      <ParameterList name="region: box">
        <Parameter name="low coordinate" type="Array(double)" value="{2, 3, 5}"/>
        <Parameter name="high coordinate" type="Array(double)" value="{4, 5, 8}"/>
      </ParameterList>
    </ParameterList>
    <ParameterList name="MIDDLE SECTION">
      <ParameterList name="region: box">
        <Parameter name="low coordinate" type="Array(double)" value="{2, 3, 3}"/>
        <Parameter name="high coordinate" type="Array(double)" value="{4, 5, 5}"/>
      </ParameterList>
    </ParameterList>
    <ParameterList name="BOTTOM SECTION">
      <ParameterList name="region: box">
        <Parameter name="low coordinate" type="Array(double)" value="{2, 3, 0}"/>
        <Parameter name="high coordinate" type="Array(double)" value="{4, 5, 3}"/>
      </ParameterList>
    </ParameterList>
    <ParameterList name="INFLOW SURFACE">
      <ParameterList name="region: labeled set">
        <Parameter name="label"  type="string" value="sideset_2"/>
        <Parameter name="file"   type="string" value="F_area_mesh.exo"/>
        <Parameter name="format" type="string" value="Exodus II"/>
        <Parameter name="entity" type="string" value="face"/>
      </ParameterList>
    </ParameterList>
    <ParameterList name="OUTFLOW PLANE">
      <ParameterList name="region: plane">
        <Parameter name="point" type="Array(double)" value="{0.5, 0.5, 0.5}"/>
        <Parameter name="normal" type="Array(double)" value="{0, 0, 1}"/>
      </ParameterList>
    </ParameterList>
    <ParameterList name="BLOODY SAND">
      <ParameterList name="region: color function">
        <Parameter name="file" type="string" value="F_area_col.txt"/>
        <Parameter name="value" type="int" value="25"/>
      </ParameterList>
    </ParameterList>
    <ParameterList name="FLUX PLANE">
      <ParameterList name="region: polygon">
        <Parameter name="number of points" type="int" value="5"/>
        <Parameter name="points" type="Array(double)" value="{-0.5, -0.5, -0.5,
                                                               0.5, -0.5, -0.5,
                                                               0.8, 0.0, 0.0,
                                                               0.5,  0.5, 0.5,
                                                              -0.5, 0.5, 0.5}"/>
       </ParameterList>
    </ParameterList>
  </ParameterList>
</ParameterList>

In this example, TOP SECTION, MIDDLE SECTION and BOTTOM SECTION are three box-shaped volumetric regions. INFLOW SURFACE is a surface region defined in an Exodus II-formatted labeled set file and OUTFLOW PLANE is a planar region. BLOODY SAND is a volumetric region defined by the value 25 in color function file.

All

No parameters required.

region_all-spec

• “empty” [bool] True This is simply here to avoid issues with empty lists.

Example:

<ParameterList name="domain">  <!-- parent list -->
  <ParameterList name="region: all">
  </ParameterList>
</ParameterList>

Box

List region: box defines a region bounded by coordinate-aligned planes. Boxes are allowed to be of zero thickness in only one direction in which case they are equivalent to planes.

region-box-spec

• “low coordinate” [Array(double)] Location of the boundary point with the lowest coordinates.

• “high coordinate” [Array(double)] Location of the boundary points with the highest coordinates.

Example:

<ParameterList name="WELL">  <!-- parent list -->
  <ParameterList name="region: box">
    <Parameter name="low coordinate" type="Array(double)" value="{-5.0,-5.0, -5.0}"/>
    <Parameter name="high coordinate" type="Array(double)" value="{5.0, 5.0,  5.0}"/>
  </ParameterList>
</ParameterList>

Plane

List region: plane defines a plane using a point lying on the plane and normal to the plane.

region-plane-spec

• “normal” [Array(double)] Normal to the plane.

• “point” [Array(double)] Point in space.

Example:

<ParameterList name="TOP_SECTION"> <!-- parent list -->
  <ParameterList name="region: plane">
    <Parameter name="point" type="Array(double)" value="{2, 3, 5}"/>
    <Parameter name="normal" type="Array(double)" value="{1, 1, 0}"/>
    <ParameterList name="expert parameters">
      <Parameter name="tolerance" type="double" value="1.0e-05"/>
    </ParameterList>
  </ParameterList>
</ParameterList>

Halfspace

List region: halfspace defines a halfspace determined by a plane and outward normal to the plane.

region_half_space-spec

• “normal” [Array(double)] Normal to the plane.

• “point” [Array(double)] Point in space.

Example:

<ParameterList name="TOP_SECTION"> <!-- parent list -->
  <ParameterList name="region: halfspace">
    <Parameter name="point" type="Array(double)" value="{2, 3, 5}"/>
    <Parameter name="normal" type="Array(double)" value="{1, 1, 0}"/>
    <ParameterList name="expert parameters">
      <Parameter name="tolerance" type="double" value="1.0e-05"/>
    </ParameterList>
  </ParameterList>
</ParameterList>

Labeled Set

The list region: labeled set defines a named set of mesh entities existing in an input mesh file. This is the same file that contains the computational mesh. The name of the entity set is given by label. For example, a mesh file in the Exodus II format can be processed to tag cells, faces and/or nodes with specific labels, using a variety of external tools. Regions based on such sets are assigned a user-defined label for Amanzi, which may or may not correspond to the original label in the exodus file. Note that the file used to express this labeled set may be in any Amanzi-supported mesh format (the mesh format is specified in the parameters for this option). The entity parameter may be necessary to specify a unique set. For example, an Exodus file requires cell, face or node as well as a label (which is an integer). The resulting region will have the dimensionality associated with the entities in the indicated set.

region-labeled-set-spec

• “label” [string] Set per label defined in the mesh file.

• “file” [string] File name.

• “format” [string] Currently, we only support mesh files in the “Exodus II” format.

• “entity” [string] Type of the mesh object (cell, face, etc).

Example:

<ParameterList name="AQUIFER">
  <ParameterList name="region: labeled set">
    <Parameter name="entity" type="string" value="cell"/>
    <Parameter name="file" type="string" value="porflow4_4.exo"/>
    <Parameter name="format" type="string" value="Exodus II"/>
    <Parameter name="label" type="string" value="1"/>
  </ParameterList>
</ParameterList>

Function Color

The list region: color function defines a region based a specified integer color, value, in a structured color function file, file. The format of the color function file is given below in the “Tabulated function file format” section. As shown in the file, the color values may be specified at the nodes or cells of the color function grid. A computational cell is assigned the ‘color’ of the data grid cell containing its cell centroid (cell-based colors) or the data grid nearest its cell-centroid (node-based colors). Computational cells sets are then built from all cells with the specified color Value.

In order to avoid, gaps and overlaps in specifying materials, it is strongly recommended that regions be defined using a single color function file.

region-color-function-spec

• “file” [string] File name containing color function.

• “value” [int] Color that defines the set in the tabulated function file.

Example:

<ParameterList name="SOIL_TOP">
  <ParameterList name="region: color function">
    <Parameter name="file" type="string" value="geology_resamp_2D.tf3"/>
    <Parameter name="value" type="int" value="1"/>
  </ParameterList>
</ParameterList>

Polygon

The list region: polygon defines a polygonal region on which mesh faces and nodes can be queried. NOTE that one cannot ask for cells in a polygonal surface region. In 2D, the polygonal region is a line and is specified by 2 points. In 3D, the polygonal region is specified by an arbitrary number of points. In both cases the point coordinates are given as a linear array. The polygon can be non-convex.

This provides a set of faces with a normal for computing flux.

The polygonal surface region can be queried for a normal. In 2D, the normal is defined as [Vy,-Vx] where [Vx,Vy] is the vector from point 1 to point 2. In 3D, the normal of the polygon is defined by the order in which points are specified.

region_polygon-spec

• “number of points” [int] Number of polygon points.

• “points” [Array(double)] Point coordinates in a linear array.

Example:

<ParameterList name="XY_PENTAGON">
  <ParameterList name="region: polygon">
    <Parameter name="number of points" type="int" value="5"/>
    <Parameter name="points" type="Array(double)" value="{-0.5, -0.5, -0.5,
                                                           0.5, -0.5, -0.5,
                                                           0.8, 0.0, 0.0,
                                                           0.5,  0.5, 0.5,
                                                          -0.5, 0.5, 0.5}"/>
    <ParameterList name="expert parameters">
      <Parameter name="tolerance" type="double" value="1.0e-3"/>
    </ParameterList>
  </ParameterList>
</ParameterList>

Cylinder

List region: cylinder defines an infinite cylinder determined by a symmetry axis, point on this axis and radius.

region_cylinder-spec

• “axis” [Array(double)] symmetry axis

• “point” [Array(double)] point on a symmetry axis

• “radius” [double] cylinder radius

Example:

<ParameterList name="TOP_SECTION"> <!-- parent list -->
  <ParameterList name="region: cylinder">
    <Parameter name="normal" type="Array(double)" value="{1, 1, 0}"/>
    <Parameter name="axis" type="Array(double)" value="{0, 0, 1}"/>
    <ParameterList name="expert parameters">
      <Parameter name="tolerance" type="double" value="1.0e-05"/>
    </ParameterList>
  </ParameterList>
</ParameterList>

Point

List region: point defines a point in space. This region consists of cells containing this point.

region-point-spec

• “coordinate” [Array(double)] Location of point in space.

Example:

<ParameterList name="DOWN_WIND150"> <!-- parent list defining the name -->
  <ParameterList name="region: point">
    <Parameter name="coordinate" type="Array(double)" value="{-150.0, 0.0, 0.0}"/>
  </ParameterList>
</ParameterList>

Logical

The list region: logical defines logical operations on regions allow for more advanced region definitions. At this time the logical region allows for logical operations on a list of regions. union and intersection are self-evident. In the case of subtraction, subtraction is performed from the first region in the list. The complement is a special case in that it is the only case that operates on single region, and returns the complement to it within the domain ENTIRE_DOMAIN. Currently, multi-region booleans are not supported in the same expression.

region-logical-spec

• “operation” [string] defines operation on the list of regions. One of: “union”, “intersect”, “subtract”, “complement”

• “regions” [Array(string)] specifies the list of involved regions.

Example:

<ParameterList name="LOWER_LAYERs">
  <ParameterList name="region: logical">
    <Parameter name="operation" type="string" value="union"/>
    <Parameter name="regions" type="Array(string)" value="{Middle1, Middle2, Bottom}"/>
  </ParameterList>
</ParameterList>

Boundary

List region: boundary defines a set of all boundary faces. Using this definition, faces located on the domain boundary are extracted.

region-boundary-spec

• “entity” [string] Type of the mesh object. Unclear whether this is used or can be other things than “face”?

Example:

<ParameterList name="DOMAIN_BOUNDARY"> <!-- parent list names the region -->
  <ParameterList name="region: boundary">
    <Parameter name="entity" type="string" value="face"/>
  </ParameterList>
</ParameterList>

Enumerated set

List region: enumerated set defines a set of mesh entities via the list of input global ids. Note that global ids are not defined correctly when parallel mesh is created on a fly.

region-enumerated-spec

• “entity” [string] Type of the mesh object. One of: “cell”, “face”, “edge”, “node”

• “entity gids” [Array(int)] List of the global IDs of the entities.

Example:

<ParameterList name="WELL"> <!-- parent list -->
  <ParameterList name="region: enumerated set">
    <Parameter name="entity" type="string" value="face"/>
    <Parameter name="entity gids" type="Array(int)" value="{1, 12, 23, 34}"/>
  </ParameterList>
</ParameterList>

Box volume fractions

List region: box volume fraction defines a region bounded by a box not aligned with coordinate axes. Boxes are allowed to be of zero thickness in only one direction in which case they are equivalent to rectangles on a plane or segments on a line.

region-box-volume-fractions-spec

• “corner coordinate” [Array(double)] Location of one box corner.

• “opposite corner coordinate” [Array(double)] Location of the opposite box corner.

• “normals” [Array(double)] Normals to sides in a linear array. Default is columns of the identity matrix. The normals may be scaled arbitrarily but must be orthogonal to one another and form the right coordinate frame.

Example:

<ParameterList name="BASIN">  <!-- parent list -->
  <ParameterList name="region: box volume fractions">
    <Parameter name="corner coordinate" type="Array(double)" value="{-1.0,-1.0, 1.0}"/>
    <Parameter name="opposite corner coordinate" type="Array(double)" value="{1.0, 1.0, 1.0}"/>
    <Parameter name="normals" type="Array(double)" value="{1.0, 0.0, 0.0
                                                           0.0, 2.0, 0.0,
                                                           0.0, 0.0, 3.0}"/>
  </ParameterList>
</ParameterList>

This example defines a degenerate box, a square on a surface z=1.

Line Segment

List region: line segment desribes a region defined by a line segment. This region is a set of cells which intersect with a line segment. The line segment is allowed to intersect with one or more cells. Zero length line segments are allowed. The line segment is defined by its ends points.

region-line-segment-spec

• “end coordinate” [Array(double)] Location of one end of a line segment.

• “opposite end coordinate” [Array(double)] Location of the opposite end of a line segment.

Example:

<ParameterList name="WELL"> <!-- parent list -->
   <ParameterList name="region: line segment">
     <Parameter name="end coordinate" type="Array(double)" value="{497542.44, 5393755.77, 0.0}"/>
     <Parameter name="opposite end coordinate" type="Array(double)" value="{497542.44, 5393755.77, 100.0}"/>
   </ParameterList>
 </ParameterList>

Notes and example

• Surface files contain labeled triangulated face sets. The user is responsible for ensuring that the intersections with other surfaces in the problem, including the boundaries, are exact (i.e. that surface intersections are watertight where applicable), and that the surfaces are contained within the computational domain. If nodes in the surface fall outside the domain, the elements they define are ignored.

  Examples of surface files are given in the Exodus II file format here.

• Region names must NOT be repeated.

<ParameterList>  <!-- parent list -->
<ParameterList name="regions">
  <ParameterList name="_TOP SECTION">
    <ParameterList name="region: box">
      <Parameter name="low coordinate" type="Array(double)" value="{2, 3, 5}"/>
      <Parameter name="high coordinate" type="Array(double)" value="{4, 5, 8}"/>
    </ParameterList>
  </ParameterList>
  <ParameterList name="_MIDDLE SECTION">
    <ParameterList name="region: box">
      <Parameter name="low coordinate" type="Array(double)" value="{2, 3, 3}"/>
      <Parameter name="high coordinate" type="Array(double)" value="{4, 5, 5}"/>
    </ParameterList>
  </ParameterList>
  <ParameterList name="_BOTTOM SECTION">
    <ParameterList name="region: box">
      <Parameter name="low coordinate" type="Array(double)" value="{2, 3, 0}"/>
      <Parameter name="high coordinate" type="Array(double)" value="{4, 5, 3}"/>
    </ParameterList>
  </ParameterList>
  <ParameterList name="_INFLOW SURFACE">
    <ParameterList name="region: labeled set">
      <Parameter name="label"  type="string" value="_SIDESET2"/>
      <Parameter name="file"   type="string" value="_MESH.exo"/>
      <Parameter name="format" type="string" value="Exodus II"/>
      <Parameter name="entity" type="string" value="face"/>
    </ParameterList>
  </ParameterList>
  <ParameterList name="_OUTFLOW PLANE">
    <ParameterList name="region: plane">
      <Parameter name="point" type="Array(double)" value="{0.5, 0.5, 0.5}"/>
      <Parameter name="normal" type="Array(double)" value="{0, 0, 1}"/>
    </ParameterList>
  </ParameterList>
  <ParameterList name="_BLOODY SAND">
    <ParameterList name="region: color function">
      <Parameter name="file" type="string" value="_FAREA_COLOR.txt"/>
      <Parameter name="value" type="int" value="25"/>
    </ParameterList>
  </ParameterList>
  <ParameterList name="_FLUX PLANE">
    <ParameterList name="region: polygon">
      <Parameter name="number of points" type="int" value="5"/>
      <Parameter name="points" type="Array(double)" value="{-0.5, -0.5, -0.5,
                                                             0.5, -0.5, -0.5,
                                                             0.8, 0.0, 0.0,
                                                             0.5,  0.5, 0.5,
                                                            -0.5, 0.5, 0.5}"/>
    </ParameterList>
  </ParameterList>
  <ParameterList name="_ENTIRE MESH">
    <ParameterList name="region: all"/>
  </ParameterList>
</ParameterList>
</ParameterList>

In this example, _TOP SECTION, _MIDDLE SECTION and _BOTTOM SECTION are three box-shaped volumetric regions. _INFLOW SURFACE is a surface region defined in an Exodus II-formatted labeled set file and _OUTFLOW PLANE is a planar region. _BLOODY SAND is a volumetric region defined by the value 25 in color function file.

DATABASES

Thermodynamic database

The thermodynamic database has a few sections, some of them are optional.

Primary species

The “primary species” section is a list of primary species, one sublist for a species. Each sublist is named after the species and contains the following parameters:

primary_species-spec

• “ion size parameter” [double] is an empirical parameter that provides agreement between measured activity coefficients and ionic strength. In theory, it is the diameter of the hydrated ion.

• “charge” [int] is the ion charge. The net charge of an ion is non-zero since the total number of electrons is unequal to the total number of protons.

• “gram molecular weight” [double] is amount of a molecular substance whose weight, in grams, is numerically equal to the molecular weight of that substance.

<ParameterList name="thermodynamic database">
  <ParameterList name="primary species">
    <ParameterList name="H+">
      <Parameter name="ion size parameter" type="double" value="9.0"/>
      <Parameter name="charge" type="int" value="1"/>
      <Parameter name="gram molecular weight" type="double" value="1.0079"/>
    </ParameterList>
    <ParameterList name="Ca++">
      <Parameter name="ion size parameter" type="double" value="6.0"/>
      <Parameter name="charge" type="int" value="2"/>
      <Parameter name="gram molecular weight" type="double" value="40.078"/>
    </ParameterList>
  </ParameterList>
</ParameterList>

Here is the short list of species that could be used for models. Each line has four data fields: name of a species, ion size parameter, charge, and atomic mass [u].

Al+++      9.0     3.0    26.9815
HCO3-      4.0    -1.0    61.0171
HPO4--     4.0    -2.0    95.9793
Cl-        3.0    -1.0    35.4527
CO2(aq)    3.0     0.0    44.01
Cs137      2.5     1.0   132.9054
F-         3.5    -1.0    18.9984
Fe++       6.0     2.0    55.847
K+         3.0     1.0    39.0983
Mg++       8.0     2.0    24.30
Na+        4.0     1.0    22.9898
N2(aq)     3.0     0.0    28.0135
NO3-       3.0    -1.0    62.0049
O2(aq)     3.0     0.0    31.9988
Pb_210     1.0     0.0   210.00
Pu_238     1.0     0.0   238.00
Ra_226     1.0     0.0   226.00
SiO2(aq)   3.0     0.0    60.0843
SO4--      4.0    -2.0    96.0636
Sr90       5.0     2.0    87.6200
Tc_99      1.0     0.0    99.00
Th_230     1.0     0.0   230.00
Tritium    9.0     0.0     1.01
U_234      1.0     0.0   234.00
UO2++      4.5     2.0   270.028
Zn++       6.0     2.0    65.39

Aqueous equilibrium complexes

The “aqueous equilibrium complexes” section is a list of aqueous equiliblium reactions. Each sublist is named after a secondary species and contains the following parameters:

aqueous_complex-spec

• “ion size parameter” [double] is an empirical parameter that provides agreement between measured activity coefficients and ionic strength. In theory, it is the diameter of the hydrated ion.

• “charge” [int] is the ion charge. The net charge of an ion is non-zero since the total number of electrons is unequal to the total number of protons.

• “gram molecular weight” [double] is amount of a molecular substance whose weight, in grams, is numerically equal to the molecular weight of that substance.

• “reaction” [string] is the equilibrium reaction equation.

• “equilibrium constant” [double] is the value of reaction quotient at chemical equilibrium. The reaction quotient is defined mathematically as the ratio of the activities of the product species to that of the reactant species. The stoichiometric coefficients are taken into account as exponents of activities.

In a non-isothermal simulation, the equlibrium constant is a function of temperature. In such a case we have one additional parameter “temperature” and parameter “equilibrium constant” becomes a list that defines a piecewice linear function logK(T). This is the tabular function with linear forms.

aqueous_complex_temprerature-spec

• “T” [Array(double)] is the array of temperature points.

• “Keq” [Array(double)] is the matching array of equilibrium constant values.

<ParameterList name="thermodynamic database">
  <ParameterList name="aqueous equilibrium complexes">
    <ParameterList name="OH-">
      <Parameter name="ion size parameter" type="double" value="3.5"/>
      <Parameter name="charge" type="int" value="-1"/>
      <Parameter name="gram molecular weight" type="double" value="17.0073"/>
      <Parameter name="reaction" type="string" value="1.0 H2O  -1.0 H+"/>
      <Parameter name="equilibrium constant" type="double" value="13.9951"/>
      <Parameter name="temperature" type="double" value="298.15"/>
      <ParameterList name="equilibrium constant">
        <Parameter name="T" type="Array(double)" value="{273.15,  298.15,  333.15,  373.15,  423.15,  473.15,  523.15,  573.15}"/>
        <Parameter name="Keq" type="Array(double)" value="{14.9398, 13.9951, 13.0272, 12.2551, 11.6308, 11.2836, 11.1675, 11.3002}"/>
      </ParameterList>
    <ParameterList name="HCO3-">
      <Parameter name="ion size parameter" type="double" value="4.0"/>
      <Parameter name="charge" type="int" value="-1"/>
      <Parameter name="gram molecular weight" type="double" value="61.0171"/>
      <Parameter name="reaction" type="string" value="1.0 H2O  -1.0 H+  1.0 CO2(aq)"/>
      <Parameter name="equilibrium constant" type="double" value="6.3447"/>
      <Parameter name="temperature" type="double" value="298.15"/>
    </ParameterList>
  </ParameterList>
</ParameterList>

Below are a few examples of aqueous reaction. Each line has five fields: reaction equation, logarithm of the equlibrium constant at 25C, ion size parameter, ion charge, and the gram molecular weight.

AlHPO4+    =  1.0 Al+++ 1.0 HPO4--    -7.4       4.0     1.0   122.961
CaCl+      =  1.0 Ca++  1.0 Cl-        0.6956    4.0     1.0    75.5307
CaCl2(aq)  =  1.0 Ca++  2.0 Cl-        0.6436    3.0     0.0   110.9834
CaCO3(aq)  =  1.0 Ca+2  1.0 CO3-2     -3.151     3.5     0.0    17.0073
CaHCO3+    =  1.0 Ca++  1.0 HCO3-     -1.0467    4.0     1.0   101.0951
CaHPO4(aq) =  1.0 Ca++  1.0 HPO4--    -2.7400    3.0     0.0   136.0573
CO3--      = -1.0 H+    1.0 HCO3-     10.3288    4.5    -2.0    60.0092
FeCl+      =  1.0 Fe++  1.0 Cl-        0.1605    4.0     1.0    91.2997
FeCl2(aq)  =  1.0 Fe++  2.0 Cl-        2.4541    3.0     0.0   126.752
FeCl4--    =  1.0 Fe++  4.0 Cl-       -1.9       4.0    -2.0   197.658
FeHCO3+    =  1.0 Fe++  1.0 HCO3-     -2.72      4.0     1.0   116.864
FeHPO4(aq) =  1.0 Fe++  1.0 HPO4--    -3.6       3.0     0.0   151.826
FeF+       =  1.0 Fe++  1.0 F-        -1.36      4.0     1.0    74.8454
FeSO4(aq)  =  1.0 Fe++  1.0 SO4--     -2.2       3.0     0.0   151.911
H2PO4-     =  1.0 H+    1.0 HPO4--    -7.2054    4.0    -1.0    96.9872
H3PO4(aq)  =  2.0 H+    1.0 HPO4--    -9.3751    3.0     0.0    97.9952
H2SO4(aq)  =  2.0 H+    1.0 SO4--      1.0209    3.0     0.0    98.0795
HCl(aq)    =  1.0 H+    1.0 Cl-       -0.67      3.0     0.0    36.4606
HNO3(aq)   =  1.0 H+    1.0 NO3-       1.3025    3.0     0.0    63.0129
HSO4-      =  1.0 H+    1.0 SO4--     -1.9791    4.0    -1.0    97.0715
KCl(aq)    =  1.0 K+    1.0 Cl-        1.4946    3.0     0.0    74.551
KHPO4-     =  1.0 K+    1.0 HPO4--    -0.78      4.0    -1.0   135.078
KSO4-      =  1.0 K+    1.0 SO4--     -0.8796    4.0    -1.0   135.162
NaHCO3(aq) =  1.0 Na+   1.0 HCO3-     -0.1541    3.0     0.0    84.0069
NaCl(aq)   =  1.0 Na+   1.0 Cl-        0.777     3.0     0.0    58.4425
NaF(aq)    =  1.0 Na+   1.0 F-         0.9976    3.0     0.0    41.9882
NaHPO4-    =  1.0 Na+   1.0 HPO4--    -0.92      4.0    -1.0   118.969
NaNO3(aq)  =  1.0 Na+   1.0 NO3-       1.044     3.0     0.0    84.9947
NaSO4-     =  1.0 Na+   1.0 SO4--     -0.82      4.0    -1.0   119.053
MgCO3(aq)  =  1.0 Mg+2  1.0 CO3-2     -2.928     3.5     0.0    17.0073
OH-        =  1.0 H2O  -1.0 H+        13.9951    3.5    -1.0    17.00730
P2O7----   = -1.0 H2O   2.0 HPO4--     3.7463    4.0    -4.0   173.9433
PO4---     = -1.0 H+    1.0 HPO4--    12.3218    4.0    -3.0    94.9714
UO2Cl+     =  1.0 Cl-   1.0 UO2++     -0.1572    4.0     1.0   305.48
UO2Cl2(aq) =  2.0 Cl-   1.0 UO2++      1.1253    3.0     0.0   340.933
UO2F+      =  1.0 F-    1.0 UO2++     -5.0502    4.0     1.0   289.026
UO2F2(aq)  =  2.0 F-    1.0 UO2++     -8.5403    3.0     0.0   308.024
UO2F3-     =  3.0 F-    1.0 UO2++    -10.7806    4.0    -1.0   327.023
UO2F4--    =  4.0 F-    1.0 UO2++    -11.5407    4.0    -2.0   346.021
UO2HPO4(aq)= 1.0 HPO4-- 1.0 UO2++     -8.4398    3.0     0.0   366.007
UO2NO3+    =  1.0 NO3-  1.0 UO2++     -0.2805    4.0     1.0   332.033
UO2SO4(aq) =  1.0 SO4-- 1.0 UO2++     -3.0703    3.0     0.0   366.091
UO2(SO4)2-- = 2.0 SO4-- 1.0 UO2++     -3.9806    4.0    -2.0   462.155

Al2(OH)2++++  = -2.0 H+    2.0 Al+++   2.0 H2O         7.6902    5.5    4.0    87.9778
Al3(OH)4(5+)  = -4.0 H+    3.0 Al+++   4.0 H2O        13.8803    6.0    5.0   148.9740
Al(OH)2+      =  2.0 H2O   1.0 Al+++  -2.0 H+         10.5945    4.0    1.0    60.9962
Al(OH)3(aq)   =  3.0 H2O   1.0 Al+++  -3.0 H+         16.1577    3.0    0.0    78.0034
AlH2PO4++     =  1.0 Al+++ 1.0 H+      1.0 HPO4--     -3.1       4.5    2.0   123.969
AlO2-         =  2.0 H2O   1.0 Al+++  -4.0 H+         22.8833    4.0   -1.0    58.9803
AlOH++        =  1.0 H2O   1.0 Al+++  -1.0 H+          4.9571    4.5    2.0    43.9889
CaCO3(aq)     = -1.0 H+    1.0 Ca++    1.0 HCO3-       7.0017    3.0    0.0   100.0872
CaH2PO4+      =  1.0 Ca++  1.0 H+      1.0 HPO4--     -1.4000    4.0    1.0   137.0652
CaP2O7--      = -1.0 H2O   1.0 Ca++    2.0 HPO4--     -3.0537    4.0   -2.0   214.0213
CaPO4-        = -1.0 H+    1.0 Ca++    1.0 HPO4--      5.8618    4.0   -1.0   135.0494
CaOH+         = -1.0 H+    1.0 Ca++    1.0 H2O        12.8500    4.0    1.0    57.0853
CO2(aq)       = -1.0 H2O   1.0 H+      1.0 HCO3-      -6.3447    3.0    0.0    44.0098
H2P2O7--      = -1.0 H2O   2.0 H+      2.0 HPO4--    -12.0709    4.0   -2.0   175.9592
H2S(aq)       =  2.0 H+    1.0 SO4--  -2.0 O2(aq)    131.329     3.0    0.0    34.0819
H3P2O7-       = -1.0 H2O   2.0 HPO4--  3.0 H+        -14.4165    4.0   -1.0   176.9671
H4P2O7(aq)    = -1.0 H2O   2.0 HPO4--  4.0 H+        -15.9263    3.0    0.0   177.9751
HAlO2(aq)     =  2.0 H2O   1.0 Al+++  -3.0 H+         16.4329    3.0    0.0    59.9883
HCO3-         =  1.0 H2O  -1.0 H+      1.0 CO2(aq)    6.34470    4.0   -1.0    61.01710
HO2-          =  1.0 H2O  -1.0 H+      0.5 O2(aq)     28.302     4.0   -1.0    33.0067
HP2O7---      = -1.0 H2O   1.0 H+      2.0 HPO4--     -5.4498    4.0   -3.0   174.9513
HS-           =  1.0 H+    1.0 SO4--  -2.0 O2(aq)    138.317     3.5   -1.0    33.0739
Fe2(OH)2++++  =  1.0 H2O   2.0 Fe++    0.5 O2(aq)    -14.0299    5.5    4.0   145.709
FeCO3(aq)     =  1.0 Fe++ -1.0 H+      1.0 HCO3-       5.5988    3.0    0.0   115.856
FeH2PO4+      =  1.0 Fe++  1.0 H+      1.0 HPO4--     -2.7       4.0    1.0   152.834
Fe(OH)2(aq)   =  2.0 H2O   1.0 Fe++   -2.0 H+         20.6       3.0    0.0    89.8617
Fe(OH)3-      =  3.0 H2O   1.0 Fe++   -3.0 H+         31.0       4.0   -1.0   106.869
Fe(OH)4--     =  4.0 H2O   1.0 Fe++   -4.0 H+         46.0       4.0   -2.0   123.876
FeOH+         =  1.0 H2O   1.0 Fe++   -1.0 H+          9.5       4.0    1.0    72.8543
FeOH++        =  0.5 H2O   1.0 Fe++    0.25 O2(aq)    -6.3       4.5    2.0    72.8543
FePO4-        =  1.0 Fe++ -1.0 H+      1.0 HPO4--      4.3918    4.0   -1.0   150.818
KHSO4(aq)     =  1.0 K+    1.0 H+      1.0 SO4--      -0.8136    3.0    0.0   136.17
KOH(aq)       =  1.0 H2O   1.0 K+     -1.0 H+         14.46      3.0    0.0    56.1056
KP2O7---      = -1.0 H2O   1.0 K+      2.0 HPO4--      1.4286    4.0   -3.0   213.042
MgOH+         =  1.0 H2O  -1.0 H+      1.0 Mg++      11.78510    4.0    1.0    41.3123
NaCO3-        = -1.0 H+    1.0 HCO3-   1.0 Na+         9.8144    4.0   -1.0    82.9990
NaOH(aq)      =  1.0 H2O   1.0 Na+    -1.0 H+         14.7948    3.0    0.0    39.9971
NH3(aq)       =  1.5 H2O   0.5 N2(aq) -0.75 O2(aq)    58.2305    3.0    0.0    17.0306
UO2CO3(aq)    = -1.0 H+    1.0 HCO3-   1.0 UO2++       0.6634    3.0    0.0   330.037
UO2(CO3)2--   = -2.0 H+    2.0 HCO3-   1.0 UO2++       3.7467    4.0   -2.0   390.046
UO2(CO3)3---- = -3.0 H+    3.0 HCO3-   1.0 UO2++       9.4302    4.0   -4.0   450.055
UO2H2PO4+     =  1.0 H+    1.0 HPO4--  1.0 UO2++     -11.6719    4.0    1.0   367.015
UO2H3PO4++    =  2.0 H+    1.0 HPO4--  1.0 UO2++     -11.3119    4.5    2.0   368.023
UO2OH+        =  1.0 H2O  -1.0 H+      1.0 UO2++       5.2073    4.0    1.0   287.035
UO2PO4-       = -1.0 H+    1.0 HPO4--  1.0 UO2++      -2.0798    4.0   -1.0   364.999
UO2(OH)2(aq)  =  2.0 H2O  -2.0 H+      1.0 UO2++      10.3146    3.0    0.0   304.042
UO2(OH)3-     =  3.0 H2O  -3.0 H+      1.0 UO2++      19.2218    4.0   -1.0   321.05
UO2(OH)4--    =  4.0 H2O  -4.0 H+      1.0 UO2++      33.0291    4.0   -2.0   338.057
(UO2)2OH+++   =  1.0 H2O  -1.0 H+      2.0 UO2++       2.7072    5.0    3.0   557.063
(UO2)2(OH)2++ =  2.0 H2O  -2.0 H+      2.0 UO2++       5.6346    4.5    2.0   574.07
(UO2)3(OH)4++ =  4.0 H2O  -4.0 H+      3.0 UO2++      11.929     4.5    2.0   878.112
(UO2)3(OH)5+  =  5.0 H2O  -5.0 H+      3.0 UO2++      15.5862    4.0    1.0   895.12
(UO2)3(OH)7-  =  7.0 H2O  -7.0 H+      3.0 UO2++      31.0508    4.0   -1.0   929.135
(UO2)4(OH)7+  =  7.0 H2O  -7.0 H+      4.0 UO2++      21.9508    4.0    1.0   1199.16
UO2(H2PO4)(H3PO4)+ = 3.0 H+ 2.0 HPO4-- 1.0 UO2++     -22.7537    4.0    1.0   465.01
UO2(H2PO4)2(aq) =    2.0 H+ 2.0 HPO4-- 1.0 UO2++     -21.7437    3.0    0.0   464.002
Zn(OH)2(aq)   =  2.0 H2O  -2.0 H+      1.0 Zn++       17.3282    3.0    0.0    99.4047
Zn(OH)3-      =  3.0 H2O  -3.0 H+      1.0 Zn++       28.8369    4.0   -1.0   116.41200
Zn(OH)4--     =  4.0 H2O  -4.0 H+      1.0 Zn++       41.6052    4.0   -2.0   133.41940
ZnOH+         =  1.0 H2O  -1.0 H+      1.0 Zn++        8.9600    4.0    1.0    82.39730

Ca2UO2(CO3)3(aq) =  2.0 Ca++ -3.0 H+     3.0 HCO3-   1.0 UO2++        0.2864   4.0    0.0   530.215
CaUO2(CO3)3--    =  1.0 Ca++ -3.0 H+     3.0 HCO3-   1.0 UO2++        3.8064   4.0   -2.0   530.215
CH4(aq)          =  1.0 H2O   1.0 H+     1.0 HCO3-  -2.0 O2(aq)     144.141    3.0    0.0     0.0
NaAlO2(aq)       =  2.0 H2O   1.0 Na+    1.0 Al+++  -4.0 H+          23.6266   3.0    0.0    81.9701
NaHP2O7--        = -1.0 H2O   1.0 Na+    1.0 H+      2.0 HPO4--      -6.8498   4.0   -2.0   197.941
NaHSiO3(aq)      =  1.0 H2O   1.0 Na+   -1.0 H+      1.0 SiO2(aq)     8.304    3.0    0.0   100.081
Fe+++            = -0.5 H2O   1.0 Fe++   1.0 H+      0.25 O2(aq)     -8.49     9.0    3.0    55.847
Fe3(OH)4(5+)     =  2.5 H2O   3.0 Fe++  -1.0 H+      0.75 O2(aq)    -19.1699   6.0    5.0   235.57
Fe(OH)2+         =  1.5 H2O   1.0 Fe++  -1.0 H+      0.25 O2(aq)     -2.82     4.0    1.0    89.8617
Fe(OH)3(aq)      =  2.5 H2O   1.0 Fe++  -2.0 H+      0.25 O2(aq)      3.51     3.0    0.0   106.869
Fe(OH)4-         =  3.5 H2O   1.0 Fe++  -3.0 H+      0.25 O2(aq)     13.11     4.0   -1.0   123.876
FeCO3+           = -0.5 H2O   1.0 Fe++   1.0 HCO3-   0.25 O2(aq)     -7.8812   4.0    1.0   115.856
MgHCO3+          =  1.0 H2O  -1.0 H+     1.0 CO2(aq) 1.0 Mg++         5.309    4.0    1.0    85.3221
N3-              =  0.5 H2O  -1.0 H+     1.5 N2(aq) -0.25 O2(aq)     77.7234   4.0   -1.0    42.0202
NH4+             =  1.5 H2O   1.0 H+     0.5 N2(aq) -0.75 O2(aq)     48.9895   2.5    1.0    18.0385
U+++             = -0.5 H2O   1.0 H+     1.0 UO2++  -0.75 O2(aq)     64.8028   5.0    3.0   238.029
U++++            = -1.0 H2O   2.0 H+     1.0 UO2++  -0.5 O2(aq)      33.949    5.5    4.0   238.029
UO2+             =  0.5 H2O  -1.0 H+     1.0 UO2++  -0.25 O2(aq)     20.0169   4.0    1.0   270.028
UO2OSi(OH)3+     =  2.0 H2O  -1.0 H+     1.0 SiO2(aq) 1.0 UO2++       2.4810   9.0    1.0   365.135
(UO2)2CO3(OH)3-  =  3.0 H2O  -4.0 H+     1.0 HCO3-   2.0 UO2++       11.2229   4.0   -1.0   651.087

Fe(SO4)2- = -0.5 H2O  1.0 Fe++  1.0 H+      2.0 SO4--   0.25 O2(aq)   -11.7037   4.0   -1.0   247.974
FeCl++    = -0.5 H2O  1.0 Fe++  1.0 H+      1.0 Cl-     0.25 O2(aq)    -7.6792   4.5    2.0    91.2997
FeCl2+    = -0.5 H2O  1.0 Fe++  1.0 H+      2.0 Cl-     0.25 O2(aq)   -10.62     4.0    1.0   126.752
FeCl4-    = -0.5 H2O  1.0 Fe++  1.0 H+      4.0 Cl-     0.25 O2(aq)    -7.7      4.0   -1.0   197.658
FeF++     = -0.5 H2O  1.0 Fe++  1.0 H+      1.0 F-      0.25 O2(aq)   -12.6265   4.5    2.0    74.8454
FeF2+     = -0.5 H2O  1.0 Fe++  1.0 H+      2.0 F-      0.25 O2(aq)   -16.8398   4.0    1.0    93.8438
FeH2PO4++ = -0.5 H2O  1.0 Fe++  2.0 H+      1.0 HPO4--  0.25 O2(aq)   -12.66     4.5    2.0   152.834
FeHPO4+   = -0.5 H2O  1.0 Fe++  1.0 H+      1.0 HPO4--  0.25 O2(aq)   -18.67     4.0    1.0   151.826
FeNO3++   = -0.5 H2O  1.0 Fe++  1.0 H+      1.0 NO3-    0.25 O2(aq)    -9.49     4.5    2.0   117.852
FeSO4+    = -0.5 H2O  1.0 Fe++  1.0 H+      1.0 SO4--   0.25 O2(aq)   -10.4176   4.0    1.0   151.911
MgUO2(CO3)3-- = 3.0 H2O -6.0 H+ 3.0 CO2(aq) 1.0 Mg++    1.00 UO2++     23.9105   3.0   -2.0   500.0
NH4SO4-   =  1.5 H2O  1.0 H+    0.5 N2(aq)  1.0 SO4--  -0.75 O2(aq)    57.2905   4.0   -1.0   114.102

Isotherms

The “isotherms” section is the list of sorption processes that relate the sorbed concentration at the solid surface to the aqueous concentration in contact with the solid at constant temperature. It is a function of the free ion primary species concentrations, not total conentrations. A sorption isotherm may represent equilibrium or kinetic processes depending on the data used to fit the isotherm. Each sublist has two parameters:

isotherms-spec

• “model” [string] specifies the model name. The available options are “linear”, `”langmuir”, and “freundlich”.

• “parameters” [Array(double)] is the list of model parameters. The distribution coefficient K is in the first position.

<ParameterList name="isotherms">
  <ParameterList name="A">
    <Parameter name="model" type="string" value="linear"/>
    <Parameter name="parameters" type="Array(double)" value="{10.0}"/>
  </ParameterList>
  <ParameterList name="B">
    <Parameter name="model" type="string" value="langmuir"/>
    <Parameter name="parameters" type="Array(double)" value="{30.0, 0.1}"/>
  </ParameterList>
  <ParameterList name="C">
    <Parameter name="model" type="string" value="freundlich"/>
    <Parameter name="parameters" type="Array(double)" value="{1.5, 1.25}"/>
  </ParameterList>
</ParameterList>

A few examples are given below. Each line has three fields: primary species name, adsorption isotherm model, and parameters. The number of parameters and their meaning depends on the model; although the first one is always the distribution coefficient.

Pu_238   linear    461168.4
U_234    linear    329406.0
Th_230   linear   1482327.0
Ra_226   linear     41175.75
Pb_210   linear   3294060.0
Tc_99    linear       988.218

NOTE: The parameters provided here are global. The state field isotherm_kd overwrites any global data given here.

General kinetics

The “general kinetics” section describes kinetic (slow) reactions in the form

3 A(aq) + 2 B(aq) <-> C(aq) + 3 D(aq) + 4 E(aq)

where a number (stoichiometic coefficient) is followed by a species name. Implementation of the generalized kinetic formulation is based of publication “Multicomponent reactive transport modeling in variably saturated porous media using a generalized formulation for kinetically controlled reactions” by K. Ulrich Mayer (Water Res. Research, Vol.38:9, 1174, doi:10.1029/2001WR000862, 2002).

The list of parameters for each reaction includes

general_kinetics-spec

• “reactants” [string] is the left-hand side of the above equation.

• “products” [string] is the right-hand side of the above equation.

• “forward rate” [double] is the forward reaction rate.

• “backward rate” [double] is the reverse reaction rate.

• “reaction orders (reactants/products)” [Array(double)] is the list of reaction orders with respect to activies of dissolved species. The list includes orders for reactants (first) and products (second).

<ParameterList name="general kinetics">
  <ParameterList name="general_0">
    <Parameter name="reactants" type="string" value="1.0 Tritium"/>
    <Parameter name="products" type="string" value=""/>
    <Parameter name="forward rate" type="double" value="1.78577e-09"/>
    <Parameter name="backward rate" type="double" value="0.0"/>
    <Parameter name="reaction orders (reactants/products)" type="Array(double)" value="{3.0, 2.0, 0.0, 0.0, 0.0}"/>
  </ParameterList>
</ParameterList>

This example describes decay of Tritium.

Mineral thermodynamics and kinetics

The class for mineral reaction, should be written with the mineral as the reactant:

Calcite = 1.0 Ca++ + 1.0 HCO3- -1.0 H+

The “mineral kinetics” section includes at the moment both thermodynamic and kinetic data. The list of parameters for each reaction includes

mineral-spec

• “rate model” [string] is the model name for reaction rate [mol/s]. Available option is TST.

• “equilibrium constant” [double] is logarithm of the equilibrium constant.

• “rate constant” [double] is log10 of the reaction rate constant.

• “modifiers” [string] is the list of pairs: species name and value of modyfying exponent, so that the string has always even number of words.

• “gram molecular weight” [double] is amount of a molecular substance whose weight, in grams, is numerically equal to the molecular weight of that substance.

• “reaction” [string] is the mineral reaction equation.

• “molar volume” [double] is the molar volume, [m^3 / mol].

• “specific surface area” [double] is the specific mineral surface area.

The reaction rate is the dissolution rate for the mineral, so it is positive for dissolution and negative for precipitation. We assume that the sublist name is the mineral name. We also assume that the mineral reaction includes only primary species and that the modifying species are only primary species.

<ParameterList name="mineral kinetics">
  <ParameterList name="Kaolinite">
    <Parameter name="rate model" type="string" value="TST"/>
    <Parameter name="rate constant" type="double" value="-8.967"/>
    <Parameter name="modifiers" type="string" value="H+  7.77000E-01"/>
    <Parameter name="gram molecular weight" type="double" value="258.16"/>
    <Parameter name="reaction" type="string" value="5.0 H2O  -6.0 H+  2.0 Al+++  2.0 SiO2(aq)"/>
    <Parameter name="equilibrium constant" type="double" value="7.57"/>
    <Parameter name="molar volume" type="double" value="9.952e-05"/>
    <Parameter name="specific surface area" type="double" value="100.0"/>
  </ParameterList>
</ParameterList>

A few examples of mineral reactions is beflow. Each line in has five columns: mineral reaction, logarithm of equilibrium constant, gram molecular weight [g/mol], molar volume [cm^3/mol], and specific surface area [cm^2 mineral / cm^3 bulk].

Halite       = 1.0 Na+   1.0 Cl-                                         1.58550   58.4425   27.0150  1.0
Gypsum       = 2.0 H2O   1.0 SO4-2   1.0 Ca++                           -4.581    172.1722   74.21216 1.0
Calcite      =-1.0 H+    1.0 HCO3-   1.0 Ca++                            1.8487   100.087    36.934   1.0
Gibbsite     = 3.0 H2O   1.0 Al+++  -3.0 H+                              7.756     78.0036   31.956   1.0
Schoepite    = 3.0 H2O  -2.0 H+      1.0 UO2++                           4.8333   322.058    66.08    1.0
Basaluminite =15.0 H2O -10.0 H+      4.0 Al+++    1.0 SO4--             22.2511   464.140   218.934   1.0
Ferrihydrite = 2.5 H2O   1.0 Fe++   -2.0 H+       0.25 O2(aq)           -3.594    106.869    23.99    1.0
Jurbanite    = 6.0 H2O   1.0 Al+++  -1.0 H+       1.0 SO4--             -3.23     230.064   126.0     1.0
Kaolinite    = 5.0 H2O  -6.0 H+      2.0 Al+++    2.0 SiO2(aq)           7.570    258.160    99.520   1.0
Soddyite     = 4.0 H2O  -4.0 H+      1.0 SiO2(aq) 2.0 UO2++              0.392    668.169   131.27    1.0
K-Feldspar   = 2.0 H2O   1.0 K+      1.0 Al+++   -4.0 H+  3.0 SiO2(aq)  -0.2753   278.332   108.87    1.0
Polyhalite   = 2.0 H2O   1.0 Mg++    2.0 Ca++     2.0 K+  4.0 SO4-2    -13.7440   218.1     100.9722  1.0

Ion exchange sites

The “ion exchange sites” section is the list of ion exchange reactions. Each sublist is named after the exchanger site and has the following parameters:

ion_exchange-spec

• “location” [string] is the mineral name.

• “charge” [int] is the exchanger charge.

<ParameterList name="ion exchange sites">
  <ParameterList name="X-">
    <Parameter name="location" type="string" value="Bulk"/>
    <Parameter name="charge" type="int" value="-1"/>
  </ParameterList>
</ParameterList>

Ion exchange complexes

The class for ion exchange complexation reaction

NaX <===> Na+ + X-

The “ion exchange complexes” is the list of ion exchange complexation reactions. We assume that the reactions are always written as one complex equals one primary species plus one exchange site. Each complex is defined by the following parameters:

ion_exchange_complex-spec

• “reaction” [string] is the exchange reaction.

• “equilibrium constant” [double] is the logarithm of the equilibrium constant.

<ParameterList name="ion exchange complexes">
  <ParameterList name="Na+X">
    <Parameter name="reaction" type="string" value="1.0 Na+  1.0 X-"/>
    <Parameter name="equilibrium constant" type="double" value="1.0"/>
  </ParameterList>
  <ParameterList name="Ca++X">
    <Parameter name="reaction" type="string" value="1.0 Ca++  2.0 X-"/>
    <Parameter name="equilibrium constant" type="double" value="0.316228"/>
  </ParameterList>
</ParameterList>

A few additional examples, reaction equation and the equilibrium coefficient:

Al+++X = 1.0 Al+++ 3.0 X-    1.71133
Ca0.5X = 0.5 Ca++  1.0 X-   -0.99
H+X    = 1.0 H+    1.0 X-    0.0251189
Mg++X  = 1.0 Mg++  2.0 X-    0.1666

Surface complex sites

The “surface complex sites” section is the list of surface sites and related data. Each site has one parameter “density”.

<ParameterList name="surface complex sites">
  <ParameterList name=">davis_OH">
    <Parameter name="density" type="double" value="0.156199"/>
  </ParameterList>
</ParameterList>

NOTE: Each instance of this class should contain a single unique surface site (e.g. >FeOH) and ALL surface complexes associated with that site!

Surface complexes

The “surface complexes” is the list of surface complexation reactions. Each reaction is defined by the following parameters:

surface_complex-spec

• “reaction” [string] is a surface complexation reaction involing the complex site and primary species.

• “charge” [int] is the charge of the complex.

• “equilibrium constant” [double] is the logarithm of the reaction equilibrium coefficeint.

<ParameterList name="surface complexes">
  <ParameterList name="(>davis_O)UO2+">
    <Parameter name="charge" type="int" value="1"/>
    <Parameter name="reaction" type="string" value="1.0 >davis_OH  -1.0 H+  1.0 UO2++"/>
    <Parameter name="equilibrium constant" type="double" value="-0.444"/>
  </ParameterList>
</ParameterList>

A few additional examples. Each line has three fields: reaction, logarithm of equailibrium coefficient, and charge.

>FeOH2+_s     = 1.0 >FeOH_s   1.0 H+                        -7.18   1.0
>FeO-_w       = 1.0 >FeOH_w  -1.0 H+                         8.82  -1.0
>FeOHUO2++    = 1.0 >FeOH     1.0 UO2++                     -6.63   2.0
>SiO-         =-1.0 H+        1.0 >SiOH                      0.0
>SiOH2+       = 1.0 H+        1.0 >SiOH                      0.0
>AlOUO2+      = 1.0 >AlOH    -1.0 H+   1.0 UO2++            -3.13   1.0
>FeOHZn+_s    = 1.0 >FeOH_s  -1.0 H+   1.0 Zn++             -0.66   1.0
>SiOUO3H3++   = 1.0 >SiOH     1.0 H2O  1.0 UO2++             5.18   2.0
>UO2++        = 1.0 UO2++    -1.0 Ca++ 1.0 >Ca++            -5.12   0.0
>SiOUO3H2+    = 1.0 >SiOH     1.0 H2O  -1.0 H+  1.0 UO2++    5.18   1.0
>SiOUO3H      = 1.0 >SiOH     1.0 H2O  -2.0 H+  1.0 UO2++    5.18   0.0
>SiOUO3-      = 1.0 >SiOH     1.0 H2O  -3.0 H+  1.0 UO2++   12.35  -1.0
>SiOUO2(OH)2- = 1.0 >SiOH     2.0 H2O  -3.0 H+  1.0 UO2++   12.35  -1.0
>FeOHUO3      = 1.0 >FeOH     1.0 H2O  -2.0 H+  1.0 UO2++    3.05   0.0

Radioactive decay

The “radiaoctive decay” section is the list of decay reactions for aqueous and sorbed species. It not deal with decay of solid phase. Currently, it does not deal with decay of solid phase.

radioactive_decay-spec

• “reactant” [string] is the species name.

• “product” [string] is the species name.

• “half life” [double] is the half-life time of decay.

<ParameterList name="radioactive decay">
  <ParameterList name="complex_0">
    <Parameter name="reactant" type="string" value="A"/>
    <Parameter name="product" type="string" value="B"/>
    <Parameter name="half life" type="double" value="3.9447e+08"/>
  </ParameterList>
  <ParameterList name="complex_1">
    <Parameter name="reactant" type="string" value="B"/>
    <Parameter name="product" type="string" value="C"/>
    <Parameter name="half life" type="double" value="7.8894e+08"/>
  </ParameterList>
  <ParameterList name="complex_2">
    <Parameter name="reactant" type="string" value="C"/>
    <Parameter name="product" type="string" value=""/>
    <Parameter name="half life" type="double" value="1.57788e+08"/>
  </ParameterList>
</ParameterList>

A few additionale examples:

Cs137  -->  1.0 Cs137    half_life 30.2 years
Pb_210 -->               half_life 22.2 years
Pu_238 -->  1.0 U_234    half_life 87.7 years
Ra_226 -->  1.0 Pb_210   half_life 1.6e3 years
Th_230 -->  1.0 Ra_226   half_life 7.54e4 years
U_234  -->  1.0 Th_230   half_life 2.45e5 years
Tc_99  -->               half_life 2.111e5 years
Sr90   -->  1.0 Sr90     half_life 28.8 years

OUTPUT DATA

Amanzi uses a few ways to communicate simulation data to the user that includes a short file with observations and full-scale visualization files.

Observation file

A user may request any number of specific observations from Amanzi. Each labeled observation data quantity involves a field quantity, a model, a region from which it will extract its source data, and a list of discrete times for its evaluation. The observations are evaluated during the simulation and returned to the calling process through one of Amanzi arguments.

observations-spec

• “observation data” [list] can accept multiple lists for named observations.

  • “observation output filename” [string] user-defined name for the file that the observations are written to. The file name can contain relative or absolute path to an existing directory only.

  • “time unit” [string] defines time unit for output data. Available options are “s”, “h”, “d”, and “y”. Default is “s”.

  • “mass unit” [string] defines mass unit for output data. Available options are “g”, “lb”, and “lb”. Default is “kg”.

  • “length unit” [string] defines length unit for output data.

    Available options are “cm”, “in”, “ft”, “yd” , “m”, and “km”. Default is “m”.

  • “concentration unit” [string] defines concentration unit for output data.

    Available options are “molar”, and “SI”. Default is “molar”.

  • “precision” [int] defines the number of significant digits. Default is 16.

  • OBSERVATION [list] user-defined label, can accept values for “variables”, “functional”, “region”, “times”, and TSPS (see below).

    • “domain name” [string] name of the domain. Typically, it is either “domain” for the matrix/subsurface or “fracture” for the fracture network.

    • “variable” [string] a list of field quantities taken from the list of available field quantities.

    • “functional” [string] the label of a function to apply to each of the variables in the variable list (Function options detailed below)

    • “region” [string] the label of a user-defined region

    • “cycles start period stop” [Array(int)] the first entry is the start cycle, the second is the cycle period, and the third is the stop cycle or -1 in which case there is no stop cycle. A visualization dump shall be written at such cycles that satisfy cycle = start + n*period, for n=0,1,2,… and cycle < stop if stop != -1.0.

    • “cycles start period stop n” [Array(int)] if multiple cycles start-period-stop parameters are needed, then use these parameters with n=0,1,2,…, and not the single “cycles start period stop” parameter.

    • “cycles” [Array(int)] an array of discrete cycles that at which a visualization dump shall be written.

    • “times start period stop” [Array(double)] the first entry is the start time, the second is the time period, and the third is the stop time or -1 in which case there is no stop time. A visualization dump shall be written at such times that satisfy time = start + n*period, for n=0,1,2,… and time < stop if stop != -1.0.

    • “times start period stop n” [Array(double) if multiple start-period-stop parameters are needed, then use this these parameters with n=0,1,2,…, and not the single “times start period stop” parameter.

    • “times” [Array(double)] an array of discrete times that at which a visualization dump shall be written.

    • “delimiter” [string] the string used to delimit columns in the observation file output, default is “,”.

    • “interpolation” [string] the string which defines the interpolation method to compute observation. Works ONLY with Line Segment region at the moment. Available options “linear” and “constant”. Default is “linear”

    • “weighting” [string] the string defined the weighting function applied to compute observation. Works ONLY with Line Segment region at the moment. Available options “flux norm” and “none”. Default is “none”. “flux norm” is the absolute value of the Darcy flux in a cell.

Available observation variables:

• volumetric water content [-] (volume water / bulk volume)

• aqueous saturation [-] (volume water / volume pore space)

• aqueous pressure [Pa]

• hydraulic head [m]

• permeability-weighted hydraulic head [m]

• drawdown [m]

• permeability-weighted drawdown [m]

• volumetric water content [-]

• gravimetric water content [-]

• water table [m]

• SOLUTE aqueous concentration [mol/m^3]

• SOLUTE gaseous concentration [mol/m^3]

• SOLUTE sorbed concentration [mol/kg]

• SOLUTE free ion concentration

• x-, y-, z- aqueous volumetric flux [m/s]

• material id [-]

• aqueous mass flow rate [kg/s] (when funtional=”integral”)

• aqueous volumetric flow rate [m^3/s] (when functional=”integral”)

• fractures aqueous volumetric flow rate [m^3/s] (when functional=”integral”)

• SOLUTE volumetric flow rate [mol/s] (when functional=”integral”)

• SOLUTE breakthrough curve [mol] (when functional=”integral”)

• pH [-]

• centroid x [m]

Observations drawdown and permeability-weighted are calculated with respect to the value registered at the first time it was requested.

The following observations are point-type obervations: “water table”, “drawdown”.

The following observations are integrated continuously in time but saved only at specified times: “SOLUTE breakthrough curve”.

The following observation functionals are currently supported. All of them operate on the variables identified.

observation_data-spec

• “observation data: point” returns the value of the field quantity at a point.

• “observation data: integral” returns the integral of the field quantity over the region specified.

• “observation data: extensive integral” returns the integral of an extensive variable over the region specified. Note that this should be used over the above Integral when the variable to be integrated is an extensive quantity, i.e. water content or flux.

• “observation data: minimum” and “observation data: maximum” returns the minimum (respectively maximum) of the field quantity over the region specified.

<ParameterList>  <!-- parent list -->
<ParameterList name="observation data">
  <Parameter name="observation output filename" type="string" value="_OUTPUT.out"/>
  <Parameter name="precision" type="int" value="10"/>
  <ParameterList name="_ANY OBSERVATION NAME">
    <Parameter name="region" type="string" value="_REGION"/>
    <Parameter name="functional" type="string" value="observation data: point"/>
    <Parameter name="variable" type="string" value="volumetric water content"/>

    <Parameter name="cycles" type="Array(int)" value="{100000, 200000, 400000, 500000}"/>
    <Parameter name="cycles start period stop" type="Array(int)" value="{0, 100, -1}"/>

    <Parameter name="times" type="Array(double)" value="{101.0, 303.0, 422.0}"/>
    <Parameter name="times start period stop 0" type="Array(double)" value="{0.0, 10.0, 100.0}"/>
    <Parameter name="times start period stop 1" type="Array(double)" value="{100.0, 25.0, -1.0}"/>
  </ParameterList>

  <ParameterList name="_ANY OBSERVATION NAME B">  <!-- another observation -->
    ...
  </ParameterList>
</ParameterList>
</ParameterList>

In this example, we collect “volumetric water content” on four selected cycles, every 100 cycles, three selected times, every 10 seconds from 0 to 100, and every 25 seconds after that.

Checkpoint file

A user may request periodic dumps of Amanzi checkpoint data. The user has no explicit control over the content of these files, but has the guarantee that the Amanzi run will be reproducible (with accuracy determined by machine round errors and randomness due to execution in a parallel computing environment). Therefore, output controls for checkpoint data are limited to file name generation and writing frequency, by numerical cycle number.

checkpoint_file-spec

• “checkpoint data” [list] can accept a file name base [string] and cycle data [list] used to generate the file base name or directory base name that is used in writing checkpoint data.

  • “file name base” [string] (“checkpoint”)

  • “file name digits” [int] Default value is 5.

  • “cycles start period stop” [Array(int)] the first entry is the start cycle, the second is the cycle period, and the third is the stop cycle or -1 in which case there is no stop cycle. A visualization dump shall be written at such cycles that satisfy cycle = start + n*period, for n=0,1,2,… and cycle < stop if stop != -1.0.

  • “cycles start period stop n” [Array(int)] if multiple cycles start period stop parameters are needed, then use these parameters with n=0,1,2,…, and not the single “cycles start period stop” parameter.

  • “cycles” [Array(int)] an array of discrete cycles that at which a visualization dump shall be written.

  • “times start period stop” [Array(double)] the first entry is the start time, the second is the time period, and the third is the stop time or -1 in which case there is no stop time. A visualization dump shall be written at such times that satisfy time = start + n*period, for n=0,1,2,… and time < stop if stop != -1.0.

  • “times start period stop n” [Array(double)] if multiple start period stop parameters are needed, then use this these parameters with n=0,1,2,…, and not the single “times start period stop” parameter.

  • “times” [Array(double)] an array of discrete times that at which a visualization dump shall be written.

<ParameterList>  <!-- parent list -->
<ParameterList name="checkpoint data">
  <Parameter name="file name base" type="string" value="_CHKPOINT"/>
  <Parameter name="file name digits" type="int" value="5"/>

  <Parameter name="cycles start period stop" type="Array(int)" value="{0, 100, -1}"/>
  <Parameter name="cycles" type="Array(int)" value="{999, 1001}"/>

  <Parameter name="times start period stop 0" type="Array(double)" value="{0.0, 10.0, 100.0}"/>
  <Parameter name="times start period stop 1" type="Array(double)" value="{100.0, 25.0, -1.0}"/>
  <Parameter name="times" type="Array(double)" value="{101.0, 303.0, 422.0}"/>
</ParameterList>
</ParameterList>

In this example, checkpoint data files are written when the cycle number is a multiple of 100.

Visualization file

A user may request periodic writes of field data for the purposes of visualization. The user will specify explicitly what is to be included in the file at each snapshot. Visualization files can only be written at intervals corresponding to the numerical time step values or intervals corresponding to the cycle number; writes are controlled by time step cycle number.

vis-spec

• “visualization data” [list] can accept a file name base [string] and cycle data [list] that is used to generate the file base name or directory base name that is used in writing visualization data. It can also accept a set of lists to specify which field quantities to write. The file name can contain relative or absolute path to an existing directory only.

  • “file name base” [string] (“amanzi_vis”)

  • “file format” [string] (“XDMF”) Amanzi supports two types of visualization files. XDMF is the default and preferred method, but does not correctly handle general polyhedra. Serial, 3D general polyhedral support is supported by the “SILO” option. This will eventually be extended to parallel, 2/3D support, but this is not yet implemented.

  • “cycles start period stop” [Array(int)] the first entry is the start cycle, the second is the cycle period, and the third is the stop cycle or -1 in which case there is no stop cycle. A visualization dump shall be written at such cycles that satisfy cycle = start + n*period, for n=0,1,2,… and cycle < stop if stop != -1.0.

  • “cycles start period stop n” [Array(int)] if multiple cycles start-period-stop parameters are needed, then use these parameters with n=0,1,2,…, and not the single “cycles start period stop” parameter.

  • “cycles” [Array(int)] an array of discrete cycles that at which a visualization dump shall be written.

  • “times start period stop” [Array(double)] the first entry is the start time, the second is the time period, and the third is the stop time or -1 in which case there is no stop time. A visualization dump shall be written at such times that satisfy time = start + n*period, for n=0,1,2,… and time < stop if stop != -1.0.

  • “times start period stop n” [Array(double) if multiple start-period-stop parameters are needed, then use this these parameters with n=0,1,2,…, and not the single “times start period stop” parameter.

  • “times” [Array(double)] an array of discrete times that at which a visualization dump shall be written.

  • “dynamic mesh” [bool] (false) write mesh data for every visualization dump, this facilitates visualizing deforming meshes.

  • “write regions” [list] contains auxiliary fields with region ids to write into the visualization file.

    • “REGION_NAME” [Array(string)] the user-defined field name and the list of assigned regions. The first entry in the regions array is marked with the value 1.0 in the array, the second with the value 2.0, and so forth. The code ignores entries in the regions array that are not valid regions that contain cells.

  • “write partitions” [bool] (false) if this parameter is true, then write an array into the visualization file that contains the rank number of the processor that owns a mesh cell.

  • “time units” [string] time units, e.g. “y” or “h”.

Example:

<ParameterList>  <!-- parent list -->
<ParameterList name="visualization data">
  <Parameter name="file name base" type="string" value="_PLOT"/>

  <Parameter name="cycles start period stop" type="Array(int)" value="{0, 100, -1}"/>
  <Parameter name="cycles" type="Array(int)" value="{999, 1001}"/>

  <Parameter name="times start period stop 0" type="Array(double)" value="{0.0, 10.0, 100.0}"/>
  <Parameter name="times start period stop 1" type="Array(double)" value="{100.0, 25.0, -1.0}"/>
  <Parameter name="times" type="Array(double)" value="{101.0, 303.0, 422.0}"/>

  <Parameter name="dynamic mesh" type="bool" value="false"/>

  <ParameterList name="write regions">
    <Parameter name="regions" type="Array(string)" value="{_OBS1, _OBS2, _OBS3}"/>
    <Parameter name="wells" type="Array(string)" value="{_OBS4}"/>
  </ParameterList>
</ParameterList>
</ParameterList>

Walkabout file

A user may request periodic dumps of Walkabout data. Output controls for Walkabout data are limited to file name generation and writing frequency, by numerical cycle number or time.

walkabout-spec

• “walkabout data” [list] can accept a file name base [string] and cycle data [list] used to generate the file base name or directory base name that is used in writing Walkabout data.

  • “file name base” [string] The file name can contain relative or absolute path to an existing directory only. Default is “walkabout”.

  • “file name digits” [int] specify the number of digits that should be appended to the file name for the cycle number. Default is 5.

  • “cycles start period stop” [Array(int)] the first entry is the start cycle, the second is the cycle period, and the third is the stop cycle or -1 in which case there is no stop cycle. A visualization dump shall be written at such cycles that satisfy cycle = start + n*period, for n=0,1,2,… and cycle < stop if stop != -1.0.

  • “cycles start period stop n” [Array(int)] if multiple cycles start-period-stop parameters are needed, then use these parameters with n=0,1,2,…, and not the single “cycles start period stop” parameter.

  • “cycles” [Array(int)] an array of discrete cycles that at which a visualization dump shall be written.

  • “times start period stop” [Array(double)] the first entry is the start time, the second is the time period, and the third is the stop time or -1 in which case there is no stop time. A visualization dump shall be written at such times that satisfy time = start + n*period, for n=0,1,2,… and time < stop if stop != -1.0.

  • “times start period stop n” [Array(double)] if multiple start-period-stop parameters are needed, then use this these parameters with n=0,1,2,…, and not the single “times start period stop” parameter.

  • “times” [Array(double)] an array of discrete times that at which a visualization dump shall be written.

  • “write regions” [list] contains three lists of equal size with region names, material names, and material ids to write into the output file.

    • “region names” [Array(string)] specifies names of regions.

    • “material names” [Array(int)] specifies names of materials.

    • “material ids” [Array(int)] specifies material ids.

Example:

<ParameterList>  <!-- parent list -->
<ParameterList name="walkabout data">
  <Parameter name="file name base" type="string" value="_WALKABOUT"/>
  <Parameter name="file name digits" type="int" value="5"/>

  <Parameter name="cycles start period stop" type="Array(int)" value="{0, 100, -1}"/>
  <Parameter name="cycles" type="Array(int)" value="{999, 1001}"/>

  <Parameter name="times start period stop 0" type="Array(double)" value="{0.0, 10.0, 100.0}"/>
  <Parameter name="times start period stop 1" type="Array(double)" value="{100.0, 25.0, -1.0}"/>
  <Parameter name="times" type="Array(double)" value="{101.0, 303.0, 422.0}"/>

  <ParameterList name="write regions">
    <Parameter name="region names" type="Array(string)" value="{_REGION1, _REGION2}"/>
    <Parameter name="material names" type="Array(string)" value="{_MAT1, _MAT2}"/>
    <Parameter name="material ids" type="Array(int)" value="{1000, 2000}"/>
  </ParameterList>
</ParameterList>
</ParameterList>

In this example, walkabout data files are written when the cycle number is a multiple of 100.

Mesh info

A user may request to dump mesh information. Mesh information includes coordinates of cell centroids written is the order consistent with all output fields.

mesh_info-spec

• “filename” [string] is name of the HDF5 file where coordinates of the centroids are dumped.

<ParameterList>  <!-- parent list -->
<ParameterList name="mesh info">
  <Parameter name="filename" type="string" value="centroids"/>
</ParameterList>

<ParameterList name="mesh info fracture">
  <Parameter name="filename" type="string" value="centroids_fracture"/>
</ParameterList>
</ParameterList>

INPUT DATA

This section describes format and purpose of various input files. In addition it explain how to verify the input information (e.g. regions) using special sublists.

Input analysis

This list contains data collected by the input parser of a higher-level spec.

input_analysis-spec

• “used boundary condition regions” [Array(string)] provides list of boundary regions for analysis. The simulator will print number of faces and total area of these regions if verbosity level is equal to or above high.

• “used source and sink regions” [Array(string)] provides list of source and sink regions for analysis. The simulator will print number of cells and the total volume of these regions if verbosity level is equal to or above high.

• “used observation regions” [Array(string)] provides list of observation regions for analysis. The simulator will print number of faces(or cells) and the total area (or volume) of these regions if verbosity level is equal to or above high.

<ParameterList>  <!-- parent list -->
<ParameterList name="analysis">
  <Parameter name="used boundary condition regions" type="Array(string)" value="{_REG1,_REG2}"/>
  <Parameter name="used source and sink regions" type="Array(string)" value="{_REG3,_REG4}"/>
  <Parameter name="used observation regions" type="Array(string)" value="{_REG5}"/>
  <ParameterList name="verbose object">
    <Parameter name="verbosity level" type="string" value="high"/>
  </ParameterList>
</ParameterList>
</ParameterList>

Tabulated EOS format (Amanzi)

The tabulated equation of state uses a free format for integers, doubles and text string; however, the order of block is fixed and string cannot have spaces. The first three lines have the following format:

Amanzi
N1 scale shift temperature
N2 scale shift pressure

where Ni is the number of data points in either T or p directions, scale and shift are the scale factor and shift required to convert data to SI units, and field_name is either “temperature” or “pressure”

The second block has seven lines and the above format, except for N which is now the total number of data points in the T-p table.

The third block provides data for physical fields using the following order:

T p rho internal_energy enthalpy Cv Cp viscosity thermal_conductivity phase

where phase is the first letter of phase, “l” or “g”.

Example for liquid water

In this example the following units are used for 3x3 tabulated data: temperature [K], pressure [bar], density [kg/m3], internal energy [kJ/mol], enthalpy [kJ/mol], Cv [J/mol*K], Cp [J/mol*K], viscosity [Pa*s], and thermal conductivity [W/m*K]. Thus, to convert these data to SI units, we specify non-unity factors for pressure, internal energy, and enthalpy.

Amanzi
3 1.0000 0.0 temperature
3 1.0e+5 0.0 pressure

9 1.0000 0.0 density
9 1.0e+3 0.0 internal_energy
9 1000.0 0.0 enthalpy
9 1.0 0.0 cv
9 1.0 0.0 cp
9 1.0 0.0 viscosity
9 1.0 0.0 thermal_conductivity

280.15 0.5  999.879 0.530097 0.530997 75.6574 75.6793 0.00142708 0.574347  l
290.15 0.5  998.754  1.28544  1.28634 75.0821 75.4233 0.00107982 0.593025  l
300.15 0.5  996.493  2.03905  2.03995 74.3983 75.3170 0.000851   0.610549  l

280.15 1.0  999.904  0.530085 0.531887 75.6536 75.6756 0.00142701 0.574372  l
290.15 1.0  998.777  1.28540  1.28720  75.0787 75.4203 0.00107979  0.593048  l
300.15 1.0  996.515  2.03898  2.04078  74.3954 75.3145 0.000850991  0.610572  l

280.15 1.2 999.913  0.53008 0.532242 75.6521  75.6741  0.00142699  0.574382  l
290.15 1.2 998.787  1.28538 1.28754 75.0774  75.4191  0.00107978 0.593058  l
300.15 1.2 996.524  2.03895 2.04112 74.3942  75.3135  0.000850987  0.610581 l

Tabulated EOS format (FEHM)

The tabulated equation of state uses a free format for integers, doubles and text string; however, the order of block is fixed and string cannot have spaces.

Example for liquid water

nonuniform                                                Grid type
   8    6     9                                          Array dimensions
   1.000000        0.000000                               Temperature factor and offset
   1.000000        0.000000                               Pressure factor and offset
   19                                                     Saturation line closeness flag
Table temperatures --->
 -200.00   -195.00   -190.00   -185.00   -180.00   -175.00   -170.00   -165.00
Table pressures --->
 6.5074101E-02  0.1195323   0.2031954   0.3243768   0.4918777   0.7149075
Array property types --->
density
dd/dt
dd/dp
enthalpy
dh/dt
dh/dp
viscosity
dv/dt
dv/dp
Density array --->
  900.9539  2.986970  2.792525  2.623569  2.474998  2.343097  2.225062  2.118721
  901.0682  878.6923  5.241142  4.904944  4.613682  4.357906  4.130922  3.927764
  901.2436  878.8931  855.6431  8.584902  8.033060  7.558094  7.142697  6.775020
  901.4971  879.1833  855.9796  831.6210  13.32351  12.44936  11.70498  11.05858
  901.8464  879.5831  856.4427  832.1661  806.3937  19.82610  18.47106  17.33520
  902.3099  880.1130  857.0560  832.8870  807.2585  779.6672  28.56710  26.48487
dD/dT array --->
 -4.401056  -4.465137  -3.634006E-02 -3.175270E-02 -2.8047204E-02 -2.4993611E-02 -2.2437597E-02 -2.0270992E-02
 -4.404730  -4.554480  -4.633789     -6.274590E-02 -5.4703809E-02 -4.8275996E-02 -4.3014195E-02 -3.8629483E-02
 -4.400110  -4.560047  -4.748645     -4.847290     -0.1026809     -8.9036323E-02 -7.8307390E-02 -6.9622092E-02
 -4.393451  -4.551746  -4.756232     -4.995862     -5.119995      -0.1618524     -0.1390779     -0.1215700
 -4.384296  -4.540369  -4.741705     -5.004901     -5.313653      -5.472827      -0.2490898     -0.2111700
 -4.372241  -4.525384  -4.722595     -4.979755     -5.321985      -5.728076      -5.937891      -0.3788673
dD/dP array --->
  1.676756  46.57508  44.38469  41.44989  38.92557  36.72093  34.77290  33.03539
  2.097235  1.917826  44.96321  43.16014  40.24044  37.75665  35.60374  33.71167
  2.093754  2.396778  2.227240  43.98545  42.51919  39.50047  36.97467  34.81088
  2.088050  2.390179  2.769690  2.630571  43.65722  42.49656  39.24163  36.58062
  2.081284  2.380574  2.756290  3.241876  3.168733  44.03999  43.17746  39.50084
  2.072423  2.367801  2.738755  3.216566  3.853050  3.902956  45.26769  44.73188
Enthalpy array --->
 -11.08680  202.0743    207.3656  212.5940  217.7786  222.9311  228.0593  233.1689
 -11.04795  -1.441901   205.8394  211.2797  216.6258  221.9066  227.1395  232.3360
 -10.98825  -1.385264   8.315416  209.1396  214.7777  220.2795  225.6877  231.0275
 -10.90171  -1.303152   8.391649  18.23019  211.8896  217.7864  223.4902  229.0631
 -10.78201  -1.189503   8.497250  18.32486  28.35650  213.9910  220.2263  226.1903
 -10.62244  -1.037911   8.638263  18.45154  28.46319  38.76015  215.3184  222.0043
dH/dT array --->
  1.917285  1.925377  1.051970  1.041299  1.033711  1.028070  1.023779  1.020430
  1.915049  1.933634  1.946059  1.078633  1.062691  1.051370  1.042940  1.036470
  1.914508  1.930367  1.957850  1.975565  1.113985  1.091000  1.074800  1.062750
  1.913727  1.929336  1.953334  1.991954  2.016201  1.160051  1.127669  1.105069
  1.912657  1.927926  1.951436  1.985925  2.038849  2.071369  1.219936  1.174710
  1.911247  1.926070  1.948945  1.982493  2.030861  2.102872  2.146351  1.298700
dH/dP array --->
 0.5700204  -31.36753  -27.58228  -23.85831  -20.98058  -18.67626  -16.78886  -15.21603
 0.7135022  0.5404872  -28.02403  -25.00951  -21.72656  -19.19770  -17.17040  -15.50379
 0.7139086  0.6773382  0.5040666  -25.57932  -23.12046  -20.11376  -17.81494  -15.97750
 0.7144175  0.6781194  0.6298733  0.4559285  -23.83204  -21.78361  -18.91836  -16.75610
 0.7151018  0.6791806  0.6314830  0.5667921  0.3892052  -22.65932  -20.92485  -18.07487
 0.7160630  0.6805639  0.6335793  0.5700380  0.4816387  0.2933709  -22.00557  -20.52189
Viscosity array --->
 2.0667833E-04  9.9999997E-06  5.9198064E-06  6.2753802E-06  6.6270109E-06  6.9747775E-06  7.3187521E-06  7.6590040E-06
 2.0685646E-04  1.7128805E-04  9.9999997E-06  6.2796817E-06  6.6325765E-06  6.9813382E-06  7.3261003E-06  7.6669758E-06
 2.0713007E-04  1.7152553E-04  1.4429199E-04  6.2892664E-06  6.6436955E-06  6.9936523E-06  7.3393444E-06  7.6809429E-06
 2.0752632E-04  1.7186940E-04  1.4459844E-04  1.2321399E-04  9.9999997E-06  7.0168890E-06  7.3631581E-06  7.7051855E-06
 2.0807389E-04  1.7234446E-04  1.4502162E-04  1.2360205E-04  1.0634069E-04  7.0617830E-06  7.4066052E-06  7.7475661E-06
 2.0880280E-04  1.7297656E-04  1.4558440E-04  1.2411769E-04  1.0682854E-04  9.2474213E-05  9.9999997E-06  7.8235007E-06
dV/dT array --->
-8.2993447E-06 -7.0931769E-06  7.1518159E-08  7.0720446E-08  6.9939738E-08  6.9174121E-08  6.8422651E-08  6.7684958E-08
-8.3128207E-06 -6.2673280E-06 -5.4209727E-06  7.1008344E-08  7.0165655E-08  6.9352382E-08  6.8563757E-08  6.7796513E-08
-8.3213863E-06 -6.2838080E-06 -4.8463658E-06 -4.2460219E-06  7.0438588E-08  6.9564884E-08  6.8729058E-08  6.7924752E-08
-8.3338045E-06 -6.2927879E-06 -4.8655406E-06 -3.8467074E-06 -3.4165248E-06  6.9696760E-08  6.8829650E-08  6.7998009E-08
-8.3509904E-06 -6.3052271E-06 -4.8742413E-06 -3.8680932E-06 -3.1414711E-06 -2.8306706E-06  6.8578309E-08  6.7805061E-08
-8.3739160E-06 -6.3218395E-06 -4.8858869E-06 -3.8755861E-06 -3.1643483E-06 -2.6455236E-06 -2.4201813E-06  6.6710982E-08
dV/dP array --->
 2.6134942E-06 -1.4940065E-08  4.0690097E-08  7.1375538E-08  9.5491849E-08  1.1454858E-07  1.2968569E-07  1.4171002E-07
 3.2705682E-06  2.2667357E-06  4.9334119E-08  1.0053603E-07  1.2079731E-07  1.3665340E-07  1.4908871E-07  1.5883744E-07
 3.2700600E-06  2.8380018E-06  2.0271643E-06  1.1456318E-07  1.6409528E-07  1.7355011E-07  1.8090707E-07  1.8653024E-07
 3.2694163E-06  2.8367640E-06  2.5274601E-06  1.8711472E-06  1.8563070E-07  2.3600580E-07  2.3299235E-07  2.3078385E-07
 3.2685762E-06  2.8350339E-06  2.5246798E-06  2.3140383E-06  1.7851987E-06  2.6802255E-07  3.2196141E-07  3.0296016E-07
 3.2674695E-06  2.8328298E-06  2.5211041E-06  2.3085806E-06  2.1822507E-06  1.7638789E-06  3.6895875E-07  4.3009413E-07
Saturation line intersection array --->
  F   T   F   F   F   F   F   F
  F   F   T   F   F   F   F   F
  F   F   F   T   F   F   F   F
  F   F   F   F   T   F   F   F
  F   F   F   F   F   T   F   F
  F   F   F   F   F   F   T   F
Saturation line closeness array --->
  0  10   0   0   0   0   0   0
  0   0  11   0   0   0   0   0
  0   0   0  12   0   0   0   0
  0   0   0   0  13   0   0   0
  0   0   0   0   0  14   0   0
  0   0   0   0   0   0  15   0

    6               Number of saturation line vertices
Pressure            Temperature         dP/dT               Cell      Intersection_type
 6.5074101E-02       -195.0000           1.0891641E-02         487        6
 0.1195323           -190.0000           1.6732618E-02         541        6
 0.2031954           -185.0000           2.4236280E-02         595        6
 0.3243768           -180.0000           3.3500183E-02         649        6
 0.4918777           -175.0000           4.4605963E-02         703        6
 0.7149075           -170.0000           5.7636499E-02         757        6


Liquid properties at saturation line vertices --->
Density     dD/dT      dD/dP      Enthalpy   dH/dT      dH/dP      Viscosity       dV/dT          dV/dP
  878.6282  -4.465137   1.176809  -1.459917   1.925377  0.3308200  1.7121245E-04  -7.0931769E-06  1.3881693E-06
  855.5234  -4.633789   1.431348   8.288394   1.946059  0.3229902  1.4418318E-04  -5.4209727E-06  1.3005098E-06
  831.4067  -4.847290   1.768377   18.19324   1.975565  0.3049078  1.2306188E-04  -4.2460219E-06  1.2552354E-06
  806.0210  -5.119995   2.224947   28.31119   2.016201  0.2704890  1.0613137E-04  -3.4165248E-06  1.2496769E-06
  779.0295  -5.472827   2.858964   38.71334   2.071369  0.2098658  9.2187336E-05  -2.8306706E-06  1.2862695E-06
  749.9777  -5.937891   3.772253   49.49191   2.146351  0.1066512  8.0373306E-05  -2.4201813E-06  1.3734772E-06

Vapour properties at saturation line vertices --->
Density     dD/dT           dD/dP     Enthalpy  dH/dT     dH/dP      Viscosity      dV/dT           dV/dP
  2.986970  -3.8888931E-02  46.57508  202.0743  1.058261  -31.36753  5.5601986E-06  7.1921569E-08  -1.4940065E-08
  5.241141  -6.7239381E-02  44.96321  205.8395  1.088046  -28.02403  5.9224931E-06  7.1437718E-08   4.9334119E-08
  8.584903  -0.1103685      43.98545  209.1396  1.127609  -25.57932  6.2892664E-06  7.0885832E-08   1.1456318E-07
  13.32350  -0.1748287      43.65722  211.8897  1.179343  -23.83204  6.6661905E-06  7.0139684E-08   1.8563070E-07
  19.82610  -0.2710075      44.03999  213.9910  1.247070  -22.65932  7.0617830E-06  6.8964439E-08   2.6802255E-07
  28.56710  -0.4164467      45.26769  215.3184  1.337180  -22.00557  7.4888940E-06  6.6921345E-08   3.6895875E-07

Tabulated function file format

The following ASCII input file format supports the definition of a tabulated function defined over a grid. Several XML input Parameters refer to files in this format. The file consists of the following records (lines). Each record is on a single line, except for the DATAVAL record which may be split across multiple lines.

1. DATATYPE: An integer value: 0 for integer data, 1 for real data.

• An integer-valued file is used to define a ‘color’ function used in the definition of a region.

2. GRIDTYPE: A string that specifies the type of grid used to define the function. The format of the rest of the file is contingent upon this value. The currently supported options are uniform rectilinear grids in 1, 2 and 3-D, which are indicated by the values 1DCoRectMesh, 2DCoRectMesh and 3DCoRectMesh, respectively (names adopted from XDMF).

For the uniform rectilinear grids, the remaining records are as follows. Several records take 1, 2 or 3 values depending on the space dimension of the grid.

3. NXNYNZ: 3 (or 2, 1) integer values (NX, NY, NZ) giving the number of zones in the x, y and z coordinate directions, respectively.

4. CORNER1: 3 (or 2, 1) floating point values (X1, Y1, Z1) giving the coordinate of the first corner of the domain.

5. CORNER2: 3 (or 2, 1) floating point values (X2, Y2, Z2) giving the coordinate of the second corner of the domain. The grid points r_{i,j,k} = (x_i, y_j, z_j) are defined as:

x_i = X1 + i*(X2-X1)/NX, 0 <= i <= NX

y_j = Y1 + j*(Y2-Y1)/NY, 0 <= j <= NY

z_k = Z1 + k*(Z2-Z1)/NZ, 0 <= k <= NZ

The (i,j,k) grid cell is defined by the corner grid points r_{i-1,j-1,k-1} and r_{i,j,k}, for 1 <= i <= NX, 1 <= j <= NY, 1 <= k <= NZ. Note that the corner points are any pair of opposite corner points; the ordering of grid points and cells starts at CORNER1 and ends at CORNER2.

6. DATALOC: An integer value: 0 for cell-based data, 1 for point-based data.

7. DATACOL: An integer (N) giving the number of “columns” in the data. This is the number of values per grid cell/point. N=1 for a scalar valued function; N>1 for a N-vector valued function.

• only a single column is currently supported.

8. DATAVAL: The values of the function on the cells/points of the grid. The values should appear in Fortran array order were the values stored in the Fortran array A(N,NX,NY,NZ) (A(N,0:NX,0:NY,0:NZ) for point-based data). That is, the column index varies fastest, x grid index next fastest, etc.

Example

As an example, consider the following integer-valued function in 2-D:

          +-----+-----+-----+ (2.0,3.0)
          |     |     |     |
          |  2  |  1  |  1  |
          |     |     |     |
          +-----+-----+-----+
          |     |     |     |
          |  5  |  1  |  2  |
          |     |     |     |
(0.0,0.0) +-----+-----+-----+

The corresponding input file would be:

0
2DCoRectMesh
3 2
0.0 0.0
2.0 3.0
0
1
5 1 2 2 1 1