Changes
Page history
move release Wiki to open-source repository
authored
Feb 28, 2020
by
Johannes Keller
Hide whitespace changes
Inline
Side-by-side
howtos/howto_forward.md
0 → 100644
View page @
bb817b1b
# 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
](
../tutorial/input_file
)
.
## Table of Contents ##
[[
_TOC_
]]
# 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 readibility 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, on 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
](
../tutorial/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 whitespace. 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 computaion 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 everyting,
`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!
