diff --git a/scripts/spack/configs/versions.yaml b/scripts/spack/configs/versions.yaml index 7ebda5d8..4ef79642 100644 --- a/scripts/spack/configs/versions.yaml +++ b/scripts/spack/configs/versions.yaml @@ -20,7 +20,7 @@ packages: - spec: "@2.26.0" mfem: require: - - spec: "@4.8.0.1" + - spec: "@4.9.0" petsc: require: - spec: "@3.21.6" diff --git a/src/tests/CMakeLists.txt b/src/tests/CMakeLists.txt index 150f1dfa..1c1f8777 100644 --- a/src/tests/CMakeLists.txt +++ b/src/tests/CMakeLists.txt @@ -64,6 +64,52 @@ foreach( test ${tribol_tests} ) endforeach() +#------------------------------------------------------------------------------ +# Add MPI tests +#------------------------------------------------------------------------------ +if ( TRIBOL_USE_MPI ) + + set( mpi_tests + tribol_par_sparse_mat.cpp + tribol_par_vector.cpp + ) + + set(mpi_test_depends tribol gtest) + + foreach( test ${mpi_tests} ) + + get_filename_component( test_name ${test} NAME_WE ) + + blt_add_executable( + NAME ${test_name}_test + SOURCES ${test} + OUTPUT_DIR ${TEST_OUTPUT_DIRECTORY} + DEPENDS_ON ${mpi_test_depends} ${tribol_device_depends} + FOLDER tribol/tests ) + + blt_add_test( NAME ${test_name} + COMMAND ${test_name}_test + NUM_MPI_TASKS 2 ) + + if (ENABLE_CUDA) + set_target_properties(${test_name}_test PROPERTIES CUDA_SEPARABLE_COMPILATION On) + endif() + + if (ENABLE_HIP) + target_compile_options(${test_name}_test PRIVATE -fgpu-rdc) + target_compile_options(${test_name}_test PRIVATE + $<$: + -ggdb + > + ) + target_link_options(${test_name}_test PRIVATE -fgpu-rdc) + target_link_options(${test_name}_test PRIVATE --hip-link) + endif() + + endforeach() + +endif() + #------------------------------------------------------------------------------ # Add redecomp tests #------------------------------------------------------------------------------ diff --git a/src/tests/tribol_par_sparse_mat.cpp b/src/tests/tribol_par_sparse_mat.cpp new file mode 100644 index 00000000..cac3cf88 --- /dev/null +++ b/src/tests/tribol_par_sparse_mat.cpp @@ -0,0 +1,405 @@ +// 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_MPI +#include +#endif + +#include "mfem.hpp" + +#include "tribol/utils/ParVector.hpp" +#include "tribol/utils/ParSparseMat.hpp" + +class ParSparseMatTest : public ::testing::Test { + protected: + mfem::Array GetRowStarts( MPI_Comm comm, HYPRE_BigInt global_size ) + { + int rank, num_procs; + MPI_Comm_rank( comm, &rank ); + MPI_Comm_size( comm, &num_procs ); + + int local_size = global_size / num_procs; + int remainder = global_size % num_procs; + if ( rank < remainder ) { + local_size++; + } + + mfem::Array row_starts( num_procs + 1 ); + std::vector local_sizes( num_procs ); + MPI_Allgather( &local_size, 1, MPI_INT, local_sizes.data(), 1, MPI_INT, comm ); + row_starts[0] = 0; + for ( int i = 0; i < num_procs; ++i ) { + row_starts[i + 1] = row_starts[i] + local_sizes[i]; + } + if ( HYPRE_AssumedPartitionCheck() ) { + auto total_dofs = row_starts[num_procs]; + row_starts.SetSize( 3 ); + row_starts[0] = row_starts[rank]; + row_starts[1] = row_starts[rank + 1]; + row_starts[2] = total_dofs; + } + + return row_starts; + } +}; + +// Test Construction +TEST_F( ParSparseMatTest, Construction ) +{ + int rank; + MPI_Comm_rank( MPI_COMM_WORLD, &rank ); + int num_procs; + MPI_Comm_size( MPI_COMM_WORLD, &num_procs ); + constexpr int size = 10; + int local_size = size / num_procs + ( rank < ( size % num_procs ) ? 1 : 0 ); + if ( rank == 0 ) std::cout << "Testing Construction..." << std::endl; + + auto row_starts_array = GetRowStarts( MPI_COMM_WORLD, size ); + + // 1. From mfem::HypreParMatrix* + mfem::HypreParMatrix* m1 = + tribol::ParSparseMat::diagonalMatrix( MPI_COMM_WORLD, size, row_starts_array, 1.0 ).release(); + tribol::ParSparseMat psm1( m1 ); + EXPECT_EQ( psm1.Height(), local_size ); + + // 2. From unique_ptr + auto m2 = std::unique_ptr( + tribol::ParSparseMat::diagonalMatrix( MPI_COMM_WORLD, size, row_starts_array, 2.0 ).release() ); + tribol::ParSparseMat psm2( std::move( m2 ) ); + EXPECT_EQ( psm2.Height(), local_size ); + + // 3. From SparseMatrix rvalue + mfem::SparseMatrix diag( local_size ); + for ( int i = 0; i < local_size; ++i ) diag.Set( i, i, 3.0 ); + diag.Finalize(); + + tribol::ParSparseMat psm3( MPI_COMM_WORLD, (HYPRE_BigInt)size, row_starts_array.GetData(), std::move( diag ) ); + EXPECT_EQ( psm3.Height(), local_size ); + + mfem::Vector x( local_size ), y( local_size ); + x = 1.0; + psm3->Mult( x, y ); + EXPECT_NEAR( y.Max(), 3.0, 1e-12 ); +} + +// Test View +TEST_F( ParSparseMatTest, View ) +{ + int rank; + MPI_Comm_rank( MPI_COMM_WORLD, &rank ); + if ( rank == 0 ) std::cout << "Testing View..." << std::endl; + + auto row_starts = GetRowStarts( MPI_COMM_WORLD, 10 ); + tribol::ParSparseMat A = tribol::ParSparseMat::diagonalMatrix( MPI_COMM_WORLD, 10, row_starts, 2.0 ); + + // Construct View + tribol::ParSparseMatView view( &A.get() ); + + EXPECT_EQ( view.Height(), A.Height() ); + + // Operate on View + tribol::ParSparseMat B = view * 2.0; + mfem::Vector x( A.Width() ), y( A.Height() ); + x = 1.0; + B->Mult( x, y ); + EXPECT_NEAR( y.Max(), 4.0, 1e-12 ); +} + +// Test Addition +TEST_F( ParSparseMatTest, Addition ) +{ + int rank; + MPI_Comm_rank( MPI_COMM_WORLD, &rank ); + if ( rank == 0 ) std::cout << "Testing Addition..." << std::endl; + + auto row_starts = GetRowStarts( MPI_COMM_WORLD, 10 ); + tribol::ParSparseMat A = tribol::ParSparseMat::diagonalMatrix( MPI_COMM_WORLD, 10, row_starts, 2.0 ); + tribol::ParSparseMat B = tribol::ParSparseMat::diagonalMatrix( MPI_COMM_WORLD, 10, row_starts, 3.0 ); + + // A + B + tribol::ParSparseMat C = A + B; + mfem::Vector x( A.Width() ), y( A.Height() ); + x = 1.0; + C->Mult( x, y ); + // Result should be (2+3)*1 = 5 + EXPECT_NEAR( y.Max(), 5.0, 1e-12 ); + EXPECT_NEAR( y.Min(), 5.0, 1e-12 ); + + // A += B + A += B; + A->Mult( x, y ); + EXPECT_NEAR( y.Max(), 5.0, 1e-12 ); +} + +// Test Subtraction +TEST_F( ParSparseMatTest, Subtraction ) +{ + int rank; + MPI_Comm_rank( MPI_COMM_WORLD, &rank ); + if ( rank == 0 ) std::cout << "Testing Subtraction..." << std::endl; + + auto row_starts = GetRowStarts( MPI_COMM_WORLD, 10 ); + tribol::ParSparseMat A = tribol::ParSparseMat::diagonalMatrix( MPI_COMM_WORLD, 10, row_starts, 5.0 ); + tribol::ParSparseMat B = tribol::ParSparseMat::diagonalMatrix( MPI_COMM_WORLD, 10, row_starts, 2.0 ); + + // A - B + tribol::ParSparseMat C = A - B; + mfem::Vector x( A.Width() ), y( A.Height() ); + x = 1.0; + C->Mult( x, y ); + // Result should be (5-2)*1 = 3 + EXPECT_NEAR( y.Max(), 3.0, 1e-12 ); + + // A -= B + A -= B; + A->Mult( x, y ); + EXPECT_NEAR( y.Max(), 3.0, 1e-12 ); +} + +// Test Scalar Multiplication +TEST_F( ParSparseMatTest, ScalarMult ) +{ + int rank; + MPI_Comm_rank( MPI_COMM_WORLD, &rank ); + if ( rank == 0 ) std::cout << "Testing Scalar Multiplication..." << std::endl; + + auto row_starts = GetRowStarts( MPI_COMM_WORLD, 10 ); + tribol::ParSparseMat A = tribol::ParSparseMat::diagonalMatrix( MPI_COMM_WORLD, 10, row_starts, 2.0 ); + + // A * s + tribol::ParSparseMat B = A * 3.0; + mfem::Vector x( A.Width() ), y( A.Height() ); + x = 1.0; + B->Mult( x, y ); + EXPECT_NEAR( y.Max(), 6.0, 1e-12 ); + + // s * A + tribol::ParSparseMat C = 4.0 * A; + C->Mult( x, y ); + EXPECT_NEAR( y.Max(), 8.0, 1e-12 ); +} + +// Test Matrix Multiplication +TEST_F( ParSparseMatTest, MatrixMult ) +{ + int rank; + MPI_Comm_rank( MPI_COMM_WORLD, &rank ); + if ( rank == 0 ) std::cout << "Testing Matrix Multiplication..." << std::endl; + + auto row_starts = GetRowStarts( MPI_COMM_WORLD, 10 ); + tribol::ParSparseMat A = tribol::ParSparseMat::diagonalMatrix( MPI_COMM_WORLD, 10, row_starts, 2.0 ); + tribol::ParSparseMat B = tribol::ParSparseMat::diagonalMatrix( MPI_COMM_WORLD, 10, row_starts, 3.0 ); + + // A * B + tribol::ParSparseMat C = A * B; + mfem::Vector x( A.Width() ), y( A.Height() ); + x = 1.0; + C->Mult( x, y ); + // Result should be (2*3)*1 = 6 + EXPECT_NEAR( y.Max(), 6.0, 1e-12 ); + + // A *= B + A *= B; + A->Mult( x, y ); + EXPECT_NEAR( y.Max(), 6.0, 1e-12 ); +} + +// Test Matrix-Vector Multiplication +TEST_F( ParSparseMatTest, MatVecMult ) +{ + int rank; + MPI_Comm_rank( MPI_COMM_WORLD, &rank ); + if ( rank == 0 ) std::cout << "Testing Matrix-Vector Multiplication..." << std::endl; + + auto row_starts = GetRowStarts( MPI_COMM_WORLD, 10 ); + tribol::ParSparseMat A = tribol::ParSparseMat::diagonalMatrix( MPI_COMM_WORLD, 10, row_starts, 2.0 ); + mfem::HypreParVector x_hypre( A.get(), 1 ); + x_hypre = 1.0; + tribol::ParVectorView x( &x_hypre ); + + // y = A * x + tribol::ParVector y = A * x; + EXPECT_NEAR( y.Max(), 2.0, 1e-12 ); +} + +// Test Vector-Matrix Multiplication +TEST_F( ParSparseMatTest, VecMatMult ) +{ + int rank; + MPI_Comm_rank( MPI_COMM_WORLD, &rank ); + if ( rank == 0 ) std::cout << "Testing Vector-Matrix Multiplication..." << std::endl; + + auto row_starts = GetRowStarts( MPI_COMM_WORLD, 10 ); + tribol::ParSparseMat A = tribol::ParSparseMat::diagonalMatrix( MPI_COMM_WORLD, 10, row_starts, 3.0 ); + mfem::HypreParVector x_hypre( A.get(), 0 ); + x_hypre = 1.0; + tribol::ParVectorView x( &x_hypre ); + + // y = x^T * A + tribol::ParVector y = x * A; + EXPECT_NEAR( y.Max(), 3.0, 1e-12 ); + EXPECT_NEAR( y.Min(), 3.0, 1e-12 ); +} + +// Test Elimination +TEST_F( ParSparseMatTest, Elimination ) +{ + int rank; + MPI_Comm_rank( MPI_COMM_WORLD, &rank ); + if ( rank == 0 ) std::cout << "Testing Elimination..." << std::endl; + + auto row_starts = GetRowStarts( MPI_COMM_WORLD, 10 ); + tribol::ParSparseMat A = tribol::ParSparseMat::diagonalMatrix( MPI_COMM_WORLD, 10, row_starts, 3.0 ); + + // Eliminate row 0 (globally) + // Determine if I own row 0 + mfem::Array rows_to_elim; + int row_starts_idx = HYPRE_AssumedPartitionCheck() ? 0 : rank; + if ( row_starts[row_starts_idx] == 0 ) { + rows_to_elim.Append( 0 ); + } + A.EliminateRows( rows_to_elim ); + + // Check if row 0 is identity (or zero with diagonal 1) + // Diagonal matrix means we can just check multiplication + mfem::HypreParVector x_hypre( A.get(), 1 ); + x_hypre = 1.0; + tribol::ParVectorView x( &x_hypre ); + tribol::ParVector y = A * x; // y = A * x + + // if rank owns row 0, the result for that row should be 0.0 * x[0] = 0.0 (since diag is 0.0) + // other rows should be 3.0 + + // local row 0 on rank 0 is global row 0 + if ( rank == 0 ) { + EXPECT_NEAR( y[0], 0.0, 1e-12 ); + for ( int i = 1; i < y.Size(); ++i ) { + EXPECT_NEAR( y[i], 3.0, 1e-12 ); + } + } else { + for ( int i = 0; i < y.Size(); ++i ) { + EXPECT_NEAR( y[i], 3.0, 1e-12 ); + } + } + + A = tribol::ParSparseMat::diagonalMatrix( MPI_COMM_WORLD, 10, row_starts, 3.0 ); + int num_procs; + MPI_Comm_size( MPI_COMM_WORLD, &num_procs ); + + // Eliminate last local col + auto last_local_col = A.Width() - 1; + mfem::Array cols_to_elim( { last_local_col } ); + + tribol::ParSparseMat Ae = A.EliminateCols( cols_to_elim ); + + // Now check A * e_last = 0 + // Create vector with 1 at last_local_col, 0 elsewhere + mfem::HypreParVector x_hypre_last( A.get(), 1 ); + x_hypre_last = 0.0; + x_hypre_last[last_local_col] = 1.0; + tribol::ParVectorView x_last( &x_hypre_last ); + tribol::ParVector y_last = A * x_last; + EXPECT_NEAR( y_last[last_local_col], 0.0, 1e-12 ); + + // Check Ae * e_last = original value + tribol::ParVector ye = Ae * x_last; + double expected_val = 3.0; + + EXPECT_NEAR( ye[last_local_col], expected_val, 1e-12 ); +} + +// Test Transpose and Square +TEST_F( ParSparseMatTest, TransposeSquare ) +{ + int rank; + MPI_Comm_rank( MPI_COMM_WORLD, &rank ); + if ( rank == 0 ) std::cout << "Testing Transpose and Square..." << std::endl; + + auto row_starts = GetRowStarts( MPI_COMM_WORLD, 10 ); + tribol::ParSparseMat A = tribol::ParSparseMat::diagonalMatrix( MPI_COMM_WORLD, 10, row_starts, 2.0 ); + + // Transpose (Diagonal matrix is symmetric) + tribol::ParSparseMat At = A.transpose(); + mfem::Vector x( A.Width() ), y( A.Height() ); + x = 1.0; + At->Mult( x, y ); + EXPECT_NEAR( y.Max(), 2.0, 1e-12 ); + + // Square + tribol::ParSparseMat A2 = A.square(); + A2->Mult( x, y ); + EXPECT_NEAR( y.Max(), 4.0, 1e-12 ); +} + +// Test RAP +TEST_F( ParSparseMatTest, RAP ) +{ + int rank; + MPI_Comm_rank( MPI_COMM_WORLD, &rank ); + if ( rank == 0 ) std::cout << "Testing RAP..." << std::endl; + + // Use Identity for P to simplify testing: P^T * A * P = I * A * I = A + auto row_starts = GetRowStarts( MPI_COMM_WORLD, 10 ); + tribol::ParSparseMat A = tribol::ParSparseMat::diagonalMatrix( MPI_COMM_WORLD, 10, row_starts, 5.0 ); + tribol::ParSparseMat P = tribol::ParSparseMat::diagonalMatrix( MPI_COMM_WORLD, 10, row_starts, 1.0 ); + tribol::ParSparseMat R = tribol::ParSparseMat::diagonalMatrix( MPI_COMM_WORLD, 10, row_starts, 1.0 ); + + // RAP(P) + tribol::ParSparseMat Res1 = A.RAP( P ); + mfem::Vector x( A.Width() ), y( A.Height() ); + x = 1.0; + Res1->Mult( x, y ); + EXPECT_NEAR( y.Max(), 5.0, 1e-12 ); + + // RAP(R, A, P) + tribol::ParSparseMat Res2 = tribol::ParSparseMat::RAP( R, A, P ); + Res2->Mult( x, y ); + EXPECT_NEAR( y.Max(), 5.0, 1e-12 ); +} + +// Test Accessors +TEST_F( ParSparseMatTest, Accessors ) +{ + int rank; + MPI_Comm_rank( MPI_COMM_WORLD, &rank ); + int num_procs; + MPI_Comm_size( MPI_COMM_WORLD, &num_procs ); + constexpr int size = 10; + int local_size = size / num_procs + ( rank < ( size % num_procs ) ? 1 : 0 ); + if ( rank == 0 ) std::cout << "Testing Accessors..." << std::endl; + + auto row_starts = GetRowStarts( MPI_COMM_WORLD, size ); + tribol::ParSparseMat A = tribol::ParSparseMat::diagonalMatrix( MPI_COMM_WORLD, size, row_starts, 1.0 ); + + // get() + EXPECT_EQ( A.Height(), local_size ); + + // operator-> + EXPECT_EQ( A->Height(), local_size ); +} + +//------------------------------------------------------------------------------ +#include "axom/slic/core/SimpleLogger.hpp" + +int main( int argc, char* argv[] ) +{ + int result = 0; + + MPI_Init( &argc, &argv ); + + ::testing::InitGoogleTest( &argc, argv ); + + axom::slic::SimpleLogger logger; + + result = RUN_ALL_TESTS(); + + MPI_Finalize(); + + return result; +} diff --git a/src/tests/tribol_par_vector.cpp b/src/tests/tribol_par_vector.cpp new file mode 100644 index 00000000..9ca74ca1 --- /dev/null +++ b/src/tests/tribol_par_vector.cpp @@ -0,0 +1,305 @@ +// 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_MPI +#include +#endif + +#include "mfem.hpp" + +#include "tribol/utils/ParVector.hpp" +#include "tribol/utils/ParSparseMat.hpp" + +class ParVectorTest : public ::testing::Test { + protected: + mfem::Array GetRowStarts( MPI_Comm comm, HYPRE_BigInt global_size ) + { + int rank, num_procs; + MPI_Comm_rank( comm, &rank ); + MPI_Comm_size( comm, &num_procs ); + + int local_size = global_size / num_procs; + int remainder = global_size % num_procs; + if ( rank < remainder ) { + local_size++; + } + + mfem::Array row_starts( num_procs + 1 ); + std::vector local_sizes( num_procs ); + MPI_Allgather( &local_size, 1, MPI_INT, local_sizes.data(), 1, MPI_INT, comm ); + row_starts[0] = 0; + for ( int i = 0; i < num_procs; ++i ) { + row_starts[i + 1] = row_starts[i] + local_sizes[i]; + } + if ( HYPRE_AssumedPartitionCheck() ) { + auto total_dofs = row_starts[num_procs]; + row_starts.SetSize( 3 ); + row_starts[0] = row_starts[rank]; + row_starts[1] = row_starts[rank + 1]; + row_starts[2] = total_dofs; + } + + return row_starts; + } +}; + +// Test Construction +TEST_F( ParVectorTest, Construction ) +{ + int rank; + MPI_Comm_rank( MPI_COMM_WORLD, &rank ); + int num_procs; + MPI_Comm_size( MPI_COMM_WORLD, &num_procs ); + constexpr int size = 10; + int local_size = size / num_procs + ( rank < ( size % num_procs ) ? 1 : 0 ); + if ( rank == 0 ) std::cout << "Testing Construction..." << std::endl; + + auto row_starts_array = GetRowStarts( MPI_COMM_WORLD, size ); + + // 1. From mfem::HypreParVector* + mfem::HypreParVector* v1 = new mfem::HypreParVector( MPI_COMM_WORLD, size, row_starts_array.GetData() ); + tribol::ParVector pv1( v1 ); + EXPECT_EQ( pv1.Size(), local_size ); + + // 2. From unique_ptr + auto v2 = std::make_unique( MPI_COMM_WORLD, size, row_starts_array.GetData() ); + tribol::ParVector pv2( std::move( v2 ) ); + EXPECT_EQ( pv2.Size(), local_size ); + + // 3. Template constructor + tribol::ParVector pv3( MPI_COMM_WORLD, size, row_starts_array.GetData() ); + EXPECT_EQ( pv3.Size(), local_size ); +} + +// Test View +TEST_F( ParVectorTest, View ) +{ + int rank; + MPI_Comm_rank( MPI_COMM_WORLD, &rank ); + if ( rank == 0 ) std::cout << "Testing View..." << std::endl; + + auto row_starts = GetRowStarts( MPI_COMM_WORLD, 10 ); + tribol::ParVector v( MPI_COMM_WORLD, 10, row_starts.GetData() ); + v.Fill( 1.0 ); + + // Construct View + tribol::ParVectorView view( &v.get() ); + + EXPECT_EQ( view.Size(), v.Size() ); + EXPECT_NEAR( view.Max(), 1.0, 1e-12 ); + + // Operate on View + tribol::ParVector v2 = view * 2.0; + EXPECT_NEAR( v2.Max(), 2.0, 1e-12 ); +} + +// Test Accessors +TEST_F( ParVectorTest, Accessors ) +{ + int rank; + MPI_Comm_rank( MPI_COMM_WORLD, &rank ); + if ( rank == 0 ) std::cout << "Testing Accessors..." << std::endl; + + auto row_starts = GetRowStarts( MPI_COMM_WORLD, 10 ); + tribol::ParVector v( MPI_COMM_WORLD, 10, row_starts.GetData() ); + v.Fill( 0.0 ); + if ( v.Size() > 0 ) { + v[0] = 5.0; + EXPECT_NEAR( v[0], 5.0, 1e-12 ); + EXPECT_NEAR( ( *v.operator->() )[0], 5.0, 1e-12 ); + } + + v.Fill( 3.0 ); + EXPECT_NEAR( v.Max(), 3.0, 1e-12 ); + EXPECT_NEAR( v.Min(), 3.0, 1e-12 ); +} + +// Test Addition and Subtraction +TEST_F( ParVectorTest, AddSub ) +{ + int rank; + MPI_Comm_rank( MPI_COMM_WORLD, &rank ); + if ( rank == 0 ) std::cout << "Testing Addition and Subtraction..." << std::endl; + + auto row_starts = GetRowStarts( MPI_COMM_WORLD, 10 ); + tribol::ParVector v1( MPI_COMM_WORLD, 10, row_starts.GetData() ); + tribol::ParVector v2( MPI_COMM_WORLD, 10, row_starts.GetData() ); + v1.Fill( 2.0 ); + v2.Fill( 3.0 ); + + // v1 + v2 + tribol::ParVector v3 = v1 + v2; + EXPECT_NEAR( v3.Max(), 5.0, 1e-12 ); + + // v1 += v2 + v1 += v2; + EXPECT_NEAR( v1.Max(), 5.0, 1e-12 ); + + // v1 - v2 + tribol::ParVector v4 = v1 - v2; + EXPECT_NEAR( v4.Max(), 2.0, 1e-12 ); + + // v1 -= v2 + v1 -= v2; + EXPECT_NEAR( v1.Max(), 2.0, 1e-12 ); +} + +// Test Scalar Multiplication +TEST_F( ParVectorTest, ScalarMult ) +{ + int rank; + MPI_Comm_rank( MPI_COMM_WORLD, &rank ); + if ( rank == 0 ) std::cout << "Testing Scalar Multiplication..." << std::endl; + + auto row_starts = GetRowStarts( MPI_COMM_WORLD, 10 ); + tribol::ParVector v( MPI_COMM_WORLD, 10, row_starts.GetData() ); + v.Fill( 2.0 ); + + // v * s + tribol::ParVector v2 = v * 3.0; + EXPECT_NEAR( v2.Max(), 6.0, 1e-12 ); + + // s * v + tribol::ParVector v3 = 4.0 * v; + EXPECT_NEAR( v3.Max(), 8.0, 1e-12 ); + + // v *= s + v *= 5.0; + EXPECT_NEAR( v.Max(), 10.0, 1e-12 ); +} + +// Test Component-wise Multiplication and Division +TEST_F( ParVectorTest, ComponentWise ) +{ + int rank; + MPI_Comm_rank( MPI_COMM_WORLD, &rank ); + if ( rank == 0 ) std::cout << "Testing Component-wise operations..." << std::endl; + + auto row_starts = GetRowStarts( MPI_COMM_WORLD, 10 ); + tribol::ParVector v1( MPI_COMM_WORLD, 10, row_starts.GetData() ); + tribol::ParVector v2( MPI_COMM_WORLD, 10, row_starts.GetData() ); + v1.Fill( 2.0 ); + v2.Fill( 4.0 ); + + // multiply + tribol::ParVector v3 = v1.multiply( v2 ); + EXPECT_NEAR( v3.Max(), 8.0, 1e-12 ); + + // multiply in-place + v1.multiplyInPlace( v2 ); + EXPECT_NEAR( v1.Max(), 8.0, 1e-12 ); + + // divide + tribol::ParVector v4 = v1.divide( v2 ); + EXPECT_NEAR( v4.Max(), 2.0, 1e-12 ); + + // divide in-place + v1.divideInPlace( v2 ); + EXPECT_NEAR( v1.Max(), 2.0, 1e-12 ); +} + +// Test Move and Release +TEST_F( ParVectorTest, MoveAndRelease ) +{ + int rank; + MPI_Comm_rank( MPI_COMM_WORLD, &rank ); + if ( rank == 0 ) std::cout << "Testing Move and Release..." << std::endl; + + auto row_starts = GetRowStarts( MPI_COMM_WORLD, 10 ); + tribol::ParVector v1( MPI_COMM_WORLD, 10, row_starts.GetData() ); + v1.Fill( 7.0 ); + + // Move constructor + tribol::ParVector v2( std::move( v1 ) ); + EXPECT_NEAR( v2.Max(), 7.0, 1e-12 ); + EXPECT_EQ( v1.operator->(), nullptr ); + + // Move assignment + tribol::ParVector v3( MPI_COMM_WORLD, 10, row_starts.GetData() ); + v3 = std::move( v2 ); + EXPECT_NEAR( v3.Max(), 7.0, 1e-12 ); + EXPECT_EQ( v2.operator->(), nullptr ); + + // Release + mfem::HypreParVector* raw = v3.release(); + EXPECT_NE( raw, nullptr ); + EXPECT_NEAR( raw->Max(), 7.0, 1e-12 ); + delete raw; +} + +// Test Fill +TEST_F( ParVectorTest, Fill ) +{ + int rank; + MPI_Comm_rank( MPI_COMM_WORLD, &rank ); + if ( rank == 0 ) std::cout << "Testing Fill..." << std::endl; + + auto row_starts = GetRowStarts( MPI_COMM_WORLD, 10 ); + tribol::ParVector v( MPI_COMM_WORLD, 10, row_starts.GetData() ); + + v.Fill( 1.0 ); + EXPECT_NEAR( v.Max(), 1.0, 1e-12 ); + EXPECT_NEAR( v.Min(), 1.0, 1e-12 ); + + v.Fill( 2.5 ); + EXPECT_NEAR( v.Max(), 2.5, 1e-12 ); + EXPECT_NEAR( v.Min(), 2.5, 1e-12 ); +} + +// Test Copy +TEST_F( ParVectorTest, Copy ) +{ + int rank; + MPI_Comm_rank( MPI_COMM_WORLD, &rank ); + if ( rank == 0 ) std::cout << "Testing Copy..." << std::endl; + + auto row_starts = GetRowStarts( MPI_COMM_WORLD, 10 ); + tribol::ParVector v1( MPI_COMM_WORLD, 10, row_starts.GetData() ); + v1.Fill( 3.0 ); + + // Copy constructor + tribol::ParVector v2( v1 ); + EXPECT_NEAR( v2.Max(), 3.0, 1e-12 ); + + // Verify it's a deep copy + v1.Fill( 4.0 ); + EXPECT_NEAR( v1.Max(), 4.0, 1e-12 ); + EXPECT_NEAR( v2.Max(), 3.0, 1e-12 ); + + // Copy assignment + tribol::ParVector v3( MPI_COMM_WORLD, 10, row_starts.GetData() ); + v3.Fill( 5.0 ); + v3 = v2; + EXPECT_NEAR( v3.Max(), 3.0, 1e-12 ); + + // Verify deep copy for assignment + v2.Fill( 6.0 ); + EXPECT_NEAR( v2.Max(), 6.0, 1e-12 ); + EXPECT_NEAR( v3.Max(), 3.0, 1e-12 ); +} + +//------------------------------------------------------------------------------ +#include "axom/slic/core/SimpleLogger.hpp" + +int main( int argc, char* argv[] ) +{ + int result = 0; + + MPI_Init( &argc, &argv ); + + ::testing::InitGoogleTest( &argc, argv ); + + axom::slic::SimpleLogger logger; + + result = RUN_ALL_TESTS(); + + MPI_Finalize(); + + return result; +} diff --git a/src/tribol/CMakeLists.txt b/src/tribol/CMakeLists.txt index 5ec833f6..6b753287 100644 --- a/src/tribol/CMakeLists.txt +++ b/src/tribol/CMakeLists.txt @@ -49,6 +49,8 @@ set(tribol_headers utils/ContactPlaneOutput.hpp utils/DataManager.hpp utils/Math.hpp + utils/ParVector.hpp + utils/ParSparseMat.hpp utils/TestUtils.hpp ) @@ -80,6 +82,8 @@ set(tribol_sources utils/ContactPlaneOutput.cpp utils/Math.cpp + utils/ParVector.cpp + utils/ParSparseMat.cpp utils/TestUtils.cpp ) diff --git a/src/tribol/interface/mfem_tribol.cpp b/src/tribol/interface/mfem_tribol.cpp index 68874175..9f29bcd9 100644 --- a/src/tribol/interface/mfem_tribol.cpp +++ b/src/tribol/interface/mfem_tribol.cpp @@ -335,23 +335,28 @@ std::unique_ptr getMfemBlockJacobian( IndexT cs_id ) cs_id ) ); // creates a block Jacobian on the parent mesh/parent-linked boundary submesh based on the element Jacobians stored in // the coupling scheme's method data + const std::vector> all_info{ + { 0, BlockSpace::MORTAR }, { 0, BlockSpace::NONMORTAR }, { 1, BlockSpace::LAGRANGE_MULTIPLIER } }; if ( cs->isEnzymeEnabled() ) { - auto dfdx = cs->getMfemJacobianData()->GetMfemDfDxFullJacobian( *cs->getMethodData() ); - auto dfdn = cs->getMfemJacobianData()->GetMfemDfDnJacobian( *cs->getDfDnMethodData() ); - auto dndx = cs->getMfemJacobianData()->GetMfemDnDxJacobian( *cs->getDnDxMethodData() ); - auto dfdx_nconst = std::unique_ptr( - mfem::ParMult( &static_cast( dfdn->GetBlock( 0, 0 ) ), - &static_cast( dndx->GetBlock( 0, 0 ) ) ) ); - dfdx->SetBlock( 0, 0, - mfem::ParAdd( dfdx_nconst.get(), &static_cast( dfdx->GetBlock( 0, 0 ) ) ) ); - dfdx_nconst = std::unique_ptr( - mfem::ParMult( &static_cast( dfdn->GetBlock( 1, 0 ) ), - &static_cast( dndx->GetBlock( 0, 0 ) ) ) ); - dfdx->SetBlock( 1, 0, - mfem::ParAdd( dfdx_nconst.get(), &static_cast( dfdx->GetBlock( 1, 0 ) ) ) ); + auto dfdx = cs->getMfemJacobianData()->GetMfemBlockJacobian( *cs->getMethodData(), all_info, all_info ); + const std::vector> nonmortar_info{ { 0, BlockSpace::NONMORTAR } }; + auto dfdn = cs->getMfemJacobianData()->GetMfemBlockJacobian( *cs->getDfDnMethodData(), all_info, nonmortar_info ); + auto dndx = + cs->getMfemJacobianData()->GetMfemBlockJacobian( *cs->getDnDxMethodData(), nonmortar_info, nonmortar_info ); + + auto block_00 = ( ParSparseMatView( &static_cast( dfdn->GetBlock( 0, 0 ) ) ) * + &static_cast( dndx->GetBlock( 0, 0 ) ) ) + + &static_cast( dfdx->GetBlock( 0, 0 ) ); + dfdx->SetBlock( 0, 0, block_00.release() ); + + auto block_10 = ( ParSparseMatView( &static_cast( dfdn->GetBlock( 1, 0 ) ) ) * + &static_cast( dndx->GetBlock( 0, 0 ) ) ) + + &static_cast( dfdx->GetBlock( 1, 0 ) ); + dfdx->SetBlock( 1, 0, block_10.release() ); + return dfdx; } else { - return cs->getMfemJacobianData()->GetMfemBlockJacobian( cs->getMethodData() ); + return cs->getMfemJacobianData()->GetMfemBlockJacobian( *cs->getMethodData(), all_info, all_info ); } } diff --git a/src/tribol/mesh/MfemData.cpp b/src/tribol/mesh/MfemData.cpp index 3e552c3d..9aa7dce1 100644 --- a/src/tribol/mesh/MfemData.cpp +++ b/src/tribol/mesh/MfemData.cpp @@ -9,6 +9,8 @@ #ifdef BUILD_REDECOMP +#include + #include "axom/slic.hpp" #include "shared/infrastructure/Profiling.hpp" @@ -857,7 +859,7 @@ MfemJacobianData::MfemJacobianData( const MfemMeshData& parent_data, const MfemS mfem::Vector submesh_parent_data( submesh2parent_vdof_list_.Size() ); 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( + submesh_parent_vdof_xfer_ = std::make_unique( TRIBOL_COMM_WORLD, 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() ); @@ -913,386 +915,144 @@ void MfemJacobianData::UpdateJacobianXfer() update_data_ = std::make_unique( parent_data_, submesh_data_ ); } -std::unique_ptr MfemJacobianData::GetMfemBlockJacobian( const MethodData* method_data ) const +std::unique_ptr MfemJacobianData::GetMfemBlockJacobian( + const MethodData& method_data, const std::vector>& row_info, + const std::vector>& col_info ) const { - // 0 = displacement DOFs, 1 = lagrange multiplier DOFs - // (0,0) block is empty (for now using SINGLE_MORTAR with approximate tangent) - // (1,1) block is a diagonal matrix with ones on the diagonal of submesh nodes without a Lagrange multiplier DOF - // (0,1) and (1,0) are symmetric (for now using SINGLE_MORTAR with approximate tangent) - const auto& elem_map_1 = parent_data_.GetElemMap1(); - const auto& elem_map_2 = parent_data_.GetElemMap2(); - // empty data structures are needed even when no meshes are on rank since TransferToParallelSparse() needs to be - // called on all ranks (even those without data) - auto mortar_elems = ArrayT( 0, 0 ); - auto nonmortar_elems = ArrayT( 0, 0 ); - auto lm_elems = ArrayT( 0, 0 ); - auto elem_J_1_ptr = std::make_unique>( 0, 0 ); - auto elem_J_2_ptr = std::make_unique>( 0, 0 ); - const ArrayT* elem_J_1 = elem_J_1_ptr.get(); - const ArrayT* elem_J_2 = elem_J_2_ptr.get(); - // this means both of the meshes exist - if ( method_data != nullptr && !elem_map_1.empty() && !elem_map_2.empty() ) { - mortar_elems = method_data->getBlockJElementIds()[static_cast( BlockSpace::MORTAR )]; - for ( auto& mortar_elem : mortar_elems ) { - mortar_elem = elem_map_1[static_cast( mortar_elem )]; - } - nonmortar_elems = method_data->getBlockJElementIds()[static_cast( BlockSpace::NONMORTAR )]; - for ( auto& nonmortar_elem : nonmortar_elems ) { - nonmortar_elem = elem_map_2[static_cast( nonmortar_elem )]; - } - lm_elems = method_data->getBlockJElementIds()[static_cast( BlockSpace::LAGRANGE_MULTIPLIER )]; - for ( auto& lm_elem : lm_elems ) { - lm_elem = elem_map_2[static_cast( lm_elem )]; - } - // get (1,0) block - elem_J_1 = &method_data->getBlockJ()( static_cast( BlockSpace::LAGRANGE_MULTIPLIER ), - static_cast( BlockSpace::MORTAR ) ); - elem_J_2 = &method_data->getBlockJ()( static_cast( BlockSpace::LAGRANGE_MULTIPLIER ), - static_cast( BlockSpace::NONMORTAR ) ); - } - // move to submesh level - auto submesh_J = - GetUpdateData().submesh_redecomp_xfer_10_->TransferToParallelSparse( lm_elems, mortar_elems, *elem_J_1 ); - submesh_J += - GetUpdateData().submesh_redecomp_xfer_10_->TransferToParallelSparse( lm_elems, nonmortar_elems, *elem_J_2 ); - submesh_J.Finalize(); - - // transform J values from submesh to (global) parent mesh - mfem::Array J( submesh_J.NumNonZeroElems() ); - // This copy is needed to convert mfem::SparseMatrix int J values to the HYPRE_BigInt values the mfem::HypreParMatrix - // constructor needs - auto* J_int = submesh_J.GetJ(); - for ( int i{ 0 }; i < J.Size(); ++i ) { - J[i] = J_int[i]; - } - auto submesh_vector_fes = parent_data_.GetSubmeshFESpace(); - auto mpi = redecomp::MPIUtility( submesh_vector_fes.GetComm() ); - auto submesh_dof_offsets = ArrayT( mpi.NRanks() + 1, mpi.NRanks() + 1 ); - // we need the dof offsets of each rank. check if mfem stores this or if we - // need to create it. - if ( HYPRE_AssumedPartitionCheck() ) { - submesh_dof_offsets[mpi.MyRank() + 1] = submesh_vector_fes.GetDofOffsets()[1]; - mpi.Allreduce( &submesh_dof_offsets, MPI_SUM ); - } else { - for ( int i{ 0 }; i < mpi.NRanks(); ++i ) { - submesh_dof_offsets[i] = submesh_vector_fes.GetDofOffsets()[i]; - } - } - // the submesh to parent vdof map only exists for vdofs on rank, so J values - // not on rank will need to be transferred to the rank that the vdof exists on - // to query the map. the steps are laid out below. - - // step 1) query J values on rank for their parent vdof and package J values - // not on rank to send - auto send_J_by_rank = redecomp::MPIArray( &mpi ); - auto J_idx = redecomp::MPIArray( &mpi ); - auto est_num_J = submesh_J.NumNonZeroElems() / mpi.NRanks(); - for ( int r{}; r < mpi.NRanks(); ++r ) { - if ( r == mpi.MyRank() ) { - send_J_by_rank[r].shrink(); - J_idx[r].shrink(); - } else { - send_J_by_rank[r].reserve( est_num_J ); - J_idx[r].reserve( est_num_J ); - } - } - for ( int j{}; j < submesh_J.NumNonZeroElems(); ++j ) { - if ( J[j] >= submesh_dof_offsets[mpi.MyRank()] && J[j] < submesh_dof_offsets[mpi.MyRank() + 1] ) { - J[j] = submesh2parent_vdof_list_[J[j] - submesh_dof_offsets[mpi.MyRank()]]; - } else { - for ( int r{}; r < mpi.NRanks(); ++r ) { - if ( J[j] >= submesh_dof_offsets[r] && J[j] < submesh_dof_offsets[r + 1] ) { - send_J_by_rank[r].push_back( J[j] - submesh_dof_offsets[r] ); - J_idx[r].push_back( j ); - break; - } - } - } + // Determine block structure + int max_row_block = 0; + for ( auto info : row_info ) { + if ( info.first > max_row_block ) max_row_block = info.first; } - // step 2) sends the J values to the ranks that own them - auto recv_J_by_rank = redecomp::MPIArray( &mpi ); - recv_J_by_rank.SendRecvArrayEach( send_J_by_rank ); - // step 3) query the on-rank map to recover J values - for ( int r{}; r < mpi.NRanks(); ++r ) { - for ( auto& recv_J : recv_J_by_rank[r] ) { - recv_J = submesh2parent_vdof_list_[recv_J]; - } - } - // step 4) send the updated parent J values back and update the J vector - send_J_by_rank.SendRecvArrayEach( recv_J_by_rank ); - for ( int r{}; r < mpi.NRanks(); ++r ) { - for ( int j{}; j < send_J_by_rank[r].size(); ++j ) { - J[J_idx[r][j]] = send_J_by_rank[r][j]; - } + int max_col_block = 0; + for ( auto info : col_info ) { + if ( info.first > max_col_block ) max_col_block = info.first; } - // create block operator - auto block_J = std::make_unique( block_offsets_ ); - block_J->owns_blocks = 1; + SLIC_ERROR_ROOT_IF( max_row_block > GetUpdateData().submesh_redecomp_xfer_.shape()[0] || + max_col_block > GetUpdateData().submesh_redecomp_xfer_.shape()[1], + axom::fmt::format( "No transfer object for row {0} and col {1}", max_row_block, max_col_block ) ); - // fill block operator - auto& submesh_fes = submesh_data_.GetSubmeshFESpace(); - auto& parent_trial_fes = *parent_data_.GetParentCoords().ParFESpace(); - // NOTE: we don't call MatrixTransfer::ConvertToHypreParMatrix() because the - // trial space is on the parent mesh, not the submesh - auto J_full = std::make_unique( mpi.MPIComm(), submesh_fes.GetVSize(), - submesh_fes.GlobalVSize(), parent_trial_fes.GlobalVSize(), - submesh_J.GetI(), J.GetData(), submesh_J.GetData(), - submesh_fes.GetDofOffsets(), parent_trial_fes.GetDofOffsets() ); - auto J_true = std::unique_ptr( - mfem::RAP( submesh_fes.Dof_TrueDof_Matrix(), J_full.get(), parent_trial_fes.Dof_TrueDof_Matrix() ) ); - - // Create ones on diagonal of eliminated mortar tdofs, i.e. inactive dofs (CSR sparse matrix -> HypreParMatrix) - // I vector - mfem::Array rows( submesh_fes.GetTrueVSize() + 1 ); - rows = 0; - auto mortar_tdofs_ct = 0; - for ( int i{ 0 }; i < submesh_fes.GetTrueVSize(); ++i ) { - if ( mortar_tdofs_ct < mortar_tdof_list_.Size() && mortar_tdof_list_[mortar_tdofs_ct] == i ) { - ++mortar_tdofs_ct; - } - rows[i + 1] = mortar_tdofs_ct; - } - // J vector - mfem::Array mortar_tdofs( mortar_tdof_list_ ); - // data vector - mfem::Vector ones( mortar_tdofs_ct ); - ones = 1.0; - mfem::SparseMatrix inactive_sm( rows.GetData(), mortar_tdofs.GetData(), ones.GetData(), submesh_fes.GetTrueVSize(), - submesh_fes.GetTrueVSize(), false, false, true ); - auto inactive_hpm = std::make_unique( J_true->GetComm(), J_true->GetGlobalNumRows(), - J_true->GetRowStarts(), &inactive_sm ); - // Have the mfem::HypreParMatrix manage the data pointers - rows.GetMemory().ClearOwnerFlags(); - mortar_tdofs.GetMemory().ClearOwnerFlags(); - ones.GetMemory().ClearOwnerFlags(); - inactive_sm.GetMemoryI().ClearOwnerFlags(); - inactive_sm.GetMemoryJ().ClearOwnerFlags(); - inactive_sm.GetMemoryData().ClearOwnerFlags(); - - block_J->SetBlock( 0, 1, J_true->Transpose() ); - block_J->SetBlock( 1, 0, J_true.release() ); - block_J->SetBlock( 1, 1, inactive_hpm.release() ); - - return block_J; -} + const mfem::Array& row_offsets = ( max_row_block == 0 ) ? disp_offsets_ : block_offsets_; + const mfem::Array& col_offsets = ( max_col_block == 0 ) ? disp_offsets_ : block_offsets_; -// TODO: Merge with GetMfemBlockJacobian() to avoid code duplication -std::unique_ptr MfemJacobianData::GetMfemDfDxFullJacobian( const MethodData& method_data ) const -{ - // create block operator - auto block_J = std::make_unique( block_offsets_ ); + auto block_J = std::make_unique( row_offsets, col_offsets ); block_J->owns_blocks = 1; - // these are Tribol element ids - auto mortar_elems = method_data.getBlockJElementIds()[static_cast( BlockSpace::MORTAR )]; - // convert them to redecomp element ids - const auto& elem_map_1 = parent_data_.GetElemMap1(); - for ( auto& mortar_elem : mortar_elems ) { - mortar_elem = elem_map_1[static_cast( mortar_elem )]; - } + // Map unique (r_blk, c_blk) -> list of (row_space, col_space) pairs + std::map, std::vector>> block_contribs; - // these are Tribol element ids - auto nonmortar_elems = method_data.getBlockJElementIds()[static_cast( BlockSpace::NONMORTAR )]; - // convert them to redecomp element ids - const auto& elem_map_2 = parent_data_.GetElemMap2(); - for ( auto& nonmortar_elem : nonmortar_elems ) { - nonmortar_elem = elem_map_2[static_cast( nonmortar_elem )]; + for ( const auto& r_pair : row_info ) { + for ( const auto& c_pair : col_info ) { + block_contribs[{ r_pair.first, c_pair.first }].push_back( { r_pair.second, c_pair.second } ); + } } - // these are Tribol element ids - auto lm_elems = method_data.getBlockJElementIds()[static_cast( BlockSpace::LAGRANGE_MULTIPLIER )]; - // convert them to redecomp element ids - for ( auto& lm_elem : lm_elems ) { - lm_elem = elem_map_2[static_cast( lm_elem )]; - } + // 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() }; - // transfer (0, 0) block (residual dof rows, displacement dof cols) - auto submesh_J = GetUpdateData().submesh_redecomp_xfer_00_->TransferToParallelSparse( - mortar_elems, mortar_elems, - method_data.getBlockJ()( static_cast( BlockSpace::MORTAR ), static_cast( BlockSpace::MORTAR ) ) ); - submesh_J += GetUpdateData().submesh_redecomp_xfer_00_->TransferToParallelSparse( - mortar_elems, nonmortar_elems, - method_data.getBlockJ()( static_cast( BlockSpace::MORTAR ), static_cast( BlockSpace::NONMORTAR ) ) ); - submesh_J += GetUpdateData().submesh_redecomp_xfer_00_->TransferToParallelSparse( - nonmortar_elems, mortar_elems, - method_data.getBlockJ()( static_cast( BlockSpace::NONMORTAR ), static_cast( BlockSpace::MORTAR ) ) ); - submesh_J += GetUpdateData().submesh_redecomp_xfer_00_->TransferToParallelSparse( - nonmortar_elems, nonmortar_elems, - method_data.getBlockJ()( static_cast( BlockSpace::NONMORTAR ), static_cast( BlockSpace::NONMORTAR ) ) ); - submesh_J.Finalize(); - auto submesh_J_hypre = GetUpdateData().submesh_redecomp_xfer_00_->ConvertToHypreParMatrix( submesh_J, false ); - // Matrix returned by mfem::RAP copies all existing data and owns its data - auto parent_J_hypre = - std::unique_ptr( mfem::RAP( submesh_J_hypre.get(), submesh_parent_vdof_xfer_.get() ) ); - block_J->SetBlock( - 0, 0, mfem::RAP( parent_J_hypre.get(), parent_data_.GetParentCoords().ParFESpace()->Dof_TrueDof_Matrix() ) ); - - // transfer (0, 1) block (residual dof rows, lagrange multiplier dof cols) - submesh_J = GetUpdateData().submesh_redecomp_xfer_01_->TransferToParallelSparse( - mortar_elems, lm_elems, - method_data.getBlockJ()( static_cast( BlockSpace::MORTAR ), - static_cast( BlockSpace::LAGRANGE_MULTIPLIER ) ) ); - submesh_J += GetUpdateData().submesh_redecomp_xfer_01_->TransferToParallelSparse( - nonmortar_elems, lm_elems, - method_data.getBlockJ()( static_cast( BlockSpace::NONMORTAR ), - static_cast( BlockSpace::LAGRANGE_MULTIPLIER ) ) ); - submesh_J.Finalize(); - submesh_J_hypre = GetUpdateData().submesh_redecomp_xfer_01_->ConvertToHypreParMatrix( submesh_J, false ); - // Matrix returned by mfem::ParMult copies row and column starts since last arg is true. All other data is copied and - // owned by the new matrix. - parent_J_hypre = std::unique_ptr( - mfem::ParMult( std::unique_ptr( submesh_parent_vdof_xfer_->Transpose() ).get(), - submesh_J_hypre.get(), true ) ); - block_J->SetBlock( 0, 1, - mfem::RAP( parent_data_.GetParentCoords().ParFESpace()->Dof_TrueDof_Matrix(), parent_J_hypre.get(), - submesh_data_.GetSubmeshFESpace().Dof_TrueDof_Matrix() ) ); - - // transfer (1, 0) block (gap dof rows, displacement dof cols) - submesh_J = GetUpdateData().submesh_redecomp_xfer_10_->TransferToParallelSparse( - lm_elems, mortar_elems, - method_data.getBlockJ()( static_cast( BlockSpace::LAGRANGE_MULTIPLIER ), - static_cast( BlockSpace::MORTAR ) ) ); - submesh_J += GetUpdateData().submesh_redecomp_xfer_10_->TransferToParallelSparse( - lm_elems, nonmortar_elems, - method_data.getBlockJ()( static_cast( BlockSpace::LAGRANGE_MULTIPLIER ), - static_cast( BlockSpace::NONMORTAR ) ) ); - submesh_J.Finalize(); - submesh_J_hypre = GetUpdateData().submesh_redecomp_xfer_10_->ConvertToHypreParMatrix( submesh_J, false ); - // Matrix returned by mfem::ParMult copies row and column starts since last arg is true. All other data is copied and - // owned by the new matrix. - parent_J_hypre = std::unique_ptr( std::unique_ptr( - mfem::ParMult( submesh_J_hypre.get(), submesh_parent_vdof_xfer_.get(), true ) ) ); - block_J->SetBlock( 1, 0, - mfem::RAP( submesh_data_.GetSubmeshFESpace().Dof_TrueDof_Matrix(), parent_J_hypre.get(), - parent_data_.GetParentCoords().ParFESpace()->Dof_TrueDof_Matrix() ) ); - - // Create ones on diagonal of eliminated mortar tdofs, i.e. inactive dofs (CSR sparse matrix -> HypreParMatrix) - // I vector - auto& submesh_fes = submesh_data_.GetSubmeshFESpace(); - mfem::Array rows( submesh_fes.GetTrueVSize() + 1 ); - rows = 0; - auto mortar_tdofs_ct = 0; - for ( int i{ 0 }; i < submesh_fes.GetTrueVSize(); ++i ) { - if ( mortar_tdofs_ct < mortar_tdof_list_.Size() && mortar_tdof_list_[mortar_tdofs_ct] == i ) { - ++mortar_tdofs_ct; - } - rows[i + 1] = mortar_tdofs_ct; - } - // J vector - mfem::Array mortar_tdofs( mortar_tdof_list_ ); - // data vector - mfem::Vector ones( mortar_tdofs_ct ); - ones = 1.0; - mfem::SparseMatrix inactive_sm( rows.GetData(), mortar_tdofs.GetData(), ones.GetData(), submesh_fes.GetTrueVSize(), - submesh_fes.GetTrueVSize(), false, false, true ); - auto inactive_hpm = std::make_unique( TRIBOL_COMM_WORLD, submesh_fes.GlobalTrueVSize(), - submesh_fes.GetTrueDofOffsets(), &inactive_sm ); - // Have the mfem::HypreParMatrix manage the data pointers - rows.GetMemory().SetHostPtrOwner( false ); - mortar_tdofs.GetMemory().SetHostPtrOwner( false ); - ones.GetMemory().SetHostPtrOwner( false ); - inactive_sm.SetDataOwner( false ); - inactive_hpm->SetOwnerFlags( 3, inactive_hpm->OwnsOffd(), inactive_hpm->OwnsColMap() ); - block_J->SetBlock( 1, 1, inactive_hpm.release() ); + // Iterate over unique blocks + for ( const auto& entry : block_contribs ) { + int r_blk = entry.first.first; + int c_blk = entry.first.second; + const auto& contribs = entry.second; - return block_J; -} + std::unique_ptr submesh_J; -// TODO: Merge with GetMfemBlockJacobian() to avoid code duplication -std::unique_ptr MfemJacobianData::GetMfemDfDnJacobian( const MethodData& method_data ) const -{ - // create block operator - auto block_J = std::make_unique( block_offsets_, disp_offsets_ ); - block_J->owns_blocks = 1; + for ( const auto& pair : contribs ) { + BlockSpace rs = pair.first; + BlockSpace cs = pair.second; - // these are Tribol element ids - auto mortar_elems = method_data.getBlockJElementIds()[static_cast( BlockSpace::MORTAR )]; - // convert them to redecomp element ids - const auto& elem_map_1 = parent_data_.GetElemMap1(); - for ( auto& mortar_elem : mortar_elems ) { - mortar_elem = elem_map_1[static_cast( mortar_elem )]; - } + // Get block from method_data + const auto& J_block = method_data.getBlockJ()( static_cast( rs ), static_cast( cs ) ); - // these are Tribol element ids - auto nonmortar_elems = method_data.getBlockJElementIds()[static_cast( BlockSpace::NONMORTAR )]; - // convert them to redecomp element ids - const auto& elem_map_2 = parent_data_.GetElemMap2(); - for ( auto& nonmortar_elem : nonmortar_elems ) { - nonmortar_elem = elem_map_2[static_cast( nonmortar_elem )]; - } + // Map element IDs to redecomp IDs + const auto& row_elem_ids_tribol = method_data.getBlockJElementIds()[static_cast( rs )]; + ArrayT row_redecomp_ids; + row_redecomp_ids.reserve( row_elem_ids_tribol.size() ); + for ( auto id : row_elem_ids_tribol ) { + row_redecomp_ids.push_back( ( *elem_map_by_space[static_cast( rs )] )[static_cast( id )] ); + } - // these are Tribol element ids - auto lm_elems = method_data.getBlockJElementIds()[static_cast( BlockSpace::LAGRANGE_MULTIPLIER )]; - // convert them to redecomp element ids - for ( auto& lm_elem : lm_elems ) { - lm_elem = elem_map_2[static_cast( lm_elem )]; - } + const auto& col_elem_ids_tribol = method_data.getBlockJElementIds()[static_cast( cs )]; + ArrayT col_redecomp_ids; + col_redecomp_ids.reserve( col_elem_ids_tribol.size() ); + for ( auto id : col_elem_ids_tribol ) { + col_redecomp_ids.push_back( ( *elem_map_by_space[static_cast( cs )] )[static_cast( id )] ); + } - // transfer (0, 0) block (residual dof rows, displacement dof cols) - auto submesh_J = GetUpdateData().submesh_redecomp_xfer_00_->TransferToParallelSparse( - mortar_elems, nonmortar_elems, - method_data.getBlockJ()( static_cast( BlockSpace::MORTAR ), static_cast( BlockSpace::NONMORTAR ) ) ); - submesh_J += GetUpdateData().submesh_redecomp_xfer_00_->TransferToParallelSparse( - nonmortar_elems, nonmortar_elems, - method_data.getBlockJ()( static_cast( BlockSpace::NONMORTAR ), static_cast( BlockSpace::NONMORTAR ) ) ); - submesh_J.Finalize(); - auto submesh_J_hypre = GetUpdateData().submesh_redecomp_xfer_00_->ConvertToHypreParMatrix( submesh_J, false ); - // Matrix returned by mfem::RAP copies all existing data and owns its data - auto parent_J_hypre = - std::unique_ptr( mfem::RAP( submesh_J_hypre.get(), submesh_parent_vdof_xfer_.get() ) ); - block_J->SetBlock( - 0, 0, mfem::RAP( parent_J_hypre.get(), parent_data_.GetParentCoords().ParFESpace()->Dof_TrueDof_Matrix() ) ); - - // transfer (1, 0) block (gap dof rows, displacement dof cols) - submesh_J = GetUpdateData().submesh_redecomp_xfer_10_->TransferToParallelSparse( - lm_elems, nonmortar_elems, - method_data.getBlockJ()( static_cast( BlockSpace::LAGRANGE_MULTIPLIER ), - static_cast( BlockSpace::NONMORTAR ) ) ); - submesh_J.Finalize(); - submesh_J_hypre = GetUpdateData().submesh_redecomp_xfer_10_->ConvertToHypreParMatrix( submesh_J, false ); - // Matrix returned by mfem::ParMult copies row and column starts since last arg is true. All other data is copied and - // owned by the new matrix. - parent_J_hypre = std::unique_ptr( - mfem::ParMult( submesh_J_hypre.get(), submesh_parent_vdof_xfer_.get(), true ) ); - block_J->SetBlock( 1, 0, - mfem::RAP( submesh_data_.GetSubmeshFESpace().Dof_TrueDof_Matrix(), parent_J_hypre.get(), - parent_data_.GetParentCoords().ParFESpace()->Dof_TrueDof_Matrix() ) ); + // Pick transfer object + redecomp::MatrixTransfer* xfer = GetUpdateData().submesh_redecomp_xfer_( r_blk, c_blk ).get(); - return block_J; -} + // No transfer object for LAGRANGE_MULTIPLER, LAGRANGE_MULTIPLIER block + if ( xfer != nullptr ) { + auto J_contrib = xfer->TransferToParallelSparse( row_redecomp_ids, col_redecomp_ids, J_block ); + if ( !submesh_J ) { + submesh_J = std::make_unique( std::move( J_contrib ) ); + } else { + ( *submesh_J ) += J_contrib; + } + } + } -// TODO: Merge with GetMfemBlockJacobian() to avoid code duplication -std::unique_ptr MfemJacobianData::GetMfemDnDxJacobian( const MethodData& method_data ) const -{ - // create block operator - auto block_J = std::make_unique( disp_offsets_, disp_offsets_ ); - block_J->owns_blocks = 1; + if ( submesh_J ) { + submesh_J->Finalize(); + + // Pick xfer again for conversion + redecomp::MatrixTransfer* xfer = GetUpdateData().submesh_redecomp_xfer_( r_blk, c_blk ).get(); + + auto submesh_J_hypre = xfer->ConvertToHypreParMatrix( *submesh_J, false ); + + mfem::HypreParMatrix* block_mat = nullptr; + + if ( r_blk == 0 && c_blk == 0 ) { + ParSparseMatView submesh_J_view( submesh_J_hypre.get() ); + auto parent_J = submesh_J_view.RAP( *submesh_parent_vdof_xfer_ ); + ParSparseMatView parent_P( parent_data_.GetParentCoords().ParFESpace()->Dof_TrueDof_Matrix() ); + block_mat = parent_J.RAP( parent_P ).release(); + } else if ( r_blk == 0 && c_blk == 1 ) { + auto parent_J = submesh_parent_vdof_xfer_->transpose() * submesh_J_hypre.get(); + block_mat = ParSparseMat::RAP( parent_data_.GetParentCoords().ParFESpace()->Dof_TrueDof_Matrix(), parent_J, + submesh_data_.GetSubmeshFESpace().Dof_TrueDof_Matrix() ) + .release(); + } else if ( r_blk == 1 && c_blk == 0 ) { + ParSparseMatView submesh_J_view( submesh_J_hypre.get() ); + auto parent_J = submesh_J_view * ( *submesh_parent_vdof_xfer_ ); + ParSparseMatView submesh_P( submesh_data_.GetSubmeshFESpace().Dof_TrueDof_Matrix() ); + ParSparseMatView parent_P( parent_data_.GetParentCoords().ParFESpace()->Dof_TrueDof_Matrix() ); + block_mat = ParSparseMat::RAP( submesh_P, parent_J, parent_P ).release(); + } - // these are Tribol element ids - auto nonmortar_elems = method_data.getBlockJElementIds()[static_cast( BlockSpace::NONMORTAR )]; - // convert them to redecomp element ids - const auto& elem_map_2 = parent_data_.GetElemMap2(); - for ( auto& nonmortar_elem : nonmortar_elems ) { - nonmortar_elem = elem_map_2[static_cast( nonmortar_elem )]; + block_J->SetBlock( r_blk, c_blk, block_mat ); + } } - // transfer (0, 0) block (residual dof rows, displacement dof cols) - auto submesh_J = GetUpdateData().submesh_redecomp_xfer_00_->TransferToParallelSparse( - nonmortar_elems, nonmortar_elems, - method_data.getBlockJ()( static_cast( BlockSpace::NONMORTAR ), static_cast( BlockSpace::NONMORTAR ) ) ); - submesh_J.Finalize(); - auto submesh_J_hypre = GetUpdateData().submesh_redecomp_xfer_00_->ConvertToHypreParMatrix( submesh_J, false ); - // Matrix returned by mfem::RAP copies all existing data and owns its data - auto parent_J_hypre = - std::unique_ptr( mfem::RAP( submesh_J_hypre.get(), submesh_parent_vdof_xfer_.get() ) ); - block_J->SetBlock( - 0, 0, mfem::RAP( parent_J_hypre.get(), parent_data_.GetParentCoords().ParFESpace()->Dof_TrueDof_Matrix() ) ); + // Handle Inactive DOFs for (1, 1) + bool has_11 = false; + for ( auto rb : row_info ) + if ( rb.first == 1 ) has_11 = true; + bool col_has_1 = false; + for ( auto cb : col_info ) + if ( cb.first == 1 ) col_has_1 = true; + has_11 = has_11 && col_has_1; + + if ( has_11 ) { + auto& submesh_fes_full = submesh_data_.GetSubmeshFESpace(); + ParSparseMat inactive_hpm_full = + ParSparseMat::diagonalMatrix( TRIBOL_COMM_WORLD, 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() ); + } + } return block_J; } MfemJacobianData::UpdateData::UpdateData( const MfemMeshData& parent_data, const MfemSubmeshData& submesh_data ) + : submesh_redecomp_xfer_( 2, 2 ) { auto dual_submesh_fes = &submesh_data.GetSubmeshFESpace(); auto primal_submesh_fes = &parent_data.GetSubmeshFESpace(); @@ -1301,13 +1061,13 @@ MfemJacobianData::UpdateData::UpdateData( const MfemMeshData& parent_data, const primal_submesh_fes = parent_data.GetLORMeshFESpace(); } // create a matrix transfer operator for moving data from redecomp to the submesh - submesh_redecomp_xfer_00_ = std::make_unique( + submesh_redecomp_xfer_( 0, 0 ) = std::make_unique( *primal_submesh_fes, *primal_submesh_fes, *parent_data.GetRedecompResponse().FESpace(), *parent_data.GetRedecompResponse().FESpace() ); - submesh_redecomp_xfer_01_ = std::make_unique( *primal_submesh_fes, *dual_submesh_fes, - *parent_data.GetRedecompResponse().FESpace(), - *submesh_data.GetRedecompGap().FESpace() ); - submesh_redecomp_xfer_10_ = std::make_unique( + submesh_redecomp_xfer_( 0, 1 ) = std::make_unique( + *primal_submesh_fes, *dual_submesh_fes, *parent_data.GetRedecompResponse().FESpace(), + *submesh_data.GetRedecompGap().FESpace() ); + submesh_redecomp_xfer_( 1, 0 ) = std::make_unique( *dual_submesh_fes, *primal_submesh_fes, *submesh_data.GetRedecompGap().FESpace(), *parent_data.GetRedecompResponse().FESpace() ); } diff --git a/src/tribol/mesh/MfemData.hpp b/src/tribol/mesh/MfemData.hpp index de2846d8..17798cce 100644 --- a/src/tribol/mesh/MfemData.hpp +++ b/src/tribol/mesh/MfemData.hpp @@ -12,6 +12,7 @@ #ifdef BUILD_REDECOMP #include +#include #include #include "mfem.hpp" @@ -22,6 +23,7 @@ #include "tribol/common/BasicTypes.hpp" #include "tribol/common/Parameters.hpp" #include "tribol/mesh/MethodCouplingData.hpp" +#include "tribol/utils/ParSparseMat.hpp" namespace tribol { @@ -1649,37 +1651,19 @@ class MfemJacobianData { void UpdateJacobianXfer(); /** - * @brief Returns symmetric, off-diagonal Jacobian contributions as an mfem::BlockOperator + * @brief Returns a Jacobian as an mfem::BlockOperator * * @param method_data Method data holding element Jacobians + * @param row_info List of {block_row_index, BlockSpace} pairs. Since a single block row in the output matrix might + * aggregate DOFs from multiple Tribol spaces (e.g. Mortar and NonMortar spaces might both map to the Displacement + * block 0), this vector defines the mapping from each Tribol space to its corresponding block row index. + * @param col_info List of {block_col_index, BlockSpace} pairs. Similar to row_info, this defines the mapping from + * each Tribol space to its corresponding block column index. * @return std::unique_ptr */ - std::unique_ptr GetMfemBlockJacobian( const MethodData* method_data ) const; - - /** - * @brief Returns full, potentially non-symmetric derivative of the force w.r.t. nodal coordinates as an - * mfem::BlockOperator - * - * @param method_data Method data holding element Jacobians - * @return std::unique_ptr - */ - std::unique_ptr GetMfemDfDxFullJacobian( const MethodData& method_data ) const; - - /** - * @brief Returns the derivative of the force w.r.t. the normal direction as an mfem::BlockOperator - * - * @param method_data Method data holding element Jacobians - * @return std::unique_ptr - */ - std::unique_ptr GetMfemDfDnJacobian( const MethodData& method_data ) const; - - /** - * @brief Returns the derivative of the normal direction w.r.t. the nodal coordinates as an mfem::BlockOperator - * - * @param method_data Method data holding element Jacobians - * @return std::unique_ptr - */ - std::unique_ptr GetMfemDnDxJacobian( const MethodData& method_data ) const; + std::unique_ptr GetMfemBlockJacobian( + const MethodData& method_data, const std::vector>& row_info, + const std::vector>& col_info ) const; private: /** @@ -1696,19 +1680,11 @@ class MfemJacobianData { UpdateData( const MfemMeshData& parent_data, const MfemSubmeshData& submesh_data ); /** - * @brief Redecomp to parent-linked boundary submesh transfer operator, (displacement, displacement) block - */ - std::unique_ptr submesh_redecomp_xfer_00_; - - /** - * @brief Redecomp to parent-linked boundary submesh transfer operator, (displacement, pressure) block - */ - std::unique_ptr submesh_redecomp_xfer_01_; - - /** - * @brief Redecomp to parent-linked boundary submesh transfer operator, (pressure, displacement) block + * @brief Redecomp to parent-linked boundary submesh transfer operators + * + * @note Indexed by (row_block, col_block) */ - std::unique_ptr submesh_redecomp_xfer_10_; + Array2D> submesh_redecomp_xfer_; }; /** @@ -1755,7 +1731,7 @@ class MfemJacobianData { /** * @brief Submesh to parent transfer operator */ - std::unique_ptr submesh_parent_vdof_xfer_; + std::unique_ptr submesh_parent_vdof_xfer_; /** * @brief List of submesh true dofs that only exist on the mortar surface diff --git a/src/tribol/utils/ParSparseMat.cpp b/src/tribol/utils/ParSparseMat.cpp new file mode 100644 index 00000000..d2669014 --- /dev/null +++ b/src/tribol/utils/ParSparseMat.cpp @@ -0,0 +1,239 @@ +// 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/utils/ParSparseMat.hpp" +#include + +namespace tribol { + +// ParSparseMatView implementations + +ParSparseMat operator+( const ParSparseMatView& lhs, const ParSparseMatView& rhs ) +{ + mfem::HypreParMatrix* result = mfem::Add( 1.0, *lhs.mat_, 1.0, *rhs.mat_ ); + return ParSparseMat( result ); +} + +ParSparseMat operator-( const ParSparseMatView& lhs, const ParSparseMatView& rhs ) +{ + mfem::HypreParMatrix* result = mfem::Add( 1.0, *lhs.mat_, -1.0, *rhs.mat_ ); + return ParSparseMat( result ); +} + +ParSparseMat ParSparseMatView::operator*( double s ) const +{ + mfem::HypreParMatrix* result = mfem::Add( s, *mat_, 0.0, *mat_ ); + return ParSparseMat( result ); +} + +ParSparseMat operator*( const ParSparseMatView& lhs, const ParSparseMatView& rhs ) +{ + mfem::HypreParMatrix* result = mfem::ParMult( lhs.mat_, rhs.mat_ ); + result->CopyRowStarts(); + result->CopyColStarts(); + return ParSparseMat( result ); +} + +ParVector ParSparseMatView::operator*( const ParVectorView& x ) const +{ + ParVector y( *mat_ ); + mat_->Mult( const_cast( x.get() ), y.get() ); + return y; +} + +ParSparseMat ParSparseMatView::transpose() const { return ParSparseMat( mat_->Transpose() ); } + +ParSparseMat ParSparseMatView::square() const { return *this * *this; } + +ParSparseMat ParSparseMatView::RAP( const ParSparseMatView& P ) const +{ + return ParSparseMat( mfem::RAP( mat_, P.mat_ ) ); +} + +ParSparseMat ParSparseMatView::RAP( const ParSparseMatView& A, const ParSparseMatView& P ) +{ + return ParSparseMat( mfem::RAP( A.mat_, P.mat_ ) ); +} + +ParSparseMat ParSparseMatView::RAP( const ParSparseMatView& R, const ParSparseMatView& A, const ParSparseMatView& P ) +{ + return ParSparseMat( mfem::RAP( R.mat_, A.mat_, P.mat_ ) ); +} + +void ParSparseMatView::EliminateRows( const mfem::Array& rows ) { mat_->EliminateRows( rows ); } + +ParSparseMat ParSparseMatView::EliminateCols( const mfem::Array& cols ) +{ + return ParSparseMat( mat_->EliminateCols( cols ) ); +} + +ParSparseMat operator*( double s, const ParSparseMatView& mat ) { return mat * s; } + +ParVector operator*( const ParVectorView& x, const ParSparseMatView& mat ) +{ + ParVector y( *mat.mat_, 1 ); + mat.mat_->MultTranspose( const_cast( x.get() ), y.get() ); + return y; +} + +// ParSparseMat implementations + +ParSparseMat::ParSparseMat( mfem::HypreParMatrix* mat ) : ParSparseMatView( mat ), owned_mat_( mat ) {} + +ParSparseMat::ParSparseMat( std::unique_ptr mat ) + : ParSparseMatView( mat.get() ), owned_mat_( std::move( mat ) ) +{ +} + +ParSparseMat::ParSparseMat( MPI_Comm comm, HYPRE_BigInt glob_size, HYPRE_BigInt* row_starts, mfem::SparseMatrix&& diag ) + : ParSparseMatView( nullptr ) +{ + // ParSparseMat is host only now. Make sure CSR data is copied on host in the constructor. + HYPRE_MemoryLocation old_hypre_mem_location; + HYPRE_GetMemoryLocation( &old_hypre_mem_location ); + HYPRE_SetMemoryLocation( HYPRE_MEMORY_HOST ); + owned_mat_ = std::make_unique( comm, glob_size, row_starts, &diag ); + mat_ = owned_mat_.get(); + diag.GetMemoryI().ClearOwnerFlags(); + diag.GetMemoryJ().ClearOwnerFlags(); + diag.GetMemoryData().ClearOwnerFlags(); + auto mfem_owned_arrays = 3; + owned_mat_->SetOwnerFlags( mfem_owned_arrays, owned_mat_->OwnsOffd(), owned_mat_->OwnsColMap() ); + // Return hypre's memory location to what it was before + HYPRE_SetMemoryLocation( old_hypre_mem_location ); +} + +ParSparseMat::ParSparseMat( ParSparseMat&& other ) noexcept + : ParSparseMatView( other.owned_mat_.get() ), owned_mat_( std::move( other.owned_mat_ ) ) +{ + other.mat_ = nullptr; +} + +ParSparseMat& ParSparseMat::operator=( ParSparseMat&& other ) noexcept +{ + if ( this != &other ) { + owned_mat_ = std::move( other.owned_mat_ ); + mat_ = owned_mat_.get(); + other.mat_ = nullptr; + } + return *this; +} + +mfem::HypreParMatrix* ParSparseMat::release() +{ + mat_ = nullptr; + return owned_mat_.release(); +} + +ParSparseMat& ParSparseMat::operator+=( const ParSparseMatView& other ) +{ + *this = *this + other; + return *this; +} + +ParSparseMat& ParSparseMat::operator-=( const ParSparseMatView& other ) +{ + *this = *this - other; + return *this; +} + +ParSparseMat& ParSparseMat::operator*=( const ParSparseMatView& other ) +{ + *this = *this * other; + return *this; +} + +ParSparseMat ParSparseMat::diagonalMatrix( MPI_Comm comm, HYPRE_BigInt global_size, + const mfem::Array& row_starts, double diag_val, + const mfem::Array& ordered_rows, bool skip_rows ) +{ + int num_local_rows = 0; + if ( HYPRE_AssumedPartitionCheck() ) { + num_local_rows = static_cast( row_starts[1] - row_starts[0] ); + } else { + int rank; + MPI_Comm_rank( comm, &rank ); + num_local_rows = static_cast( row_starts[rank + 1] - row_starts[rank] ); + } + + // NOTE: mfem::HypreParMatrix(MPI_Comm, HYPRE_BigInt, HYPRE_BigInt*, SparseMatrix*) does not take ownership + // of the provided CSR arrays (see MFEM docs). To avoid dangling pointers and allocator mismatches, build + // the ParCSR data using the HypreParMatrix constructor that takes ownership of raw arrays allocated with + // new[]. + + const int num_ordered_rows = ordered_rows.Size(); + if ( num_local_rows < 0 ) { + num_local_rows = 0; + } + + // Count selected diagonal entries in a first pass (do not rely on ordered_rows being unique/in-range). + HYPRE_Int num_diag_entries = 0; + int ordered_idx = 0; + for ( int i = 0; i < num_local_rows; ++i ) { + while ( ordered_idx < num_ordered_rows && ordered_rows[ordered_idx] < i ) { + ++ordered_idx; + } + const bool is_in_list = ( ordered_idx < num_ordered_rows && ordered_rows[ordered_idx] == i ); + const bool add_entry = skip_rows ? !is_in_list : is_in_list; + if ( add_entry ) { + ++num_diag_entries; + } + } + + auto* diag_i = new HYPRE_Int[num_local_rows + 1]; + auto* diag_j = ( num_diag_entries > 0 ) ? new HYPRE_Int[num_diag_entries] : nullptr; + auto* diag_data = ( num_diag_entries > 0 ) ? new mfem::real_t[num_diag_entries] : nullptr; + + // No off-diagonal entries for a purely diagonal matrix. + auto* offd_i = new HYPRE_Int[num_local_rows + 1]; + auto* offd_j = static_cast( nullptr ); + auto* offd_data = static_cast( nullptr ); + auto* offd_col_map = static_cast( nullptr ); + + diag_i[0] = 0; + for ( int i = 0; i < num_local_rows + 1; ++i ) { + offd_i[i] = 0; + } + + HYPRE_Int diag_entry_ct = 0; + ordered_idx = 0; + for ( int i = 0; i < num_local_rows; ++i ) { + while ( ordered_idx < num_ordered_rows && ordered_rows[ordered_idx] < i ) { + ++ordered_idx; + } + const bool is_in_list = ( ordered_idx < num_ordered_rows && ordered_rows[ordered_idx] == i ); + const bool add_entry = skip_rows ? !is_in_list : is_in_list; + + if ( add_entry ) { + diag_j[diag_entry_ct] = static_cast( i ); + diag_data[diag_entry_ct] = diag_val; + ++diag_entry_ct; + } + diag_i[i + 1] = diag_entry_ct; + } + + // copy row_starts to a new array + mfem::Array row_starts_copy = row_starts; + auto diag_hpm = std::make_unique( comm, global_size, global_size, row_starts_copy.GetData(), + row_starts_copy.GetData(), diag_i, diag_j, diag_data, offd_i, + offd_j, offd_data, 0, offd_col_map, true ); + diag_hpm->CopyRowStarts(); + diag_hpm->CopyColStarts(); + auto mfem_owned_arrays = 3; + diag_hpm->SetOwnerFlags( mfem_owned_arrays, diag_hpm->OwnsOffd(), diag_hpm->OwnsColMap() ); + return ParSparseMat( std::move( diag_hpm ) ); +} + +ParSparseMat ParSparseMat::diagonalMatrix( MPI_Comm comm, HYPRE_BigInt global_size, HYPRE_BigInt* row_starts, + double diag_val, const mfem::Array& ordered_rows, bool skip_rows ) +{ + int num_procs; + MPI_Comm_size( comm, &num_procs ); + int n_row_starts = HYPRE_AssumedPartitionCheck() ? 3 : num_procs + 1; + mfem::Array row_starts_array( row_starts, n_row_starts ); + return diagonalMatrix( comm, global_size, row_starts_array, diag_val, ordered_rows, skip_rows ); +} + +} // namespace tribol diff --git a/src/tribol/utils/ParSparseMat.hpp b/src/tribol/utils/ParSparseMat.hpp new file mode 100644 index 00000000..9f33967d --- /dev/null +++ b/src/tribol/utils/ParSparseMat.hpp @@ -0,0 +1,258 @@ +// 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) + +#ifndef SRC_TRIBOL_UTILS_PARSPARSEMAT_HPP_ +#define SRC_TRIBOL_UTILS_PARSPARSEMAT_HPP_ + +#include "tribol/config.hpp" +#include "mfem.hpp" +#include "tribol/utils/ParVector.hpp" + +#include +#include + +namespace tribol { + +class ParSparseMat; + +/** + * @brief Non-owning view of a mfem::HypreParMatrix + * + * This class holds a raw pointer to a mfem::HypreParMatrix and provides algebraic operations. + * It does not manage the lifetime of the matrix. + */ +class ParSparseMatView { + public: + /** + * @brief Construct from a mfem::HypreParMatrix pointer + * + * @param mat Pointer to the mfem HypreParMatrix + */ + ParSparseMatView( mfem::HypreParMatrix* mat ) : mat_( mat ) {} + + virtual ~ParSparseMatView() = default; + + /** + * @brief Access the underlying mfem::HypreParMatrix + */ + mfem::HypreParMatrix& get() { return *mat_; } + + /** + * @brief Access the underlying mfem::HypreParMatrix (const) + */ + const mfem::HypreParMatrix& get() const { return *mat_; } + + /** + * @brief Access underlying matrix members via arrow operator + */ + mfem::HypreParMatrix* operator->() { return mat_; } + + /** + * @brief Access underlying matrix members via arrow operator (const) + */ + const mfem::HypreParMatrix* operator->() const { return mat_; } + + /** + * @brief Returns the number of local rows + */ + int Height() const { return mat_->Height(); } + + /** + * @brief Returns the number of local columns + */ + int Width() const { return mat_->Width(); } + + /** + * @brief Matrix addition: returns A + B + */ + friend ParSparseMat operator+( const ParSparseMatView& lhs, const ParSparseMatView& rhs ); + + /** + * @brief Matrix subtraction: returns A - B + */ + friend ParSparseMat operator-( const ParSparseMatView& lhs, const ParSparseMatView& rhs ); + + /** + * @brief Matrix scalar multiplication: returns s * A + */ + ParSparseMat operator*( double s ) const; + + /** + * @brief Matrix multiplication: returns A * B + */ + friend ParSparseMat operator*( const ParSparseMatView& lhs, const ParSparseMatView& rhs ); + + /** + * @brief Matrix-vector multiplication: returns y = A * x + */ + ParVector operator*( const ParVectorView& x ) const; + + /** + * @brief Returns the transpose of the matrix + */ + ParSparseMat transpose() const; + + /** + * @brief Returns the square of the matrix (A * A) + */ + ParSparseMat square() const; + + /** + * @brief Returns P^T * A * P + */ + ParSparseMat RAP( const ParSparseMatView& P ) const; + + /** + * @brief Returns P^T * A * P + */ + static ParSparseMat RAP( const ParSparseMatView& A, const ParSparseMatView& P ); + + /** + * @brief Returns R * A * P + */ + static ParSparseMat RAP( const ParSparseMatView& R, const ParSparseMatView& A, const ParSparseMatView& P ); + + /** + * @brief Eliminates the rows from the matrix + * + * @param rows Array of rows to eliminate + */ + void EliminateRows( const mfem::Array& rows ); + + /** + * @brief Eliminates the columns from the matrix + * + * @param cols Array of columns to eliminate + */ + ParSparseMat EliminateCols( const mfem::Array& cols ); + + /** + * @brief Scalar-Matrix multiplication: returns s * A + */ + friend ParSparseMat operator*( double s, const ParSparseMatView& mat ); + + /** + * @brief Vector-Matrix multiplication: returns y = x^T * A (computed as A^T * x) + */ + friend ParVector operator*( const ParVectorView& x, const ParSparseMatView& mat ); + + protected: + mfem::HypreParMatrix* mat_; +}; + +/** + * @brief Wrapper class for mfem::HypreParMatrix to provide convenience operators + * + * This class owns a mfem::HypreParMatrix via a unique_ptr and adds support for + * algebraic operations like addition and scalar multiplication. + */ +class ParSparseMat : public ParSparseMatView { + public: + /** + * @brief Construct from a mfem::HypreParMatrix pointer and take ownership + * + * @param mat Pointer to the mfem HypreParMatrix + */ + explicit ParSparseMat( mfem::HypreParMatrix* mat ); + + /** + * @brief Construct from a mfem::HypreParMatrix pointer and take ownership + * + * @param mat Pointer to the mfem HypreParMatrix + */ + explicit ParSparseMat( std::unique_ptr mat ); + + /** + * @brief Construct from MPI communicator, global size, row_starts, and mfem::SparseMatrix rvalue + * + * @param comm MPI communicator + * @param glob_size Global number of rows (and columns) + * @param row_starts Global row 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 glob_size, HYPRE_BigInt* row_starts, mfem::SparseMatrix&& diag ); + + /// Template constructor forwarding arguments to mfem::HypreParMatrix constructor + template + explicit ParSparseMat( Args&&... args ) + : ParSparseMatView( nullptr ), + owned_mat_( std::make_unique( std::forward( args )... ) ) + { + mat_ = owned_mat_.get(); + } + + /// Move constructor + ParSparseMat( ParSparseMat&& other ) noexcept; + + /// Move assignment + ParSparseMat& operator=( ParSparseMat&& other ) noexcept; + + // Disable copy constructor and assignment + ParSparseMat( const ParSparseMat& ) = delete; + ParSparseMat& operator=( const ParSparseMat& ) = delete; + + /** + * @brief Access and release ownership of the HypreParMatrix pointer. The caller is now resposible for releasing the + * memory. + */ + mfem::HypreParMatrix* release(); + + /** + * @brief Matrix in-place addition: A += B + */ + ParSparseMat& operator+=( const ParSparseMatView& other ); + + /** + * @brief Matrix in-place subtraction: A -= B + */ + ParSparseMat& operator-=( const ParSparseMatView& other ); + + /** + * @brief Matrix in-place multiplication: A *= B + */ + ParSparseMat& operator*=( const ParSparseMatView& other ); + + /** + * @brief Returns a diagonal matrix with the given diagonal value + * + * @param comm MPI communicator + * @param global_size Global size of the matrix (rows and columns) + * @param row_starts Row partitioning (global offsets) + * @param diag_val Value for the diagonal entries + * @param ordered_rows Sorted array of local row indices. Defaults to empty. + * @param skip_rows If true (default), ordered_rows are skipped (zero entries). If false, ordered_rows are the only + * entries. + * @return ParSparseMat The constructed diagonal matrix + */ + static ParSparseMat diagonalMatrix( MPI_Comm comm, HYPRE_BigInt global_size, + const mfem::Array& row_starts, double diag_val, + const mfem::Array& ordered_rows = mfem::Array(), + bool skip_rows = true ); + + /** + * @brief Returns a diagonal matrix with the given diagonal value + * + * @param comm MPI communicator + * @param global_size Global size of the matrix (rows and columns) + * @param row_starts Row partitioning (global offsets) + * @param diag_val Value for the diagonal entries + * @param ordered_rows Sorted array of local row indices. Defaults to empty. + * @param skip_rows If true (default), ordered_rows are skipped (zero entries). If false, ordered_rows are the only + * entries. + * @return ParSparseMat The constructed diagonal matrix + */ + static ParSparseMat diagonalMatrix( MPI_Comm comm, HYPRE_BigInt global_size, HYPRE_BigInt* row_starts, + double diag_val, const mfem::Array& ordered_rows = mfem::Array(), + bool skip_rows = true ); + + private: + std::unique_ptr owned_mat_; +}; + +} // namespace tribol + +#endif /* SRC_TRIBOL_UTILS_PARSPARSEMAT_HPP_ */ diff --git a/src/tribol/utils/ParVector.cpp b/src/tribol/utils/ParVector.cpp new file mode 100644 index 00000000..3f548046 --- /dev/null +++ b/src/tribol/utils/ParVector.cpp @@ -0,0 +1,139 @@ +// 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/utils/ParVector.hpp" + +#include "axom/slic.hpp" + +#include "tribol/common/BasicTypes.hpp" + +namespace tribol { + +ParVector::ParVector( mfem::HypreParVector* vec ) : ParVectorView( vec ), owned_vec_( vec ) {} + +ParVector::ParVector( std::unique_ptr vec ) + : ParVectorView( vec.get() ), owned_vec_( std::move( vec ) ) +{ +} + +ParVector::ParVector( ParVector&& other ) noexcept + : ParVectorView( other.owned_vec_.get() ), owned_vec_( std::move( other.owned_vec_ ) ) +{ + other.vec_ = nullptr; +} + +ParVector& ParVector::operator=( ParVector&& other ) noexcept +{ + if ( this != &other ) { + owned_vec_ = std::move( other.owned_vec_ ); + vec_ = owned_vec_.get(); + other.vec_ = nullptr; + } + return *this; +} + +ParVector::ParVector( const ParVector& other ) + : ParVectorView( nullptr ), owned_vec_( std::make_unique( *other.vec_ ) ) +{ + vec_ = owned_vec_.get(); +} + +ParVector& ParVector::operator=( const ParVector& other ) +{ + if ( this != &other ) { + owned_vec_ = std::make_unique( *other.vec_ ); + vec_ = owned_vec_.get(); + } + return *this; +} + +mfem::HypreParVector* ParVector::release() +{ + vec_ = nullptr; + return owned_vec_.release(); +} + +ParVector operator+( const ParVectorView& lhs, const ParVectorView& rhs ) +{ + ParVector result( new mfem::HypreParVector( lhs.get() ) ); + result.get().Add( 1.0, rhs.get() ); + return result; +} + +ParVector operator-( const ParVectorView& lhs, const ParVectorView& rhs ) +{ + ParVector result( new mfem::HypreParVector( lhs.get() ) ); + result.get().Add( -1.0, rhs.get() ); + return result; +} + +ParVector ParVectorView::operator*( double s ) const +{ + ParVector result( new mfem::HypreParVector( *vec_ ) ); + result.get() *= s; + return result; +} + +ParVector operator*( double s, const ParVectorView& vec ) { return vec * s; } + +ParVector& ParVector::operator+=( const ParVectorView& other ) +{ + vec_->Add( 1.0, other.get() ); + return *this; +} + +ParVector& ParVector::operator-=( const ParVectorView& other ) +{ + vec_->Add( -1.0, other.get() ); + return *this; +} + +ParVector& ParVector::operator*=( double s ) +{ + *vec_ *= s; + return *this; +} + +ParVector ParVectorView::multiply( const ParVectorView& other ) const +{ + ParVector result( new mfem::HypreParVector( *vec_ ) ); + result.multiplyInPlace( other ); + return result; +} + +ParVector ParVectorView::divide( const ParVectorView& other ) const +{ + ParVector result( new mfem::HypreParVector( *vec_ ) ); + result.divideInPlace( other ); + return result; +} + +ParVector& ParVector::multiplyInPlace( const ParVectorView& other ) +{ + SLIC_ASSERT( vec_->Size() == other.get().Size() ); + int n = vec_->Size(); + if ( n > 0 ) { + bool use_device = vec_->UseDevice() || other.get().UseDevice(); + auto d_vec = vec_->ReadWrite( use_device ); + auto d_other = other.get().Read( use_device ); + mfem::forall_switch( use_device, n, [=] TRIBOL_HOST_DEVICE( int i ) { d_vec[i] *= d_other[i]; } ); + } + return *this; +} + +ParVector& ParVector::divideInPlace( const ParVectorView& other ) +{ + SLIC_ASSERT( vec_->Size() == other.get().Size() ); + int n = vec_->Size(); + if ( n > 0 ) { + bool use_device = vec_->UseDevice() || other.get().UseDevice(); + auto d_vec = vec_->ReadWrite( use_device ); + auto d_other = other.get().Read( use_device ); + mfem::forall_switch( use_device, n, [=] TRIBOL_HOST_DEVICE( int i ) { d_vec[i] /= d_other[i]; } ); + } + return *this; +} + +} // namespace tribol diff --git a/src/tribol/utils/ParVector.hpp b/src/tribol/utils/ParVector.hpp new file mode 100644 index 00000000..84cd3a67 --- /dev/null +++ b/src/tribol/utils/ParVector.hpp @@ -0,0 +1,204 @@ +// 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) + +#ifndef SRC_TRIBOL_UTILS_PARVECTOR_HPP_ +#define SRC_TRIBOL_UTILS_PARVECTOR_HPP_ + +#include "tribol/config.hpp" +#include "mfem.hpp" + +#include +#include + +namespace tribol { + +class ParVector; + +/** + * @brief Non-owning view of a mfem::HypreParVector + * + * This class holds a raw pointer to a mfem::HypreParVector and provides algebraic operations. + * It does not manage the lifetime of the vector. + */ +class ParVectorView { + public: + /** + * @brief Construct from a mfem::HypreParVector pointer + * + * @param vec Pointer to the mfem HypreParVector + */ + ParVectorView( mfem::HypreParVector* vec ) : vec_( vec ) {} + + virtual ~ParVectorView() = default; + + /** + * @brief Access the underlying mfem::HypreParVector + */ + mfem::HypreParVector& get() { return *vec_; } + + /** + * @brief Access the underlying mfem::HypreParVector (const) + */ + const mfem::HypreParVector& get() const { return *vec_; } + + /** + * @brief Access underlying vector members via arrow operator + */ + mfem::HypreParVector* operator->() { return vec_; } + + /** + * @brief Access underlying vector members via arrow operator (const) + */ + const mfem::HypreParVector* operator->() const { return vec_; } + + /** + * @brief Access local vector entry + */ + mfem::real_t& operator[]( int i ) { return ( *vec_ )[i]; } + + /** + * @brief Access local vector entry (const) + */ + const mfem::real_t& operator[]( int i ) const { return ( *vec_ )[i]; } + + /** + * @brief Sets all entries of the vector to the given value + * + * @param val Value to set + */ + void Fill( double val ) { *vec_ = val; } + + /** + * @brief Returns the local size of the vector + */ + int Size() const { return vec_->Size(); } + + /** + * @brief Returns the maximum value in the vector + */ + mfem::real_t Max() const { return vec_->Max(); } + + /** + * @brief Returns the minimum value in the vector + */ + mfem::real_t Min() const { return vec_->Min(); } + + /** + * @brief Component-wise multiplication: returns z[i] = x[i] * y[i] + */ + ParVector multiply( const ParVectorView& other ) const; + + /** + * @brief Component-wise division: returns z[i] = x[i] / y[i] + */ + ParVector divide( const ParVectorView& other ) const; + + /** + * @brief Vector addition: returns x + y + */ + friend ParVector operator+( const ParVectorView& lhs, const ParVectorView& rhs ); + + /** + * @brief Vector subtraction: returns x - y + */ + friend ParVector operator-( const ParVectorView& lhs, const ParVectorView& rhs ); + + /** + * @brief Vector scalar multiplication: returns s * x + */ + ParVector operator*( double s ) const; + + /** + * @brief Scalar-Vector multiplication: returns s * x + */ + friend ParVector operator*( double s, const ParVectorView& vec ); + + protected: + mfem::HypreParVector* vec_; +}; + +/** + * @brief Wrapper class for mfem::HypreParVector to provide convenience operators + * + * This class owns a mfem::HypreParVector via a unique_ptr and adds support for + * algebraic operations. + */ +class ParVector : public ParVectorView { + public: + /** + * @brief Construct from a mfem::HypreParVector pointer and take ownership + * + * @param vec Pointer to the mfem HypreParVector + */ + explicit ParVector( mfem::HypreParVector* vec ); + + /** + * @brief Construct from a mfem::HypreParVector pointer and take ownership + * + * @param vec Pointer to the mfem HypreParVector + */ + explicit ParVector( std::unique_ptr vec ); + + /// Template constructor forwarding arguments to mfem::HypreParVector constructor + template + explicit ParVector( Args&&... args ) + : ParVectorView( nullptr ), owned_vec_( std::make_unique( std::forward( args )... ) ) + { + vec_ = owned_vec_.get(); + } + + /// Move constructor + ParVector( ParVector&& other ) noexcept; + + /// Move assignment + ParVector& operator=( ParVector&& other ) noexcept; + + /// Copy constructor + ParVector( const ParVector& other ); + + /// Copy constructor (non-const; prevents the template constructor from trying to match this case) + ParVector( ParVector& other ) : ParVector( static_cast( other ) ) {} + + /// Copy assignment + ParVector& operator=( const ParVector& other ); + + /** + * @brief Access and release ownership of the HypreParVector pointer. The caller is now responsible for releasing the + * memory. + */ + mfem::HypreParVector* release(); + + /** + * @brief Vector in-place addition: x += y + */ + ParVector& operator+=( const ParVectorView& other ); + + /** + * @brief Vector in-place subtraction: x -= y + */ + ParVector& operator-=( const ParVectorView& other ); + + /** + * @brief Vector in-place multiplication: x *= s + */ + ParVector& operator*=( double s ); + + /** + * @brief Component-wise in-place multiplication: x[i] *= y[i] + */ + ParVector& multiplyInPlace( const ParVectorView& other ); + + /** + * @brief Component-wise in-place division: x[i] /= y[i] + */ + ParVector& divideInPlace( const ParVectorView& other ); + + private: + std::unique_ptr owned_vec_; +}; + +} // namespace tribol + +#endif /* SRC_TRIBOL_UTILS_PARVECTOR_HPP_ */