parallel_grid/
parallel_grid.rs

1use mpi::{
2    collective::SystemOperation, environment::Universe, topology::Communicator,
3    traits::CommunicatorCollectives,
4};
5use ndelement::types::ReferenceCellType;
6use ndgrid::{
7    grid::local_grid::SingleElementGridBuilder,
8    traits::{Builder, Entity, Grid, ParallelBuilder, ParallelGrid},
9    types::{GraphPartitioner, Ownership},
10};
11
12/// Creating a distributed parallel grid
13fn main() {
14    // The SingleElementGridBuilder is used to create the mesh
15    let mut b = SingleElementGridBuilder::<f64>::new(2, (ReferenceCellType::Quadrilateral, 1));
16
17    let universe: Universe = mpi::initialize().unwrap();
18    let comm = universe.world();
19    let rank = comm.rank();
20
21    // Add points and cells to the builder on process 0
22    let n = 10;
23    let grid = if rank == 0 {
24        let mut i = 0;
25        for y in 0..n {
26            for x in 0..n {
27                b.add_point(i, &[x as f64 / (n - 1) as f64, y as f64 / (n - 1) as f64]);
28                i += 1;
29            }
30        }
31
32        let mut i = 0;
33        for y in 0..n - 1 {
34            for x in 0..n - 1 {
35                let sw = n * y + x;
36                b.add_cell(i, &[sw, sw + 1, sw + n, sw + n + 1]);
37                i += 1;
38            }
39        }
40
41        // Distribute the grid
42        // In this example, we use Scotch to partition the grid into pieces to be handles by each process
43        b.create_parallel_grid_root(&comm, GraphPartitioner::Scotch)
44    } else {
45        // receice the grid
46        b.create_parallel_grid(&comm, 0)
47    };
48
49    // Check that owned cells are sorted ahead of ghost cells
50    let cell_count_owned = grid
51        .local_grid()
52        .entity_iter(ReferenceCellType::Quadrilateral)
53        .filter(|entity| entity.is_owned())
54        .count();
55
56    // Now check that the first `cell_count_owned` entities are actually owned.
57    for cell in grid
58        .local_grid()
59        .entity_iter(ReferenceCellType::Quadrilateral)
60        .take(cell_count_owned)
61    {
62        assert!(cell.is_owned())
63    }
64
65    // Now make sure that the indices of the global cells are in consecutive order
66    let mut cell_global_count = grid.cell_layout().local_range().0;
67
68    for cell in grid
69        .local_grid()
70        .entity_iter(ReferenceCellType::Quadrilateral)
71        .take(cell_count_owned)
72    {
73        assert_eq!(cell.global_index(), cell_global_count);
74        cell_global_count += 1;
75    }
76
77    // Get the global indices
78    let global_vertices = grid
79        .local_grid()
80        .entity_iter(ReferenceCellType::Point)
81        .filter(|e| matches!(e.ownership(), Ownership::Owned))
82        .map(|e| e.global_index())
83        .collect::<Vec<_>>();
84
85    let nvertices = global_vertices.len();
86
87    let global_cells = grid
88        .local_grid()
89        .entity_iter(ReferenceCellType::Quadrilateral)
90        .filter(|e| matches!(e.ownership(), Ownership::Owned))
91        .map(|e| e.global_index())
92        .collect::<Vec<_>>();
93
94    let ncells = global_cells.len();
95
96    let mut total_cells: usize = 0;
97    let mut total_vertices: usize = 0;
98
99    comm.all_reduce_into(&ncells, &mut total_cells, SystemOperation::sum());
100    comm.all_reduce_into(&nvertices, &mut total_vertices, SystemOperation::sum());
101
102    // Check that the total number of cells and vertices are correct
103    assert_eq!(total_cells, (n - 1) * (n - 1));
104    assert_eq!(total_vertices, n * n);
105}
parallel_grid/parallel_grid.rs

parallel_grid/
parallel_grid.rs
