CVHM2 element single proc - giorgk/ichnos GitHub Wiki
The Central Valley Hydrologic Model 2 (CVHM2) is an update to CVHM which is a modflow based integrated surface - groundwater model that simulates groundwater flow, stream aquifer interaction, conjuctive water use and more etc.
Here we use the flow output information to built a particle tracking simulation using Ichnos.
For the domain outline we read the discretization shapefile (CVHM2_DIS.shp) which can be downloaded from the CVHM2 repository. In a GIS software we convert the projection system to EPSG:3310, keep only the active cells and join them to form an outline which we can read it as follows:
cvhm_outline = shaperead('path/to/cvhm_outline_3310.shp');
We split the shapefile into polygons and print them to a file
[Xs, Ys] = polysplit(cvhm_outline.X, cvhm_outline.Y);
fid = fopen('cvhm2_domain.ich','w');
for ii = 1:length(Xs)
fprintf(fid, '%d %d\n',[length(Xs{ii,1}) ispolycw(Xs{ii,1},Ys{ii,1})]);
fprintf(fid, '%.3f %.3f\n', [Xs{ii,1}' Ys{ii,1}']');
end
fclose(fid);
The CVHM2 is defined on a grid that is rotated about 30 degrees. Therefore we have to convert the grid into an equivalent mesh. The discretization gis layer DIS shapefile is in a different projection system than 3310 therefore we first converted it to the EPSG:3310. The following code reads the converted shapefile and keeps only the active cells
dis_shp = shaperead('D:/BoxFolder/Box/CVHMII/gis_data/CVHM2_DIS_3310.shp');
idx = find([dis_shp.active]' == 1);
dis_shp_act = dis_shp(idx,:);
To convert the grid into a mesh we need to identify the unique locations of the nodes of the grid assign ids and use those ids to form the elements.
The following lengthy snippet does that and creates 2 variables the mesh nodes msh_nd and mesh elements msh_el
msh_nd = [];
msh_el = [];
for ii = 1:length(dis_shp_act)
% Here we know that the elements are all quadrilateral so we will loop
% from 1 to 4
el = []; % This will hold the element ids if the ii element
for j = 1:4
nd = [dis_shp_act(ii,1).X(j) dis_shp_act(ii,1).Y(j)];
if isempty(msh_nd)
msh_nd = [msh_nd;nd];
el = [el 1];
else
% calculate the distance with the nodes already in the list
dst = sqrt(sum((nd - msh_nd).^2, 2));
idx = find(dst < 1);
if isempty(idx)
% This is a new node and we append it in the list
msh_nd = [msh_nd;nd];
el = [el size(msh_nd,1)];
else
el = [el idx];
end
end
end
% last we append the indices of this element
msh_el = [msh_el;el];
end
In Matlab its not easy to plot a quadrilateral mesh but we can plot the mesh nodes only.
plot(msh_nd(:,1),msh_nd(:,2),'.')
The mesh elements must be oriented in a counter clockwise manner. The following code checks the orientation of the elements and reverse the indices order if needed.
for ii = 1:size(msh_el,1)
if ispolycw(msh_nd(msh_el(ii,:),1),msh_nd(msh_el(ii,:),2))
msh_el(ii,:) = msh_el(ii,[1 4 3 2]);
end
end
Typicaly the elevation data in CVHM2 are defined in the DIS file. The DIS file does not contain projection information therefore we will use the provided DIS shapefile which also contain the elevation information as well as the top and bottom active layer in each cell. The following snippet extracts cell centroids the elevations and the active layer flag
top_lay_act = [dis_shp.up_act]';
bot_lay_act = [dis_shp.low_act]';
Lay_XY = zeros(length(dis_shp),2);
for ii = 1:length(dis_shp)
Lay_XY(ii,1) = mean(dis_shp(ii).X(1:4));
Lay_XY(ii,2) = mean(dis_shp(ii).Y(1:4));
end
Lay_elev = zeros(length(dis_shp),14);
for ii = 1:13
Lay_elev(:,ii) = [dis_shp.(['lay' num2str(ii) 't_m'])]';
end
Lay_elev(:,14) = [dis_shp.lay13b_m]';
Next we will gather the top and bottom elevation into a separate matrix from the appropriate layer. The first line identifies the cells that have valid layer information, which are all the active layers, and the next two lines extract the top and bottom layer values.
id_lay = find(top_lay_act ~=0 & bot_lay_act ~= 0);
top_lay_elev = Lay_elev(sub2ind(size(Lay_elev), id_lay, top_lay_act(id_lay)));
bot_lay_elev = Lay_elev(sub2ind(size(Lay_elev), id_lay, bot_lay_act(id_lay)));
Then we can create interpolants and interpolate the cell center elevations to the mesh nodes.
Ftop = scatteredInterpolant(Lay_XY(id_lay,1),Lay_XY(id_lay,2),top_lay_elev,'linear','nearest');
Fbot = scatteredInterpolant(Lay_XY(id_lay,1),Lay_XY(id_lay,2),bot_lay_elev,'linear','nearest');
msh_top = Ftop(msh_nd(:,1),msh_nd(:,2));
msh_bot = Fbot(msh_nd(:,1),msh_nd(:,2));
For the top and bottom elevation we will use the MESH2D option. We generated previously the mesh nodes and mesh elements. Here we append the top and bottom elevation and print into one file as the top and bottom are sharing the same nodes.
fid = fopen('cvhm2_top_bot.ich','w');
fprintf(fid,'MESH2D\n');
fprintf(fid, '%d %d %.2f\n',[size(msh_nd,1) size(msh_el,1) 2500]);
fprintf(fid, '%.3f %.3f %.3f %.3f\n', [msh_nd msh_top msh_bot]');
fprintf(fid, '%d %d %d %d\n',msh_el');
fclose(fid);
The radius parameter was set 2500 because we know that all elements are rectangles with 1 mile side. Therefore within a 2.5km radius from any point within the domain there will always be at least one element within the search distance.
In this example will use the MESH2D support structure with ELEMENT interpolation. To this end we need to prepare three files that define the mesh.
We will use the same mesh as in the top and bottom elevations
fid = fopen("cvhm2_nodes.ich",'w');
fprintf(fid, '%.3f %.3f\n', msh_nd');
fclose(fid);
Similarly we already have the mesh data available from previous snippets
fid = fopen("cvhm2_elements.ich",'w');
fprintf(fid, '%d %d %d %d 0\n', msh_el');
fclose(fid);
While in modflow files the elevation is defined in the DIS file, the discretization shapefile that we have used previously already contains the layer elevation at the cell centers. In a similar way we can interpolate to the mesh nodes
msh_elev = zeros(size(msh_nd,1),14);
for ii = 1:size(Lay_elev,2)
Felev = scatteredInterpolant(Lay_XY(:,1),Lay_XY(:,2), Lay_elev(:,ii),'linear','nearest');
msh_elev(:,ii) = Felev(msh_nd(:,1),msh_nd(:,2));
end
It's always a good idea to check if there is any layer crossing. This can happend during interpolation when the layers are very thin. A quick way to check for errors is the following
any(diff(msh_elev,1,2)>0)
This should be all false.
Now we can print the mesh elevation file
fid = fopen("cvhm2_elevations.ich",'w');
fprintf(fid, [repmat('%.3f ',1,13) '%.3f\n'], msh_elev');
fclose(fid);
Next we will print the velocity. In Modflow the velocity is printed as the volumetric flow that cross the grid cells. However to keep this document short we will convert the volumetric flows to velocity by dividing with the cross area and assign the velocity to the center of the cell grid and we will use the ELEMENT interpolation.
The Modflow velocity is printed into the binary Cell by Cell output file. For CVHM2 the *.cbc is about 10 GB and can be generated after the model is run or it can be downloaded from the CVHM repository (Output.zip). To bring the cbc data into Matlab we have prepared a script that reads the entire CBC and then we saved each component of the CBC flow into a separate *.mat file so that we can easily bring the data into matlab. The following code reads the entire cbc and saves the FLOW velocity terms (Note that the first line is going to take a few hours for the CVHM2)
CBC = readCVHMcbc('path\to\cbcf.out');
% Save each flow term
for ii = 1:length(CBC)
FLOWRIGHTFACE{ii,1} = CBC{ii,1}{3,2};
end
save('path\to\FLOWRIGHTFACE.mat','FLOWRIGHTFACE','-v7.3')
for ii = 1:length(CBC)
FLOWFRONTFACE{ii,1} = CBC{ii,1}{4,2};
end
save('path\to\FLOWFRONTFACE.mat','FLOWFRONTFACE','-v7.3')
for ii = 1:length(CBC)
FLOWLOWERFACE{ii,1} = CBC{ii,1}{5,2};
end
save('path\to\FLOWLOWERFACE.mat','FLOWLOWERFACE','-v7.3')
Once we have those mat files we can retrieve the velocity information very quickly as:
load('path\to\FLOWFRONTFACE.mat')
load('path\to\FLOWRIGHTFACE.mat')
load('path\to\FLOWLOWERFACE.mat')
First let's reshape the velocity data so that the rows correspond to grid cells and the columns to time steps
% extract auxiliary variables
RC_act = [[dis_shp_act.ROW]' [dis_shp_act.COLUMN_]'];
RC_ind = sub2ind([441 98], [dis_shp_act.ROW]', [dis_shp_act.COLUMN_]');
% Initialize the variables
for j = 1:13
QX{j,1} = zeros(size(RC_act,1),length(FLOWRIGHTFACE));
QY{j,1} = zeros(size(RC_act,1),length(FLOWRIGHTFACE));
QZ{j,1} = zeros(size(RC_act,1),length(FLOWRIGHTFACE));
end
for ii = 1:length(FLOWRIGHTFACE)
for j = 1:13
tmpx = FLOWRIGHTFACE{ii,1}(:,:,j);
QX{j,1}(:,ii) = tmpx(RC_ind);
tmpy = FLOWFRONTFACE{ii,1}(:,:,j);
QY{j,1}(:,ii) = tmpy(RC_ind);
tmpz = FLOWFRONTFACE{ii,1}(:,:,j);
QZ{j,1}(:,ii) = tmpz(RC_ind);
end
end
The units of QX QY QZ are m^3/day.
To convert them to m/day we calculate the thickness of each layer.
layThick = zeros(length(RC_ind),13);
for ii = 1:13
if ii == 13
layThick(:,ii) = [dis_shp_act.lay13t_m]' - [dis_shp_act.lay13b_m]';
else
layThick(:,ii) = [dis_shp_act.(['lay' num2str(ii) 't_m'])]' - [dis_shp_act.(['lay' num2str(ii+1) 't_m'])]';
end
end
Now we can calculate the velocity by dividing with the face area. The For the VX, VY the area is 1 mile times the thickness and for the vertical velocity the face area is 1 square mile
for ii = 1:13
VX{ii,1} = QX{ii,1}./(layThick(:,ii)*1609.34);
VY{ii,1} = QY{ii,1}./(layThick(:,ii)*1609.34);
VZ{ii,1} = QZ{ii,1}./(1609.34*1609.34);
end
The velocities that we just computed have direction perpendicular to cell face. However the cells are rotated with respect the horizontal axis.
The VX VY (black) components of velocity results in the velocity vector Vxy (orange). With Ichnos we have to calculate the components Vx` and VY` (green)
First we calculate the grid rotation. The first two lines below identify the left and right lower cells of the Modflow grid.
Then we create a vector v1 from the second node of the left lower cell to the first node of the right lower cell and calculate the anlge between this vector and the horizontal axis v2. The calculated rotation is about 30.0652.
id1 = find([dis_shp.ROW]' == 441 & [dis_shp.COLUMN_]' == 1);
id2 = find([dis_shp.ROW]' == 441 & [dis_shp.COLUMN_]' == 98);
v1 = [dis_shp(id2,1).X(1) - dis_shp(id1,1).X(2) dis_shp(id2,1).Y(1) - dis_shp(id1,1).Y(2)];
v2 = [100000 0];
rot = acosd(sum(v1.*v2)/(sqrt(sum(v1.^2))*sqrt(sum(v2.^2))));
Now that we know the rotation angle we can calculate the unit vectors of the rotated grid projected on the cartesian axis. bx by are the cartesian unit vectors and bxrot byrot are the unit vectors of the rotated system
bx = [1 0];
by = [0 1];
bxrot(1) = bx(1)*cosd(rot) - bx(2)*sind(rot);
bxrot(2) = bx(1)*sind(rot) + bx(2)*cosd(rot);
byrot(1) = by(1)*cosd(rot) - by(2)*sind(rot);
byrot(2) = by(1)*sind(rot) + by(2)*cosd(rot);
We can plot those vectors:
figure()
clf
plot([0 bx(1)], [0 bx(2)],'b')
hold on
plot([0 bxrot(1)], [0 bxrot(2)],'r')
plot([0 by(1)], [0 by(2)],'b')
plot([0 byrot(1)], [0 byrot(2)],'r')
axis equal
grid on
To understand the calculations needed for the rotation we will pick the velocity of one cell e.g. 10000 and layer 9
id = 10000;
lid = 9;
vx = VX{lid,1}(id,1)
vy = VY{lid,1}(id,1)
vz = VZ{lid,1}(id,1)
The values vx=-0.0252, vy=-0.0104, vz=-3.758e-4 show the magnitude of the velocity aligned with the rotated grid.
To plot them we need to express them using the unit rotated vectors where vx_m vy_m are the velocity vectors with magnitude and direction defined on the rotated gridded system and a is the resultant vector.
vx_m = bxrot*vx;
vy_m = byrot*vy;
a = vx_m + vy_m;
Last we project the resultant vector a onto the unit cartesian vectors to calculate the x,y vecolity components to pass to ichnos
a1 = bx*dot(a,bx)/dot(bx,bx);
a2 = by*dot(a,by)/dot(by,by);
Let's plot the above vectors:
figure()
clf
plot([0 vx_m(1)], [0 vx_m(2)],'r')
hold on
plot([0 vy_m(1)], [0 vy_m(2)],'r')
plot([0 a(1)], [0 a(2)],'--r')
plot([0 a1(1)],[0 a1(2)],'g')
plot([0 a2(1)],[0 a2(2)],'g')
grid on
axis equal
The red lines correspond to the Modflow velocities defined as perpendicular to the cell faces. The green vectors are defined so that their resultant velocity vector has same direction and magnitute as the original Modflow velocity.
Next we are going to repeat the above calculations in vectorized manner for all cells.
clear VXproj VYproj VZproj
for ii = 1:13
vx_m1 = bxrot(1)*VX{ii,1};
vx_m2 = bxrot(2)*VX{ii,1};
vy_m1 = byrot(1)*VY{ii,1};
vy_m2 = byrot(2)*VY{ii,1};
VXproj{ii,1} = vx_m1 + vy_m1;
VYproj{ii,1} = vx_m2 + vy_m2;
VZproj{ii,1} = -VZ{ii,1};
end
For the vertical velocity we reverse the direction because in Modflow positive is the flow between a layer and the layer below
For this demonstration we will average the velocities of the last decade and print them as steady state.
VAVxyz = [];
for ii = 1:13
VAVxyz = [VAVxyz; ...
[mean(VXproj{ii,1}(:,end-119:end),2) ...
mean(VYproj{ii,1}(:,end-119:end),2) ...
mean(VZproj{ii,1}(:,end-119:end),2)]];
end
- Print steady state as ASCII
fid = fopen('path\to\cvhm2_velSS_000.vel','w');
fprintf(fid,'%.5e %.5e %.5e\n',VAVxyz');
fclose(fid);
- Print steady state as HDF5
fname = fullfile('G:\UCDAVIS\IchnosPaper\CVHM2\cvhm2_velSS_000.h5');
h5create(fname,'/VXYZ',[size(VAVxyz,1) 3], 'Datatype','single');
h5write(fname, '/VXYZ', VAVxyz);
For completness we provide the code snippets to prepare and print the transient state velocity which is very similar with the difference being that each component is printed into a different file or dataset.
VXtrs = [];
VYtrs = [];
VZtrs = [];
for ii = 1:13
VXtrs = [VXtrs;VXproj{ii,1}];
VYtrs = [VYtrs;VXproj{ii,1}];
VZtrs = [VZtrs;VXproj{ii,1}];
end
- Print transient state ASCII
fid = fopen('G:\UCDAVIS\IchnosPaper\CVHM2\cvhm2_velTR_VX_000.vel','w');
fprintf(fid,[repmat('%.5e ',1,701) '%.5e\n'], VXtrs');
fclose(fid);
fid = fopen('G:\UCDAVIS\IchnosPaper\CVHM2\cvhm2_velTR_VY_000.vel','w');
fprintf(fid,[repmat('%.5e ',1,701) '%.5e\n'], VYtrs');
fclose(fid);
fid = fopen('G:\UCDAVIS\IchnosPaper\CVHM2\cvhm2_velTR_VZ_000.vel','w');
fprintf(fid,[repmat('%.5e ',1,701) '%.5e\n'], VZtrs');
fclose(fid);
- Print transient state HDF5
fnameTR = fullfile('G:\UCDAVIS\IchnosPaper\CVHM2\cvhm2_velTR_000.h5');
h5create(fnameTR,'/VX',[size(VXtrs,1) size(VXtrs,2)], 'Datatype','single');
h5write(fnameTR, '/VX', VXtrs);
h5create(fnameTR,'/VY',[size(VYtrs,1) size(VYtrs,2)], 'Datatype','single');
h5write(fnameTR, '/VY', VYtrs);
h5create(fnameTR,'/VZ',[size(VZtrs,1) size(VZtrs,2)], 'Datatype','single');
h5write(fnameTR, '/VZ', VZtrs);
In the following we list the main configuration options
[Velocity]
XYZType = MESH2D
Type = DETRM
ConfigFile = path\to\velocityConfig.ini
[Domain]
Outline = path\to\cvhm2_domain.ich
TopFile = path\to\cvhm2_top_bot.ich
BottomFile =
ProcessorPolys =
[StepConfig]
Method = RK4
Direction = 1
StepSize = 50
StepSizeTime = 100000
nSteps = 5
nStepsTime = 0
minExitStepSize = 1
MaxIterationsPerStreamline = 3000
MaxProcessorExchanges = 50
AgeLimit = 100
Stuckiter = 10
AttractFile =
AttractRadius = 30
[InputOutput]
PartilceFile = path\to\cvhm2_particles.ich
WellFile =
OutputFile = path\to\out_prefix_
PrintH5 = 0
PrintASCII = 1
ParticlesInParallel = 5000
GatherOneFile = 0
[Other]
Nrealizations = 1
nThreads = 1
RunAsThread = 0
Version = 0.5.06
OutFreq = 1
Make sure that the version number is the same as the output of the ichnos.exe -v
Here we list the relevant options for hte steady state simulation
[Velocity]
Prefix = path\to\cvhm2_velSS_
LeadingZeros = 3
Suffix = .vel
Type = DETRM
[MESH2D]
NodeFile = path\to\cvhm2_nodes.ich
Meshfile = path\to\cvhm2_elements.ich
ElevationFile = path\to\cvhm2_elevations.ich
FaceIdFile =
Nlayers = 13
INTERP = ELEMENT
[Porosity]
Value = 0.1
To run particle tracking we need initial positions for particles. Here we will generate 1000 particles scattered randomly in Central Valley while their starting point will be 5 m below the top boundary.
cv_pp = polyshape(cvhm_outline.X, cvhm_outline.Y);
cv_pp_buff = cv_pp.polybuffer(-5000);
particles = zeros(1000,3);
for ii = 1:size(particles,1)
while 1
xrnd = cvhm_outline.BoundingBox(1,1) + (cvhm_outline.BoundingBox(2,1)-cvhm_outline.BoundingBox(1,1))*rand;
yrnd = cvhm_outline.BoundingBox(1,2) + (cvhm_outline.BoundingBox(2,2)-cvhm_outline.BoundingBox(1,2))*rand;
if cv_pp_buff.isinterior(xrnd,yrnd)
particles(ii,1:2) = [xrnd yrnd];
break;
end
end
end
particles(:,3) = Ftop(particles(:,1), particles(:,2))-5;
Print the particle file
fid = fopen('cvhm2_particles.ich','w');
fprintf(fid, '%d %d %.3f %.3f %.3f 0\n',[ones(1000,1) [1:1000]' particles]');
fclose(fid);
To run the simulation make sure that all input files are in the same directory or the correct locations are specified. Note that you can use relative locations with respect the directory that is returned after running pwd on powershell or terminal.
The simulation can be run as
> path\to\ichnos.exe -c cvhm2_MESH2D_ELEM_SS_config.iniThe output file can be ASCII or HDF5. Within Matlab both can be read with the following command
S = readICHNOStraj('path\to\cvhm2_MESH2D_SS_ELEM_iter_0000.traj');
To visualize the streamlines we can do the following
clf
plot(cv_pp)
hold on
for ii = 1:1000
plot(S(ii,1).p(:,1), S(ii,1).p(:,2),'b')
end
axis equal off
