C2VSim MESH2D Element Multi processor

C2VSim multi processor simulation with multiple subdomains using MESH2D Element interpolation

First make sure to setup the C2VSim simulation (see here).

Set the paths

c2vsim_path = fullfile('path','to','c2vsimfg_version101')
input_files_path = fullfile('input_files01');

Domain

The domain input files are identical to single processing and they are discussed here.

Mesh data

Nodes

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);
ND_3310(:,3) = ND(:,3);

Mesh

MSH = readIWFM_Elements(fullfile(c2vsim_path,'Preprocessor',"C2VSimFG_Elements.dat"), 32537, 142);

Elevation

STRAT = readIWFM_Stratigraphy(fullfile(c2vsim_path,'Preprocessor',"C2VSimFG_Stratigraphy.dat"), size(ND,1), 4, 105);
ft2m = 0.3048;
Lay(:,1) = STRAT(:,2)*ft2m;
cnt = 3;
for ii = 2:5
    Lay(:,ii) = Lay(:,ii-1) - sum(STRAT(:,cnt:cnt+1),2)*ft2m;
    cnt = cnt+2;
end

Define subdomains

For this example we will use the same subdomains as in the CLOUD case. See here

Splitting the mesh

In multi processor simulation each subdomain will load a part of the entire mesh. Therefore the mesh indices will be modified so that they point to the nodes on the local subdomain.

Assign processor ids to the mesh elements. For that we will use another function from the gwtools repo that calculates the barycenters of the elements

First let's split the Nodes and the Elevations which is stragntforward. This will replace the 6th column from subregion index to processor index

BC_EL = calcIWFM_BaryCenters(ND_3310,MSH);
for ii = 1:length(proc_poly)
    idx = find(inpolygon(BC_EL(:,2), BC_EL(:,3), proc_poly(ii,1).x, proc_poly(ii,1).y));
    MSH(idx,6) = ii;
end

Plot for visual inspection

hold on
for ii = 1:length(proc_poly)
    idx = find(MSH(:,6) == ii);
    plot(BC_EL(idx,2), BC_EL(idx,3),'.')
end
plot(c2vsim_outline.X, c2vsim_outline.Y,'k')
axis equal off

The processor files will containt data from the extended polygons. For this simulation we set the buffer zone equal to 10,000.

The following lengthy code snippet identifies the mesh elements within the extended polygon and makes a unique list of nodes that are used by the elements. Then reassign the mesh ids so that they correspond to the local node file.

clear SplitData
for ii = 1:length(proc_poly)
    % Find the elements within the extended polygons
    idx = find(inpolygon(BC_EL(:,2), BC_EL(:,3), ext_poly(ii,1).x, ext_poly(ii,1).y));
    SplitData(ii,1).ei_id = idx;
    % Find the node ids that the mesh elements are made up
    id_nd = MSH(idx,2:5);
    id_nd = unique(id_nd);
    id_nd(id_nd == 0,:) = [];
    SplitData(ii,1).nd_id = id_nd;
    % Print the node file
    ND_proc = ND_3310(id_nd,:);
    Lay_proc = Lay(id_nd,:);
    SplitData(ii,1).ND = ND_proc;
    SplitData(ii,1).Lay = Lay_proc;

    % Rearrange the mesh indices
    MSH_proc = nan(length(idx),5);
    for iel = 1:length(idx)
        for ind = 1:4
            id = MSH(idx(iel), ind+1);
            if id == 0
                MSH_proc(iel, ind) = 0;
                continue
            end
            k = find(ND_proc(:,3) == id);
            MSH_proc(iel, ind) = k;
        end
        MSH_proc(iel, 5) = MSH(idx(iel), 6);
    end
    SplitData(ii,1).MSH = MSH_proc;
end

Print the Node, Elevation and Mesh files

for ii = 1:length(proc_poly)
    fid = fopen(fullfile(input_files_path,['c2vsim_Nodes_000' num2str(ii-1) '.ich']),'w');
    fprintf(fid,'%.3f %.3f\n', SplitData(ii,1).ND');
    fclose(fid);

    fid = fopen(fullfile(input_files_path,['c2vsim_ELEV_000' num2str(ii-1) '.ich']),'w');
    fprintf(fid,'%.3f %.3f %.3f %.3f %.3f\n', SplitData(ii,1).Lay');
    fclose(fid);

    fid = fopen(fullfile(input_files_path,['c2vsim_MSH_000' num2str(ii-1) '.ich']),'w');
    fprintf(fid,'%d %d %d %d %d\n', SplitData(ii,1).MSH');
    fclose(fid);
end

Velocity

Steady state

Run the code snippets of this section to calculate the Vx_av, Vy_av, Vz_av variables.

Each element is associated with one velocity therefore we will use the splitted element indices.

mult = 1000000;
for ii = 1:4
    Vx_av_proc = Vx_av(SplitData(ii,1).ei_id,:);
    Vy_av_proc = Vy_av(SplitData(ii,1).ei_id,:);
    Vz_av_proc = Vz_av(SplitData(ii,1).ei_id,:);
    
    nvel_p = size(Vx_av_proc,1)*size(Vx_av_proc,2);
    mult = 1000000;
    VEL_DATA = [
        reshape(Vx_av_proc, nvel_p, 1)*mult ... % VX
        reshape(Vy_av_proc, nvel_p, 1)*mult ... % VY
        reshape(Vz_av_proc, nvel_p, 1)*mult ... % VZ
    ];
    fid = fopen(fullfile(input_files_path,['c2vsim_VelELEM_SSM_000' num2str(ii-1) '.ich']),'w');
    fprintf(fid,'%.3f %.3f %.3f\n', VEL_DATA');
    fclose(fid);
end

Run the model using the main configuration file and the velocity configuration file.
Here we execute the program with 4 processors equal to the number of subdomains.

C:\"Program Files\Microsoft MPI"\Bin\mpiexec.exe -n 4 ichnos.exe -c c2vsim_MESH2D_ELEM_SSM_config.ini