The outline can be defined by a number of polygons. The polygons can describe areas inside or outside (e.g.) holes of the domain. The format of the file is the following

n1 INOUT
X1 Y1
X2 Y2
. .
. .
Xn1 Yn1
n2 INOUT
X1 Y1
X2 Y2
. .
. .
Xn2 Yn2
. .
. .

where n1 n2 etc is the number of vertices of the polygon to follow and INOUT is 1 if the polygon describes area of the domain or 0 if it's a hole. Then n lines follow with the X and Y coordinates.

This pattern is repeated for all the polygons that the outline consist of. It is highly recommended to keep the outline as simple as possible.

The top and bottom function can be a constant scalar, a cloud of points, or a mesh. It is possible to define top and bottom as separate functions or as one function. The latter will be faster since the location interpolation is executed once for each point.

If top and bottom are separate functions then two files are needed with the following format

TYPE
[Radius Power] or [Npoints Nelem Radius]
X1 Y1 V1
X2 Y2 V2
.  .  .
.  .  .
id11 id12 id13 id14
id21 id22 id23 id24
.  .  .  .
.  .  .  .

where TYPE : is either CLOUD, MESH2D or GRID (NOT IMPLEMENTED YET) Depending on the TYPE the next line should be one of the following options

• CLOUD: [Radius Power] that set the radius and power of the IDW interpolation.
• MESH2D: [Npoints Nelem Radius] that set the number of nodes and elements of the mesh. Radius in this case is a search distance. The search distance must be set so that at least one element barycenter can be found from any point in the domain
• ⚠️ NOT IMPLEMENTED YET. GRID: [TopGridFilename BottomGridFilename]. The format of the grid files follows the gridinterp specifications which are listed here.

For CLOUD or MESH2D the next lines list the coordinates and the values. X Y V: the coordinates of the points and the elevation value. If top and bottom is defined on the same file then the bottom file can be omitted. In that case, this file section should include an extra column for the bottom.

X1 Y1 T1 B1
X2 Y2 T2 B2
.  .  .  .
.  .  .  .

For the MESH2D case only, after the coordinate and elevations the element ids should be listed. The elements indices are based on 1. The elements must be either quadrilaterals or triangles. When the elements are triangles the last index must be set as 0.

The elements must be listed in counter clockwise orientation

id11 id12 id13 id14
id21 id22 id23 0
.  .  .  .
.  .  .  .

The format for both files is the same

Nproc
ID1 n1
X1 Y1
X2 Y2
. .
. .
Xn1 Yn1
ID2 n2
X1 Y1
X2 Y2
. .
. .
Xn2 Yn2
. .
. .

where Nproc is the number of processors to use in the simulation, ID is the processor id which spans from 0 - Nproc-1. n is the number of vertices to follow

The attractors are described as vertical lines. The file format is

X1 Y1 T1 B1
X2 Y2 T2 B2
.  .  .  .
.  .  .  .

where X Y are the horizontal coordinates and T and B is the start and end of the line. This is used for well type attractors where T and B corresponds to the top and bottom of the screen length respectively.

The format of the particle file is the following

Eid Sid X Y Z RT
.   .   . . .

where

Eid is the entity id

Sid is the streamline id

x, y, z are the coordinates of the initial particle location. The entity id is a way to group particles.

RT is the release time. This is ignored for Steady state simulations.

The combination of Eid-Sid must define a single particle. For example all particles released from a certain elevation or subregion can have the same entity id Eid. The Sid is a unique number within the Eid group. Therefore, during particle tracking the pair Eid-Sid must unique. If not the algorithm will never be able to catch it and the results will be at least weird.

Lines that start with # are comments. Empty lines are not allowed.

The well file has the following format:

Npart Nlay rad
Eid X Y T B RT
.   . . . .

where the first line consist of three parameters.

Npart is the number of particles to release per well

Nlay is the number of layers

rad is the distance from the well where the particles will be released.

Eid is the entity id. It is assumed that each well is a different entity and the streamline's ids Sid will span from 0 - Npart-1

X Y T B are the coordinates of the well and the top and bottom elevation of the well screen.

RT is the release time. In steady state this is ignored.

Based on the above parameters Ichnos distributes the particles around the well in a layered/spiral kind of way. First calculates the number of particles per layer. It does so by dividing Npart/Nlay. However, these numbers are integers and the results will be an integer too. You cannot have 5.5 particles per layer. Therefore, the actual number of particles to be released will be (Npart/Nlay) * Nlay. For example if Npart=100 and Nlay=30 then the actual number of particles to be released per well is (100/30)*30 = 90 not 100.

Besides the total number of particles released, these two parameters define the distribution of particles around the well screen. The figure below shows the distribution of 100 particles for 20, 25 and 30 layers (from left to right). (For 20 and 25 the actual total number of particles is 100 but not for Nlay=30). Based on the Nlay={20,25,30} parameter, the code calculates 5, 4 and 3 particles per layer respectively. Each layer is then slightly rotated. A higher number of layers results in a distribution where two consecutive layers are rotated by only a small amount.

The format of the velocity files depends on the support structure and whether the velocity is steady state or transient and the file type e.g. ASCII or HDF5.

This is used when the velocity of the medium is defined on a cloud of points.

In Steady state we print the point coordinates and the velocities all in one file.

• ASCII format

X Y Z PROC DIAM RATIO VX VY VZ 

where:

X Y Z : are the coordinates of the point

PROC : is the id of the processor that owns the point. Note that each file list a number of points around the processor that are owned by other processors. However, each processor should be aware if actually owns the point or not.

DIAM : In numerical simulations each velocity point corresponds to a cell or an element. This is the horizontal maximum distance between the element vertices.

RATIO : Is the ratio between DIAM and the average thickness of the element.

VX VY VZ : are the components of the velocity.

• HDF5 format

For HDF5 format the file should contain 3 datasets as follows:

XYZDR : is a dataset with size np x 5 of type float. np is the number of point cloud and 5 is the number of columns that correspond to X Y Z DIAM RATIO.

PROC : is a dataset with size np x 1 of type integer with the ids of the processor that owns the point.

VXYZ : is a dataset with size np x 3 of type float. The 3 columns correspond to VX VY VZ.

• ASCII format

In the transient state the velocity data are printed into 4 files.

The name of this file should follow the naming convention Prefix_XYZ_#####.ext (see more here). The format of the file is the following:

X Y Z PROC DIAM RATIO

See above for the variable explanations

The names of these files should follow the naming convention Prefix_VX_#####.ext, Prefix_VY_#####.ext and Prefix_VZ_#####.ext.

Each file contains a 2D array where the rows correspond to the points of the XYZ file and the columns to the time steps. e.g:

V11 V12 V13 ...
V21 V22 V23 ...

where V11 is the velocity component of the point at the first row of XYZ file and the first time step.

• HDF5 format

In HDF5 format all the data are written in one file that follow the naming style of the steady state.

In transient state instead of a single VXYZ dataset there are 3 datasets with the following names VX VY VZ. Each dataset consist of a matrix np x nt where np is the number of points with known velocity and nt is the number of time steps. The type of the data should be float.

It is possible to supply a graph file so that the IDW interpolation does not use the interpolated points based on the distance. The Graph file shows how the cloud points are connected. Then the algorithm will search for the closest point and use only the connected nodes during the interpolation.

The format of the graph file is as follows

X Y Z Ncell Nvel C1 C2 C3 ... V1 V2

where
X Y Z are the coordinates of the point/element
Ncell is the number of points/elements that are connected with this row point/element.
Nvel is the number of velocity points that are associated with this point/element.
C1 C2 ... are the ids of the points/elements that this point/element is connected to.
V1 V2 ... are the ids of the velocity points that are associated with this point/element.

While it may seem redundant to print the coordinates in both files the velocity and the graph. However, this is because it is possible that each row of the graph file it may not necessarily represent a point but an element. In that case the coordinates correspond to the barycenter of the element and each element may contain more than one velocity points. e.g for higher than linear order elements.

The C1 C2 ... ids point to the rows of this graph file, while the V1 V2 ... point to the ids of the velocity file ie. VXYZ or VX,VY,VZ.

In the case of MESH2D support structure type there are 3 files that have to be prepared which define the mesh. Those files provide information about the mesh nodes, the mesh elements and the mesh layer elevation.

X Y

where X and Y are the nodes coordinates. These must use the same coordinate system as every other data that uses coordinates.

EL1 EL2 EL3 . . .

where EL# is the elevation of the nodes. The file should contain Nlayers+1 columns.

ID1 ID2 ID3 ID4 PROC

where ID# are the indices of the mesh. The indices should be based on 1 (not 0). The elements must be listed in counter clockwise orientation PROC is the processor id that owns the element. For single processor use 0.

The velocity on MESH2D uses the same format as in the CLOUD case. However we do not need to print any XYZ information.

For steady state the velocity format is

VX VY VZ 

In transient state the velocity is printed into 3 files where each file contains the X, Y and Z component of velocity using the format:

V11 V12 V13 ...
V21 V22 V23 ...

• ASCII For steady state the velocity name should follow the naming convention prefix_#####.ext while for transient state the component velocities are printed int o3 files as prefix_VX_#####.ext, prefix_VY_#####.ext, prefix_VZ_#####.ext
• HDF For HDF the velocity is printed into one file with the name convention prefix_#####.h5 using the VXYZ dataset for steady state or VX, VY, VZ.

Under MESH2D support structure ichnos provides three interpolation options. Depending the interpolation option the information printed into the files is used differently and therefore should be printed accordingly

• ELEMENT For element interpolation (the INTERP value of the MESH2D section is set to ELEMENT) the rows of the velocity file correspond to the rows of the Mesh file.
• NODE For node interpolation (the INTERP value of the MESH2D section is set to NODE) the rows of the velocity file correspond to the rows of the Node file.
• FACE For face interpolation (the INTERP value of the MESH2D section is set to FACE) the rows of the velocity file correspond to the rows of the Face file which is described below.

The face file is needed only when the INTERP value of the MESH2D section is set to FACE. This file contains the indices that indicate from which line to read the velocity and the direction.

Consider the following example.

Let's suppose that the top left element is described in the mesh file as 1 62 63 2. Based on this numbering:

• face 1 is 1-62 associated with velocity V4
• face 2 is 62-63 associated with velocity V2
• face 3 is 63-2 associated with velocity V3
• face 4 is 2-1 associated with velocity V1

Lets also suppose that the velocities V1, V2, V3, V4 are printed in the 1, 2, 3 and 4 row of the velocity file. Then the face ids for this element will be

-4 2 3 -1

The sign of the id is negative when the direction of flow is towards the element e.g. V1 and V4 otherwise is positive. If V3 has negative value which means that the flow is towards the element then we can still use the id 3 with positive sign as long as the velocity file prints the velocity with the negative sign.

In multi layer simulations the velocity is printed sequentially layer by layer as follows:

V1_1_1 V1_2_1 V1_3_1 ...
V2_1_1 V2_2_1 V2_3_1 ...
.
.
.
Vk_1_1 Vn_2_1 Vn_3_1 ...
V1_1_2 V1_2_2 V1_3_2 ...
V2_1_2 V2_2_2 V2_3_2 ...
.
.
.
Vk_1_2 Vn_2_2 Vn_3_2 ...

where Vn_t_m is the velocity of the n^th element/node/face, of the t^th timestep and m^th layer starting from the top to bottom. For ELEMENT and FACE interpolation the number of layers is equal to the velocity layers. For NODE interpolation there is an extra layer of velocity nodes. For the FACE interpolation first the horizontal face velocities from layer 1 to n are printed followed by the vertical velocities from layers 1 to (n+1) (see example).

The time step file is a sequence of time units. If the time units are days and the velocity field is based on a monthly simulation then the time step file will have the following format:

0
31
61
92
123
151
182
.
.
.

The main output file contains the streamline trajectories. The output file name has the following convention

prefix_ireal_####_iter_####_proc_0000.traj 

where

• prefix: is the OutputFile option under the section InputOutput.

• _ireal_####: For non-stochastic runs this is always _ireal_0000. Otherwise, the number indicates the realization number.

• _iter_####: Is the iteration number. This is always 0 if the total number of particles is lower than ParticlesInParallel of the InputOutput section. Otherwise, the particles will be printed into several files. The number of iterations depends on the number of total particles and the ParticlesInParallel option.

• _iproc_####: is the processor id. For single processor simulations this is always 0000

The format of the file is the following.

pid Eid Sid X Y Z VX VY VZ T
:
-9 Eid Sid EXIT_REASON

where

• pid: is the particle id that coincides with the number of steps
• Eid: is the Entity id
• Sid: is the streamline id
• X, Y, Z: are the coordinates of the particle's position at pid step
• VX, VY, VZ: are the components of the particle's Velocity at pid step
• T: is the time when the particle was found at the pid position. For steady state the initial time is 0. In backward tracking the time is negative and in forward tracking the time is positive. In steady state the time coincides with the age of the particle. In transient state T is the time. To convert the time to age you have to subtract the initial time.

At the moment the output file of the gather mode is slightly different. This will change in future version but at this point the output of gather mode contains the following data:

Eid Sid X Y Z V T
:
-9 Eid Sid EXIT_REASON

where V is the magnitude of the velocity.

The -9 pid denotes that this line holds information about the reason the particle terminated.

This is a list of reasons to terminate the tracing

• CHANGE_PROCESSOR: The tracing will continue to another processor
• EXIT_TOP: The particle exits from the top boundary of the domain
• EXIT_BOTTOM: The particle exits from the bottom boundary of the domain
• EXIT_SIDE: The particle exits from the lateral boundaries of the domain
• INIT_OUT: The initial particle position was outside the domain
• STUCK: The particle stuck in the flow field
• NO_EXIT: That's weird, isn't it? (luckily we haven't encountered that yet)
• MAX_INNER_ITER: The number of steps in the streamline exceed the limit
• FIRST_POINT_GHOST: The initial particle position was inside the domain but not inside this processor
• FAR_AWAY: For some reason the particle was found during the particle tracking far away the domain. (This is a good indication of a bug or something weird in the input files. Maybe there is a mismatch in the coordination system?)
• MAX_AGE: The particle stopped because the age reached the age limit
• ATTRACT: The particle was captured by an attractor feature
• NO_REASON: The particle stopped just because it could (we also have never encountered this).

Test expansible lists