Merge pull request #39 from aperijake/nodal_integration

include/IoMesh.h

@@ -12,7 +12,7 @@ namespace aperi {
 
 struct IoMeshParameters {
     bool upward_connectivity = true;        // create upward connectivity/adjacency in the mesh
-    bool aura_option = true;                // create aura ghosting around each MPI rank
+    bool aura_option = false;               // create aura ghosting around each MPI rank
     std::string parallel_io = "pnetcdf";    // method to use for parallel io. One of mpiio, mpiposix, or pnetcdf
     std::string decomp_method = "rcb";      // decomposition method.  One of: linear, rcb, rib, hsfc, block, cyclic, random, kway, geom_kway, metis_sfc
     std::string mesh_type = "exodusII";     // mesh type. One of: exodusii, generated

include/MeshLabelerProcessor.h

@@ -149,6 +149,8 @@ class MeshLabelerProcessor {
 
                 // Get the minimum id
                 uint64_t minimum_id = m_bulk_data->identifier(nodes[0]);
+                uint64_t *active_field_data = stk::mesh::field_data(*m_active_field, nodes[0]);
+                active_field_data[0] = 0;  // Set active value to 0 for the first node
                 size_t minimum_index = 0;
                 for (size_t i = 1; i < num_nodes; ++i) {
                     uint64_t id = m_bulk_data->identifier(nodes[i]);
@@ -157,12 +159,12 @@ class MeshLabelerProcessor {
                         minimum_index = i;
                     }
                     // Set active value to 0
-                    uint64_t *active_field_data = stk::mesh::field_data(*m_active_field, nodes[i]);
+                    active_field_data = stk::mesh::field_data(*m_active_field, nodes[i]);
                     active_field_data[0] = 0;
                 }
 
                 // Set the active value to 1 for the minimum id
-                uint64_t *active_field_data = stk::mesh::field_data(*m_active_field, nodes[minimum_index]);
+                active_field_data = stk::mesh::field_data(*m_active_field, nodes[minimum_index]);
                 active_field_data[0] = 1;
             }
         }
@@ -191,7 +193,12 @@ class MeshLabelerProcessor {
                     }
                 }
                 if (num_active_nodes != 1) {
-                    std::string message = "Nodal integration requires exactly one active node per element. Found " + std::to_string(num_active_nodes) + " active nodes.";
+                    size_t element_id = m_bulk_data->identifier(element);
+                    std::string message = "Nodal integration requires exactly one active node per element. Found " + std::to_string(num_active_nodes) + " active nodes in element " + std::to_string(element_id) + ". Nodes: \n";
+                    for (size_t i = 0; i < num_nodes; ++i) {
+                        uint64_t *active_field_data = stk::mesh::field_data(*m_active_field, nodes[i]);
+                        message += " ID: " + std::to_string(m_bulk_data->identifier(nodes[i])) + ", active value: " + std::to_string(active_field_data[0]) + "\n";
+                    }
                     throw std::runtime_error(message);
                 }
             }

include/NeighborSearchProcessor.h

@@ -8,6 +8,7 @@
 #include <stk_mesh/base/BulkData.hpp>
 #include <stk_mesh/base/Field.hpp>
 #include <stk_mesh/base/FieldBLAS.hpp>
+#include <stk_mesh/base/FieldParallel.hpp>
 #include <stk_mesh/base/ForEachEntity.hpp>
 #include <stk_mesh/base/GetEntities.hpp>
 #include <stk_mesh/base/GetNgpField.hpp>
@@ -291,6 +292,11 @@ class NeighborSearchProcessor {
         ngp_kernel_radius_field.clear_sync_state();
         ngp_kernel_radius_field.modify_on_device();
         ngp_kernel_radius_field.sync_to_host();
+
+        // Get the parallel max kernel radius
+        stk::mesh::parallel_max(*m_bulk_data, {m_kernel_radius_field});
+        ngp_kernel_radius_field.modify_on_host();
+        ngp_kernel_radius_field.sync_to_device();
     }
 
     // Create local entities on host and copy to device

test/performance_tests/gtests/SolverPerformanceTest.cpp

@@ -108,7 +108,7 @@ class SolverPerformanceTest : public SolverTest {
             size_t num_elem_z_next = num_elem_z.back() * refinement_factors[2];
 
             // If the number of elements is too large (greater than 36 million), stop the refinement
-            size_t num_elem_next = (num_elem_x_next) * (num_elem_y_next) * (num_elem_z_next) * 6;
+            size_t num_elem_next = (num_elem_x_next) * (num_elem_y_next) * (num_elem_z_next)*6;
             if (num_elem_next > 36e6) {
                 aperi::CoutP0() << "Number of elements is too large. Stopping the refinement. Number of elements: " << num_elem_next << std::endl;
                 num_refinements = i;
@@ -159,7 +159,10 @@ class SolverPerformanceTest : public SolverTest {
             aperi::CoutP0() << "  Runtime (s) / increment:      " << std::scientific << average_increment_runtime << std::endl;
 
             // Name of the benchmark: Taylor Impact: m_num_procs processors, cpu/gpu, num_elem_x x num_elem_y x num_elem_z elements, runtime per increment"
-            std::string name = "Taylor Impact: " + std::to_string(m_num_procs) + " processors, " + (using_gpu ? "gpu" : "cpu") + ", " + std::to_string(num_elem_x[i]) + " x " + std::to_string(num_elem_y[i]) + " x " + std::to_string(num_elem_z[i]) + " elements, runtime per increment";
+            std::string name = "Taylor Impact";
+            name += (reproducing_kernel ? ", Reproducing Kernel:" : ":") + std::to_string(m_num_procs) + " processors, ";
+            name += (using_gpu ? "gpu, " : "cpu, ") + std::to_string(num_elem_x[i]) + " x " + std::to_string(num_elem_y[i]) + " x " + std::to_string(num_elem_z[i]);
+            name += " elements, runtime per increment";
             // Unit of the benchmark
             std::string unit = "milliseconds";
             // Value of the benchmark