Skip to content
Giorgos Kourakos edited this page Dec 17, 2024 · 45 revisions

Welcome to the Ichnos wiki!

Here you can find information about how to prepare the input files, run the code and reading the results.

To run the code one has to create a number of input files that describe the domain boundaries and the velocity field. There are also several options for each simulation. All the input files and simulation options are defined under two main files:

Run Ichnos

Ichnos is a console application therefore you can run it from a terminal as: (Under WinOS this can be the windows PowerShell)

ichnos.exe -c config.ini

The following is a list of input arguments to Ichnos executable.

  • -v : prints the version and exits.
  • -h : prints the list of options and exits.
  • -c config.ini : Runs the model using the parameters described in the config.ini file.
  • -g -n Nproc -i Niter : Runs the model in gather mode which assembles the streamlines from the various files. This is used as a second step in multicore simulations where each core traces only part of domain. Nproc is the number of processors used in particle tracking and Niter is the number of iterations. (See more)

Running Ichnos in Parallel

Ichnos provides two types of parallelizations

1) Single domain, multi-core/thread.

In this mode the particles are divided into N groups, where N is the number of cores/threads, and those groups are traced concurrently.

2) Multi subdomain, multi-core/thread

In this mode the domain is split into subdomains. Here the user is responsible to run the simulations with as many cores as the number of subdomains. Each subdomain traces the particles inside its boundaries. Once the particle moves out of the boundaries of the subdomain then the particle tracking continues to an adjacent subdomain. In this mode the particles are split according to their relationship with the subdomains. For example, if the domain is split into 32 subdomains but the initial positions of particles are all fall within one subdomain then the core associated with this subdomain will perform all the tracing and the remaining 31 cores would be idle.

In both cases the running command is identical Under windows a multicore simulation would start as follows:

C:\"Program Files\Microsoft MPI"\Bin\mpiexec.exe -n Nproc ichnos.exe -c config.ini

while under linux the command would be similar:

mpirun  -n Nproc ichnos.exe -c config.ini

where Nproc is an integer with the number of processors.

Getting started

To get started with Ichnos check the Getting Started which contains a few guides about how to prepare the input files, run the tracking simulation and read the results.

Configuration files

The following sections describe the available options of main and velocity configuration files.

The configuration files are ASCII. While any extension works, it is recommended to use the extension *.ini because the configuration file uses the *.ini format and most editors know how to display the files in a properly coloured readable format.

There are two important configuration files. The Model or Main configuration file which is explained next and the velocity configuration file which sets all the data to describe the velocity field.

Main Configuration file

The main configuration file sets the parameters for the particle tracking options, the description of the domain, the input output files etc.

The file is divided into several sections. The sections can be listed in any order. To set a section use the following format

[sectionName]

Lines that start with # are ignored and can be used for comments.

ℹ️ Not all the options are mandatory but most of them are. Options that expect strings as input e.g. filenames can be empty. However, options that take scalar values and they are not required, they should either not be listed at all or if they are listed they should have dummy values. ℹ️

Section Velocity

The velocity section of the main configuration file contains the following options

Options Values Description
XYZType Defines the type of the structure that supports the velocity. Typically, the velocity is defined on points or elements or faces etc. At this point there are 2 types avaliable:
CLOUD This represents a cloud of points where the velocity is defined, without any further information about how the points are connected.
MESH2D This corresponds to velocities calculated from numerical simulations that use an extruded mesh in the z direction
Type Defines the velocity type. It must be a keyword from the following:
DETRM This for steady and transient state deterministic velocity fields
STOCH This is a used for stochastic velocity types (Experimental)
RWPT This is used for Random walk (under development)
ConfigFile [string] This is velocity configuration file. Similarly to the main configuration file it is recommended to use the extension *.ini.

The section should look like the following:

[Velocity]
XYZType = CLOUD
Type = DETRM
ConfigFile = path\to\velocityConfig.ini

Section Domain

This section sets options about the domain of the flow field

Since Ichnos was developed with groundwater flow fields in mind, it is assumed that the flow domains are defined by a horizontal polygon along x and y. The vertical boundaries are defined as top and bottom 2D functions.

Options Values Description
Outline [string] The file name that contains the information of the horizontal polygon. The format of the file is described here
TopFile [string] or [double] This is the file name that contains the information for the top boundary of the domain. The format is described here. If the top is constant you can supply a scalar number.
BottomFile [string] or [double] This is the file name that contains the information for the bottom boundary of the domain. The format is described here. If the Bottom is constant you can supply a scalar number. If the bottom elevation is defined on the same points as the top you can omit this file.
ProcessorPolys [string] If more than one processor is used for the simulation this file describes the boundaries of the processor domains. The format is described here

The section should look like the following:

[Domain]
Outline = path\to\outline.ich
TopFile = path\to\topfunc.ich or scalar number
BottomFile = path\to\bottomfunc.ich or scalar number
ProcessorPolys = path\to\proc_subdomains.ich

Section StepConfig

The step configuration options are the following:

Options Values Description
Method This sets the method to calculate the next step. Valid options are:
Euler The simplest and fastest but not so accurate Euler method. Requires 1 velocity calculation per step but the step size has to be quite small. This is recommended for transient state velocity fields where the step size is very small anyway
RK2 2nd order Runge-Kutta. Requires 2 velocity calculations per step
RK4 4nd order Runge-Kutta. Requires 4 velocity calculations per step. Its more time-consuming but the step can be a bit larger
RK45 Adaptive Runge-Kutta. Requires 6 velocity calculations per step. Adapts the step according to the error. If the error is small can take larger steps. If the error is large reduces the step and reevaluates. It's more advance but requires more parameters to configure
RAL Ralston's method is a 2nd order method similar to RK2 with a minimum local error bound
PECE Predict-Evaluate-Correct-Evaluate or predictor corrector method it is an iterative method calculate the next step correcting iteratively the previous estimation. Requires k iterations where k is the order of the method a user defined parameter.
Direction [double] Defines the direction of tracing. Negative values indicate backward tracing and positive values (including 0) indicate forward particle tracking
StepSize [double] The length of the step size. In case of RK45 this is only the initial step
StepSizeTime [double] The length of the step size in time units
nSteps [integer] The number of steps to take within the area of influence of the closest point
nStepsTime [integer] The number of steps to divide the time steps. This is meaningful in transient state velocity fields
minExitStepSize [double] When the particles are about to exit the domain it's a good practice to reduce the step. This is particularly important for the RK45 where the step can get quite large. This set the maximum step size when the particle is about to leave the domain.

See also the Step Control section for additional discussion regarding the step size configuration.

The section should look like the following:

[StepConfig]
Method = Euler
Direction = -1
StepSize = 50
StepSizeTime = 50
nSteps = 4
nStepsTime = 2
minExitStepSize = 0.1

Section AdaptStep

This is a list of parameters to further configure the step size when the adaptive method RK45 is selected. These are all ignored for the other methods

Options Values Description
MaxStepSize [double] This is the maximum step size. Even if the error is very small the step will never get greater than this
MinStepSize [double] This is the minimum step size. Even if the error is greater than the tolerance the step will not get any smaller. If the error is greater than tolerance and the step size is already minimum the algorithm will print a warning and continue by accepting the new point regardless of the error
IncreaseRateChange [double] When the error is small and the step size can increase this is the maximum allowable increase with respect to the previous step. This allows the step size to grow smoothly
LimitUpperDecreaseStep [double]
Tolerance [double] If the error is less than the tolerance the step size increases. If the error is higher than the tolerance the step size decreases

The section should look like the following:

[AdaptStep]
MaxStepSize = 15
MinStepSize = 0.1
IncreaseRateChange = 1.5
LimitUpperDecreaseStep = 0.15
Tolerance = 0.1

Section PECE

The following options are used only with the Predict-Evaluate-Correct-Evaluate (PECE) method

Options Values Description
Order [integer] It's the order of PECE method or the number of iterations
Tolerance [double] If the predictions of two consecutive iterations have a difference less that Tolerance, PECE stops the iterations

The section should look like this

[PECE]
Order = 6
Tolerance = 0.2

Section StoppingCriteria

The following is a list of stopping criteria. Some are reasonable reasons to terminate particles and some are used to avoid the algorithm run forever when the particles are stuck in the flow field.

Options Values Description
MaxIterationsPerStreamline [integer] Maximum number of steps per streamline
MaxProcessorExchanges [integer] Maximum number of times a particle is allowed to change processors
AgeLimit [double] If the particle exceed this age limit the particle tracking stops. Negative means no limit
Stuckiter [integer] Each time step the code calculates the bounding box of the streamline. The bounding box should always expand. If for a consecutive number of iterations there is no expansion of the bounding box the particle tracking stops
AttractFile [string] If we want the particles to stop near some features of the domain such as wells we can define a list of attractors which have the property to attract the particles and terminate their paths. Currently, only wells can be defined as attractors. The format for the attractors file is described here. This should be left blank if there are no attractors
AttractRadius [double] If the particle approach any attractor closer that this distance it will terminate particle tracing

The section should look like the following:

[StoppingCriteria]
MaxIterationsPerStreamline = 1000
MaxProcessorExchanges = 50
AgeLimit = 100
Stuckiter = 10
AttractFile = path\to\atractfile.ich
AttractRadius = 30

Section InputOutput

This is a list of options that define the input particles and the output files

Options Values Description
PartilceFile [string] A file with the initial position of the particles. The format is described here
WellFile [string] Instead of setting the particles we can set the locations of the well screen and let the code generate the initial particle locations around the well screen. The format of the file is described here
OutputFile [string] This is a prefix for the output files
PrintH5 [boolean] Print the output files as HDF5
PrintASCII [boolean] Print the output files as ASCII
Pprecision [integer] Specifies the precision of particles position. This is the number of decimal digits.
Vprecision [integer] Specifies the precision of particles velocities. The velocity is printed as scientific notation #.#####e+#
Tprecision [integer] Specifies the precision of particles time. The time is printed as scientific notation #.#####e+#
ParticlesInParallel [integer] Lets suppose that you want to run 1 million particles. Running them at once will produce a huge output file. Instead, one can set how many particles to run simultaneously (in practice this is not simultaneously) and print one output file for a group of N particles
GatherOneFile [boolean] Obsolete option

The section should look like the following:

[InputOutput]
PartilceFile = path\to\particles.ich
WellFile = path\to\wells.ich
OutputFile = path\to\out_prefix_
PrintH5 = 0
PrintASCII = 1
ParticlesInParallel = 5000
GatherOneFile = 0

Section Other

The parameters that do not fit to any category are specified under this section.

Options Values Description
Version [string] This is the version number of the executable. As the code evolves each version may have different parameter options. To make sure that there is no mismatch between input files and executables the code test if the version specified here is the correct one. You can find out about the version of the executable by running ichnos.exe -v
Nrealizations [integer] This is the number of realizations. It is used in stochastic simulations
nThreads [integer] The number of threads if you are running in multithreaded mode
RunAsThread [boolean] Run the code as Multithreaded. This is used in multi-core simulations when the particle work load is split into several processors which do not share the memory
OutFreq [double] This is used in RunAsThread mode and indicates how frequently you would like the console to report the progress
PrintVex [boolean] Prints to console the position and velocity at each step of tracing as VEX code. This is used for debuging the code and should always be 0 or not present at all which defaults to 0.

The section should look like the following:

[Other]
Nrealizations = 1
nThreads = 1
RunAsThread = 0
Version = 0.4.10
OutFreq = 1

Velocity configuration file

The velocity configuration file depends on the type of the velocity.

Similar to the main configuration file this file is split into sections

Section Velocity

The first 3 options of the table below define the name of the file that the velocity is printed. Depending on the type of velocity field the naming convention and the format of the files are different.

Steady state

  • ASCII & HDF format

For steady state, the velocity files may contain the XYZ coordinates, additional point information and the 3 components of the velocity. The filename should strictly follow the naming convention: Prefix_#####.ext e.g. examp1_v2_new_0000.vel or example1__004.ich. The format of the velocity file for the steady state is given here

Transient state

  • ASCII format

For transient state velocity fields the point coordinates and the velocities are printed into different files. The file with the XYZ data should have the format Prefix_XYZ_#####.ext while the velocity is printed into 3 different files with filename format Prefix_VX_#####.ext, Prefix_VY_#####.ext and Prefix_VZ_#####.ext for the x, y, z component of velocity respectively.

Face interpolation (Steady and Transient state)

In the case of FACE interpolation under the MESH2D option the velocity does not require three components because the direction is dictated by the face orientation. The FACE velocity in this case is split into horizontal and vertical component. Therefore the horizontal velocity is printed into the Prefix_VHOR_#####.ext file and the vertical into Prefix_VERT_#####.ext

  • HDF5 format

In HDF5 format all the data are printed into one file under different dataset names. Therefore, the filename of the velocity should be similar to the steady state case Prefix_#####.ext. The following datasets names are valid:

Dataset Name Type Applicable for Description
XYZDR [float] CLOUD Steady and Transient state The X,Y,Z coordinates, the Diameter and ratio parameter. The size must be equal to the cloud points
PROC [integer] CLOUD Steady and Transient state The Processor id that each XYZDR point belongs to. The size must be equal to the cloud points
VXYZ [float] CLOUD,MESH2D Steady state The X, Y, Z component of velocity. The size must be equal to the cloud points or mesh elements or mesh nodes times the number of layers, depending the interpolation type
VX [float] CLOUD,MESH2D Transient state (expect FACE) The X component of velocity. The size must be equal to the cloud points or mesh elements or mesh nodes times the number of layers, depending the interpolation type
VY [float] CLOUD,MESH2D Transient state (expect FACE) The Y component of velocity. The size must be equal to the cloud points or mesh elements or mesh nodes times the number of layers, depending the interpolation type
VZ [float] CLOUD,MESH2D Transient state (expect FACE) The Z component of velocity. The size must be equal to the cloud points or mesh elements or mesh nodes times the number of layers, depending the interpolation type
VHOR [float] MESH2D Steady and Transient state (FACE ONLY) The horizontal component of velocity. The size must be equal to the mesh faces times the number of layers
VVER [float] MESH2D Steady and Transient state (FACE ONLY) The vertical component of velocity. The size must be equal to the mesh faces times the number of layers + 1

The remaining options are explained in the table below:

Options Values Description
Prefix [string] This is the prefix of the velocity files. It can also be the full path if the files live in a different folder. The prefix in the above examples are examp1_v2_new_ and example1__ respectively.
LeadingZeros [integer] The second part of the velocity name is the processor id. For single processor simulations this should be always 0. This parameter controls the length of the number. For the examples above this should be set to 4 and 3 respectively
Suffix [string] This is the extension of the file. It can be anything such as .txt, .dat, .vel, .ich etc. For HDF5 the suffix should always be .h5 (lower case)
Type Defines the velocity type. It must be a keyword from the following: (Note that this parameter is also defined on the model configuration file. Both files should have the same type)
DETRM This for steady and transient state deterministic velocity fields
STOCH This is a used for stochastic velocity types (Experimental)
RWPT This is used for Random walk (under development)
TimeStepFile [string] The name of the file which contains the time steps. Leave empty for Steady state
TimeInterp Interpolation time between the time steps. Valid options are :
NEAREST During a time step the velocity is constant. Between t1 and t2 the velocity is set equal to v1
LINEAR The velocity is linearly interpolated between the time steps
RepeatTime [double] This is a double number that corresponds to the units of the time step file. In transient velocity fields the simulation range is limited. In C2VSim for example is around 40 years. Given that the groundwater velocities are rather slow, 40 years are almost never enough for the particles to exit the domain. This parameter controls what happen after the particle has exceeded the last time step. The last time step in backward tracing is the first step. Setting this parameter to zero will continue tracing the particles using constant velocity equal to the last step. Any other number will force tracing to loop N times. For example if the time units are days, setting this parameter to 365 the tracing will continue by repeating the last year of the velocity field
Multiplier [double] The velocities will be multiplied by this value before they are used. Sometimes codes multiply the groundwater velocity with large numbers to increase the precision. This multiplier is used to reverse that. For example if the velocity is 0.00004335701, then one need to store 12 digits. If we multiply this by 100000 the number will be 4.335701 where we have the same precision with fewer digits to store in the ASCII files. This value is multiplied with the velocity therefore if the printed values have been multiplied by 1000 then the multiplier should be 0.001

The section should look like the following:

[Velocity]
Prefix = path\to\Vel_05_15_new_
LeadingZeros = 4
Suffix = .ich
Type = DETRM
TimeStepFile = c2vsim_STEP_monthly.ich
TimeInterp = NEAREST
RepeatTime = 365
Multiplier = 1

Section Cloud

If the XYZType is CLOUD then the cloud section list a few parameters for cloud velocity fields.

Options Values Description
Scale [double] Not a very useful parameter. Maybe it will be removed in future versions. Until then, it should always be set to 1
Power [double] The velocities are defined as a cloud of 3D points. This parameter is the power of IDW
InitDiameter [double] This is the search radius of the IDW interpolation. However, during particle tracking this changes depending on the density of velocity points found around the particle. The algorithm needs a starting value. This value should be set so that even in the coarser areas at least 3-4 velocity points can be found
InitRatio [double] It is the ratio between the elements diameter and z extend. This should be close to an average ratio. It is used only as a starting value and it will be ignored afterwards
GraphPrefix [string] It is possible to use a graph with the cloud connectivity. This has the advantage that the IDW interpolation will always use fewer and the correct points from the cloud. The downside is that requires the preparation of a rather complicated file especially in multi-processor simulations. This is the prefix of the graph file. The naming convention is path/to/graph/file_######.grph where the number of zero padding must be equal to the one of the velocity file
Threshold [double] If the particle position is closer than this threshold to a velocity point then the velocity of the particle is set equal to the velocity point

The section should look like the following:

[Cloud]
Scale = 1
Power = 3
InitDiameter = 3000
InitRatio = 70
GraphPrefix =
Threshold = 0.1

Section Porosity

The porosity can be set as a single value or interpolation cloud. If the porosity is variable and defined on the same points/elements that the velocity is defined it is better to bake the porosity in the velocity field and set the porosity to 1

Options Values Description
Value [double] or [string] A scalar value or the name of the file that describes the porosity

The section should look like the following:

[Porosity]
Value = 0.1

Section Other

These are few general options

Options Values Description
OwnerThreshold [double] In multiprocessor simulation it is possible that a particle is moving in some overlapping areas between the processors. This parameter is a percent indicating when the running processor should stop tracing the particle. If the processor owns the particle with a percent less that OwnerThreshold then the tracing will stop for the processor and will resume to the next processor during the next particle exchange. It's good idea to avoid frequent exchanges
FrequencyStat DEPRECATED Every FrequencyStat velocity calculations the program will print some stats about how long it takes to do velocity calculations. This is useful in very large velocity fields for debugging purposes

The section should look like the following:

[General]
OwnerThreshold = 0.15
FrequencyStat = 1000

Section MESH2D

This section is needed if the MESH2D support structure is used. The parameters of MESH2D section are the following:

Options Values Description
NodeFile [string] The name of the file that contains the node coordinates. The coordinates must use the same coordinate system as the rest of the file. The format is given here
MeshFile [string] The name of the file that contains the indices of the mesh elements and the processor they belong to. The format is given here
ElevationFile [string] The name of the file that contains the layer elevations at the nodes of the previous file. The format is given here
Nlayers [integer] The number of layers
FaceIdFile [string] The name of the file with the face ids. This is used only for the FACE type interpolation
INTERP [ELEMENT] or [NODE] or [FACE] This is the type of interpolation
[MESH2D]
NodeFile = c2vsim_Nodes.ich
MeshFile = c2vsim_MSH.ich
ElevationFile = c2vsim_ELEV.ich
FaceIdFile = 
Nlayers = 4
INTERP = ELEMENT