Global field gather, global Silo zonal export. Needs cleanup.

CMakeLists.txt

@@ -178,7 +178,7 @@ else()
     set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wextra -Wno-sign-compare -Wshadow \
         -Wno-unknown-pragmas -Wno-unused-parameter")
 endif()
-set(CMAKE_CXX_FLAGS_DEBUG "-O0 -g -fpermissive")
+set(CMAKE_CXX_FLAGS_DEBUG "-O0 -g -fpermissive -fsanitize=address")
 set(CMAKE_CXX_FLAGS_RELEASE "-O3 -g -DNDEBUG -fpermissive")
 set(CMAKE_CXX_FLAGS_RELEASEASSERT "-O3 -g -fpermissive")
 

include/entity.h

@@ -83,6 +83,7 @@ public:
     const reference_element&    cell_refelem(size_t) const;
     const reference_element&    cell_refelem(const physical_element&) const;
     size_t                      cell_global_index(size_t) const;
+    size_t                      cell_local_index_by_gmsh(size_t) const;
     size_t                      cell_global_index_by_gmsh(size_t) const;
     size_t                      cell_local_dof_offset(size_t) const;
     size_t                      cell_global_dof_offset(size_t) const;
@@ -107,6 +108,7 @@ public:
 
     int                         material_tag(void) const;
     int                         gmsh_tag(void) const;
+    int                         gmsh_dim(void) const;
     int                         gmsh_elem_type(void) const;
     entity_ordering             current_elem_ordering(void) const;
     constexpr size_t            num_faces_per_elem(void) const { return 4; }

include/gmsh_io.h

@@ -88,13 +88,28 @@ struct dof_range
     int                     tag;
     int                     dim;
     size_t                  base;
-    size_t  length;
+    std::vector<size_t>     dof_offsets;
+
+    //dof_range(dof_range&&) = default;
+
+    dof_range(const entity& e, size_t b)
+    {
+        tag = e.gmsh_tag();
+        dim = e.gmsh_dim();
+        base = b;
+        dof_offsets.resize( e.num_dofs() );
+        auto dofs_per_cell = e.num_dofs()/e.num_cells();
+        assert( dofs_per_cell*e.num_cells() == e.num_dofs() );
+        for (size_t iT = 0; iT < e.num_cells(); iT++)
+            for (size_t iD = 0; iD < dofs_per_cell; iD++)
+                dof_offsets[iT*dofs_per_cell+iD] = e.cell_global_dof_offset(iT)+iD;
+    }
 };
 
 inline std::ostream& operator<<(std::ostream& os, const dof_range& dr)
 {
     os << "Tag " << dr.tag << ", dim " << dr.dim << ": (";
-    os << dr.base << ", " << dr.length << ")";
+    os << dr.base << ", " << dr.dof_offsets.size() << ")";
     return os;
 }
 #endif
@@ -126,9 +141,12 @@ class model
 
     /* Map from 3D entity tag to partition */
     std::map<int, int>                          partition_inverse_map;
-
-    /* Global dof range, only rank0 has this */
+public:
+    /* RANK 0 only: keeps global (process-wide) DOF numbering and positioning */
     std::map<int, dof_range>                    global_dof_ranges;
+private:
+    /* RANK 0 only: copy of remote entities for postpro */
+    std::map<int, std::vector<entity>>          remote_entities;
 #endif /* USE_MPI */
 
     using entofs_pair = std::pair<size_t, size_t>;
@@ -165,6 +183,42 @@ public:
     void partition_mesh(int);
     void populate_from_gmsh(void);
 
+#ifdef USE_MPI
+    const std::map<int, dof_range>& dof_ranges(void) const
+    {
+        return global_dof_ranges;
+    }
+
+    size_t global_num_dofs() const
+    {
+        size_t ret = 0;
+        for (auto& dr : global_dof_ranges)
+            ret += dr.second.dof_offsets.size();
+        return ret;
+    }
+
+    size_t global_num_cells() const
+    {
+        size_t ret = 0;
+        for (auto& e : entities)
+            ret += e.num_cells();
+        
+        for (auto& re : remote_entities)
+            for (auto& e : re.second)
+                ret += e.num_cells();
+        
+        return ret;
+    }
+
+    int entity_rank(int tag) const
+    {
+        return partition_inverse_map.at(tag)-1;
+    }
+
+    entity&         entity_global_lookup(int);
+    const entity&   entity_global_lookup(int) const;
+#endif /* USE_MPI */
+
     const element_connectivity<element_key>& connectivity() const
     {
         return conn;

include/maxwell_interface.h

@@ -18,6 +18,11 @@
 #include "param_loader.h"
 #endif
 
+#ifdef USE_MPI
+#include <mpi/mpi.h>
+#include "gmsh_mpi.h"
+#endif /* USE_MPI */
+
 namespace maxwell {
 
 #ifndef HIDE_THIS_FROM_NVCC
@@ -94,6 +99,13 @@ struct field
     void    resize(size_t);
 };
 
+#ifdef USE_MPI
+void send_field(field& f, int dst, int tag, MPI_Comm comm);
+void send_field(field& f, int dst, int tag, size_t offset, size_t length, MPI_Comm comm);
+void receive_field(field& f, int src, int tag, MPI_Comm comm);
+void receive_field(field& f, int src, int tag, size_t offset, size_t length, MPI_Comm comm);
+#endif /* USE_MPI */
+
 #ifdef ENABLE_GPU_SOLVER
 struct pinned_field
 {
@@ -314,20 +326,28 @@ void init_from_model(const model&, solver_state&);
 void init_matparams(const model&, solver_state&, const parameter_loader&);
 //void apply_operator(solver_state&, const field&, field&);
 void export_fields_to_silo(const model&, const solver_state&, const parameter_loader&);
+void export_fields_to_silo_2(const model&, solver_state&, const parameter_loader&);
 void timestep(solver_state&, const parameter_loader&, time_integrator_type);
 void prepare_sources(const model&, maxwell::solver_state&, maxwell::parameter_loader&);
 void do_sources(const model&, maxwell::solver_state&, maxwell::parameter_loader&);
 void swap(maxwell::solver_state&, const parameter_loader&);
+#ifdef USE_MPI
+void gather_field(const model&, maxwell::field&, int, MPI_Comm);
+#endif /* USE_MPI */
 
 #ifdef ENABLE_GPU_SOLVER
 void init_from_model(const model&, solver_state_gpu&);
 void init_matparams(const model&, solver_state_gpu&, const parameter_loader&);
 //void apply_operator(solver_state_gpu&, const field_gpu&, field_gpu&);
 void export_fields_to_silo(const model&, const solver_state_gpu&, const parameter_loader&);
+void export_fields_to_silo_2(const model&, solver_state_gpu&, const parameter_loader&);
 void timestep(solver_state_gpu&, const parameter_loader&, time_integrator_type);
 void prepare_sources(const model&, maxwell::solver_state_gpu&, maxwell::parameter_loader&);
 void do_sources(const model&, maxwell::solver_state_gpu&, maxwell::parameter_loader&);
 void swap(maxwell::solver_state_gpu&, const parameter_loader&);
+#ifdef USE_MPI
+void gather_field(const model&, maxwell::field_gpu&, int, MPI_Comm);
+#endif /* USE_MPI */
 #endif /* ENABLE_GPU_SOLVER */
 
 void init_boundary_conditions(const model&, const parameter_loader&, vecxd&);

src/entity.cpp

@@ -273,6 +273,13 @@ entity::cell_global_index(size_t iT) const
     return m_index_base + iT;
 }
 
+size_t
+entity::cell_local_index_by_gmsh(size_t iT) const
+{
+    assert(iT < physical_cells.size());
+    return physical_cells[iT].original_position();
+}
+
 size_t
 entity::cell_global_index_by_gmsh(size_t iT) const
 {
@@ -517,6 +524,12 @@ entity::gmsh_tag(void) const
     return tag;
 }
 
+int
+entity::gmsh_dim(void) const
+{
+    return dim;
+}
+
 int
 entity::gmsh_elem_type(void) const
 {

src/gmsh_io.cpp

@@ -235,22 +235,26 @@ model::import_gmsh_entities_rank0(void)
 #ifdef USE_MPI
             if (num_partitions > 1)
             {
+                /* If there are multiple partitions, decide to which one this
+                 * entity belongs to. If the partition is remote, send entity
+                 * to the remote process. The entity is saved also locally
+                 * to have all the data for postprocessing later. */
                 auto mp = my_partition();
                 auto tag_partition = partition_inverse_map.at(tag);
                 bool entity_is_remote = ( tag_partition != mp );
             
                 if (entity_is_remote)
                 {
-                    dof_range dr;
-                    dr.dim = dim;               dr.tag = tag;
-                    dr.base = global_dof_base;  dr.length = e.num_dofs();
-                    global_dof_ranges[tag] = dr;
+                    dof_range dr(e, global_dof_base);
+                    global_dof_ranges.emplace(tag, std::move(dr));
                     global_dof_base += e.num_dofs();
 
                     assert(tag_partition > 0);
                     int rank = tag_partition - 1;
                     e.mpi_send(rank, MPI_COMM_WORLD);
-                    std::cout << "Sending entity " << tag << " to " << rank << std::endl;
+                    std::cout << "Sending entity " << tag << " to " << rank << " " << e.num_dofs() << std::endl;
+                    /* Save locally the entity */
+                    remote_entities[tag_partition].push_back( std::move(e) );
                     continue;
                 }
             }
@@ -273,19 +277,37 @@ model::import_gmsh_entities_rank0(void)
             index_base += e.num_cells();
 
 #ifdef USE_MPI
-            dof_range dr;
-            dr.dim = dim;               dr.tag = tag;
-            dr.base = global_dof_base;  dr.length = e.num_dofs();
-            global_dof_ranges[tag] = dr;
+            dof_range dr(e, global_dof_base);
+            global_dof_ranges.emplace(tag, std::move(dr));
             global_dof_base += e.num_dofs();
 #endif /* USE_MPI */
 
             update_connectivity(e, entity_index);
-
             entities.push_back( std::move(e) );
             entity_index++;
         }
     }
+
+#ifdef USE_MPI
+    /* Compute DOF offset in the remote entities saved locally. */
+    assert( IMPLIES(num_partitions == 1, remote_entities.size() == 0) );
+    for (auto& re : remote_entities)
+    {
+        entity_index = 0;
+        dof_base = 0;
+        flux_base = 0;
+        index_base = 0;
+        for (auto& e : re.second)
+        {
+            e.base(dof_base, flux_base, index_base);
+            e.number(entity_index);
+            dof_base += e.num_dofs();
+            flux_base += e.num_fluxes();
+            index_base += e.num_cells();
+            entity_index++;
+        }
+    }
+#endif
 }
 
 #ifdef USE_MPI
@@ -301,9 +323,9 @@ model::import_gmsh_entities_rankN(void)
 
     for (auto& rtag : partition_map.at(mp))
     {
-        std::cout << "Receiving " << rtag << ", rank " << mp-1 << std::endl;
         entity e;
         e.mpi_recv(0, MPI_COMM_WORLD);
+        std::cout << "Receiving " << rtag << ", rank " << mp-1 << " " << e.num_dofs() << std::endl;
 
         for (size_t i = 0; i < e.num_cells(); i++)
         {
@@ -561,6 +583,40 @@ model::partition_mesh(int parts)
     gmm::partition(parts);
     num_partitions = parts;
 }
+
+entity&
+model::entity_global_lookup(int tag)
+{
+    /* Look for tag between local entities */
+    for (auto& e : entities)
+        if (e.gmsh_tag() == tag)
+            return e;
+
+    /* Try between remote entities */
+    int partition = partition_inverse_map.at(tag);
+    for (auto& e : remote_entities[partition])
+        if (e.gmsh_tag() == tag)
+            return e;
+    
+    throw std::invalid_argument("tag not found");
+}
+
+const entity&
+model::entity_global_lookup(int tag) const
+{
+    /* Look for tag between local entities */
+    for (const auto& e : entities)
+        if (e.gmsh_tag() == tag)
+            return e;
+
+    /* Try between remote entities */
+    int partition = partition_inverse_map.at(tag);
+    for (const auto& e : remote_entities.at(partition))
+        if (e.gmsh_tag() == tag)
+            return e;
+    
+    throw std::invalid_argument("tag not found");
+}
 #endif /* USE_MPI */
 
 std::vector<entity>::const_iterator

src/maxwell_cpu.cpp

@@ -590,6 +590,114 @@ swap(maxwell::solver_state& state, const parameter_loader&)
     std::swap(state.emf_curr, state.emf_next);
 }
 
+#ifdef USE_MPI
+void
+gather_field(const model& mod, maxwell::field& in, maxwell::field& out,
+    int root, MPI_Comm comm)
+{
+    int rank;
+    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+
+    if (rank == root)
+    {
+        auto& gdr = mod.dof_ranges();
+        for (auto& dr : gdr)
+        {
+            auto tag = dr.second.tag;
+            if (mod.entity_rank(tag) == root)
+                continue;
+
+            auto base = dr.second.base;
+            auto len = dr.second.dof_offsets.size();
+            //std::cout << "GR: " << tag << " " << len << std::endl;
+            receive_field(out, mod.entity_rank(tag), tag, base, len, comm);
+        }
+
+        for (auto& e : mod)
+        {
+            auto in_base = e.dof_base();
+            auto in_len = e.num_dofs();
+            auto out_base = gdr.at(e.gmsh_tag()).base;
+            auto out_len = gdr.at(e.gmsh_tag()).dof_offsets.size();
+            assert(in_len == out_len);
+            out.Ex.segment(out_base, out_len) = in.Ex.segment(in_base, in_len);
+            out.Ey.segment(out_base, out_len) = in.Ey.segment(in_base, in_len);
+            out.Ez.segment(out_base, out_len) = in.Ez.segment(in_base, in_len);
+            out.Hx.segment(out_base, out_len) = in.Hx.segment(in_base, in_len);
+            out.Hy.segment(out_base, out_len) = in.Hy.segment(in_base, in_len);
+            out.Hz.segment(out_base, out_len) = in.Hz.segment(in_base, in_len);
+        }
+    }
+    else
+    {
+        for (auto& e : mod)
+        {
+            //std::cout << "GS: " << e.gmsh_tag() << " " << e.num_dofs() << std::endl;
+            send_field(in, root, e.gmsh_tag(), e.dof_base(), e.num_dofs(), comm);
+        }
+    }
+}
+#endif
+
+void
+export_fields_to_silo_2(const model& mod, maxwell::solver_state& state,
+    const maxwell::parameter_loader& mpl)
+{
+    int rank;
+    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+
+    if (rank == 0)
+    {
+        size_t gdnum = mod.global_num_dofs();
+        maxwell::field f;
+        f.resize(gdnum);
+        gather_field(mod, state.emf_next, f, 0, MPI_COMM_WORLD);
+    
+        std::stringstream ss;
+        ss << mpl.sim_name() << "/timestep_" << state.curr_timestep << ".silo";
+    
+        silo sdb;
+        sdb.create_db(ss.str());
+        sdb.import_mesh_from_gmsh();
+        sdb.write_mesh(state.curr_time, state.curr_timestep);
+
+        std::vector<double> out_Ex( mod.global_num_cells() );
+        std::vector<double> out_Ey( mod.global_num_cells() );
+        std::vector<double> out_Ez( mod.global_num_cells() );
+
+        size_t nc = 0;
+        for (auto& dr : mod.dof_ranges())
+        {
+            auto etag = dr.first;
+            std::cout << etag << std::endl;
+            auto& e = mod.entity_global_lookup(etag);
+            for (size_t iT = 0; iT < e.num_cells(); iT++)
+            {
+                auto& pe = e.cell(iT);
+                auto& re = e.cell_refelem(pe);
+                vecxd phi_bar = re.basis_functions({1./3., 1./3., 1./3.});
+                auto num_bf = re.num_basis_functions();
+                auto ofs = dr.second.base + num_bf*iT;
+                auto gi = nc+e.cell_local_index_by_gmsh(iT);
+                out_Ex[gi] = f.Ex.segment(ofs, num_bf).dot(phi_bar);
+                out_Ey[gi] = f.Ey.segment(ofs, num_bf).dot(phi_bar);
+                out_Ez[gi] = f.Ez.segment(ofs, num_bf).dot(phi_bar);
+            }
+            nc += e.num_cells();
+        }
+
+        sdb.write_zonal_variable("Ex", out_Ex);
+        sdb.write_zonal_variable("Ey", out_Ey);
+        sdb.write_zonal_variable("Ez", out_Ez);
+        sdb.write_vector_expression("E", "{Ex, Ey, Ez}");    
+    }
+    else
+    {
+        maxwell::field dummy;
+        gather_field(mod, state.emf_next, dummy, 0, MPI_COMM_WORLD);
+    }
+}
+
 void
 export_fields_to_silo(const model& mod, const maxwell::solver_state& state,
     const maxwell::parameter_loader& mpl)

src/maxwell_gpu.cpp

@@ -555,6 +555,16 @@ swap(solver_state_gpu& state, const parameter_loader& mpl)
     std::swap(state.emf_curr, state.emf_next);
 }
 
+void
+gather_field(const model& mod, const maxwell::field_gpu& in, maxwell::field_gpu& f,
+    int root, MPI_Comm comm)
+{}
+
+void
+export_fields_to_silo_2(const model& mod, maxwell::solver_state_gpu& state,
+    const maxwell::parameter_loader& mpl)
+{}
+
 void
 export_fields_to_silo(const model& mod, const maxwell::solver_state_gpu& state,
     const maxwell::parameter_loader& mpl)

src/maxwell_interface.cpp

@@ -25,6 +25,67 @@ field::resize(size_t p_num_dofs)
     num_dofs = p_num_dofs;
 }
 
+#ifdef USE_MPI
+void send_field(field& f, int dst, int tag, MPI_Comm comm)
+{
+    MPI_Send(&f.num_dofs, 1, MPI_UNSIGNED_LONG_LONG, dst, tag, comm);
+    MPI_Send(f.Ex.data(), f.num_dofs, MPI_DOUBLE, dst, tag, comm);
+    MPI_Send(f.Ey.data(), f.num_dofs, MPI_DOUBLE, dst, tag, comm);
+    MPI_Send(f.Ez.data(), f.num_dofs, MPI_DOUBLE, dst, tag, comm);
+    MPI_Send(f.Hx.data(), f.num_dofs, MPI_DOUBLE, dst, tag, comm);
+    MPI_Send(f.Hy.data(), f.num_dofs, MPI_DOUBLE, dst, tag, comm);
+    MPI_Send(f.Hz.data(), f.num_dofs, MPI_DOUBLE, dst, tag, comm);
+}
+
+void send_field(field& f, int dst, int tag, size_t offset, size_t length, MPI_Comm comm)
+{
+    MPI_Send(&length, 1, MPI_UNSIGNED_LONG_LONG, dst, tag, comm);
+    assert(offset+length <= f.num_dofs);
+    
+    MPI_Send(f.Ex.data()+offset, length, MPI_DOUBLE, dst, tag, comm);
+    MPI_Send(f.Ey.data()+offset, length, MPI_DOUBLE, dst, tag, comm);
+    MPI_Send(f.Ez.data()+offset, length, MPI_DOUBLE, dst, tag, comm);
+    MPI_Send(f.Hx.data()+offset, length, MPI_DOUBLE, dst, tag, comm);
+    MPI_Send(f.Hy.data()+offset, length, MPI_DOUBLE, dst, tag, comm);
+    MPI_Send(f.Hz.data()+offset, length, MPI_DOUBLE, dst, tag, comm);
+}
+
+void receive_field(field& f, int src, int tag, MPI_Comm comm)
+{
+    MPI_Status status;
+    size_t ndofs;
+    MPI_Recv(&ndofs, 1, MPI_UNSIGNED_LONG_LONG, src, tag, comm, &status);
+    if (ndofs != f.num_dofs)
+        f.resize(ndofs);
+
+    MPI_Recv(f.Ex.data(), ndofs, MPI_DOUBLE, src, tag, comm, &status);
+    MPI_Recv(f.Ey.data(), ndofs, MPI_DOUBLE, src, tag, comm, &status);
+    MPI_Recv(f.Ez.data(), ndofs, MPI_DOUBLE, src, tag, comm, &status);
+    MPI_Recv(f.Hx.data(), ndofs, MPI_DOUBLE, src, tag, comm, &status);
+    MPI_Recv(f.Hy.data(), ndofs, MPI_DOUBLE, src, tag, comm, &status);
+    MPI_Recv(f.Hz.data(), ndofs, MPI_DOUBLE, src, tag, comm, &status);
+}
+
+void receive_field(field& f, int src, int tag, size_t offset, size_t length, MPI_Comm comm)
+{
+    MPI_Status status;
+    size_t ndofs;
+    MPI_Recv(&ndofs, 1, MPI_UNSIGNED_LONG_LONG, src, tag, comm, &status);
+    if ( (ndofs != length) or (offset+length > f.num_dofs) )
+    {
+        std::cout << "Unexpected size in receiving field" << std::endl;
+        MPI_Abort(MPI_COMM_WORLD, 1);
+    }
+
+    MPI_Recv(f.Ex.data()+offset, ndofs, MPI_DOUBLE, src, tag, comm, &status);
+    MPI_Recv(f.Ey.data()+offset, ndofs, MPI_DOUBLE, src, tag, comm, &status);
+    MPI_Recv(f.Ez.data()+offset, ndofs, MPI_DOUBLE, src, tag, comm, &status);
+    MPI_Recv(f.Hx.data()+offset, ndofs, MPI_DOUBLE, src, tag, comm, &status);
+    MPI_Recv(f.Hy.data()+offset, ndofs, MPI_DOUBLE, src, tag, comm, &status);
+    MPI_Recv(f.Hz.data()+offset, ndofs, MPI_DOUBLE, src, tag, comm, &status);
+}
+#endif /* USE_MPI */
+
 #ifdef ENABLE_GPU_SOLVER
 
 pinned_field::pinned_field()

src/maxwell_solver.cpp

@@ -221,12 +221,12 @@ void test_it(const model& mod, State& state, maxwell::parameter_loader& mpl)
         do_sources(mod, state, mpl);
         
 #ifdef USE_MPI
-        int rank;
-        MPI_Comm_rank(MPI_COMM_WORLD, &rank);
-        if (rank == 0)
+        //int rank;
+        //MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+        //if (rank == 0)
 #endif /* USE_MPI */
             if ( (silo_output_rate != 0) and (i%silo_output_rate == 0) )
-                export_fields_to_silo(mod, state, mpl);
+                export_fields_to_silo_2(mod, state, mpl);
 
         swap(state, mpl);