Howto: Forward computation

This is a step-by-step introduction, how to write a small SHEMAT-Suite input file for an easy forward computation.

For a almost-exhaustive list of inputs to SHEMAT-Suite, see Input File.

Table of Contents

• Howto: Forward computation
  • Table of Contents
• Make a SHEMAT-Suite model directory and input file
• Write input file entries
  • Grid specifications
  • Time step control
  • Solver control
  • Boundary conditions
  • Initial values
  • Properties
  • Geometry
  • Output
  • Make shemade.job, compile SHEMAT-Suite and execute

Make a SHEMAT-Suite model directory and input file

First, we make a new model directory. Your terminal should now be in your home directory, or in a directory of your choice, where you want to generate the new directory. Best practice: make the model directory outside the SHEMAT-Suite code directory in order not to interfere with the git repository. to make the model directory, type

mkdir theis_model

then move into the model directory typing

cd theis_model

and create the input file with the editor of your choice (emacs is the editor of my choice)

emacs THEIS

It is not strictly required, but naming the input files in capitals will nicely distinguish it from other files in the model directory (and there will be other files!).

Write input file entries

Now, we want to write entries in the input file. For that, you should use the text editor of your choice (emacs is definitely a nice choice!).

In general you could write anything to the input file that does not contain a hashtag (#) and SHEMAT-Suite will not care. Only hashtags followed by a correct keyword will be read in, so we can not make too many mistakes.

First, let's start by inserting a header that will not be read by SHEMAT-Suite, but it will help us, when we read the input file again.

! # =====>>>>>-----------------------------------------------------<<<<<===== # !
! # =====>>>>>   RUN INFO                                          <<<<<===== # !
! # =====>>>>>-----------------------------------------------------<<<<<===== # !

You can copy these code examples from here and paste them in your THEIS file.

OK, now we can insert the first keyword (after adding one or two newlines for readability reasons), which looks as follows:

# title
theis

These lines specify the internal title of the simulation. It makes sense to give it the same name as for the input file. But we want to look at the keyword a little more closely. In this simplest case, the input consists of three parts

1. The hashtag #
2. The keyword
3. The input

The hashtag is very important: Without it, the keyword or input will not be read at all. It has to be at the beginning of the line (for example adding a ! before # will comment it out). If a hashtag is found, SHEMAT-Suite will start looking for keywords in the line behind the hashtag. There is a certain number of keywords that will be read by SHEMAT-Suite, one of those is title, between the hashtag and the keyword, there has to be a space. Once the keyword is read, SHEMAT-Suite will look for the input that fits the keyword. The form of the input is thus dependent on the keyword. For title, the program looks for one word, in our case receiving theis.

Now, we include some more keywords that set general input for the simulation that we are planning:

# linfo
2 3 2 3

The keyword linfo looks for four integer numbers, these numbers specify, how much documentation the program will print. In this tutorial, the information about the inputs will not be complete, for more information, you can always look up the keyword in Input File.

# runmode
0

# active head

# PROPS = const
# USER = none

Via runmode one can choose different versions of the program, the input 0 gives simple forward computation. Under active the computed physical variables (in this case the hydraulic head) are specified. This is a special keyword, because it looks for the input in the same line as the hashtag. The same holds for PROPS and USER which give compilation information that has to fit to the SHEMAT-Suite executable. Depending on these choices, SHEMAT-Suite will access different parameter-computation (PROPS) or user defined modules.

Now, we entered the first few inputs, but they were very general. In the following we will start setting the stage for the simulation by specifying the grid.

Grid specifications

Now, we will specify the grid. The grid consists of rectangular cells of different sizes in x-, y-, and z-direction.

! # =====>>>>>-----------------------------------------------------<<<<<===== # !
! # =====>>>>>   GRID                                              <<<<<===== # !
! # =====>>>>>-----------------------------------------------------<<<<<===== # !

# grid
69 69 1

# delx
29*0.15625d0 6*1.666667d0 3*3.333333d0 4*5.d0 27*10.d0

# dely
27*0.15625d0 6*1.666667d0 3*3.333333d0 5*5.d0 29*10.d0

# delz
20.d0

The input under # grid specifies the number of cells: 69 cells in x and y direction, only one cell in z-direction. Thus, we are talking about a horizontal, two-dimensional model.

After specifying the number of grid cells, we specify their size. Specifying the size of the grid cells is done with keywords like # delx. The input consists of an array of float numbers.

Now we will look a little closer at the input under # delx.

# delx
29*0.15625d0 6*1.666667d0 3*3.333333d0 4*5.d0 27*10.d0

The most general input for # delx would contain a number of entries (69 entries in our case) separated by white-space. To save space, it is possible to group equal entries together, so for example 4*5.d0 is equivalent to 5.d0 5.d0 5.d0 5.d0. The second peculiarity of the input is the d0. This is due to the fact that SHEMAT-Suite is written in Fortran. In Fortran, there are different kinds of floating point numbers, for example with single precision and with double precision. The d0 is a way of expressing that the entry should be read as double precision floating point number. Note: When large or small numbers are given, the d is equivalent to the better known e for specifying the exponent of ten in a floating point number (for example 1.0d3 would be equivalent to 1000.0d0, i.e. one thousand.

We have specified the grid, now we will specify the time stepping control of the simulation.

Time step control

The time step control is minimal for this simple example, since it computes a steady state of the simulation. Thus, we do not need to specify any special time steps - special equations are implemented in SHEMAT-Suite to compute the steady state of a simulation.

! # =====>>>>>-----------------------------------------------------<<<<<===== # !
! # =====>>>>>   TIME STEP CONTROL                                 <<<<<===== # !
! # =====>>>>>-----------------------------------------------------<<<<<===== # !
! # ************************************************************************* # !
! # *** timestep control 0: steady state   timestep control 1: transient  *** # !
! # *** tunit: 31557600 year                                 *** # !
! # ***           86400 day                              *** # !
! # ***            3600 hour                                 *** # !
! # ***              60 minute                               *** # !
! # ***               1 second (default)                 *** # !
! # ************************************************************************* # !

# timestep control
0

The input for # timestep control consists of one line. The line should contain a single number, 0 for steady state and 1 for transient simulation. For transient simulations, more inputs are necessary (# timestep control would get another line of input, for example).

Solver control

After specifying some geometric inputs, this section will deal with a more abstract input, the solver parameters.

! # =====>>>>>-----------------------------------------------------<<<<<===== # !
! # =====>>>>>   SOLVER CONTROL                                    <<<<<===== # !
! # =====>>>>>-----------------------------------------------------<<<<<===== # !

# nlsolve                                  !**non-linear solver control**!
50 0

! # =====>>>>>   SOLVER FLOW

# lsolvef (linear solver control)
1.d-10 64 500

# nliterf (nonlinear iteration control)
1.0d-10 1.0d0

There are two layers of solvers, a nonlinear solver and, then inside, a linear solver. The inputs # nlsolve and #nliterf specify inputs for the outer, nonlinear solver. The two inputs, which should be known are the first inputs under each keyword. The first input under # nlsolve is the maximum number of iterations executed by the nonlinear solver (the solver will iterate until the differences between the iterations fall below a threshold or until this maximum number of iterations is reached, so a smaller number would mean: Faster, but less accurate computation). The first input under # nliterf is the threshold mentioned in the last sentence. Enlarging this number would mean faster, but less accurate computation. The other inputs under these keywords are more specialized solver parameters, which should be kept in the usual state for now.

Under # lsolvef the linear solver controls are specified. The first number is the accuracy-threshold similar to the nonlinear solver, the last number is again the maximum number of iterations (note that the linear solver should in general be more accurate than the nonlinear solver). The middle number is a specification of the kind of linear solver that should be used. This can be left as 64 for most simulations.

Boundary conditions

Boundary conditions are an important part of any simulation. Here, we need to specify hydraulic head values at the boundaries.

! # =====>>>>>-----------------------------------------------------<<<<<===== # !
! # =====>>>>>   BOUNDARY CONDITIONS                               <<<<<===== # !
! # =====>>>>>-----------------------------------------------------<<<<<===== # !
! # ************************************************************************* # !
! # *** dirichlet: bcd   neumann: bcn                                *** # !
! # *** simple = top/base/back/front/left/right                           *** # !
! # ************************************************************************* # !

! # =====>>>>>   BC HEAD

# head bcn, records=1, direction=0
1 1 1 -0.001d0 0 

# head bcd, simple=right, value=init
# head bcd, simple=back, value=init

There are two fundamentally different boundary condition specified, the first one specifies one boundary condition at a single cell (with index 1 1 1). At this cell a flow boundary condition is specified (bcn) with x direction and magnitude 0.001d0.

The second type of boundary conditions are specified with simple. This gives values to a whole boundary of the specified side of the model (the back and the right in this case). This time actual hydraulic head values are specified (Dirichlet boundary condition bcd). The values are taken from the initial conditions. Their specification will be treated in the next section.

Two more boundaries (left and front) are missing from the inputs. In this case, SHEMAT-Suite sets the boundary condition to Neumann-Boundary conditions (flow, bcn) with 0.0d0 flow, the so-called no-flow boundary condition.

Initial values

Initial variable values throughout the model domain are specified.

! # =====>>>>>-----------------------------------------------------<<<<<===== # !
! # =====>>>>>   INITIAL VALUES                                    <<<<<===== # !
! # =====>>>>>-----------------------------------------------------<<<<<===== # !
! # ************************************************************************* # !
! # *** give array with initial values or hdf5 file with initial values:  *** # !
! # *** # head init, HDF5=initial.h5                                 *** # !
! # ************************************************************************* # !

# head init
4761*20.0d0

# temp init
4761*20.0d0

The input array needs to have as many entries as there are grid cells. Initial temperatures have to be entered, even though temperature is not computed in our simulation. SHEMAT-Suite always asks for temperatures, since they may affect parameter values.

Properties

The main input that we are still missing, are the physical parameters of the fluid (mostly water) and the rock that are simulated.

! # =====>>>>>-----------------------------------------------------<<<<<===== # !
! # =====>>>>>   PROPERTIES                                        <<<<<===== # !
! # =====>>>>>-----------------------------------------------------<<<<<===== # !

! # =====>>>>>   FLUID
# fluid props
1000.0d0 5.0D-8 4218.0d0 0.65d0 1.0D-3 0.02d0              !***rhof compf cpf lamf visf sigmaf***!

The fluid properties are specified in an array of floating point numbers. Most important for the head computation are the first two: density (rhof) and compressibility (compf).

! # =====>>>>>   ROCKS/UNITS
! # ************************************************************************* # !
! # *** poro | perm-anisotropy in x | perm-anisotropy in y | perm |       *** # !
! # *** compressilbility | tc-anisotropy in x | tc-anisotropy in y | tc | *** # !
! # *** heat production | volumetric thermal capacity | dispersivity      *** # !
! # *** electrical conductivity | coupling coefficient                    *** # !
! # ************************************************************************* # !

# units
0.1d0  1.d0  1.d0 1d-12 1.d-8  1.d0  1.d0  3.35d0  0.0e-06  2300000.d0  10d0  0.d0  0.d0   !Unit1

There is a whole number of rock parameters, but not all are important for head computations. The first five are most important: Porosity, permeability anisotropy parameters (if equal to one, the permeability will be isotropic), permeability in z-direction and compressibility.

Geometry

The geometry becomes important, when more than one parameter unit is specified. In our case, there is only one set of rock parameters, thus the input is an array of ones.

! # =====>>>>>-----------------------------------------------------<<<<<===== # !
! # =====>>>>>   GEOMETRY                                          <<<<<===== # !
! # =====>>>>>-----------------------------------------------------<<<<<===== # !

# uindex
4761*1

Output

Finally, an input that specifies the kind of output written by SHEMAT-Suite. For a simple first run, the vtk-output is sufficient. It can be visualized with the Paraview visualization program.

! # =====>>>>>-----------------------------------------------------<<<<<===== # !
! # =====>>>>>   OUTPUT                                            <<<<<===== # !
! # =====>>>>>-----------------------------------------------------<<<<<===== # !

# file output: vtk

Under # file output the input vtk is given in the same line as the keyword. This is an exception, since most inputs are given in the line after the keyword.

Make shemade.job, compile SHEMAT-Suite and execute

Create a file shemade.job and write the name of the input file in it

THEIS

This is everything, shemade.job should contain. Especially empty lines will lead to errors.

If not done already, compile SHEMAT-Suite and move the executable to model directory.

Now you may execute this first small SHEMAT-Suite model! Congratulations!