diff --git a/src/examples/CMakeLists.txt b/src/examples/CMakeLists.txt index a8cbbbe2..f44ae861 100644 --- a/src/examples/CMakeLists.txt +++ b/src/examples/CMakeLists.txt @@ -101,6 +101,7 @@ if ( BUILD_REDECOMP ) set( examples mfem_mortar_lm_patch.cpp mfem_common_plane.cpp + jacobian_transfer_comparison.cpp ) foreach( example ${examples} ) diff --git a/src/examples/jacobian_transfer_comparison.cpp b/src/examples/jacobian_transfer_comparison.cpp new file mode 100644 index 00000000..fc5e210c --- /dev/null +++ b/src/examples/jacobian_transfer_comparison.cpp @@ -0,0 +1,195 @@ +// Copyright (c) 2017-2025, Lawrence Livermore National Security, LLC and +// other Tribol Project Developers. See the top-level LICENSE file for details. +// +// SPDX-License-Identifier: (MIT) + +#include +#include +#include + +#include "tribol/config.hpp" + +#ifdef TRIBOL_USE_UMPIRE +#include "umpire/ResourceManager.hpp" +#endif + +#include "mfem.hpp" + +#include "axom/core.hpp" +#include "axom/slic.hpp" + +#include "shared/mesh/MeshBuilder.hpp" + +#include "tribol/common/Parameters.hpp" +#include "tribol/interface/tribol.hpp" +#include "tribol/interface/mfem_tribol.hpp" +#include "tribol/mesh/CouplingScheme.hpp" +#include "tribol/mesh/MfemData.hpp" +#include "tribol/mesh/CouplingScheme.hpp" + +int main( int argc, char** argv ) +{ + MPI_Init( &argc, &argv ); + int rank, n_ranks; + MPI_Comm_rank( MPI_COMM_WORLD, &rank ); + MPI_Comm_size( MPI_COMM_WORLD, &n_ranks ); + +#ifdef TRIBOL_USE_UMPIRE + umpire::ResourceManager::getInstance(); +#endif + + axom::slic::SimpleLogger logger; + axom::slic::setIsRoot( rank == 0 ); + + int ref_levels = 2; + axom::CLI::App app{ "jacobian_transfer_comparison" }; + app.add_option( "-r,--refine", ref_levels, "Number of times to refine the mesh uniformly." )->capture_default_str(); + CLI11_PARSE( app, argc, argv ); + + // Setup two-block mesh + int nel_per_dir = std::pow( 2, ref_levels ); + auto elem_type = mfem::Element::HEXAHEDRON; + auto mortar_attrs = std::set( { 4 } ); + auto nonmortar_attrs = std::set( { 5 } ); + + mfem::ParMesh mesh = shared::ParMeshBuilder( + MPI_COMM_WORLD, + shared::MeshBuilder::Unify( { shared::MeshBuilder::CubeMesh( nel_per_dir, nel_per_dir, nel_per_dir, elem_type ) + .updateBdrAttrib( 4, 7 ) + .updateBdrAttrib( 6, 4 ), + shared::MeshBuilder::CubeMesh( nel_per_dir, nel_per_dir, nel_per_dir, elem_type ) + .translate( { 0.0, 0.0, 0.99 } ) + .updateBdrAttrib( 8, 5 ) } ) ); + + int dim = mesh.SpaceDimension(); + mfem::H1_FECollection fec( 1, dim ); + mfem::ParFiniteElementSpace fespace( &mesh, &fec, dim ); + mfem::ParGridFunction coords( &fespace ); + mesh.GetNodes( coords ); + + // Register Tribol coupling scheme + int cs_id = 0; + tribol::registerMfemCouplingScheme( cs_id, 0, 1, mesh, coords, mortar_attrs, nonmortar_attrs, + tribol::SURFACE_TO_SURFACE, tribol::NO_SLIDING, tribol::SINGLE_MORTAR, + tribol::FRICTIONLESS, tribol::LAGRANGE_MULTIPLIER, tribol::BINNING_GRID ); + tribol::setMPIComm( cs_id, MPI_COMM_WORLD ); + tribol::updateMfemParallelDecomposition(); + + // Setup MFEM Jacobian data + auto& cs_manager = tribol::CouplingSchemeManager::getInstance(); + auto* cs = cs_manager.findData( cs_id ); + if ( !cs->hasMfemJacobianData() ) { + cs->setMfemJacobianData( std::make_unique( + *cs->getMfemMeshData(), *cs->getMfemSubmeshData(), cs->getContactMethod() ) ); + } + auto* jac_data = cs->getMfemJacobianData(); + jac_data->UpdateJacobianXfer(); + + // Synthesize Jacobian data for comparison + int ne1 = cs->getMfemMeshData()->GetMesh1NE(); + int num_nodes_per_elem = 4; // Quad faces + int num_dofs_per_elem = num_nodes_per_elem * dim; + + // Populate ComputedElementData + tribol::ComputedElementData new_data; + new_data.row_space = tribol::BlockSpace::MORTAR; + new_data.col_space = tribol::BlockSpace::MORTAR; + new_data.jacobian_data.resize( ne1 * num_dofs_per_elem * num_dofs_per_elem ); + new_data.jacobian_offsets.resize( ne1 ); + new_data.row_elem_ids.resize( ne1 ); + new_data.col_elem_ids.resize( ne1 ); + + // Setup MethodData for comparison + if ( cs->getMethodData() == nullptr ) { + cs->allocateMethodData(); + } + auto* method_data = cs->getMethodData(); + SLIC_ASSERT( method_data != nullptr ); + + tribol::ArrayT spaces( { tribol::BlockSpace::MORTAR } ); + method_data->reserveBlockJ( std::move( spaces ), ne1 ); + + for ( int e = 0; e < ne1; ++e ) { + new_data.row_elem_ids[e] = e; + new_data.col_elem_ids[e] = e; + new_data.jacobian_offsets[e] = e * num_dofs_per_elem * num_dofs_per_elem; + + tribol::StackArray, 9> blockJ; + blockJ[0] = tribol::DeviceArray2D( num_dofs_per_elem, num_dofs_per_elem ); + + for ( int i = 0; i < num_dofs_per_elem; ++i ) { + for ( int j = 0; j < num_dofs_per_elem; ++j ) { + double val = static_cast( e + i + j ); + blockJ[0]( i, j ) = val; + new_data.jacobian_data[new_data.jacobian_offsets[e] + i + j * num_dofs_per_elem] = val; + } + } + tribol::ArrayT ids( { e } ); + method_data->storeElemBlockJ( std::move( ids ), blockJ ); + } + + // Time and assemble with old method + auto start_old = std::chrono::high_resolution_clock::now(); + + auto xfer = cs->getMfemJacobianData()->GetMfemBlockJacobian( *method_data, { { 0, tribol::BlockSpace::MORTAR } }, + { { 0, tribol::BlockSpace::MORTAR } } ); + + auto end_old = std::chrono::high_resolution_clock::now(); + + // Time and assemble with new method + std::vector contribs_vec; + if ( ne1 > 0 ) { + contribs_vec.push_back( std::move( new_data ) ); + } + + auto start_new = std::chrono::high_resolution_clock::now(); + auto par_J_new = jac_data->GetMfemJacobian( contribs_vec ); + auto end_new = std::chrono::high_resolution_clock::now(); + + // Verify results + auto* old_hypre = dynamic_cast( &xfer->GetBlock( 0, 0 ) ); + auto* new_hypre = &par_J_new.get(); + + // Check difference + shared::ParSparseMat diff_psm = shared::ParSparseMatView( old_hypre ) - shared::ParSparseMatView( new_hypre ); + + // Verify match by checking max error + double max_err = 0.0; + HYPRE_ParCSRMatrix diff_csr = diff_psm.get(); + hypre_ParCSRMatrix* diff_parcsr = (hypre_ParCSRMatrix*)diff_csr; + hypre_CSRMatrix* diag = hypre_ParCSRMatrixDiag( diff_parcsr ); + double* data = hypre_CSRMatrixData( diag ); + int num_nonzeros = hypre_CSRMatrixNumNonzeros( diag ); + for ( int i = 0; i < num_nonzeros; ++i ) { + max_err = std::max( max_err, std::abs( data[i] ) ); + } + hypre_CSRMatrix* offd = hypre_ParCSRMatrixOffd( diff_parcsr ); + data = hypre_CSRMatrixData( offd ); + num_nonzeros = hypre_CSRMatrixNumNonzeros( offd ); + for ( int i = 0; i < num_nonzeros; ++i ) { + max_err = std::max( max_err, std::abs( data[i] ) ); + } + + double global_max_err = 0.0; + MPI_Allreduce( &max_err, &global_max_err, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD ); + + if ( rank == 0 ) { + std::cout << "Old method time: " + << std::chrono::duration_cast( end_old - start_old ).count() << " us" + << std::endl; + std::cout << "New method time: " + << std::chrono::duration_cast( end_new - start_new ).count() << " us" + << std::endl; + std::cout << "Matrix difference max err: " << global_max_err << std::endl; + } + + if ( global_max_err > 1e-12 ) { + SLIC_ERROR_ROOT( "Matrices do not match!" ); + } else { + SLIC_INFO_ROOT( "Verification successful: Matrices match." ); + } + + tribol::finalize(); + MPI_Finalize(); + return 0; +} diff --git a/src/redecomp/transfer/MatrixTransfer.cpp b/src/redecomp/transfer/MatrixTransfer.cpp index e9a05c2d..cf5b29bc 100644 --- a/src/redecomp/transfer/MatrixTransfer.cpp +++ b/src/redecomp/transfer/MatrixTransfer.cpp @@ -11,6 +11,7 @@ #include "redecomp/RedecompMesh.hpp" #include "shared/math/ParSparseMat.hpp" +#include "shared/infrastructure/Profiling.hpp" namespace redecomp { @@ -52,6 +53,7 @@ mfem::SparseMatrix MatrixTransfer::TransferToParallelSparse( const axom::Array& trial_elem_idx, const axom::Array& src_elem_mat ) const { + TRIBOL_MARK_SCOPE( "TransferToParallelSparse_Method1" ); // TODO (EBC): we need a SparseMatrix-like data structure that allows HYPRE_BigInt on the columns auto parentJ = mfem::SparseMatrix( parent_test_fes_.GetVSize(), parent_trial_fes_.GlobalVSize() ); @@ -77,6 +79,7 @@ mfem::SparseMatrix MatrixTransfer::TransferToParallelSparse( const axom::Array( redecomp_test_fes_.GetMesh() ); auto trial_redecomp = dynamic_cast( redecomp_trial_fes_.GetMesh() ); + TRIBOL_MARK_BEGIN( "BuildCommunicationData" ); // List of entries in src_elem_mat that belong on each parent test space rank. // This is needed so we know which rank to send entries in src_elem_mat to. auto send_array_ids = buildSendArrayIDs( test_elem_idx ); @@ -98,7 +101,9 @@ mfem::SparseMatrix MatrixTransfer::TransferToParallelSparse( const axom::Array>(), @@ -140,10 +145,204 @@ mfem::SparseMatrix MatrixTransfer::TransferToParallelSparse( const axom::Array& test_elem_idx, + const axom::Array& trial_elem_idx, + const axom::Array& src_elem_mat_data, + const axom::Array& src_elem_mat_offsets, + bool parallel_assemble ) const +{ + TRIBOL_MARK_FUNCTION; + // verify inputs + SLIC_ERROR_IF( test_elem_idx.size() != trial_elem_idx.size() || test_elem_idx.size() != src_elem_mat_offsets.size(), + "Element index arrays and element Jacobian offsets array must be the same size." ); + for ( int i{ 0 }; i < test_elem_idx.size(); ++i ) { + auto test_e = test_elem_idx[i]; + auto trial_e = trial_elem_idx[i]; + + SLIC_ERROR_IF( test_e < 0, "Invalid primary index value." ); + SLIC_ERROR_IF( trial_e < 0, "Invalid secondary index value." ); + + auto n_test_elem_vdofs = redecomp_test_fes_.GetFE( test_e )->GetDof() * redecomp_test_fes_.GetVDim(); + auto n_trial_elem_vdofs = redecomp_trial_fes_.GetFE( trial_e )->GetDof() * redecomp_trial_fes_.GetVDim(); + auto expected_size = n_test_elem_vdofs * n_trial_elem_vdofs; + + // Check that we don't go out of bounds of the data array + SLIC_ERROR_IF( src_elem_mat_offsets[i] < 0 || src_elem_mat_offsets[i] + expected_size > src_elem_mat_data.size(), + "Matrix offset and size exceeds data array bounds." ); + } + + auto test_redecomp = dynamic_cast( redecomp_test_fes_.GetMesh() ); + auto trial_redecomp = dynamic_cast( redecomp_trial_fes_.GetMesh() ); + + TRIBOL_MARK_BEGIN( "BuildCommunicationData" ); + // List of entries in src_elem_mat that belong on each parent test space rank. + // This is needed so we know which rank to send entries in src_elem_mat to. + auto send_array_ids = buildSendArrayIDs( test_elem_idx ); + // Number of matrix entries to be sent to each parent test space rank. This + // is used to size the array of element matrix values to be sent to other ranks. + auto send_num_mat_entries = buildSendNumMatEntries( test_elem_idx, trial_elem_idx ); + // Number of test and trial vdofs received from test space redecomp ranks. + // This is used to determine the beginning and end of each element matrix + // received as a single array of values and the beginning and end of vdof indices + // received as a single array of values. + auto recv_mat_sizes = buildRecvMatSizes( test_elem_idx, trial_elem_idx ); + // List of test element offsets received from test space redecomp ranks. The + // offset is used with the parent to redecomp map (and the redecomp rank + // received from) to determine the parent element ID. Parent element ID is + // used to determine the test vldofs of the element matrix entries received + // from redecomp ranks. + auto recv_test_elem_offsets = buildRecvTestElemOffsets( *test_redecomp, test_elem_idx ); + // List of trial element global vdofs corresponding to the element matrix + // entries received from redecomp ranks. The second column of recv_mat_sizes + // determines the offset for each trial element. + auto recv_trial_elem_dofs = buildRecvTrialElemDofs( *trial_redecomp, test_elem_idx, trial_elem_idx ); + TRIBOL_MARK_END( "BuildCommunicationData" ); + + TRIBOL_MARK_BEGIN( "SetupTriplets" ); + // Intermediate storage for triplets + struct Triplet { + int row; + HYPRE_BigInt col; + double val; + }; + std::vector triplets; + // Estimate capacity to reduce reallocations + int total_recv_entries = 0; + for ( int src = 0; src < getMPIUtility().NRanks(); ++src ) { + for ( int e = 0; e < recv_mat_sizes[src].shape()[0]; ++e ) { + total_recv_entries += recv_mat_sizes[src]( e, 0 ) * recv_mat_sizes[src]( e, 1 ); + } + } + triplets.reserve( total_recv_entries ); + TRIBOL_MARK_END( "SetupTriplets" ); + + TRIBOL_MARK_BEGIN( "SendRecvEach" ); + // aggregate dense matrix values, send and assemble + getMPIUtility().SendRecvEach( + type>(), + [&send_array_ids, &send_num_mat_entries, &src_elem_mat_data, &src_elem_mat_offsets, &test_elem_idx, + &trial_elem_idx, this]( axom::IndexType dst ) { + auto send_vals = axom::Array( 0, send_num_mat_entries[dst] ); + + for ( auto src_array_idx : send_array_ids[dst] ) { + // Calculate size from FE spaces + auto test_e = test_elem_idx[src_array_idx]; + auto trial_e = trial_elem_idx[src_array_idx]; + auto n_test_elem_vdofs = redecomp_test_fes_.GetFE( test_e )->GetDof() * redecomp_test_fes_.GetVDim(); + auto n_trial_elem_vdofs = redecomp_trial_fes_.GetFE( trial_e )->GetDof() * redecomp_trial_fes_.GetVDim(); + auto size = n_test_elem_vdofs * n_trial_elem_vdofs; + + send_vals.append( + axom::ArrayView( &src_elem_mat_data[src_elem_mat_offsets[src_array_idx]], size ) ); + } + + return send_vals; + }, + [this, test_redecomp, &triplets, &recv_mat_sizes, &recv_trial_elem_dofs, &recv_test_elem_offsets]( + axom::Array&& send_vals, axom::IndexType src ) { + if ( recv_trial_elem_dofs[src].empty() ) { + return; + } + auto trial_dof_ct = 0; + auto dof_ct = 0; + // element loop + for ( int e{ 0 }; e < recv_test_elem_offsets[src].size(); ++e ) { + auto test_elem_id = test_redecomp->getParentToRedecompElems().first[src][recv_test_elem_offsets[src][e]]; + auto test_elem_dofs = mfem::Array(); + parent_test_fes_.GetElementVDofs( test_elem_id, test_elem_dofs ); + auto trial_elem_dofs = + mfem::Array( &recv_trial_elem_dofs[src][trial_dof_ct], recv_mat_sizes[src]( e, 1 ) ); + // trial loop + for ( int j{ 0 }; j < trial_elem_dofs.Size(); ++j ) { + // test loop + for ( int i{ 0 }; i < test_elem_dofs.Size(); ++i ) { + triplets.push_back( { test_elem_dofs[i], trial_elem_dofs[j], + // send_vals comes from mfem::SparseMatrix (column major) + send_vals[dof_ct + i + j * test_elem_dofs.Size()] } ); + } + } + trial_dof_ct += recv_mat_sizes[src]( e, 1 ); + dof_ct += recv_mat_sizes[src]( e, 0 ) * recv_mat_sizes[src]( e, 1 ); + } + } ); + TRIBOL_MARK_END( "SendRecvEach" ); + + TRIBOL_MARK_BEGIN( "SortTriplets" ); + // Sort triplets by row then column + std::sort( triplets.begin(), triplets.end(), []( const Triplet& a, const Triplet& b ) { + if ( a.row != b.row ) return a.row < b.row; + return a.col < b.col; + } ); + TRIBOL_MARK_END( "SortTriplets" ); + + TRIBOL_MARK_BEGIN( "CSRConstruction" ); + // Count non-zeros and merge duplicates + int num_unique_nonzeros = 0; + if ( !triplets.empty() ) { + num_unique_nonzeros = 1; + for ( size_t i = 1; i < triplets.size(); ++i ) { + if ( triplets[i].row != triplets[i - 1].row || triplets[i].col != triplets[i - 1].col ) { + num_unique_nonzeros++; + } + } + } + + auto num_rows = parent_test_fes_.GetVSize(); + int* I_ptr = new int[num_rows + 1]; + HYPRE_BigInt* J_ptr = new HYPRE_BigInt[num_unique_nonzeros]; + double* data_ptr = new double[num_unique_nonzeros]; + + // Initialize I_ptr with zeros + for ( int i = 0; i <= num_rows; ++i ) { + I_ptr[i] = 0; + } + + if ( !triplets.empty() ) { + int unique_idx = 0; + J_ptr[0] = triplets[0].col; + data_ptr[0] = triplets[0].val; + I_ptr[triplets[0].row + 1]++; + + for ( size_t i = 1; i < triplets.size(); ++i ) { + if ( triplets[i].row == triplets[i - 1].row && triplets[i].col == triplets[i - 1].col ) { + data_ptr[unique_idx] += triplets[i].val; + } else { + unique_idx++; + J_ptr[unique_idx] = triplets[i].col; + data_ptr[unique_idx] = triplets[i].val; + I_ptr[triplets[i].row + 1]++; + } + } + } + + // Transform I_ptr from counts to offsets (prefix sum) + for ( int i = 0; i < num_rows; ++i ) { + I_ptr[i + 1] += I_ptr[i]; + } + TRIBOL_MARK_END( "CSRConstruction" ); + + // Construct rectangular HypreParMatrix + shared::ParSparseMat J_full( getMPIUtility().MPIComm(), num_rows, parent_test_fes_.GlobalVSize(), + parent_trial_fes_.GlobalVSize(), I_ptr, J_ptr, data_ptr, + parent_test_fes_.GetDofOffsets(), parent_trial_fes_.GetDofOffsets() ); + + if ( !parallel_assemble ) { + return J_full; + } else { + auto P_test = parent_test_fes_.Dof_TrueDof_Matrix(); + P_test->HostRead(); + auto P_trial = parent_trial_fes_.Dof_TrueDof_Matrix(); + P_trial->HostRead(); + auto J_true = shared::ParSparseMat::RAP( P_test, J_full, P_trial ); + return J_true; + } +} + shared::ParSparseMat MatrixTransfer::ConvertToParSparseMat( mfem::SparseMatrix&& sparse, bool parallel_assemble ) const { SLIC_ERROR_IF( sparse.Height() != parent_test_fes_.GetVSize(), diff --git a/src/redecomp/transfer/MatrixTransfer.hpp b/src/redecomp/transfer/MatrixTransfer.hpp index d73b5bd7..ad4c6a1e 100644 --- a/src/redecomp/transfer/MatrixTransfer.hpp +++ b/src/redecomp/transfer/MatrixTransfer.hpp @@ -86,6 +86,26 @@ class MatrixTransfer { const axom::Array& trial_elem_idx, const axom::Array& src_elem_mat ) const; + /** + * @brief Transfers element RedecompMesh matrices to parent mfem::ParMesh + * + * @param test_elem_idx List of element IDs on the redecomp test space + * @param trial_elem_idx List of element IDs on the redecomp trial space + * @param src_elem_mat_data Flattened array of element-level dense matrices from the redecomp mesh + * @param src_elem_mat_offsets Offsets into src_elem_mat_data for each element + * @param parallel_assemble Performs parallel assembly (transforms to tdofs) + * on the HypreParMatrix if true, returns ldofs otherwise + * @return mfem::HypreParMatrix on the parent mesh (ldofs on the rows, global + * ldofs on the columns) in rectangular format + * + * @note This method constructs the parallel matrix directly, bypassing mfem::SparseMatrix. + */ + shared::ParSparseMat TransferToParallel( const axom::Array& test_elem_idx, + const axom::Array& trial_elem_idx, + const axom::Array& src_elem_mat_data, + const axom::Array& src_elem_mat_offsets, + bool parallel_assemble = true ) const; + /** * @brief Converts SparseMatrix from TransferToParallelSparse to HypreParMatrix * diff --git a/src/shared/math/ParSparseMat.cpp b/src/shared/math/ParSparseMat.cpp index aa0655f7..b0747c15 100644 --- a/src/shared/math/ParSparseMat.cpp +++ b/src/shared/math/ParSparseMat.cpp @@ -128,6 +128,23 @@ ParSparseMat::ParSparseMat( MPI_Comm comm, HYPRE_BigInt glob_size, HYPRE_BigInt* owned_mat_->SetOwnerFlags( mfem_owned_host_flag, owned_mat_->OwnsOffd(), owned_mat_->OwnsColMap() ); } +ParSparseMat::ParSparseMat( MPI_Comm comm, HYPRE_BigInt global_num_rows, HYPRE_BigInt global_num_cols, + HYPRE_BigInt* row_starts, HYPRE_BigInt* col_starts, mfem::SparseMatrix&& diag ) + : ParSparseMatView( nullptr ) +{ + owned_mat_.reset( createHypreParMatrix( [&]() { + return new mfem::HypreParMatrix( comm, global_num_rows, global_num_cols, row_starts, col_starts, &diag ); + } ) ); + mat_ = owned_mat_.get(); + diag.GetMemoryI().ClearOwnerFlags(); + diag.GetMemoryJ().ClearOwnerFlags(); + diag.GetMemoryData().ClearOwnerFlags(); + // The mfem::Memory in mfem::SparseMatrix allocates using operator new [], so mark the diag memory as owned by MFEM so + // it can be deleted correctly + constexpr int mfem_owned_host_flag = 3; + owned_mat_->SetOwnerFlags( mfem_owned_host_flag, owned_mat_->OwnsOffd(), owned_mat_->OwnsColMap() ); +} + ParSparseMat::ParSparseMat( ParSparseMat&& other ) noexcept : ParSparseMatView( other.owned_mat_.get() ), owned_mat_( std::move( other.owned_mat_ ) ) { diff --git a/src/shared/math/ParSparseMat.hpp b/src/shared/math/ParSparseMat.hpp index 93a59f68..0b8171be 100644 --- a/src/shared/math/ParSparseMat.hpp +++ b/src/shared/math/ParSparseMat.hpp @@ -246,6 +246,21 @@ class ParSparseMat : public ParSparseMatView { */ ParSparseMat( MPI_Comm comm, HYPRE_BigInt glob_size, HYPRE_BigInt* row_starts, mfem::SparseMatrix&& diag ); + /** + * @brief Construct from MPI communicator, global size, row/column starts, and mfem::SparseMatrix rvalue + * + * @param comm MPI communicator + * @param global_num_rows Global number of rows + * @param global_num_cols Global number of columns + * @param row_starts Global row partitioning + * @param col_starts Global column partitioning + * @param diag Local diagonal block SparseMatrix (rvalue) + * + * @note The HypreParMatrix will take ownership of the I, J, and Data from diag. + */ + ParSparseMat( MPI_Comm comm, HYPRE_BigInt global_num_rows, HYPRE_BigInt global_num_cols, HYPRE_BigInt* row_starts, + HYPRE_BigInt* col_starts, mfem::SparseMatrix&& diag ); + /// Template constructor forwarding arguments to mfem::HypreParMatrix constructor template explicit ParSparseMat( Args&&... args ) diff --git a/src/tests/CMakeLists.txt b/src/tests/CMakeLists.txt index 90fe6f37..c5878a58 100644 --- a/src/tests/CMakeLists.txt +++ b/src/tests/CMakeLists.txt @@ -149,6 +149,7 @@ if ( BUILD_REDECOMP AND TRIBOL_USE_MPI ) set( combined_tests tribol_mfem_common_plane.cpp tribol_mfem_mortar_lm.cpp + tribol_mfem_jacobian.cpp tribol_proximity_check.cpp tribol_redecomp_tol.cpp ) diff --git a/src/tests/tribol_mfem_jacobian.cpp b/src/tests/tribol_mfem_jacobian.cpp new file mode 100644 index 00000000..54bb2869 --- /dev/null +++ b/src/tests/tribol_mfem_jacobian.cpp @@ -0,0 +1,161 @@ +// Copyright (c) 2017-2025, Lawrence Livermore National Security, LLC and +// other Tribol Project Developers. See the top-level LICENSE file for details. +// +// SPDX-License-Identifier: (MIT) + +#include "tribol/config.hpp" + +#include + +#ifdef TRIBOL_USE_UMPIRE +#include "umpire/ResourceManager.hpp" +#endif + +#include "mfem.hpp" + +#include "axom/slic.hpp" + +#include "shared/mesh/MeshBuilder.hpp" + +#include "tribol/common/Parameters.hpp" +#include "tribol/interface/tribol.hpp" +#include "tribol/interface/mfem_tribol.hpp" +#include "tribol/mesh/CouplingScheme.hpp" +#include "tribol/mesh/MfemData.hpp" + +class MfemJacobianTest : public testing::Test { + protected: + void SetUp() override {} + void TearDown() override {} +}; + +TEST_F( MfemJacobianTest, direct_jacobian_assembly ) +{ + int n_ranks; + MPI_Comm_size( MPI_COMM_WORLD, &n_ranks ); + + // Setup simple two-cube mesh + int ref_levels = 0; + int nel_per_dir = std::pow( 2, ref_levels ); + + auto mortar_attrs = std::set( { 4 } ); + auto nonmortar_attrs = std::set( { 5 } ); + + // clang-format off + mfem::ParMesh mesh = shared::ParMeshBuilder(MPI_COMM_WORLD, shared::MeshBuilder::Unify({ + shared::MeshBuilder::CubeMesh(nel_per_dir, nel_per_dir, nel_per_dir) + .updateBdrAttrib(3, 7) + .updateBdrAttrib(1, 3) + .updateBdrAttrib(4, 7) // Mortar + .updateBdrAttrib(5, 1) + .updateBdrAttrib(6, 4), + shared::MeshBuilder::CubeMesh(nel_per_dir, nel_per_dir, nel_per_dir) + .translate({0.0, 0.0, 0.99}) // Slight overlap + .updateBdrAttrib(1, 8) + .updateBdrAttrib(3, 7) + .updateBdrAttrib(4, 7) + .updateBdrAttrib(5, 1) // Nonmortar + .updateBdrAttrib(8, 5) + })); + // clang-format on + + int dim = mesh.SpaceDimension(); + int order = 1; + mfem::H1_FECollection fe_coll( order, dim ); + mfem::ParFiniteElementSpace par_fe_space( &mesh, &fe_coll, dim ); + mfem::ParGridFunction coords( &par_fe_space ); + mesh.GetNodes( coords ); + + // Register coupling scheme + int cs_id = 0; + int mesh1_id = 0; + int mesh2_id = 1; + tribol::registerMfemCouplingScheme( cs_id, mesh1_id, mesh2_id, mesh, coords, mortar_attrs, nonmortar_attrs, + tribol::SURFACE_TO_SURFACE, tribol::NO_SLIDING, tribol::SINGLE_MORTAR, + tribol::FRICTIONLESS, tribol::LAGRANGE_MULTIPLIER, tribol::BINNING_GRID ); + + // Build internal MfemData + tribol::updateMfemParallelDecomposition(); + + // Access CouplingScheme and MfemJacobianData + auto& cs_manager = tribol::CouplingSchemeManager::getInstance(); + auto* cs = cs_manager.findData( cs_id ); + ASSERT_NE( cs, nullptr ); + + auto* jac_data = cs->getMfemJacobianData(); + if ( jac_data == nullptr ) { + auto* mesh_data = cs->getMfemMeshData(); + auto* submesh_data = cs->getMfemSubmeshData(); + ASSERT_NE( mesh_data, nullptr ); + ASSERT_NE( submesh_data, nullptr ); + + cs->setMfemJacobianData( + std::make_unique( *mesh_data, *submesh_data, cs->getContactMethod() ) ); + jac_data = cs->getMfemJacobianData(); + } + ASSERT_NE( jac_data, nullptr ); + + jac_data->UpdateJacobianXfer(); + + // Synthesize ComputedElementData for mortar-mortar block + std::vector contributions; + + int num_dofs_per_elem = 8 * dim; // Hex element, 8 nodes + int mat_size = num_dofs_per_elem * num_dofs_per_elem; + + tribol::ComputedElementData contrib; + contrib.row_space = tribol::BlockSpace::MORTAR; + contrib.col_space = tribol::BlockSpace::MORTAR; + + auto* mesh_data = cs->getMfemMeshData(); + + if ( mesh_data && mesh_data->GetMesh1NE() > 0 ) { + contrib.row_elem_ids.push_back( 0 ); + contrib.col_elem_ids.push_back( 0 ); + + contrib.jacobian_data.resize( mat_size ); + for ( int i = 0; i < mat_size; ++i ) contrib.jacobian_data[i] = 1.0; + + contrib.jacobian_offsets.push_back( 0 ); + + contributions.push_back( contrib ); + } + + // Assemble Jacobian + auto ParJ = jac_data->GetMfemJacobian( contributions ); + + // Verify non-zeros if contributions were added + int local_contrib = contributions.size(); + int global_contrib = 0; + MPI_Allreduce( &local_contrib, &global_contrib, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD ); + + if ( global_contrib > 0 ) { + EXPECT_GT( ParJ->NNZ(), 0 ); + } else { + if ( n_ranks == 1 ) { + FAIL() << "No surface elements found on mesh 1 in serial run."; + } + } +} + +int main( int argc, char* argv[] ) +{ + int result = 0; + + MPI_Init( &argc, &argv ); + + ::testing::InitGoogleTest( &argc, argv ); + +#ifdef TRIBOL_USE_UMPIRE + umpire::ResourceManager::getInstance(); +#endif + + axom::slic::SimpleLogger logger; + + result = RUN_ALL_TESTS(); + + tribol::finalize(); + MPI_Finalize(); + + return result; +} diff --git a/src/tests/tribol_mortar_jacobian.cpp b/src/tests/tribol_mortar_jacobian.cpp index 1ac730ab..cb29da30 100644 --- a/src/tests/tribol_mortar_jacobian.cpp +++ b/src/tests/tribol_mortar_jacobian.cpp @@ -3,31 +3,22 @@ // // SPDX-License-Identifier: (MIT) -// Tribol includes -#include "tribol/interface/tribol.hpp" -#include "tribol/common/Parameters.hpp" +#include +#include +#include -// Axom includes -#include "axom/slic.hpp" +#include "gtest/gtest.h" -// MFEM includes #include "mfem.hpp" -// gtest includes -#include "gtest/gtest.h" +#include "axom/slic.hpp" -// c++ includes -#include -#include -#include +#include "tribol/interface/tribol.hpp" +#include "tribol/common/Parameters.hpp" using RealT = tribol::RealT; -/*! - * Test fixture class with some setup necessary to test - * and compute the Jacobian matrix for an implicit - * mortar method with Lagrange multipliers - */ +// Test fixture for implicit mortar method Jacobian class MortarJacTest : public ::testing::Test { public: int numNodes; @@ -44,17 +35,12 @@ class MortarJacTest : public ::testing::Test { void setupTribol( tribol::IndexT* conn1, tribol::IndexT* conn2, tribol::ContactMethod method ) { - // Note, this assumes that numNodes is the total number of - // nodes encompassing the two meshes that will be registered - // with tribol, and that the conn1 and conn2 connectivity arrays - // reflect a global, contiguous index space - - // grab coordinate data + // Get coordinates RealT* x = this->x; RealT* y = this->y; RealT* z = this->z; - // register the mesh with tribol + // Register Tribol mesh int cellType = static_cast( tribol::UNDEFINED_ELEMENT ); switch ( this->numNodesPerFace ) { case 4: { @@ -69,7 +55,6 @@ class MortarJacTest : public ::testing::Test { const int mortarMeshId = 0; const int nonmortarMeshId = 1; - // register mesh tribol::registerMesh( mortarMeshId, 1, this->numNodes, conn1, cellType, x, y, z, tribol::MemorySpace::Host ); tribol::registerMesh( nonmortarMeshId, 1, this->numNodes, conn2, cellType, x, y, z, tribol::MemorySpace::Host ); @@ -81,7 +66,7 @@ class MortarJacTest : public ::testing::Test { tribol::Array1D fy2( this->numNodes ); tribol::Array1D fz2( this->numNodes ); - // initialize force arrays + // Initialize forces for ( int i = 0; i < this->numNodes; ++i ) { fx1[i] = 0.; fy1[i] = 0.; @@ -100,18 +85,17 @@ class MortarJacTest : public ::testing::Test { pressures = tribol::ArrayT( this->numNodes, this->numNodes ); // length of total mesh to use global connectivity to index - // initialize gaps and pressures. Initialize all - // nonmortar pressures to 1.0 + // Initialize gaps and pressures for ( int i = 0; i < this->numNodes; ++i ) { gaps[i] = 0.; pressures[i] = 1.; } - // register nodal gaps and pressures + // Register gaps and pressures tribol::registerMortarGaps( nonmortarMeshId, gaps.data() ); tribol::registerMortarPressures( nonmortarMeshId, pressures.data() ); - // register coupling scheme + // Register coupling scheme const int csIndex = 0; tribol::registerCouplingScheme( csIndex, mortarMeshId, nonmortarMeshId, tribol::SURFACE_TO_SURFACE, tribol::NO_CASE, method, tribol::FRICTIONLESS, tribol::LAGRANGE_MULTIPLIER, diff --git a/src/tribol/common/Parameters.hpp b/src/tribol/common/Parameters.hpp index c944a135..d5dd2e6b 100644 --- a/src/tribol/common/Parameters.hpp +++ b/src/tribol/common/Parameters.hpp @@ -32,12 +32,12 @@ constexpr int ANY_MESH = -1; */ enum class LoggingLevel { - UNDEFINED, ///! Undefined - DEBUG, ///! Debug and higher - INFO, ///! Info and higher - WARNING, ///! Warning and higher - ERROR, ///! Errors only - NUM_LOGGING_LEVELS = ERROR + Undefined, ///! Undefined + Debug, ///! Debug and higher + Info, ///! Info and higher + Warning, ///! Warning and higher + Error, ///! Errors only + NumLoggingLevels = Error }; /*! diff --git a/src/tribol/interface/tribol.cpp b/src/tribol/interface/tribol.cpp index f31f1f69..bad3cbe3 100644 --- a/src/tribol/interface/tribol.cpp +++ b/src/tribol/interface/tribol.cpp @@ -310,9 +310,9 @@ void setLoggingLevel( IndexT cs_id, LoggingLevel log_level ) SLIC_ERROR_IF( !cs, "tribol::setLoggingLevel(): " << "invalid CouplingScheme id." ); - if ( !in_range( static_cast( log_level ), static_cast( LoggingLevel::NUM_LOGGING_LEVELS ) ) ) { + if ( !in_range( static_cast( log_level ), static_cast( LoggingLevel::NumLoggingLevels ) ) ) { SLIC_INFO_ROOT( "tribol::setLoggingLevel(): Logging level not an option; " << "using 'warning' level." ); - cs->setLoggingLevel( LoggingLevel::WARNING ); + cs->setLoggingLevel( LoggingLevel::Warning ); } else { cs->setLoggingLevel( log_level ); } diff --git a/src/tribol/mesh/CouplingScheme.cpp b/src/tribol/mesh/CouplingScheme.cpp index c7c3bdce..c98a0d22 100644 --- a/src/tribol/mesh/CouplingScheme.cpp +++ b/src/tribol/mesh/CouplingScheme.cpp @@ -340,7 +340,7 @@ CouplingScheme::CouplingScheme( IndexT cs_id, IndexT mesh_id1, IndexT mesh_id2, m_couplingSchemeInfo.cs_case_info = NO_CASE_INFO; m_couplingSchemeInfo.cs_enforcement_info = NO_ENFORCEMENT_INFO; - m_loggingLevel = LoggingLevel::UNDEFINED; + m_loggingLevel = LoggingLevel::Undefined; } // end CouplingScheme::CouplingScheme() @@ -1176,21 +1176,21 @@ bool CouplingScheme::init() void CouplingScheme::setSlicLoggingLevel() { // set slic logging level for coupling schemes that have API modified logging levels - if ( this->m_loggingLevel != LoggingLevel::UNDEFINED ) { + if ( this->m_loggingLevel != LoggingLevel::Undefined ) { switch ( this->m_loggingLevel ) { - case LoggingLevel::DEBUG: { + case LoggingLevel::Debug: { axom::slic::setLoggingMsgLevel( axom::slic::message::Debug ); break; } - case LoggingLevel::INFO: { + case LoggingLevel::Info: { axom::slic::setLoggingMsgLevel( axom::slic::message::Info ); break; } - case LoggingLevel::WARNING: { + case LoggingLevel::Warning: { axom::slic::setLoggingMsgLevel( axom::slic::message::Warning ); break; } - case LoggingLevel::ERROR: { + case LoggingLevel::Error: { axom::slic::setLoggingMsgLevel( axom::slic::message::Error ); break; } diff --git a/src/tribol/mesh/MfemData.cpp b/src/tribol/mesh/MfemData.cpp index 36a2775e..f08f4af2 100644 --- a/src/tribol/mesh/MfemData.cpp +++ b/src/tribol/mesh/MfemData.cpp @@ -860,7 +860,7 @@ MfemJacobianData::MfemJacobianData( const MfemMeshData& parent_data, const MfemS submesh_parent_data = 1.0; // This constructor copies all of the data, so don't worry about ownership of the CSR data submesh_parent_vdof_xfer_ = std::make_unique( - TRIBOL_COMM_WORLD, submesh_fes.GetVSize(), submesh_fes.GlobalVSize(), parent_fes.GlobalVSize(), + parent_fes.GetComm(), submesh_fes.GetVSize(), submesh_fes.GlobalVSize(), parent_fes.GlobalVSize(), submesh_parent_I.data(), submesh2parent_vdof_list_.GetData(), submesh_parent_data.GetData(), submesh_fes.GetDofOffsets(), parent_fes.GetDofOffsets() ); @@ -1037,9 +1037,9 @@ std::unique_ptr MfemJacobianData::GetMfemBlockJacobian( if ( has_11 ) { auto& submesh_fes_full = submesh_data_.GetSubmeshFESpace(); - shared::ParSparseMat inactive_hpm_full = - shared::ParSparseMat::diagonalMatrix( TRIBOL_COMM_WORLD, submesh_fes_full.GlobalTrueVSize(), - submesh_fes_full.GetTrueDofOffsets(), 1.0, mortar_tdof_list_, false ); + auto comm = parent_data_.GetParentCoords().ParFESpace()->GetComm(); + shared::ParSparseMat inactive_hpm_full = shared::ParSparseMat::diagonalMatrix( + comm, submesh_fes_full.GlobalTrueVSize(), submesh_fes_full.GetTrueDofOffsets(), 1.0, mortar_tdof_list_, false ); if ( block_J->IsZeroBlock( 1, 1 ) ) { block_J->SetBlock( 1, 1, inactive_hpm_full.release() ); @@ -1082,6 +1082,160 @@ const MfemJacobianData::UpdateData& MfemJacobianData::GetUpdateData() const return *update_data_; } +shared::ParSparseMat MfemJacobianData::GetMfemJacobian( const std::vector& contributions ) const +{ + std::unique_ptr par_J; + + // Maps BlockSpaces (MORTAR, NONMORTAR, LAGRANGE_MULTIPLIER) to a tribol element map + const std::vector*> elem_map_by_space{ &parent_data_.GetElemMap1(), &parent_data_.GetElemMap2(), + &parent_data_.GetElemMap2() }; + + auto comm = parent_data_.GetParentCoords().ParFESpace()->GetComm(); + + // Iterate over all possible blocks: (0,0), (0,1), (1,0), (1,1) + // 0: Displacement (MORTAR/NONMORTAR) + // 1: Pressure/Gap (LAGRANGE_MULTIPLIER) + for ( int r_blk = 0; r_blk < 2; ++r_blk ) { + for ( int c_blk = 0; c_blk < 2; ++c_blk ) { + // Check if we have a transfer operator for this block + if ( GetUpdateData().submesh_redecomp_xfer_.shape()[0] <= r_blk || + GetUpdateData().submesh_redecomp_xfer_.shape()[1] <= c_blk || + !GetUpdateData().submesh_redecomp_xfer_( r_blk, c_blk ) ) { + continue; + } + + axom::Array row_redecomp_ids; + axom::Array col_redecomp_ids; + axom::Array jacobian_data; + axom::Array jacobian_offsets; + + // Aggregate data for this block pair + for ( const auto& contrib : contributions ) { + int contrib_r_blk = ( contrib.row_space == BlockSpace::LAGRANGE_MULTIPLIER ) ? 1 : 0; + int contrib_c_blk = ( contrib.col_space == BlockSpace::LAGRANGE_MULTIPLIER ) ? 1 : 0; + + if ( contrib_r_blk == r_blk && contrib_c_blk == c_blk ) { + int current_offset = jacobian_data.size(); + row_redecomp_ids.reserve( row_redecomp_ids.size() + contrib.row_elem_ids.size() ); + for ( auto id : contrib.row_elem_ids ) { + row_redecomp_ids.push_back( + ( *elem_map_by_space[static_cast( contrib.row_space )] )[static_cast( id )] ); + } + col_redecomp_ids.reserve( col_redecomp_ids.size() + contrib.col_elem_ids.size() ); + for ( auto id : contrib.col_elem_ids ) { + col_redecomp_ids.push_back( + ( *elem_map_by_space[static_cast( contrib.col_space )] )[static_cast( id )] ); + } + // NOTE (EBC): This can be removed when Axom PR 1819 goes in + if ( contrib.jacobian_data.size() > 0 ) { + jacobian_data.append( axom::ArrayView( contrib.jacobian_data ) ); + } + jacobian_offsets.reserve( jacobian_offsets.size() + contrib.jacobian_offsets.size() ); + for ( auto offset : contrib.jacobian_offsets ) { + jacobian_offsets.push_back( current_offset + offset ); + } + } + } + + // Check globally if any rank has data for this block + int local_has_data = row_redecomp_ids.empty() ? 0 : 1; + int global_has_data = 0; + MPI_Allreduce( &local_has_data, &global_has_data, 1, MPI_INT, MPI_MAX, comm ); + + if ( global_has_data ) { + redecomp::MatrixTransfer* xfer = GetUpdateData().submesh_redecomp_xfer_( r_blk, c_blk ).get(); + auto submesh_J = + xfer->TransferToParallel( row_redecomp_ids, col_redecomp_ids, jacobian_data, jacobian_offsets, false ); + + shared::ParSparseMatView submesh_J_view( &submesh_J.get() ); + std::unique_ptr contrib_J; + + if ( r_blk == 0 && c_blk == 0 ) { + auto parent_J = submesh_J_view.RAP( *submesh_parent_vdof_xfer_ ); + shared::ParSparseMatView parent_P( parent_data_.GetParentCoords().ParFESpace()->Dof_TrueDof_Matrix() ); + contrib_J = std::make_unique( parent_J.RAP( parent_P ) ); + } else if ( r_blk == 0 && c_blk == 1 ) { + auto parent_J = submesh_parent_vdof_xfer_->transpose() * submesh_J; + contrib_J = std::make_unique( + shared::ParSparseMat::RAP( parent_data_.GetParentCoords().ParFESpace()->Dof_TrueDof_Matrix(), parent_J, + submesh_data_.GetSubmeshFESpace().Dof_TrueDof_Matrix() ) ); + } else if ( r_blk == 1 && c_blk == 0 ) { + auto parent_J = submesh_J_view * ( *submesh_parent_vdof_xfer_ ); + shared::ParSparseMatView submesh_P( submesh_data_.GetSubmeshFESpace().Dof_TrueDof_Matrix() ); + shared::ParSparseMatView parent_P( parent_data_.GetParentCoords().ParFESpace()->Dof_TrueDof_Matrix() ); + contrib_J = + std::make_unique( shared::ParSparseMat::RAP( submesh_P, parent_J, parent_P ) ); + } else { + // (1, 1) block + shared::ParSparseMatView submesh_P( submesh_data_.GetSubmeshFESpace().Dof_TrueDof_Matrix() ); + contrib_J = std::make_unique( submesh_J_view.RAP( submesh_P ) ); + } + + if ( !par_J ) { + par_J = std::move( contrib_J ); + } else { + ( *par_J ) += *contrib_J; + } + } + } + } + + if ( !par_J ) { + int target_r_blk = 0; + int target_c_blk = 0; + if ( !contributions.empty() ) { + target_r_blk = ( contributions[0].row_space == BlockSpace::LAGRANGE_MULTIPLIER ) ? 1 : 0; + target_c_blk = ( contributions[0].col_space == BlockSpace::LAGRANGE_MULTIPLIER ) ? 1 : 0; + } + + auto& row_fes = + ( target_r_blk == 0 ) ? *parent_data_.GetParentCoords().ParFESpace() : submesh_data_.GetSubmeshFESpace(); + auto& col_fes = + ( target_c_blk == 0 ) ? *parent_data_.GetParentCoords().ParFESpace() : submesh_data_.GetSubmeshFESpace(); + + if ( target_r_blk == target_c_blk ) { + return shared::ParSparseMat::diagonalMatrix( comm, row_fes.GlobalTrueVSize(), row_fes.GetTrueDofOffsets(), 0.0, + mfem::Array(), true ); + } else { + mfem::SparseMatrix empty_diag( row_fes.GetTrueVSize(), col_fes.GetTrueVSize() ); + empty_diag.Finalize(); + return shared::ParSparseMat( comm, row_fes.GlobalTrueVSize(), col_fes.GlobalTrueVSize(), + row_fes.GetTrueDofOffsets(), col_fes.GetTrueDofOffsets(), std::move( empty_diag ) ); + } + } + + return std::move( *par_J ); +} + +JacobianContributions::JacobianContributions( std::initializer_list> blocks ) +{ + for ( const auto& block : blocks ) { + ComputedElementData data; + data.row_space = block.first; + data.col_space = block.second; + contributions_.push_back( std::move( data ) ); + } +} + +void JacobianContributions::reserve( int n_pairs, int n_entries_per_pair ) +{ + for ( auto& contrib : contributions_ ) { + contrib.row_elem_ids.reserve( n_pairs ); + contrib.col_elem_ids.reserve( n_pairs ); + contrib.jacobian_data.reserve( n_pairs * n_entries_per_pair ); + contrib.jacobian_offsets.reserve( n_pairs ); + } +} + +void JacobianContributions::push_back( int block_idx, int row_elem_id, int col_elem_id, const double* data, int size ) +{ + auto& contrib = contributions_[block_idx]; + contrib.row_elem_ids.push_back( row_elem_id ); + contrib.col_elem_ids.push_back( col_elem_id ); + contrib.jacobian_offsets.push_back( contrib.jacobian_data.size() ); + contrib.jacobian_data.append( axom::ArrayView( data, size ) ); +} + } // namespace tribol #endif /* BUILD_REDECOMP */ diff --git a/src/tribol/mesh/MfemData.hpp b/src/tribol/mesh/MfemData.hpp index 0f432880..53c13378 100644 --- a/src/tribol/mesh/MfemData.hpp +++ b/src/tribol/mesh/MfemData.hpp @@ -1631,6 +1631,60 @@ class MfemSubmeshData { std::unique_ptr update_data_; }; +/** + * @brief Struct to hold computed element data for Jacobian assembly + */ +struct ComputedElementData { + BlockSpace row_space; ///< Block space for row elements + BlockSpace col_space; ///< Block space for column elements + axom::Array row_elem_ids; ///< Tribol element IDs for rows + axom::Array col_elem_ids; ///< Tribol element IDs for columns + axom::Array jacobian_data; ///< Flattened Jacobian data + axom::Array jacobian_offsets; ///< Offsets into data for each element +}; + +/** + * @brief Helper class to manage Jacobian contributions for different block spaces + */ +class JacobianContributions { + public: + /** + * @brief Construct a new JacobianContributions object + * + * @param blocks List of {row_space, col_space} pairs defining the Jacobian blocks + */ + JacobianContributions( std::initializer_list> blocks ); + + /** + * @brief Reserve memory for each block contribution + * + * @param n_pairs Number of interface pairs + * @param n_entries_per_pair Number of entries in the element Jacobian block + */ + void reserve( int n_pairs, int n_entries_per_pair ); + + /** + * @brief Add an element Jacobian contribution to a specific block + * + * @param block_idx Index of the block space pair (in the order provided to the constructor) + * @param row_elem_id Tribol element ID for the row space + * @param col_elem_id Tribol element ID for the column space + * @param data Pointer to the flattened element Jacobian data (column-major) + * @param size Number of entries in the element Jacobian data + */ + void push_back( int block_idx, int row_elem_id, int col_elem_id, const double* data, int size ); + + /** + * @brief Return the underlying contributions vector + * + * @return const std::vector& + */ + const std::vector& get() const { return contributions_; } + + private: + std::vector contributions_; +}; + /** * @brief Simplifies transfer of Jacobian matrix data between MFEM and Tribol */ @@ -1665,6 +1719,14 @@ class MfemJacobianData { const MethodData& method_data, const std::vector>& row_info, const std::vector>& col_info ) const; + /** + * @brief Returns a Jacobian as a single ParSparseMat + * + * @param contributions List of element computed data chunks + * @return ParSparseMat + */ + shared::ParSparseMat GetMfemJacobian( const std::vector& contributions ) const; + private: /** * @brief Creates and stores data that changes when the redecomp mesh is diff --git a/src/tribol/utils/TestUtils.cpp b/src/tribol/utils/TestUtils.cpp index 110fa74d..6da885f4 100644 --- a/src/tribol/utils/TestUtils.cpp +++ b/src/tribol/utils/TestUtils.cpp @@ -1125,7 +1125,7 @@ int TestMesh::tribolSetupAndUpdate( ContactMethod method, EnforcementMethod enfo setPlotCycleIncrement( csIndex, 1 ); } - setLoggingLevel( csIndex, LoggingLevel::WARNING ); + setLoggingLevel( csIndex, LoggingLevel::Warning ); if ( method == COMMON_PLANE && enforcement == PENALTY ) { PenaltyConstraintType constraint_type =