Simple Face Example

In this example we will trace particles in a 3x3 mesh grid, where each element is quadrilateral with dimensions 100x100 m.

Domain

Outline

The domain is a 300x300 m rectangle

domain_poly = [0 0;300 0;300 300; 0 300];
fid = fopen('simpleFace_outline.ich', 'w');
fprintf(fid, '%d %d\n', [size(domain_poly,1) 1]);
fprintf(fid, '%.1f %.1f\n', domain_poly');
fclose(fid);

Top and Bottom

For simplisity we will keep the top and bottom elevations constant. Therefore we can define under the Domain section the following options:

TopFile = 10
BottomFile = 0

Prepare Mesh

Nodes

[NDX, NDY] = meshgrid(0:100:300,0:100:300);
NDY = flipud(NDY);
[IX, IY] = meshgrid(1:4,1:4);
linind = sub2ind(size(IX),IY,IX);
idx = reshape(linind,4*4,1);
ND = [NDX(idx) NDY(idx)];

Print node file

fid = fopen('SimpleFace_Node.ich','w');
fprintf(fid,'%.2f %.2f\n', ND');
fclose(fid);

Mesh

msh = nan(9,4);
cnt = 1;
for ii = 1:size(linind,1)-1
    for jj = 1:size(linind,2)-1
        msh(cnt,:) = [linind(ii,jj) linind(ii+1,jj) linind(ii+1,jj+1) linind(ii,jj+1)];
        cnt = cnt + 1;
    end
end

The elements must be oriented in counter clockwise direction. To check this we do the following which should result all zeros

for ii = 1:9
    c(ii,1) = ispolycw(ND(msh(ii,:),1), ND(msh(ii,:),2));
end

Print the mesh file

fid = fopen('SimpleFace_Mesh.ich','w');
fprintf(fid,'%d %d %d %d %d\n', [msh zeros(size(msh,1),1)]');
fclose(fid);

Elevation file

The elevation file defines the layers. Here we will use only one layer

Elev = ones(size(ND,1),2);
Elev(:,1) = Elev(:,1)*10;
Elev(:,2) = Elev(:,2)*0;
fid = fopen('SimpleFace_Elev.ich','w');
fprintf(fid,['%.2f' repmat(' %.2f',1,size(Elev,2) - 1) '\n'], Elev');
fclose(fid);

Face ids

The following loop builds a list of unique face ids and assignes ids the the mesh faces.

First we initialize the face_nodes variables which will hold the indices of the node ids that compose each face. We also initialize the face_Vels variable which holds the face ids (the row in the face_nodes variable where each element face can be found). cnt_fc is a counter. Since we dont know the exact number of unique faces that we will end up.

face_nodes = nan(36,2);
face_Vels = zeros(size(msh,1),4);
cnt_fc = 1;

Now we loop through the elements ii and then through the element faces j. The ida, idb hold the node ids of the current face. When we visit the first element ii==1 we know that all faces are new. If not, we check to see if the face already exist in the face_nodes variable. The face ids can be stored as [ida idb] or [idb ida] therefore we check both ways. If the face [ida idb] is new eg. isempty(idx) = true, we add a new face, if not we add the face id in the face_Vels variable.

for ii = 1:size(msh,1)
    for j = 1:size(msh,2)
        if j == 4
            ida = msh(ii,4);
            idb = msh(ii,1);
        else
            ida = msh(ii,j);
            idb = msh(ii,j+1);
        end
        if ii == 1
            face_nodes(cnt_fc,:) = [ida idb];
            face_Vels(ii,j) = cnt_fc;
            cnt_fc = cnt_fc + 1;
        else
            idx = find(ida == face_nodes(1:cnt_fc-1,1) & idb == face_nodes(1:cnt_fc-1,2));
            if isempty(idx)
                idx = find(idb == face_nodes(1:cnt_fc-1,1) & ida == face_nodes(1:cnt_fc-1,2));
            end

            if isempty(idx)
                face_nodes(cnt_fc,:) = [ida idb];
                face_Vels(ii,j) = cnt_fc;
                cnt_fc = cnt_fc + 1;
            else
                face_Vels(ii,j) = idx;
            end
        end
    end
end
face_nodes(cnt_fc:end,:) = [];

Following the Ichnos convention for the face flow direction (see here), we reverse the signs of the first and fourth faces for each element.

face_Vels(:,1) = -1*face_Vels(:,1);
face_Vels(:,4) = -1*face_Vels(:,4);

Print the face id file

fid = fopen('SimpleFace_faceIds.ich','w');
fprintf(fid,'%d %d %d %d\n', face_Vels');
fclose(fid);

Velocity

Since this example is not a product of a simulation, we will hard code the velocity field using the following flow directions

The next variable is nx3. The first and second columns contain the nodes ids of the face with a non zero velocity and the third column is the velocity value

NonZeroFaces = [3 7 -1;...
                8 7 -0.5
                2 6 1;...
                6 7 2;...
                7 11 -0.5;...
                6 10 0.5;...
                10 11 3;...
                11 15 -0.5;...
                10 14 -4];

Notice that the flow through the face 3 - 7 is negative. For the left lower element [3 4 8 7] the face 3-7 is the 4th and a positive sign velocity sign indicates that the velocity aligns with the direction indicated by the face ids. The face id for the face 3-7 is -12 (see face_Vels(7,4) 7th elementm 4th face). Negative face id denotes inflow towards the element but here we need outflow, therefore we reverse the sign. The same face for the element above this [2 3 7 6] is the 2nd which has a positive face id = 12 (see face_Vels(4,2) 4th element 2nd face). The positive sign for this face indicates outflow but for this element we need inflow therefore we need again to reverse it. Basicaly the reverse sign works for both elements.

Next we initialize all velocities to zero and then we change the non zero ones. Again we have to search for the face ids to find the right face.

faceVelValues = zeros(size(face_nodes,1),1);
for ii = 1:size(NonZeroFaces,1)
    ida = NonZeroFaces(ii,1);
    idb = NonZeroFaces(ii,2);
    idx = find(ida == face_nodes(1:cnt_fc-1,1) & idb == face_nodes(1:cnt_fc-1,2));
    if isempty(idx)
        idx = find(idb == face_nodes(1:cnt_fc-1,1) & ida == face_nodes(1:cnt_fc-1,2));
    end
    faceVelValues(idx,1) = NonZeroFaces(ii,3);
end

Print the horizontal velocities

fid = fopen('SimpleFace_vel_VHOR_000.ich','w');
fprintf(fid,'%.2f\n', faceVelValues');
fclose(fid);

For the vertical velocities we assume that the bottom is zero and top equal to 0.5

faceVelVerVal = zeros(size(msh,1),2);
faceVelVerVal(:,1) = 0.5;
faceVelVerVal = reshape(faceVelVerVal, size(faceVelVerVal,1)*size(faceVelVerVal,2),1);

Print the vertical velocities file

fid = fopen('SimpleFace_vel_VVER_000.ich','w');
fprintf(fid,'%.2f\n', faceVelVerVal');
fclose(fid);

Particle positions

For this example we will release particles from the left lower element as follows:

xpart = linspace(5.32, 95.45,42)';
pp = [ones(length(xpart),1) (1:length(xpart))' xpart 10.32*ones(length(xpart),1) 5.78*ones(length(xpart),1)];
fid = fopen('SimpleFace_part.txt','w');
fprintf(fid,'%d %d %.2f %.2f %.2f 0\n',pp');
fclose(fid);

Main Configuration file

Velocity

[Velocity]
XYZType = MESH2D
Type = DETRM
ConfigFile = SimpleFace_vf.ini

Domain

[Domain]
Outline = simpleFace_outline.ich
TopFile = 10
BottomFile = 0
ProcessorPolys = 

Step configuration

[StepConfig]
Method = Euler
Direction = 1
StepSize = 5
StepSizeTime = 100000
nSteps = 5
nStepsTime = 0
minExitStepSize = 1

Input output

[InputOutput]
ParticleFile = SimpleFace_part.txt
WellFile = 
OutputFile = output/SimpleFace_out
PrintH5 = 1
PrintASCII = 1
ParticlesInParallel = 5000
GatherOneFile = 0

Velocity Configuration file

Velocity

[Velocity]
Prefix = SimpleFace_vel_
LeadingZeros = 3
Suffix = .ich
Type = DETRM
TimeStepFile =
TimeInterp =
RepeatTime = 0
Multiplier = 1

MESH2D

[MESH2D]
NodeFile = SimpleFace_Node.ich
MeshFile = SimpleFace_Mesh.ich
ElevationFile = SimpleFace_Elev.ich
FaceIdFile = SimpleFace_faceIds.ich
Nlayers = 1
INTERP = FACE

Porosity

[Porosity]
Value = 1

Run Ichnos

To run this simulation we have to execute the following command under a terminal (e.g. windows powershell)

ichnos.exe -c SimpleFace_config.ini

This will produce an output similar to the following

G:\UCDAVIS\ichnos\ichnos_code\cmake-build-debug\ichnos.exe -c SImpleFace_config.ini
Ichnos started at Wed Aug 28 08:42:02 2024

Ichnos will run using 1 processors
Reading input data...
--> Configuration file: SImpleFace_config.ini
        Reading XYZ data from SimpleFace_vf.ini...
        Point Set Building time: 0.0001065
--> Velocity configuration file: SimpleFace_vf.ini
        Read Velocity in : 0.0006891
Tracing particles...
-10 % |5 of 42
-20 % |9 of 42
-30 % |13 of 42
-40 % |17 of 42
-50 % |22 of 42
-60 % |26 of 42
-70 % |30 of 42
-80 % |34 of 42
-90 % |38 of 42
Total simulation Time : 0.209663
Ichnos Finished at Wed Aug 28 08:42:03 2024


Process finished with exit code 0

Post process output

Ichnos prints the streamlines into a file with extension *.traj or *.h5 or both depending on the input output flags

In Matlab both of these files can be read using the same command

S = readICHNOStraj('output\SimpleFace_out_iter_0000.traj');

or

S = readICHNOStraj('output\SimpleFace_out_iter_0000.h5');

The readICHNOStraj matlab function is part of our gwtools repository.

The following figure shows the output of S structure

Eid and Sid are the ids that we setup in the particle file. p, v and t contain the positions the velocities and the time at each particle position. ER is the termination reason. See here for a list of termination reasons. Xend, Yend are the last point of the particle tracking, Time is the total travel time and Len is the length of the streamline

Next we will plot the streamlines overalayed on the mesh. Note that for large number of streamlines e.g. hundrends to thousands with many particle points, plotting can be very tricky in matlab.

figure()
clf
hold on
for ii = 1:9
    plot(ND(msh(ii,[1:4 1]),1), ND(msh(ii,[1:4 1]),2),'b')
end
for ii = 1:length(S)
    plot(S(ii,1).p(:,1),S(ii,1).p(:,2),'r')
end

---

---