Skip to content

Commit dca8bfd

Browse files
cbritopachecocbritopacheco
committed
Update
1 parent fcb1a7c commit dca8bfd

File tree

7 files changed

+177
-133
lines changed

7 files changed

+177
-133
lines changed

‎examples/PETSc/MPI_Poisson.cpp‎

Lines changed: 1 addition & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -53,37 +53,9 @@ int main(int argc, char** argv)
5353
{
5454
PETSc::Variational::GridFunction gf(vh);
5555
gf = Cos(F::x);
56-
gf.save("Poisson." + std::to_string(world.rank()) + ".gf");
56+
// gf.save("Poisson." + std::to_string(world.rank()) + ".gf");
5757
}
5858

59-
// Vec v;
60-
// VecCreate(PETSC_COMM_WORLD, &v);
61-
62-
// size_t sz;
63-
// if (world.rank() == 0)
64-
// {
65-
// sz = 7;
66-
// }
67-
// else
68-
// {
69-
// sz = 8;
70-
// }
71-
// VecSetSizes(v, sz, 10);
72-
// VecSetFromOptions(v);
73-
74-
// std::vector<int> ghosts;
75-
// if (world.rank() == 0)
76-
// {
77-
// ghosts = { 2, 3 };
78-
// }
79-
// else
80-
// {
81-
// ghosts = { 4, 5, 6 };
82-
// }
83-
84-
// VecMPISetGhost(v, ghosts.size(), ghosts.data());
85-
86-
8759
PetscFinalize();
8860
}
8961

‎src/Rodin/MPI/Variational/P1/P1.h‎

Lines changed: 134 additions & 37 deletions
Original file line numberDiff line numberDiff line change
@@ -20,7 +20,13 @@ namespace Rodin::Variational
2020
Geometry::Mesh<Context::MPI>, P1<Range, Geometry::Mesh<Context::MPI>>>
2121
{
2222
public:
23-
using Bimap =
23+
struct GhostBimap
24+
{
25+
std::vector<Index> left;
26+
FlatMap<Index, Index> right;
27+
};
28+
29+
using IndexBimap =
2430
boost::bimap<
2531
boost::bimaps::vector_of<Index>,
2632
boost::bimaps::unordered_set_of<Index>>;
@@ -133,7 +139,8 @@ namespace Rodin::Variational
133139
: m_mesh(mesh),
134140
m_fes(mesh.getShard())
135141
{
136-
const auto& comm = mesh.getContext().getCommunicator();
142+
const auto& ctx = mesh.getContext();
143+
const auto& comm = ctx.getCommunicator();
137144
const auto& shard = mesh.getShard();
138145
const auto& halo = shard.getHalo(0);
139146
const auto& owner = shard.getOwner(0);
@@ -148,63 +155,147 @@ namespace Rodin::Variational
148155
{
149156
if (shard.isOwned(0, i))
150157
{
158+
const Index id = mesh.getGlobalIndex(0, i);
159+
const auto& dofs = m_fes.getDOFs(0, i);
160+
assert(dofs.size() == 1);
161+
const Index local = dofs[0];
162+
const Index global = dofIdx + m_offset;
163+
const auto [it, inserted] = m_local_to_global.right.insert({ global, local });
164+
assert(inserted);
165+
dofIdx++;
166+
151167
for (const Index& peer : halo.at(i))
152168
{
153169
assert(comm.rank() >= 0);
154170
assert(peer != static_cast<Index>(comm.rank()));
155-
const Index& global = mesh.getGlobalIndex(0, i);
156-
push[peer].push_back({ global, dofIdx + m_offset });
171+
push[peer].push_back({ id, global });
157172
}
158-
m_loc2glob.right.insert({ dofIdx + m_offset, i });
159-
dofIdx++;
160173
}
161174
else
162175
{
163176
pull.try_emplace(owner.at(i));
164177
}
165178
}
166179

167-
assert(m_owned == dofIdx);
168-
169-
std::vector<boost::mpi::request> requests;
180+
std::vector<boost::mpi::request> irecv;
170181
for (auto& [owner, requested] : pull)
171-
requests.push_back(comm.irecv(owner, 0, pull[owner]));
182+
irecv.push_back(comm.irecv(owner, 0, pull[owner]));
183+
184+
std::vector<boost::mpi::request> isend;
172185
for (const auto& [peer, requested] : push)
173-
requests.push_back(comm.isend(peer, 0, push[peer]));
174-
boost::mpi::wait_all(requests.begin(), requests.end());
186+
isend.push_back(comm.isend(peer, 0, push[peer]));
187+
188+
boost::mpi::wait_all(isend.begin(), isend.end());
189+
190+
boost::mpi::wait_all(irecv.begin(), irecv.end());
175191

176192
for (const auto& [owner, requested] : pull)
177193
{
178-
for (size_t i = 0; i < requested.size(); ++i)
194+
for (const auto& [id, global] : requested)
179195
{
180-
const auto local = mesh.getLocalIndex(0, requested[i].first);
181-
assert(local);
182-
const Index& global = requested[i].second;
183-
m_loc2glob.right.insert({ global, *local });
196+
const auto i = mesh.getLocalIndex(0, id);
197+
assert(i);
198+
const auto& dofs = m_fes.getDOFs(0, *i);
199+
assert(dofs.size() == 1);
200+
const auto [it, inserted] = m_local_to_global.right.insert({ global, dofs[0] });
201+
assert(inserted);
184202
}
185203
}
186204

187-
for (int r = 0; r < comm.size(); ++r) {
188-
if (comm.rank() == r) {
189-
std::cout << "=== rank " << r << " ===\n";
190-
for (auto it = m_loc2glob.left.begin(); it != m_loc2glob.left.end(); ++it) {
191-
std::cout
192-
<< " local=" << it->first
193-
<< " -> global=" << it->second
194-
<< "\n";
195-
}
196-
std::cout << std::flush;
197-
}
198-
comm.barrier(); // wait so next rank prints cleanly
199-
}
200-
201-
205+
const size_t begin = m_offset, end = m_offset + m_owned;
206+
Index offset = 0;
207+
for (size_t i = 0; i < m_fes.getSize(); ++i)
208+
{
209+
const Index global = this->getGlobalIndex(i);
210+
if (global < begin || end <= global)
211+
{
212+
auto [it, inserted] = m_ghosts.right.insert({ global, offset });
213+
assert(inserted);
214+
m_ghosts.left.push_back(offset);
215+
offset++;
216+
}
217+
}
202218
}
203219

204220
P1(const MeshType& mesh, size_t vdim)
205221
: m_mesh(mesh),
206222
m_fes(mesh.getShard(), vdim)
207-
{}
223+
{
224+
static thread_local std::vector<Index> s_send;
225+
226+
const auto& ctx = mesh.getContext();
227+
const auto& comm = ctx.getCommunicator();
228+
const auto& shard = mesh.getShard();
229+
const auto& halo = shard.getHalo(0);
230+
const auto& owner = shard.getOwner(0);
231+
232+
m_owned = halo.size();
233+
const size_t inclusive = boost::mpi::scan(comm, m_owned, std::plus<size_t>());
234+
m_offset = inclusive - m_owned;
235+
236+
FlatMap<Index, std::vector<std::pair<Index, std::vector<Index>>>> push, pull;
237+
Index dofIdx = 0;
238+
for (size_t i = 0; i < shard.getVertexCount(); ++i)
239+
{
240+
if (shard.isOwned(0, i))
241+
{
242+
const Index id = mesh.getGlobalIndex(0, i);
243+
const auto& dofs = m_fes.getDOFs(0, i);
244+
245+
s_send.clear();
246+
for (const Index& local : dofs)
247+
{
248+
const Index global = dofIdx + m_offset;
249+
s_send.push_back(global);
250+
251+
const auto [it, inserted] = m_local_to_global.right.insert({ global, local });
252+
assert(inserted);
253+
254+
dofIdx++;
255+
}
256+
257+
for (const Index& peer : halo.at(i))
258+
{
259+
assert(comm.rank() >= 0);
260+
assert(peer != static_cast<Index>(comm.rank()));
261+
push[peer].push_back({ id, s_send });
262+
}
263+
}
264+
else
265+
{
266+
pull.try_emplace(owner.at(i));
267+
}
268+
}
269+
270+
std::vector<boost::mpi::request> irecv;
271+
for (auto& [owner, requested] : pull)
272+
irecv.push_back(comm.irecv(owner, 0, pull[owner]));
273+
274+
std::vector<boost::mpi::request> isend;
275+
for (const auto& [peer, requested] : push)
276+
isend.push_back(comm.isend(peer, 0, push[peer]));
277+
278+
boost::mpi::wait_all(isend.begin(), isend.end());
279+
280+
boost::mpi::wait_all(irecv.begin(), irecv.end());
281+
282+
for (const auto& [owner, requested] : pull)
283+
{
284+
for (const auto& [id, global] : requested)
285+
{
286+
const auto i = mesh.getLocalIndex(0, id);
287+
assert(i);
288+
const auto& dofs = m_fes.getDOFs(0, *i);
289+
assert(dofs.size() >= 0);
290+
assert(dofs.size() == global.size());
291+
for (size_t k = 0; k < global.size(); k++)
292+
{
293+
const auto [it, inserted] = m_local_to_global.right.insert({ global[k], dofs[k] });
294+
assert(inserted);
295+
}
296+
}
297+
}
298+
}
208299

209300
P1(const P1& other) = default;
210301

@@ -229,7 +320,7 @@ for (int r = 0; r < comm.size(); ++r) {
229320
*/
230321
Index getGlobalIndex(Index globalIdx) const
231322
{
232-
return m_loc2glob.left.at(globalIdx).get_right();
323+
return m_local_to_global.left.at(globalIdx).get_right();
233324
}
234325

235326
/**
@@ -238,8 +329,8 @@ for (int r = 0; r < comm.size(); ++r) {
238329
*/
239330
Optional<Index> getLocalIndex(Index globalIdx) const
240331
{
241-
auto find = m_loc2glob.right.find(globalIdx);
242-
if (find == m_loc2glob.right.end())
332+
auto find = m_local_to_global.right.find(globalIdx);
333+
if (find == m_local_to_global.right.end())
243334
return std::nullopt;
244335
else
245336
return find->get_left();
@@ -361,13 +452,19 @@ for (int r = 0; r < comm.size(); ++r) {
361452
return typename FESType::InverseMapping(v);
362453
}
363454

455+
const auto& getGhosts() const
456+
{
457+
return m_ghosts;
458+
}
459+
364460
private:
365461
std::reference_wrapper<const MeshType> m_mesh;
366462
FESType m_fes;
367463

368464
size_t m_offset;
369465
size_t m_owned;
370-
Bimap m_loc2glob;
466+
IndexBimap m_local_to_global;
467+
GhostBimap m_ghosts;
371468
};
372469
}
373470

‎src/Rodin/PETSc/IO/GridFunctionPrinter.h‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,8 @@
77
#ifndef RODIN_PETSC_IO_GRIDFUNCTIONPRINTER_H
88
#define RODIN_PETSC_IO_GRIDFUNCTIONPRINTER_H
99

10+
#include <petscvec.h>
11+
1012
#include "Rodin/Context/Local.h"
1113

1214
#include "Rodin/IO/ForwardDecls.h"
@@ -16,7 +18,6 @@
1618

1719
#include "Rodin/MPI/Context/MPI.h"
1820
#include "Rodin/PETSc/Math/Vector.h"
19-
#include <petscvec.h>
2021

2122
namespace Rodin::IO
2223
{
@@ -59,11 +60,11 @@ namespace Rodin::IO
5960
const auto& shard = fes.getShard();
6061
const auto& read = gf.acquire().getArrayRead();
6162
const size_t sz = shard.getSize();
62-
VecView(gf.getData(), PETSC_VIEWER_STDOUT_WORLD);
6363
assert(read.raw);
6464
assert(read.acquired);
6565
for (size_t i = 0; i < sz; ++i)
6666
os << read.raw[i] << "\n";
67+
gf.flush();
6768
}
6869
else
6970
{

0 commit comments

Comments
 (0)