-
Notifications
You must be signed in to change notification settings - Fork 3
C2VSim Face Mesh
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');
The preparation of the mesh data is described here
This is going to be a long process!! 😩
The prerequisite steps are identical to the node case
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
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
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
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
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);