diff --git a/apps/multiphysics/src/Solvers/TLQS_solver_3D/include/tlqs_solver_3D.hpp b/apps/multiphysics/src/Solvers/TLQS_solver_3D/include/tlqs_solver_3D.hpp index 8f02d2e8d..221e58b73 100644 --- a/apps/multiphysics/src/Solvers/TLQS_solver_3D/include/tlqs_solver_3D.hpp +++ b/apps/multiphysics/src/Solvers/TLQS_solver_3D/include/tlqs_solver_3D.hpp @@ -281,7 +281,7 @@ class TLQS3D : public Solver const CArrayKokkos& F_elem, const CArrayKokkos& K_elem, const CArrayKokkos& displacement_iter, - const CArrayKokkos& r0 + MPICArrayKokkos& r0 ); // inputs: mesh.num_nodes, mesh.elems_in_node, mesh.num_nodes_in_elem, mesh.nodes_in_elem, K_elem, rk, p @@ -289,10 +289,14 @@ class TLQS3D : public Solver double get_alpha( const size_t num_nodes, const size_t num_nodes_in_elem, + const size_t num_owned_nodes, + const RaggedRightArrayKokkos& elems_in_node, const DCArrayKokkos& nodes_in_elem, const CArrayKokkos& K_elem, const double rktrk, - const CArrayKokkos& p + const CArrayKokkos& p, + MPICArrayKokkos& temporary, + const DCArrayKokkos shared_tally_owned_nodes ); void get_rkp1( @@ -301,10 +305,10 @@ class TLQS3D : public Solver const size_t num_nodes_in_elem, const DCArrayKokkos& nodes_in_elem, const CArrayKokkos& K_elem, - const CArrayKokkos& rk, + const MPICArrayKokkos& rk, const CArrayKokkos& p, const double alpha, - const CArrayKokkos& rkp1 + MPICArrayKokkos& rkp1 ); // **** Functions defined in post_process.cpp **** // @@ -371,12 +375,12 @@ class TLQS3D : public Solver ); */ // **** Functions defined in chebyshev_smoothing.cpp **** // - void apply_chebyshev_preconditioner(const CArrayKokkos& rk, - const CArrayKokkos& zkp1, - const CArrayKokkos& D_inv, - const CArrayKokkos& zk, + void apply_chebyshev_preconditioner(const MPICArrayKokkos& rk, + const MPICArrayKokkos& zkp1, + const MPICArrayKokkos& D_inv, + const MPICArrayKokkos& zk, const CArrayKokkos& delta_z, - const CArrayKokkos& temporary, + MPICArrayKokkos& temporary, const CArrayKokkos& K_elem, const size_t num_nodes, const RaggedRightArrayKokkos& elems_in_node, @@ -387,7 +391,7 @@ class TLQS3D : public Solver const int degree ); - void get_diagonal_inverse(CArrayKokkos& D_inv, + void get_diagonal_inverse(MPICArrayKokkos& D_inv, const CArrayKokkos& K_elem, const size_t num_nodes, const RaggedRightArrayKokkos& elems_in_node, @@ -397,15 +401,17 @@ class TLQS3D : public Solver void get_chebyshev_bounds(double& alpha, double& beta, - const CArrayKokkos& D_inv, + const MPICArrayKokkos& D_inv, const CArrayKokkos& K_elem, const size_t num_nodes, const RaggedRightArrayKokkos& elems_in_node, const size_t num_nodes_in_elem, const DCArrayKokkos& nodes_in_elem, - CArrayKokkos& v_scratch, - CArrayKokkos& w_scratch, - const int max_iters + MPICArrayKokkos& v_scratch, + MPICArrayKokkos& w_scratch, + const int max_iters, + const int num_owned_nodes, + const DCArrayKokkos shared_tally_owned_nodes ); }; diff --git a/apps/multiphysics/src/Solvers/TLQS_solver_3D/src/cgm_functions.cpp b/apps/multiphysics/src/Solvers/TLQS_solver_3D/src/cgm_functions.cpp index 9e6f19880..edc96855a 100644 --- a/apps/multiphysics/src/Solvers/TLQS_solver_3D/src/cgm_functions.cpp +++ b/apps/multiphysics/src/Solvers/TLQS_solver_3D/src/cgm_functions.cpp @@ -43,7 +43,7 @@ void TLQS3D::get_r0( const CArrayKokkos & F_elem, const CArrayKokkos & K_elem, const CArrayKokkos & displacement_iter, - const CArrayKokkos & r0 + MPICArrayKokkos & r0 ) { // getting r0 = (02F - 01F) - K * displacement_iter @@ -51,7 +51,6 @@ void TLQS3D::get_r0( const size_t num_elems_in_node = elems_in_node.stride(node_gid); for (size_t p = 0; p < 3; p++) { - const size_t global_dof = 3 * node_gid + p; double val = 0.0; // Sum contributions from all elements containing this node @@ -78,58 +77,81 @@ void TLQS3D::get_r0( const size_t node_gid_b = nodes_in_elem(elem_gid, b); for (size_t q = 0; q < 3; q++) { const size_t local_dof_b = 3 * b + q; - const size_t global_dof_b = 3 * node_gid_b + q; - val -= K_elem(elem_gid, local_dof, local_dof_b) * displacement_iter(global_dof_b); + val -= K_elem(elem_gid, local_dof, local_dof_b) * displacement_iter(node_gid_b, q); //std::cout << "K_ELEM: " << K_elem(elem_gid, local_dof, local_dof_b) << std::endl; } } } - r0(global_dof) = val; + r0(node_gid, p) = val; } }); Kokkos::fence(); + r0.communicate(); } // end get_r0 double TLQS3D::get_alpha( const size_t num_nodes, const size_t num_nodes_in_elem, + const size_t num_owned_nodes, + const RaggedRightArrayKokkos& elems_in_node, const DCArrayKokkos& nodes_in_elem, const CArrayKokkos& K_elem, const double rktrk, - const CArrayKokkos& p) + const CArrayKokkos& p, + MPICArrayKokkos& temporary, + const DCArrayKokkos shared_tally_owned_nodes + ) { + // Kernel 1: compute temporary = K * p + FOR_ALL(node_gid, 0, num_nodes, { + for (size_t p_dir = 0; p_dir < 3; p_dir++) { + double val = 0.0; - // denominator: p^T * K * p - // first compute Kp = K * p via assembly-free matvec - // then dot with p - double ptkp = 0.0; - double loc_ptkp = 0.0; - FOR_REDUCE_SUM(elem_gid, 0, K_elem.dims(0), loc_ptkp, { + for (size_t elem_lid = 0; elem_lid < elems_in_node.stride(node_gid); elem_lid++) { + const size_t elem_gid = elems_in_node(node_gid, elem_lid); - for (size_t a = 0; a < num_nodes_in_elem; a++) { - const size_t node_gid_a = nodes_in_elem(elem_gid, a); - for (size_t p_dir = 0; p_dir < 3; p_dir++) { - const size_t local_dof_a = 3 * a + p_dir; - const size_t global_dof_a = 3 * node_gid_a + p_dir; + size_t local_node_lid = num_nodes_in_elem; + for (size_t a = 0; a < num_nodes_in_elem; a++) { + if (nodes_in_elem(elem_gid, a) == node_gid) { + local_node_lid = a; + break; + } + } + + const size_t local_dof = 3 * local_node_lid + p_dir; - double Kp_val = 0.0; for (size_t b = 0; b < num_nodes_in_elem; b++) { const size_t node_gid_b = nodes_in_elem(elem_gid, b); for (size_t q = 0; q < 3; q++) { - const size_t local_dof_b = 3 * b + q; - const size_t global_dof_b = 3 * node_gid_b + q; - Kp_val += K_elem(elem_gid, local_dof_a, local_dof_b) * p(global_dof_b); + const size_t local_dof_b = 3 * b + q; + val += K_elem(elem_gid, local_dof, local_dof_b) * p(node_gid_b, q); } } - loc_ptkp += p(global_dof_a) * Kp_val; + } + temporary(node_gid, p_dir) = val; + } + }); + MATAR_FENCE(); + + // temporary ghost entries needed for dot product + temporary.communicate(); + + // Kernel 2: p^T * temporary over owned nodes only, then Allreduce + double ptkp = 0.0; + double loc_ptkp = 0.0; + FOR_REDUCE_SUM(node_gid, 0, (int)num_owned_nodes, loc_ptkp, { + if(shared_tally_owned_nodes(node_gid)){ + for (int j = 0; j < 3; j++) { + loc_ptkp += p(node_gid, j) * temporary(node_gid, j); } } }, ptkp); - Kokkos::fence(); - //std::cout << "PTKP: " << ptkp << std::endl; + MATAR_FENCE(); + + MPI_Allreduce(MPI_IN_PLACE, &ptkp, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); - return rktrk / (ptkp+1E-16); + return rktrk / (ptkp + 1e-16); } // end get_alpha void TLQS3D::get_rkp1( @@ -138,17 +160,16 @@ void TLQS3D::get_rkp1( const size_t num_nodes_in_elem, const DCArrayKokkos& nodes_in_elem, const CArrayKokkos& K_elem, - const CArrayKokkos& rk, + const MPICArrayKokkos& rk, const CArrayKokkos& p, const double alpha, - const CArrayKokkos& rkp1) + MPICArrayKokkos& rkp1) { // r_{k+1} = r_k - alpha * K * p FOR_ALL(node_gid, 0, num_nodes, { const size_t num_elems_in_node = elems_in_node.stride(node_gid); for (size_t p_dir = 0; p_dir < 3; p_dir++) { - const size_t global_dof = 3 * node_gid + p_dir; double Kp_val = 0.0; for (size_t elem_lid = 0; elem_lid < num_elems_in_node; elem_lid++) { @@ -169,14 +190,14 @@ void TLQS3D::get_rkp1( const size_t node_gid_b = nodes_in_elem(elem_gid, b); for (size_t q = 0; q < 3; q++) { const size_t local_dof_b = 3 * b + q; - const size_t global_dof_b = 3 * node_gid_b + q; - Kp_val += K_elem(elem_gid, local_dof, local_dof_b) * p(global_dof_b); + Kp_val += K_elem(elem_gid, local_dof, local_dof_b) * p(node_gid_b, q); } } } - rkp1(global_dof) = rk(global_dof) - alpha * Kp_val; + rkp1(node_gid, p_dir) = rk(node_gid, p_dir) - alpha * Kp_val; } }); Kokkos::fence(); + rkp1.communicate(); } // end get_rkp1 \ No newline at end of file diff --git a/apps/multiphysics/src/Solvers/TLQS_solver_3D/src/chebyshev_smoothing.cpp b/apps/multiphysics/src/Solvers/TLQS_solver_3D/src/chebyshev_smoothing.cpp index cb76438ea..f98b16a94 100644 --- a/apps/multiphysics/src/Solvers/TLQS_solver_3D/src/chebyshev_smoothing.cpp +++ b/apps/multiphysics/src/Solvers/TLQS_solver_3D/src/chebyshev_smoothing.cpp @@ -39,12 +39,12 @@ ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * @brief Applies a matrix-free Chebyshev polynomial preconditioner using a thread-safe, * node-based gathering approach (no atomic operations required). */ -void TLQS3D::apply_chebyshev_preconditioner(const CArrayKokkos& rk, - const CArrayKokkos& zkp1, - const CArrayKokkos& D_inv, - const CArrayKokkos& zk, +void TLQS3D::apply_chebyshev_preconditioner(const MPICArrayKokkos& rk, + const MPICArrayKokkos& zkp1, + const MPICArrayKokkos& D_inv, + const MPICArrayKokkos& zk, const CArrayKokkos& delta_z, - const CArrayKokkos& temporary, + MPICArrayKokkos& temporary, const CArrayKokkos& K_elem, const size_t num_nodes, const RaggedRightArrayKokkos& elems_in_node, @@ -54,22 +54,21 @@ void TLQS3D::apply_chebyshev_preconditioner(const CArrayKokkos& rk, const double beta, const int degree) { - const size_t total_dofs = 3 * num_nodes; - // Compute Chebyshev parameters based on spectral bounds const double d = (beta + alpha) / 2.0; const double c = (beta - alpha) / 2.0; // --- Step 1: Initialize the 3-term recurrence (Iteration k = 0) --- - FOR_ALL(i, 0, total_dofs, { - zk(i) = 0.0; - delta_z(i) = (1.0 / d) * D_inv(i) * rk(i); - zk(i) += delta_z(i); + FOR_ALL(i, 0, (int)num_nodes, + j, 0, 3, { + zk(i,j) = 0.0; + delta_z(i,j) = (1.0 / d) * D_inv(i,j) * rk(i,j); + zk(i,j) += delta_z(i,j); }); MATAR_FENCE(); double rho_prev = c / (2.0 * d); - + // --- Step 2: Recurrence Loop (Iteration k = 1 to degree-1) --- for (int k = 1; k < degree; ++k) { double rho_k = 1.0 / (2.0 * d / c - rho_prev); @@ -81,7 +80,6 @@ void TLQS3D::apply_chebyshev_preconditioner(const CArrayKokkos& rk, const size_t num_elems_in_node = elems_in_node.stride(node_gid); for (size_t p = 0; p < 3; p++) { - const size_t global_dof = 3 * node_gid + p; double val = 0.0; // Sum contributions from all elements containing this global node @@ -105,24 +103,25 @@ void TLQS3D::apply_chebyshev_preconditioner(const CArrayKokkos& rk, for (size_t q = 0; q < 3; q++) { const size_t local_dof_b = 3 * b + q; - const size_t global_dof_b = 3 * node_gid_b + q; - val += K_elem(elem_gid, local_dof, local_dof_b) * zk(global_dof_b); + val += K_elem(elem_gid, local_dof, local_dof_b) * zk(node_gid_b, q); } } } // Directly assign to scratch array without atomics or clearing passes - temporary(global_dof) = val; + temporary(node_gid, p) = val; } }); MATAR_FENCE(); // ----------------------------------------------------------------- + temporary.communicate(); // Perform the vector updates using MATAR 1D indexing - FOR_ALL(i, 0, total_dofs, { - delta_z(i) = rho_k * delta_z(i) + gamma_k * D_inv(i) * (rk(i) - temporary(i)); - zk(i) += delta_z(i); + FOR_ALL(i, 0, (int)num_nodes, + j, 0, 3, { + delta_z(i,j) = rho_k * delta_z(i,j) + gamma_k * D_inv(i,j) * (rk(i,j) - temporary(i,j)); + zk(i,j) += delta_z(i,j); }); MATAR_FENCE(); @@ -130,8 +129,9 @@ void TLQS3D::apply_chebyshev_preconditioner(const CArrayKokkos& rk, } // --- Step 3: Finalize Output --- - FOR_ALL(i, 0, total_dofs, { - zkp1(i) = zk(i); + FOR_ALL(i, 0, (int)num_nodes, + j, 0, 3, { + zkp1(i,j) = zk(i,j); }); MATAR_FENCE(); } @@ -146,7 +146,7 @@ void TLQS3D::apply_chebyshev_preconditioner(const CArrayKokkos& rk, * @param num_nodes_in_elem Number of nodes per element (e.g., 64 for cubic hex) * @param nodes_in_elem Array mapping element ID to its global node IDs */ -void TLQS3D::get_diagonal_inverse(CArrayKokkos& D_inv, +void TLQS3D::get_diagonal_inverse(MPICArrayKokkos& D_inv, const CArrayKokkos& K_elem, const size_t num_nodes, const RaggedRightArrayKokkos& elems_in_node, @@ -157,7 +157,6 @@ void TLQS3D::get_diagonal_inverse(CArrayKokkos& D_inv, const size_t num_elems_in_node = elems_in_node.stride(node_gid); for (size_t p_dir = 0; p_dir < 3; p_dir++) { - const size_t global_dof = 3 * node_gid + p_dir; double diag = 0.0; for (size_t elem_lid = 0; elem_lid < num_elems_in_node; elem_lid++) { @@ -179,10 +178,11 @@ void TLQS3D::get_diagonal_inverse(CArrayKokkos& D_inv, } // Invert with safety epsilon protection - D_inv(global_dof) = 1.0 / (diag + 1e-16); + D_inv(node_gid, p_dir) = 1.0 / (diag + 1e-16); } }); MATAR_FENCE(); + D_inv.communicate(); } /** @@ -191,22 +191,24 @@ void TLQS3D::get_diagonal_inverse(CArrayKokkos& D_inv, */ void TLQS3D::get_chebyshev_bounds(double& alpha, double& beta, - const CArrayKokkos& D_inv, + const MPICArrayKokkos& D_inv, const CArrayKokkos& K_elem, const size_t num_nodes, const RaggedRightArrayKokkos& elems_in_node, const size_t num_nodes_in_elem, const DCArrayKokkos& nodes_in_elem, - CArrayKokkos& v_scratch, - CArrayKokkos& w_scratch, - const int max_iters) + MPICArrayKokkos& v_scratch, + MPICArrayKokkos& w_scratch, + const int max_iters, + const int num_owned_nodes, + const DCArrayKokkos shared_tally_owned_nodes) { - const size_t total_dofs = 3 * num_nodes; double lambda_max = 0.0; // Initialize initial guess vector v to 1.0 - FOR_ALL(i, 0, total_dofs, { - v_scratch(i) = 1.0; + FOR_ALL(i, 0, (int)num_nodes, + j, 0, 3, { + v_scratch(i,j) = 1.0; }); MATAR_FENCE(); @@ -218,7 +220,6 @@ void TLQS3D::get_chebyshev_bounds(double& alpha, const size_t num_elems_in_node = elems_in_node.stride(node_gid); for (size_t p = 0; p < 3; p++) { - const size_t global_dof = 3 * node_gid + p; double val = 0.0; for (size_t elem_lid = 0; elem_lid < num_elems_in_node; elem_lid++) { @@ -239,15 +240,14 @@ void TLQS3D::get_chebyshev_bounds(double& alpha, for (size_t q = 0; q < 3; q++) { const size_t local_dof_b = 3 * b + q; - const size_t global_dof_b = 3 * node_gid_b + q; - val += K_elem(elem_gid, local_dof, local_dof_b) * v_scratch(global_dof_b); + val += K_elem(elem_gid, local_dof, local_dof_b) * v_scratch(node_gid_b, q); } } } // Fused operation: Scale the accumulated stiffness action by D_inv - w_scratch(global_dof) = val * D_inv(global_dof); + w_scratch(node_gid, p) = val * D_inv(node_gid, p); } }); MATAR_FENCE(); @@ -257,33 +257,55 @@ void TLQS3D::get_chebyshev_bounds(double& alpha, double v_dot_v = 0.0; double local_v_dot_w = 0.0; - FOR_REDUCE_SUM(i, 0, total_dofs, local_v_dot_w, { - local_v_dot_w += v_scratch(i) * w_scratch(i); + FOR_REDUCE_SUM(i, 0, (int)num_owned_nodes, local_v_dot_w, { + if(shared_tally_owned_nodes(i)){ + for (int j = 0; j < 3; j++) { + local_v_dot_w += v_scratch(i,j) * w_scratch(i,j); + } + } }, v_dot_w); + MPI_Allreduce(MPI_IN_PLACE, &v_dot_w, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + double local_v_dot_v = 0.0; - FOR_REDUCE_SUM(i, 0, total_dofs, local_v_dot_v, { - local_v_dot_v += v_scratch(i) * v_scratch(i); + FOR_REDUCE_SUM(i, 0, (int)num_owned_nodes, local_v_dot_v, { + if(shared_tally_owned_nodes(i)){ + for (int j = 0; j < 3; j++) { + local_v_dot_v += v_scratch(i,j) * v_scratch(i,j); + } + } }, v_dot_v); MATAR_FENCE(); + MPI_Allreduce(MPI_IN_PLACE, &v_dot_v, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + // Calculate current estimate of the maximum eigenvalue lambda_max = v_dot_w / (v_dot_v + 1e-16); // 3. Normalize w vector to update our guess v: v = w / ||w||_2 double w_norm2 = 0.0; double local_w_norm2 = 0.0; - FOR_REDUCE_SUM(i, 0, total_dofs, local_w_norm2, { - local_w_norm2 += w_scratch(i) * w_scratch(i); + FOR_REDUCE_SUM(i, 0, (int)num_owned_nodes, local_w_norm2, { + if(shared_tally_owned_nodes(i)){ + for (int j = 0; j < 3; j++) { + local_w_norm2 += w_scratch(i,j) * w_scratch(i,j); + } + } }, w_norm2); MATAR_FENCE(); + + MPI_Allreduce(MPI_IN_PLACE, &w_norm2, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); double inv_norm = 1.0 / (sqrt(w_norm2) + 1e-16); - FOR_ALL(i, 0, total_dofs, { - v_scratch(i) = w_scratch(i) * inv_norm; + FOR_ALL(i, 0, (int)num_owned_nodes, + j, 0, 3, { + v_scratch(i,j) = w_scratch(i,j) * inv_norm; }); MATAR_FENCE(); + + // updating stale indices for following iteration + v_scratch.communicate(); } // --- 4. Apply Heuristic Bounding Box --- diff --git a/apps/multiphysics/src/Solvers/TLQS_solver_3D/src/elem_arrays.cpp b/apps/multiphysics/src/Solvers/TLQS_solver_3D/src/elem_arrays.cpp index e87bb6598..55e25ff89 100644 --- a/apps/multiphysics/src/Solvers/TLQS_solver_3D/src/elem_arrays.cpp +++ b/apps/multiphysics/src/Solvers/TLQS_solver_3D/src/elem_arrays.cpp @@ -105,7 +105,7 @@ void TLQS3D::get_gradients( const double dpsig_k2 = inv_J[2][0]*gauss_point_grad_basis(k,0) + inv_J[2][1]*gauss_point_grad_basis(k,1) + inv_J[2][2]*gauss_point_grad_basis(k,2); for (int j = 0; j < 3; j++) { - const double u_total = displacement(node_gid, j) + displacement_step(3*node_gid + j); + const double u_total = displacement(node_gid, j) + displacement_step(node_gid, j); grad_u[j][0] += u_total * dpsig_k0; grad_u[j][1] += u_total * dpsig_k1; grad_u[j][2] += u_total * dpsig_k2; diff --git a/apps/multiphysics/src/Solvers/TLQS_solver_3D/src/tlqs_execute.cpp b/apps/multiphysics/src/Solvers/TLQS_solver_3D/src/tlqs_execute.cpp index 22280f7b3..309e1c4c9 100644 --- a/apps/multiphysics/src/Solvers/TLQS_solver_3D/src/tlqs_execute.cpp +++ b/apps/multiphysics/src/Solvers/TLQS_solver_3D/src/tlqs_execute.cpp @@ -54,6 +54,15 @@ void TLQS3D::execute(SimulationParameters_t& SimulationParamaters, swage::Mesh_t& mesh, State_t& State) { + if (mesh.num_dims != 3) { + Kokkos::abort("TLQS SOLVER ONLY SUPPORTED IN 3D!!!"); + } + // Get MPI ranks and num ranks + int rank; + int num_ranks; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + MPI_Comm_size(MPI_COMM_WORLD, &num_ranks); + // Conveinent local variables double fuzz = SimulationParamaters.DynamicOptions.fuzz; double tiny = SimulationParamaters.DynamicOptions.tiny; @@ -108,15 +117,15 @@ void TLQS3D::execute(SimulationParameters_t& SimulationParamaters, CArrayKokkos zkp1(3 * mesh.num_nodes); */ // conjugate gradient method vectors - CArrayKokkos p(3*mesh.num_nodes); - CArrayKokkos rk(3*mesh.num_nodes); - CArrayKokkos rkp1(3*mesh.num_nodes); + CArrayKokkos p(mesh.num_nodes, 3); + MPICArrayKokkos rk(mesh.num_nodes, 3); + MPICArrayKokkos rkp1(mesh.num_nodes, 3); // Anderson acceleration variables size_t window_size = 1; const size_t max_hist = (window_size > 1) ? (window_size - 1) : 1; DCArrayKokkos anderson_weights(max_hist); - CArrayKokkos curr_anderson_residual(3*mesh.num_nodes); + //CArrayKokkos curr_anderson_residual(3*mesh.num_nodes); CArrayKokkos hist_anderson_residual(3*mesh.num_nodes,window_size); hist_anderson_residual.set_values(0); CArrayKokkos hist_displacement_iter(3*mesh.num_nodes,window_size); @@ -127,6 +136,14 @@ void TLQS3D::execute(SimulationParameters_t& SimulationParamaters, const double max_weight = 50.0; // Threshold to catch exploded weights const double fine_floor = 1e-10; // Residual norm below which Anderson is unsafe + // *********************************************************** + // TEMPORARY ALLOCATION AS ANDERSON IS NOT BEING SET FOR MPI YET + // *********************************************************** + CArrayKokkos curr_anderson_residual(mesh.num_nodes,3); + // *********************************************************** + // TEMPORARY ALLOCATION AS ANDERSON IS NOT BEING SET FOR MPI YET + // *********************************************************** + // QR variables FArrayKokkos Q(3*mesh.num_nodes, max_hist ,"Q"); CArrayKokkos R(max_hist, max_hist, "R"); @@ -148,16 +165,16 @@ void TLQS3D::execute(SimulationParameters_t& SimulationParamaters, // displacement_iter_kp1: result of the CG solve for this Picard iteration G(x_k), // then overwritten with the Anderson-accelerated update x_{k+1}. // Reset to zero before every CG solve. - CArrayKokkos displacement_step(3*mesh.num_nodes); /// current load-step displacement estimate - CArrayKokkos displacement_iter_k(3*mesh.num_nodes); /// x_k (Picard iterate in) - CArrayKokkos displacement_iter_kp1(3*mesh.num_nodes); /// G(x_k) then x_{k+1} + CArrayKokkos displacement_step(mesh.num_nodes,3); /// current load-step displacement estimate + CArrayKokkos displacement_iter_k(mesh.num_nodes,3); /// x_k (Picard iterate in) + CArrayKokkos displacement_iter_kp1(mesh.num_nodes,3); /// G(x_k) then x_{k+1} // variables for chebyshev smoothing - CArrayKokkos D_inv(3 * mesh.num_nodes); - CArrayKokkos zk(3 * mesh.num_nodes); - CArrayKokkos zkp1(3 * mesh.num_nodes); - CArrayKokkos delta_z(3 * mesh.num_nodes); - CArrayKokkos temporary(3 * mesh.num_nodes); + MPICArrayKokkos D_inv(mesh.num_nodes, 3); + MPICArrayKokkos zk(mesh.num_nodes, 3); + MPICArrayKokkos zkp1(mesh.num_nodes, 3); + CArrayKokkos delta_z(mesh.num_nodes, 3); + MPICArrayKokkos temporary(mesh.num_nodes, 3); // Algebraic Multigrid variables @@ -186,6 +203,16 @@ void TLQS3D::execute(SimulationParameters_t& SimulationParamaters, CArrayKokkos P9; CArrayKokkos P10; */ + // setting comm plans to node communication plan for all MPI arrays + auto& node_communication_plan = State.node.displacement.comm_plan_; + rk.initialize_comm_plan(*node_communication_plan); + rkp1.initialize_comm_plan(*node_communication_plan); + D_inv.initialize_comm_plan(*node_communication_plan); + zk.initialize_comm_plan(*node_communication_plan); + zkp1.initialize_comm_plan(*node_communication_plan); + temporary.initialize_comm_plan(*node_communication_plan); + + // Create mesh writer MeshWriter mesh_writer; @@ -386,7 +413,7 @@ void TLQS3D::execute(SimulationParameters_t& SimulationParamaters, get_chebyshev_bounds(alpha, beta, D_inv, K_elem, mesh.num_nodes, mesh.elems_in_node, mesh.num_nodes_in_elem, mesh.nodes_in_elem, - zk, temporary, 15); // Running 15 power iterations + zk, temporary, 15, mesh.num_owned_nodes, mesh.shared_tally_owned_nodes); // Running 15 power iterations // getting r0 = (02F - 01F) - K * displacement_iter_k get_r0(mesh.num_nodes, mesh.elems_in_node, mesh.num_nodes_in_elem, mesh.nodes_in_elem, F_elem, K_elem, displacement_iter_kp1, rk); @@ -398,9 +425,10 @@ void TLQS3D::execute(SimulationParameters_t& SimulationParamaters, // z0 = M_inv * r0, p0 = z0 //get_z0(mesh.num_nodes, mesh.num_nodes_in_elem, mesh.nodes_in_elem, mesh.elems_in_node, K_elem, rk, zk); - FOR_ALL(i, 0, 3 * mesh.num_nodes, { + FOR_ALL(i, 0, (int)mesh.num_nodes, + j, 0, 3, { //zk(i) = M_inv(i) * rk(i); - p(i) = zk(i); + p(i,j) = zk(i,j); }); Kokkos::fence(); @@ -410,17 +438,24 @@ void TLQS3D::execute(SimulationParameters_t& SimulationParamaters, // r_k^T * z_k double rktzk = 0.0; double loc_rktzk = 0.0; - FOR_REDUCE_SUM(i, 0, 3*mesh.num_nodes, loc_rktzk, { - loc_rktzk += rk(i) * zk(i); + FOR_REDUCE_SUM(i, 0, (int)mesh.num_owned_nodes, loc_rktzk, { + if(mesh.shared_tally_owned_nodes(i)) { + for (int j = 0; j < 3; j++) { + loc_rktzk += rk(i,j) * zk(i,j); + } + } }, rktzk); Kokkos::fence(); + MPI_Allreduce(MPI_IN_PLACE, &rktzk, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + // alpha_k = (r_k^T * z_k) / (p_k^T * K * p_k) - double alpha_k = get_alpha(mesh.num_nodes, mesh.num_nodes_in_elem, mesh.nodes_in_elem, K_elem, rktzk, p); + double alpha_k = get_alpha(mesh.num_nodes, mesh.num_nodes_in_elem, mesh.num_owned_nodes, mesh.elems_in_node, mesh.nodes_in_elem, K_elem, rktzk, p, temporary, mesh.shared_tally_owned_nodes); // displacement_iter_kp1 = displacement_iter_kp1 + alpha_k * p_k - FOR_ALL(i, 0, 3*mesh.num_nodes, { - displacement_iter_kp1(i) += alpha_k * p(i); + FOR_ALL(i, 0, (int)mesh.num_nodes, + j, 0, 3, { + displacement_iter_kp1(i, j) += alpha_k * p(i, j); }); Kokkos::fence(); @@ -442,10 +477,15 @@ void TLQS3D::execute(SimulationParameters_t& SimulationParamaters, // check convergence on true residual norm double rkp1trkp1 = 0.0; double loc_rkp1trkp1 = 0.0; - FOR_REDUCE_SUM(i, 0, 3*mesh.num_nodes, loc_rkp1trkp1, { - loc_rkp1trkp1 += rkp1(i) * rkp1(i); + FOR_REDUCE_SUM(i, 0, (int)mesh.num_owned_nodes, loc_rkp1trkp1, { + if(mesh.shared_tally_owned_nodes(i)){ + for (int j = 0; j < 3; j++) { + loc_rkp1trkp1 += rkp1(i,j) * rkp1(i,j); + } + } }, rkp1trkp1); Kokkos::fence(); + MPI_Allreduce(MPI_IN_PLACE, &rkp1trkp1, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); double norm = sqrt(rkp1trkp1); //std::cout << "CGM iter " << cgm_iter << " residual norm: " << norm << "\n"; if (norm < 1.0/*1E-10*/) { @@ -455,23 +495,30 @@ void TLQS3D::execute(SimulationParameters_t& SimulationParamaters, // r_{k+1}^T * z_{k+1} double rkp1tzkp1 = 0.0; double loc_rkp1tzkp1 = 0.0; - FOR_REDUCE_SUM(i, 0, 3*mesh.num_nodes, loc_rkp1tzkp1, { - loc_rkp1tzkp1 += rkp1(i) * zkp1(i); + FOR_REDUCE_SUM(i, 0, (int)mesh.num_owned_nodes, loc_rkp1tzkp1, { + if(mesh.shared_tally_owned_nodes(i)){ + for (int j = 0; j < 3; j++) { + loc_rkp1tzkp1 += rkp1(i,j) * zkp1(i,j); + } + } }, rkp1tzkp1); Kokkos::fence(); + MPI_Allreduce(MPI_IN_PLACE, &rkp1tzkp1, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); // beta_k = (r_{k+1}^T * z_{k+1}) / (r_k^T * z_k) double beta_k = rkp1tzkp1 / (rktzk + 1e-16); // p_{k+1} = z_{k+1} + beta_k * p_k - FOR_ALL(i, 0, 3*mesh.num_nodes, { - p(i) = zkp1(i) + beta_k * p(i); + FOR_ALL(i, 0, (int)mesh.num_nodes, + j, 0, 3, { + p(i,j) = zkp1(i,j) + beta_k * p(i,j); }); // update rk, zk for next iteration - FOR_ALL(i, 0, 3*mesh.num_nodes, { - rk(i) = rkp1(i); - zk(i) = zkp1(i); + FOR_ALL(i, 0, (int)mesh.num_nodes, + j, 0, 3, { + rk(i,j) = rkp1(i,j); + zk(i,j) = zkp1(i,j); }); Kokkos::fence(); @@ -501,12 +548,13 @@ void TLQS3D::execute(SimulationParameters_t& SimulationParamaters, // *************************************************** // --- Step 1: compute Anderson residual f_k = G(x_k) - x_k --- - FOR_ALL(i, 0, 3*mesh.num_nodes, { - curr_anderson_residual(i) = displacement_iter_kp1(i) - displacement_iter_k(i); + FOR_ALL(i, 0, (int)mesh.num_nodes, + j, 0, 3, { + curr_anderson_residual(i,j) = displacement_iter_kp1(i,j) - displacement_iter_k(i,j); }); Kokkos::fence(); - // Compute current residual norm for safeguarding checks + /* // Compute current residual norm for safeguarding checks double safe_norm = 0.0; double safe_loc_norm = 0.0; FOR_REDUCE_SUM(i, 0, 3*mesh.num_nodes, safe_loc_norm, { @@ -559,7 +607,7 @@ void TLQS3D::execute(SimulationParameters_t& SimulationParamaters, break; } } - /* // --- Diagnostic Verification Block --- + // --- Diagnostic Verification Block --- std::vector alpha(m_diff + 1, 0.0); alpha[0] = anderson_weights.host(0); @@ -575,7 +623,7 @@ void TLQS3D::execute(SimulationParameters_t& SimulationParamaters, alpha_sum += a; } std::cout << "\nTotal Alpha Sum (Should be 1.0): " << alpha_sum << std::endl; - // ------------------------------------- */ + // ------------------------------------- if (weights_are_valid) { // --- Step 4a: Apply Accelerated Update --- @@ -601,6 +649,7 @@ void TLQS3D::execute(SimulationParameters_t& SimulationParamaters, } } + */ auto point_E = std::chrono::steady_clock::now(); auto elapsed_E = std::chrono::duration_cast(point_E - point_D).count(); //std::cout << "Time elapsed for anderson: " << elapsed_E << " ms\n"; @@ -609,8 +658,9 @@ void TLQS3D::execute(SimulationParameters_t& SimulationParamaters, // *************************************************** // update displacement step vector for convergence check and next iteration - FOR_ALL(i, 0, 3*mesh.num_nodes, { - displacement_step(i) += displacement_iter_kp1(i); + FOR_ALL(i, 0, (int)mesh.num_nodes, + j, 0, 3, { + displacement_step(i,j) += displacement_iter_kp1(i,j); }); Kokkos::fence(); @@ -630,10 +680,16 @@ void TLQS3D::execute(SimulationParameters_t& SimulationParamaters, double norm = sqrt(norm_num / norm_den); */ double norm = 0.0; double loc_norm = 0.0; - FOR_REDUCE_SUM(i, 0, 3*mesh.num_nodes, loc_norm, { - loc_norm += curr_anderson_residual(i) * curr_anderson_residual(i); + FOR_REDUCE_SUM(i, 0, (int)mesh.num_owned_nodes, loc_norm, { + if(mesh.shared_tally_owned_nodes(i)){ + for (int j = 0; j < 3; j++) { + loc_norm += curr_anderson_residual(i,j) * curr_anderson_residual(i,j); + } + } }, norm); + MPI_Allreduce(MPI_IN_PLACE, &norm, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + std::cout << "ITER: " << iter << " ANDERSON RESIDUAL NORM: " << norm << std::endl; if (norm < 1E-16 && iter > 1) { std::cout << "PICARD CONVERGED AT ITER: " << iter+1 << std::endl; @@ -641,8 +697,9 @@ void TLQS3D::execute(SimulationParameters_t& SimulationParamaters, } // Update x_k <- x_{k+1} for the next Picard iteration. - FOR_ALL(i, 0, 3*mesh.num_nodes, { - displacement_iter_k(i) = displacement_iter_kp1(i); + FOR_ALL(i, 0, (int)mesh.num_nodes, + j, 0, 3, { + displacement_iter_k(i,j) = displacement_iter_kp1(i,j); }); Kokkos::fence(); auto point_F = std::chrono::steady_clock::now(); @@ -655,7 +712,7 @@ void TLQS3D::execute(SimulationParameters_t& SimulationParamaters, // updating total displacement for next load step FOR_ALL(i, 0, (int)mesh.num_nodes, j, 0, 3, { - State.node.displacement(i,j) += displacement_step(3*i + j); + State.node.displacement(i,j) += displacement_step(i, j); State.node.coords(i,j) = State.node.coords_t0(i,j) + State.node.displacement(i,j); }); Kokkos::fence(); diff --git a/apps/multiphysics/src/boundary_conditions/displacement/cyclic_displacement_bc.hpp b/apps/multiphysics/src/boundary_conditions/displacement/cyclic_displacement_bc.hpp index 7d5ca47b9..034bd2417 100644 --- a/apps/multiphysics/src/boundary_conditions/displacement/cyclic_displacement_bc.hpp +++ b/apps/multiphysics/src/boundary_conditions/displacement/cyclic_displacement_bc.hpp @@ -89,7 +89,7 @@ static void displacement(const swage::Mesh_t& mesh, double disp_next = disp_bc_global_vars(bdy_set, 0)*sin((time_value+dt-time_start)*ang_freq); const double total_step_target = disp_next - disp_curr; - displacement_step(3 * bdy_node_gid + constrained_dir) = total_step_target; + displacement_step(bdy_node_gid, constrained_dir) = total_step_target; for (size_t elem_lid = 0; elem_lid < num_elems_in_node; elem_lid++) { const size_t elem_gid = mesh.elems_in_node(bdy_node_gid, elem_lid); diff --git a/apps/multiphysics/src/boundary_conditions/displacement/piston_displacement_bc.hpp b/apps/multiphysics/src/boundary_conditions/displacement/piston_displacement_bc.hpp index e1d791870..40a7cca33 100644 --- a/apps/multiphysics/src/boundary_conditions/displacement/piston_displacement_bc.hpp +++ b/apps/multiphysics/src/boundary_conditions/displacement/piston_displacement_bc.hpp @@ -85,7 +85,7 @@ static void displacement(const swage::Mesh_t& mesh, const double total_step_target = disp_bc_global_vars(bdy_set, 0) * (dt / (time_end - time_start)); - displacement_step(3 * bdy_node_gid + constrained_dir) = total_step_target; + displacement_step(bdy_node_gid, constrained_dir) = total_step_target; for (size_t elem_lid = 0; elem_lid < num_elems_in_node; elem_lid++) { const size_t elem_gid = mesh.elems_in_node(bdy_node_gid, elem_lid); diff --git a/apps/multiphysics/src/boundary_conditions/displacement/total_displacement_bc.hpp b/apps/multiphysics/src/boundary_conditions/displacement/total_displacement_bc.hpp index 746a453e9..e2b7b6d35 100644 --- a/apps/multiphysics/src/boundary_conditions/displacement/total_displacement_bc.hpp +++ b/apps/multiphysics/src/boundary_conditions/displacement/total_displacement_bc.hpp @@ -86,9 +86,9 @@ static void displacement(const swage::Mesh_t& mesh, const double total_step_target_1 = disp_bc_global_vars(bdy_set, 1) * (dt / (time_end - time_start)); const double total_step_target_2 = disp_bc_global_vars(bdy_set, 2) * (dt / (time_end - time_start)); - displacement_step(3 * bdy_node_gid) = total_step_target_0; - displacement_step(3 * bdy_node_gid + 1) = total_step_target_1; - displacement_step(3 * bdy_node_gid + 2) = total_step_target_2; + displacement_step(bdy_node_gid, 0) = total_step_target_0; + displacement_step(bdy_node_gid, 1) = total_step_target_1; + displacement_step(bdy_node_gid, 2) = total_step_target_2; for (size_t elem_lid = 0; elem_lid < num_elems_in_node; elem_lid++) { const size_t elem_gid = mesh.elems_in_node(bdy_node_gid, elem_lid); diff --git a/apps/multiphysics/src/common/include/mesh_io.hpp b/apps/multiphysics/src/common/include/mesh_io.hpp index 556f41146..a6652d992 100644 --- a/apps/multiphysics/src/common/include/mesh_io.hpp +++ b/apps/multiphysics/src/common/include/mesh_io.hpp @@ -5948,6 +5948,8 @@ class MeshWriter num_nodes_in_elem, num_dims, solver_id, + mpi_rank, + mpi_size, mat_den_id, mat_pres_id, mat_sie_id, @@ -5980,6 +5982,38 @@ class MeshWriter // save the graphics time graphics_times(graphics_id) = time_value; + std::vector local_mat_counts(num_mats, 0ULL); + for (size_t mi = 0; mi < num_mats; mi++) { + local_mat_counts[mi] = + static_cast(State.MaterialToMeshMaps.num_mat_elems.host(mi)); + } + std::vector gathered_mat_elems; + if (mpi_rank == 0) { + gathered_mat_elems.assign(static_cast(mpi_size) * num_mats, 0ULL); + } + MPI_Gather(local_mat_counts.data(), + static_cast(num_mats), + MPI_UNSIGNED_LONG_LONG, + mpi_rank == 0 ? gathered_mat_elems.data() : nullptr, + static_cast(num_mats), + MPI_UNSIGNED_LONG_LONG, + 0, + MPI_COMM_WORLD); + + const unsigned long long local_owned_elems_ull = static_cast(n_owned_elems); + std::vector gathered_owned_elems; + if (mpi_rank == 0) { + gathered_owned_elems.assign(static_cast(mpi_size), 0ULL); + } + MPI_Gather(&local_owned_elems_ull, + 1, + MPI_UNSIGNED_LONG_LONG, + mpi_rank == 0 ? gathered_owned_elems.data() : nullptr, + 1, + MPI_UNSIGNED_LONG_LONG, + 0, + MPI_COMM_WORLD); + // check to see if an mesh state was written bool write_mesh_state = false; if( num_elem_scalar_vars > 0 || @@ -5998,24 +6032,40 @@ class MeshWriter write_mat_pt_state = true; } - // call the vtm file writer - std::string mat_fields_name = "mat"; - write_vtm_Pn(graphics_times, - elem_fields_name, - mat_fields_name, - time_value, - graphics_id, - num_mats, - write_mesh_state, - write_mat_pt_state, - solver_id); - - // call the pvd file writer - write_pvd(graphics_times, - time_value, - graphics_id, - solver_id, - mpi_rank); + MPI_Barrier(MPI_COMM_WORLD); + + if (mpi_rank == 0) { + // call the vtm file writer + std::string mat_fields_name = "mat"; + //write_vtm_Pn(graphics_times, + // elem_fields_name, + // mat_fields_name, + // time_value, + // graphics_id, + // num_mats, + // write_mesh_state, + // write_mat_pt_state, + // solver_id); + write_vtm_Pn(graphics_times, + elem_fields_name, + mat_fields_name, + time_value, + graphics_id, + num_mats, + write_mesh_state, + write_mat_pt_state, + solver_id, + mpi_size, + gathered_owned_elems.empty() ? nullptr : gathered_owned_elems.data(), + gathered_mat_elems.empty() ? nullptr : gathered_mat_elems.data()); + + // call the pvd file writer + write_pvd(graphics_times, + time_value, + graphics_id, + solver_id, + mpi_rank); + } // increment graphics id counter @@ -6112,6 +6162,8 @@ class MeshWriter const size_t num_nodes_in_elem, const size_t num_dims, const size_t solver_id, + int mpi_rank, + int mpi_size, // field slot IDs (-1 means "not requested") const int mat_den_id, const int mat_pres_id, @@ -6134,7 +6186,7 @@ class MeshWriter // Number of Gauss points per element: (2*Pn)^num_dims // Read directly from the reference element so we stay consistent with // whatever quadrature order was set up for this run. - const size_t num_gp_per_elem = mesh.num_gauss_in_elem; + const size_t num_gp_per_elem = ref_elem.qpt_grad_basis.dims(0); const size_t num_total_gp = num_mat_elems * num_gp_per_elem; const size_t num_scalar_vars = mat_scalar_var_names.size(); @@ -6145,15 +6197,35 @@ class MeshWriter // ----------------------------------------------------------------------- FILE* fp; - char filename[128]; - int max_len = sizeof filename; + char filename[512]; + int max_len = static_cast(sizeof filename); + int str_output_len; // File lives next to the standard per-material .vtu files. - // The "_Pn" suffix distinguishes this point-cloud dataset from the - // element-averaged dataset written by write_vtu. - snprintf(filename, max_len, - "vtk/data/Fierro.solver%zu.%s.%05d.vtu", - solver_id, partname.c_str(), graphics_id); + // The "_Pn" suffix (already part of partname) distinguishes this + // point-cloud dataset from the element-averaged dataset written by + // write_vtu. Filename convention matches write_vtu exactly so the + // .vtm DataSet references (Fierro.solverN.matM_Pn.FFFFF_rRRRR.vtu) + // resolve correctly under MPI. + if (mpi_size > 1) { + str_output_len = snprintf(filename, + static_cast(max_len), + "vtk/data/Fierro.solver%zu.%s.%05d_r%04d.vtu", + solver_id, + partname.c_str(), + graphics_id, + mpi_rank); + } + else { + str_output_len = snprintf(filename, + static_cast(max_len), + "vtk/data/Fierro.solver%zu.%s.%05d.vtu", + solver_id, + partname.c_str(), + graphics_id); + } + + if (str_output_len >= max_len) { fputs("Filename length exceeded; string truncated", stderr); } fp = fopen(filename, "w"); if (!fp) { @@ -6184,7 +6256,7 @@ class MeshWriter fprintf(fp, " \n"); fprintf(fp, " \n"); - + DCArrayKokkos x_phys(num_mat_elems, num_gp_per_elem, 3); x_phys.set_values(0); FOR_ALL(elem, 0, num_mat_elems, { @@ -6274,10 +6346,6 @@ class MeshWriter // Scalar fields // ------------------------------------------------------------------- - // Helper lambda-style macro to avoid repeating the loop boilerplate. - // Written as an ordinary block; C++11 lambdas would also work if the - // codebase already uses them. - // density if (mat_den_id >= 0) { fprintf(fp, " & graphics_times, - const std::string& elem_part_name, - const std::string& mat_part_name, - double time_value, - int graphics_id, - int num_mats, - bool write_mesh_state, - bool write_mat_pt_state, - const size_t solver_id) + void write_vtm_Pn(CArray& graphics_times, + const std::string& elem_part_name, + const std::string& mat_part_name, + double time_value, + int graphics_id, + size_t num_mats_global, + bool write_mesh_state, + bool write_mat_pt_state, + const size_t solver_id, + int mpi_size, + const unsigned long long* owned_elems_by_rank, + const unsigned long long* mat_elem_counts_by_rank) { - // The _Pn files are named mat_part_name + mat_id + "_Pn" - // e.g. "mat0_Pn" → Fierro.solver0.mat0_Pn.00001.vtu - // This matches the partname constructed in write_mesh: - // std::string pn_name = "mat" + str_mat_val + "_Pn"; - for (int file_id = 0; file_id <= graphics_id; file_id++) { FILE* out[1]; - char filename[100]; - int max_len = sizeof filename; + char filename[512]; + int max_len = static_cast(sizeof filename); int str_output_len; - str_output_len = snprintf(filename, max_len, - "vtk/data/Fierro.solver%zu.%05d.vtm", solver_id, file_id); - if (str_output_len >= max_len) { - fputs("Filename length exceeded; string truncated", stderr); - } + str_output_len = + snprintf(filename, static_cast(max_len), "vtk/data/Fierro.solver%zu.%05d.vtm", solver_id, file_id); + + if (str_output_len >= max_len) { fputs("Filename length exceeded; string truncated", stderr); } out[0] = fopen(filename, "w"); + if (!out[0]) { + std::cerr << "[MeshWriter] Failed to open VTM file: " << filename << std::endl; + continue; + } fprintf(out[0], "\n"); - fprintf(out[0], "\n"); + fprintf(out[0], "\n"); fprintf(out[0], " \n"); size_t block_id = 0; - // ------------------------------------------------------------------- - // Block 0: element-averaged mesh (node + element state) - // ------------------------------------------------------------------- + // --------------------------------------------------------------- + // Block: element-averaged mesh (same schema as write_vtm) + // --------------------------------------------------------------- if (write_mesh_state) { - fprintf(out[0], " \n", - block_id); - block_id++; + int mesh_pieces = 0; + if (mpi_size > 1 && owned_elems_by_rank != nullptr) { + for (int r = 0; r < mpi_size; r++) { + if (owned_elems_by_rank[static_cast(r)] > 0ULL) { + mesh_pieces++; + } + } + } + else if (owned_elems_by_rank != nullptr && owned_elems_by_rank[0] > 0ULL) { + mesh_pieces = 1; + } + else if (owned_elems_by_rank == nullptr) { + mesh_pieces = 1; + } + + if (mesh_pieces > 0) { + fprintf(out[0], " \n", block_id); + block_id++; - fprintf(out[0], " \n"); - fprintf(out[0], " \n", - file_id, solver_id, elem_part_name.c_str(), file_id, - graphics_times(file_id)); - fprintf(out[0], " \n"); + int ds_index = 0; + if (mpi_size > 1 && owned_elems_by_rank != nullptr) { + for (int r = 0; r < mpi_size; r++) { + if (owned_elems_by_rank[static_cast(r)] == 0ULL) { + continue; + } + fprintf(out[0], + " \n", + ds_index, + r, + solver_id, + elem_part_name.c_str(), + file_id, + r); + ds_index++; + } + } + else { + fprintf(out[0], + " \n", + solver_id, + elem_part_name.c_str(), + file_id); + } - fprintf(out[0], " \n"); + fprintf(out[0], " \n"); + } } - // ------------------------------------------------------------------- - // Block 1: per-material data + // --------------------------------------------------------------- + // Block: per-material data // - // Each material gets its own named sub-block so ParaView shows a - // clean tree. Inside each sub-block there are up to two Pieces: - // • "Elements" — element-averaged material fields (matN.vtu) - // • "GaussPoints" — Gauss-point cloud (matN_Pn.vtu) - // ------------------------------------------------------------------- - if (write_mat_pt_state) { - - fprintf(out[0], " \n", - block_id); - block_id++; + // Each material that has data gets its own "MatN" sub-block + // (contiguous mat_block_idx, same reasoning as write_vtm to avoid + // duplicate sibling indices). Inside MatN, each rank that owns + // data for that material contributes TWO DataSet children: + // - "Elements" -> matN.vtu (element-averaged field) + // - "GaussPoints"-> matN_Pn.vtu (Gauss-point cloud) + // Both are direct children of MatN with unique indices, so + // ParaView won't collapse them the way the old layout did. + // --------------------------------------------------------------- + if (write_mat_pt_state && mat_elem_counts_by_rank != nullptr && num_mats_global > 0) { + int mat_pieces = 0; + for (size_t mat_id = 0; mat_id < num_mats_global; mat_id++) { + if (mpi_size > 1) { + for (int r = 0; r < mpi_size; r++) { + if (mat_elem_counts_by_rank[static_cast(r) * num_mats_global + mat_id] > 0ULL) { + mat_pieces++; + } + } + } + else if (mat_elem_counts_by_rank[mat_id] > 0ULL) { + mat_pieces++; + } + } - for (size_t mat_id = 0; mat_id < (size_t)num_mats; mat_id++) { + if (mat_pieces > 0) { + fprintf(out[0], " \n", block_id); + size_t mat_block_idx = 0; - // Open a sub-block for this material - fprintf(out[0], - " \n", - mat_id, mat_id); + for (size_t mat_id = 0; mat_id < num_mats_global; mat_id++) { + int pieces_this_mat = 0; + if (mpi_size > 1) { + for (int r = 0; r < mpi_size; r++) { + if (mat_elem_counts_by_rank[static_cast(r) * num_mats_global + mat_id] > 0ULL) { + pieces_this_mat++; + } + } + } + else if (mat_elem_counts_by_rank[mat_id] > 0ULL) { + pieces_this_mat = 1; + } + if (pieces_this_mat == 0) { + continue; + } - size_t piece_id = 0; + fprintf(out[0], " \n", mat_block_idx, mat_id); - // Piece 0: element-averaged material fields - if (write_mat_pt_state) { - fprintf(out[0], - " \n", - piece_id); - fprintf(out[0], - " \n", - file_id, solver_id, mat_part_name.c_str(), - mat_id, file_id, graphics_times(file_id)); - fprintf(out[0], " \n"); - piece_id++; - } + int ds_index = 0; + if (mpi_size > 1) { + for (int r = 0; r < mpi_size; r++) { + const unsigned long long nm = + mat_elem_counts_by_rank[static_cast(r) * num_mats_global + mat_id]; + if (nm == 0ULL) { + continue; + } - // Piece 1: Gauss-point cloud - if (write_mat_pt_state) { - fprintf(out[0], - " \n", - piece_id); - // filename pattern: mat_part_name + mat_id + "_Pn" - // e.g. Fierro.solver0.mat0_Pn.00001.vtu - fprintf(out[0], - " \n", - file_id, solver_id, mat_part_name.c_str(), - mat_id, file_id, graphics_times(file_id)); - fprintf(out[0], " \n"); - piece_id++; - } + // Elements (element-averaged material field) + fprintf(out[0], + " \n", + ds_index, + mat_id, + r, + solver_id, + mat_part_name.c_str(), + mat_id, + file_id, + r); + ds_index++; - fprintf(out[0], " \n"); // close MatN sub-block + // GaussPoints (Pn point cloud) + fprintf(out[0], + " \n", + ds_index, + mat_id, + r, + solver_id, + mat_part_name.c_str(), + mat_id, + file_id, + r); + ds_index++; + } + } + else { + // Elements + fprintf(out[0], + " \n", + ds_index, + mat_id, + solver_id, + mat_part_name.c_str(), + mat_id, + file_id); + ds_index++; - } // end for mat_id + // GaussPoints + fprintf(out[0], + " \n", + ds_index, + mat_id, + solver_id, + mat_part_name.c_str(), + mat_id, + file_id); + ds_index++; + } - fprintf(out[0], " \n"); // close Mat block + fprintf(out[0], " \n"); + mat_block_idx++; + } - } // end if any material state + fprintf(out[0], " \n"); + } + } fprintf(out[0], " \n"); fprintf(out[0], ""); @@ -7368,6 +7524,10 @@ class MeshWriter } MPI_Barrier(MPI_COMM_WORLD); + int rank; + int world_size; + mesh_io_mpi_detail::query_world_rank_size(rank,world_size); + // ----------------------------------------------------------------------- // FILE 1 – Material-point state // @@ -7411,9 +7571,11 @@ class MeshWriter FILE* out_elem_state; char filename[128]; int max_len = sizeof filename; - snprintf(filename, max_len, - "state/mat_pt_state_t_%6.4e_mat_id_%d.txt", - time_value, mat_id); + if (world_size == 1) { + snprintf(filename, max_len, "state/mat_pt_state_t_%6.4e_mat_id_%d.txt", time_value, mat_id); + } else { + snprintf(filename, max_len, "state/mat_pt_state_t_%6.4e_mat_id_%d_rank_%d.txt", time_value, mat_id, rank); + } out_elem_state = fopen(filename, "w"); if (!out_elem_state) { @@ -7506,10 +7668,6 @@ class MeshWriter } fprintf(out_elem_state, "\n"); - int rank; - int world_size; - mesh_io_mpi_detail::query_world_rank_size(rank,world_size); - // ---- Data rows -------------------------------------------------- // // One row per element. pt_id uses gp=0 since there is exactly @@ -7607,7 +7765,11 @@ class MeshWriter FILE* out_point_state; char filename[128]; int max_len = sizeof filename; - snprintf(filename, max_len, "state/node_state_t_%6.4e.txt", time_value); + if (world_size == 1) { + snprintf(filename, max_len, "state/node_state_t_%6.4e.txt", time_value); + } else { + snprintf(filename, max_len, "state/node_state_t_%6.4e_rank_%d.txt", time_value, rank); + } out_point_state = fopen(filename, "w"); if (!out_point_state) { @@ -7786,7 +7948,7 @@ class MeshWriter // Derived sizes // ----------------------------------------------------------------------- - const size_t num_gp_per_elem = mesh.num_gauss_in_elem; + const size_t num_gp_per_elem = ref_elem.qpt_grad_basis.dims(0); // ----------------------------------------------------------------------- // Gauss-point physical coordinates (isoparametric mapping) @@ -7811,6 +7973,10 @@ class MeshWriter //} MPI_Barrier(MPI_COMM_WORLD); + int rank; + int world_size; + mesh_io_mpi_detail::query_world_rank_size(rank,world_size); + // ----------------------------------------------------------------------- // FILE 1 – Material-point (Gauss-point) state // @@ -7840,7 +8006,11 @@ class MeshWriter FILE* out_elem_state; char filename[128]; int max_len = sizeof filename; - snprintf(filename, max_len, "state/mat_pt_state_t_%6.4e_mat_id_%d.txt", time_value, mat_id); + if (world_size == 1) { + snprintf(filename, max_len, "state/mat_pt_state_t_%6.4e_mat_id_%d.txt", time_value, mat_id); + } else { + snprintf(filename, max_len, "state/mat_pt_state_t_%6.4e_mat_id_%d_rank_%d.txt", time_value, mat_id, rank); + } out_elem_state = fopen(filename, "w"); if (!out_elem_state) { @@ -7925,10 +8095,6 @@ class MeshWriter } fprintf(out_elem_state, "\n"); - int rank; - int world_size; - mesh_io_mpi_detail::query_world_rank_size(rank,world_size); - // ---- Data rows -------------------------------------------------- for (size_t elem = 0; elem < State.MaterialToMeshMaps.num_mat_elems.host(mat_id); elem++) { const size_t elem_rid = @@ -8035,7 +8201,11 @@ class MeshWriter FILE* out_point_state; char filename[128]; int max_len = sizeof filename; - snprintf(filename, max_len, "state/node_state_t_%6.4e.txt", time_value); + if (world_size == 1) { + snprintf(filename, max_len, "state/node_state_t_%6.4e.txt", time_value); + } else { + snprintf(filename, max_len, "state/node_state_t_%6.4e_rank_%d.txt", time_value, rank); + } out_point_state = fopen(filename, "w"); if (!out_point_state) { diff --git a/apps/multiphysics/src/driver.cpp b/apps/multiphysics/src/driver.cpp index 63dfa04df..3f24e0f71 100644 --- a/apps/multiphysics/src/driver.cpp +++ b/apps/multiphysics/src/driver.cpp @@ -151,7 +151,9 @@ void Driver::initialize() // Partition the mesh to all ranks if(world_size != 1) { // pass through the partitioning function if not a single rank - elements::partition_mesh(initial_mesh, mesh, initial_node_coords, final_node_coords, element_communication_plan, node_communication_plan, world_size, rank); + elements::partition_mesh(initial_mesh, mesh, initial_node_coords, final_node_coords, element_communication_plan, node_communication_plan, world_size, rank); + MPI_Allreduce(MPI_IN_PLACE, &initial_mesh.num_gauss_in_elem, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + mesh.num_gauss_in_elem = initial_mesh.num_gauss_in_elem; // Verify communication plans (matches ELEMENTS decomp_example pattern) // element_communication_plan.verify_graph_communicator(); // node_communication_plan.verify_graph_communicator(); @@ -376,6 +378,9 @@ void Driver::initialize() if (State.node.vel_n0.size() > 0) { State.node.vel_n0.initialize_comm_plan(node_communication_plan); } + if (State.node.displacement.size() > 0) { + State.node.displacement.initialize_comm_plan(node_communication_plan); + } if (State.GaussPoints.shock_detector.size() > 0){ State.GaussPoints.shock_detector.initialize_comm_plan(element_communication_plan); }