Setting up input files for Node based interpolation with Mesh2D structure

The first step is to setup C2VSim to generate the required data. This is described here

Let's now set some common paths

c2vsim_path = fullfile('path','to','C2VsimV1','c2vsim-working');
input_files_path = fullfile('path','to','input','files');

Mesh data

The preparation of the mesh data is described here

Domain data

See here

Velocity data

In C2VSim the velocity is not defined on the nodes. However, for demonstration purposes we will interpolate/extrapolate the velocity that is defined in the center of the elements to the nodes.

Read the velocity

VelOut = readIWFM_Velocity(fullfile(c2vsim_path,'Results','C2VSimFG_GW_VELOUTFL.out'), 32537, 4, 1000);

Convert the coordinates of the element barycenters from feet to meter and to 3310

[lat,lon] = projinv(projcrs(26910),VelOut.ND(:,1)*0.3048, VelOut.ND(:,2)*0.3048);
[BCX_3310, BCY_3310] = projfwd(projcrs(3310),lat, lon);

Create an average velocity flow field for the last decade of the simulation

for ii = 1:size(VelOut.VX,3)
    Vx_av(:,ii) = mean(VelOut.VX(:,385:504,ii),2);
    Vy_av(:,ii) = mean(VelOut.VY(:,385:504,ii),2);
    Vz_av(:,ii) = mean(VelOut.VZ(:,385:504,ii),2);
end

Convert the velocity units from FT/MON to m/day.

Vx_av = Vx_av .* (0.3048/30.437);
Vy_av = Vy_av .* (0.3048/30.437);
Vz_av = Vz_av .* (0.3048/30.437);

Read the element nodes and convert the coordinates to 3310

ND = readIWFM_Nodes(fullfile(c2vsim_path,'Preprocessor',"C2VSimFG_Nodes.dat"));
[lat,lon] = projinv(projcrs(26910),ND(:,1), ND(:,2));
[ND_3310(:,1), ND_3310(:,2)] = projfwd(projcrs(3310),lat, lon);

For each layer interpolate/extrapolate the velocity to the mesh nodes.

Vx_nd = zeros(size(ND_3310,1),size(Vx_av,2));
Vy_nd = zeros(size(ND_3310,1),size(Vy_av,2));
Vz_nd = zeros(size(ND_3310,1),size(Vz_av,2));
for ii = 1:size(Vx_av,2)
    Fx = scatteredInterpolant(BCX_3310, BCY_3310, Vx_av(:,ii),'linear','nearest');
    Fy = scatteredInterpolant(BCX_3310, BCY_3310, Vy_av(:,ii),'linear','nearest');
    Fz = scatteredInterpolant(BCX_3310, BCY_3310, Vz_av(:,ii),'linear','nearest');
    Vx_nd(:,ii) = Fx(ND_3310(:,1), ND_3310(:,2));
    Vy_nd(:,ii) = Fy(ND_3310(:,1), ND_3310(:,2));
    Vz_nd(:,ii) = Fz(ND_3310(:,1), ND_3310(:,2));
end

Assign the velocities of the first and last layer to the first and last layer of element nodes. For the intermediate layers we will average the velocities of the layer above and below.

Vx_nd_final = zeros(size(ND_3310,1),size(Vx_av,2)+1);
Vy_nd_final = zeros(size(ND_3310,1),size(Vy_av,2)+1);
Vz_nd_final = zeros(size(ND_3310,1),size(Vz_av,2)+1);
Vx_nd_final(:,1) = Vx_nd(:,1);
Vx_nd_final(:,end) = Vx_nd(:,end);
Vy_nd_final(:,1) = Vy_nd(:,1);
Vy_nd_final(:,end) = Vy_nd(:,end);
Vz_nd_final(:,1) = Vz_nd(:,1);
Vz_nd_final(:,end) = Vz_nd(:,end);
for ii = 2:size(Vx_av,2)
    Vx_nd_final(:,ii) = (Vx_nd(:,ii-1)+Vx_nd(:,ii))/2;
    Vy_nd_final(:,ii) = (Vy_nd(:,ii-1)+Vy_nd(:,ii))/2;
    Vz_nd_final(:,ii) = (Vz_nd(:,ii-1)+Vz_nd(:,ii))/2;
end

Reshape the velocity data and write them to files

VXYZ = [reshape(Vx_nd_final, size(Vx_nd_final,1)*size(Vx_nd_final,2),1)...
    reshape(Vy_nd_final, size(Vy_nd_final,1)*size(Vy_nd_final,2),1)...
    reshape(Vz_nd_final, size(Vz_nd_final,1)*size(Vz_nd_final,2),1)
    ];

Print the velocity to file

mult = 0.000001;
fid = fopen(fullfile(input_files_path,'c2vsim_VelNODE_SS_0000.vel'),'w');
fprintf(fid,'%.5f %.5f %.5f\n', [VXYZ/mult]');
fclose(fid);