Setting up input files for Face 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

This is going to be a long process!! :weary:

Prerequisite steps

The prerequisite steps are identical to the node case

Prepare the face velocities

In IWFM the horizontal velocity is defined on the element faces as flows. Therefore, getting the horizontal velocity is straightforward

For the vertical velocity of the top layer we will use the deep percolation which is also defined on the element face.

The vertical flow exchange between the model layers is defined on the element nodes and therefore have to allocate the flows to the elements

Vertical Flows

Read the vertical flows

DP_flows = h5read(GB_bud_info.Filename,...
    [GB_bud_info.Groups(2).Name GB_bud_info.Name GB_bud_info.Groups(2).Datasets(5).Name]);
for ii = 1:3
    VertFlows{ii,1} = h5read(GB_bud_info.Filename,...
    [GB_bud_info.Groups(1+ii).Name GB_bud_info.Name GB_bud_info.Groups(1+ii).Datasets(29).Name]);
end

Distribute the vertical flows between the layers to the elements

for ii = 1:length(VertFlows)
    VertFlowsElem{ii,1} = zeros(size(MSH,1), size(VertFlows{ii,1},2));
end
for ii = 1:size(ND_3310,1)
    [elems, JJ] = find(MSH == ii);
    nd_area = zeros(length(elems),1);
    for j = 1:length(elems)
        nd_area(j,1) = ElemNodesAreas(elems(j),JJ(j));
    end
    nd_area = nd_area./sum(nd_area);
    for j = 1:length(elems)
        for k = 1:length(VertFlows)
            VertFlowsElem{k,1}(elems(j),:) = VertFlowsElem{k,1}(elems(j),:) + VertFlows{k,1}(ii,:)*nd_area(j);
        end
    end
end

Join the deep percolation and the vertical flows in one variable and convert them from cuft/month to m/day. The positive deep percolation and vertical flow point downward. Hence, we negate the velocities

Fsy = scatteredInterpolant(ND_3310(:,1),ND_3310(:,2),GW_param.PN(:,1),'linear','nearest');
VertFlowsFinal{1,1} =  -((((ft2m^3)*DP_flows)./area_elem)'./ndays)'./Fsy(BC_elem(:,1),BC_elem(:,2));
for ii = 2:4
    Fsy = scatteredInterpolant(ND_3310(:,1),ND_3310(:,2),GW_param.PN(:,ii),'linear','nearest');
    VertFlowsFinal{ii,1} =  -((((ft2m^3)*VertFlowsElem{ii-1,1})./area_elem)'./ndays)'./Fsy(BC_elem(:,1),BC_elem(:,2));
end

Face flows

Read the face ids. The ids correspond to the element ids left and right of the face. They also indicate the flow direction which is always from face(:,1) -> face(:,2).

faces = h5read(GB_bud_info.Filename,...
    [GB_bud_info.Groups(1).Name GB_bud_info.Name GB_bud_info.Groups(1).Datasets(18).Name])';

The next task is to identify which faces each element consists of and the direction of flow for each face. This is a lengthy code which can be found in this gist.

The result of the above gist is the variable MSH_faces which contains the face ids that compose each element. If the id is negative then the flow direction is toward the element.

Next we have to identify the face centers and the direction of the faces. This is another involved loop which can be found here

The face_CC_NRM variable contains the coordinates of the faces centers and the normals.

Next calculate the face areas

face_areas = zeros(size(face_nodes,1),4);
for ii = 1:size(face_nodes,1)
    dst = pdist(ND_3310(face_nodes(ii,:),:));
    xx = [0 dst dst 0];
    for jj = 1:4
        dst_nda = Lay(face_nodes(ii,1),jj) - Lay(face_nodes(ii,1),jj+1);
        dst_ndb = Lay(face_nodes(ii,2),jj) - Lay(face_nodes(ii,2),jj+1);
        yy = [0 0 dst_nda dst_ndb];
        face_areas(ii,jj) = polyarea(xx,yy);
    end
end

Read the face flows for each layer

clear FaceFlows
for ii = 1:4
    FaceFlows{ii,1} = h5read(GB_bud_info.Filename,...
    [GB_bud_info.Groups(1+ii).Name GB_bud_info.Name GB_bud_info.Groups(1+ii).Datasets(9).Name]);
end

Convert face flows from cuft/month to m/day.

for ii = 1:length(FaceFlows)
    Fsy = scatteredInterpolant(ND_3310(:,1),ND_3310(:,2),GW_param.PN(:,ii),'linear','nearest');
    FaceFlows_m_days{ii,1} = ((((ft2m^3)*FaceFlows{ii,1})./face_areas(:,ii))'./ndays)'./Fsy(face_CC_NRM(:,1), face_CC_NRM(:,2));
end

Compute a Steady State velocity field

Prepare a steady state flow field for the last decade

Vface_av = zeros(length(face_areas), length(FaceFlows_m_days));
Vvert_av = zeros(length(area_elem), length(FaceFlows_m_days)+1);
for k = 1:length(FaceFlows_m_days)
    Vface_av(:,k) = mean(FaceFlows_m_days{k,1}(:,385:504),2);
    Vvert_av(:,k) = mean(VertFlowsFinal{k,1}(:,385:504),2);
end

Write the data into file

The velocity file contains the data in the following order:

First we list the face velocities starting from the first layer up to the n layer followed by the vertical velocities for layers 1 up to n+1 layer.

Reshape the data

VDATA = [reshape(Vface_av,size(Vface_av,1)*size(Vface_av,2),1); ...
    reshape(Vvert_av,size(Vvert_av,1)*size(Vvert_av,2),1)];

and print the data

mult = 0.000001;
fid = fopen(fullfile(input_files_path,'c2vsim_VelFACE_SS_000.ich'),'w');
fprintf(fid,'%.5f\n', VDATA/mult);
fclose(fid);

For the face velocity field we need one extra file with the face ids

fid = fopen(fullfile(input_files_path,'c2vsim_FaceIds.ich'),'w');
fprintf(fid,'%d %d %d %d\n', fliplr(MSH_faces)');
fclose(fid);