phynexis-deps/mfem
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge remote-tracking branch 'origin/master' into stefanozampini/petsc-3.21

Browse Source
This commit is contained in:
Stefano Zampini
2024-04-02 12:04:25 +02:00
parent 503a2e2529 7bf2a79920
commit 040ae1609b
8 changed files with 298 additions and 34 deletions
Download Patch File Download Diff File Expand all files Collapse all files
fem/pgridfunc.cpp
+6 -2
View File
@@ -683,7 +683,7 @@ void ParGridFunction::ProjectBdrCoefficient(
gcomm.Reduce<int>(values_counter.HostReadWrite(), GroupCommunicator::Sum);
// Accumulate the values globally.
gcomm.Reduce<real_t>(values.HostReadWrite(), GroupCommunicator::Sum);
// Only the values in the master are guaranteed to be correct!
for (int i = 0; i < values.Size(); i++)
{
if (values_counter[i])
@@ -691,6 +691,8 @@ void ParGridFunction::ProjectBdrCoefficient(
(*this)(i) = values(i)/values_counter[i];
}
}
// Broadcast values to other processors to have a consistent GridFunction
gcomm.Bcast<real_t>((*this).HostReadWrite());
#ifdef MFEM_DEBUG
Array<int> ess_vdofs_marker;
@@ -737,7 +739,7 @@ void ParGridFunction::ProjectBdrCoefficientTangent(VectorCoefficient &vcoeff,
gcomm.Reduce<int>(values_counter.HostReadWrite(), GroupCommunicator::Sum);
// Accumulate the values globally.
gcomm.Reduce<real_t>(values.HostReadWrite(), GroupCommunicator::Sum);
// Only the values in the master are guaranteed to be correct!
for (int i = 0; i < values.Size(); i++)
{
if (values_counter[i])
@@ -745,6 +747,8 @@ void ParGridFunction::ProjectBdrCoefficientTangent(VectorCoefficient &vcoeff,
(*this)(i) = values(i)/values_counter[i];
}
}
// Broadcast values to other processors to have a consistent GridFunction
gcomm.Bcast<real_t>((*this).HostReadWrite());
#ifdef MFEM_DEBUG
Array<int> ess_vdofs_marker;
fem/pgridfunc.hpp
-2
View File
@@ -254,12 +254,10 @@ public:
using GridFunction::ProjectBdrCoefficient;
// Only the values in the master are guaranteed to be correct!
void ProjectBdrCoefficient(VectorCoefficient &vcoeff,
const Array<int> &attr) override
{ ProjectBdrCoefficient(NULL, &vcoeff, attr); }
// Only the values in the master are guaranteed to be correct!
void ProjectBdrCoefficient(Coefficient *coeff[],
const Array<int> &attr) override
{ ProjectBdrCoefficient(coeff, NULL, attr); }
linalg/hypre.cpp
+45 -22
View File
@@ -701,6 +701,47 @@ static void CopyCSR_J(const int nnz, const MemoryIJData &mem_csr,
}
#endif
// Method called after hypre_CSRMatrixReorder()
static void SyncBackCSR(SparseMatrix *csr, MemoryIJData &mem_csr)
{
const MemoryClass hypre_mc = GetHypreMemoryClass();
const bool data_shallow = CanShallowCopy(csr->GetMemoryData(), hypre_mc);
#if !defined(HYPRE_BIGINT) && defined(MFEM_DEBUG)
const bool J_shallow = CanShallowCopy(csr->GetMemoryJ(), hypre_mc);
MFEM_ASSERT(J_shallow == data_shallow, "unsupported state");
#endif
if (data_shallow)
{
// I is not modified
#ifndef HYPRE_BIGINT
csr->GetMemoryJ().Sync(mem_csr.J);
#else
// We use nnz = csr->GetMemoryJ().Capacity() which is the same as the
// value used in CopyConvertMemory() in CopyCSR().
CopyCSR_J(csr->GetMemoryJ().Capacity(), mem_csr, csr->GetMemoryJ());
#endif
csr->GetMemoryData().Sync(mem_csr.data);
}
}
// Method called after hypre_CSRMatrixReorder()
static void SyncBackBoolCSR(Table *bool_csr, MemoryIJData &mem_csr)
{
const MemoryClass hypre_mc = GetHypreMemoryClass();
const bool J_shallow = CanShallowCopy(bool_csr->GetJMemory(), hypre_mc);
if (J_shallow)
{
// I is not modified
#ifndef HYPRE_BIGINT
bool_csr->GetJMemory().Sync(mem_csr.J);
#else
// No need to sync the J array back to the Table
#endif
}
}
// static method
signed char HypreParMatrix::HypreCsrToMem(hypre_CSRMatrix *h_mat,
MemoryType h_mat_mt,
@@ -798,12 +839,7 @@ HypreParMatrix::HypreParMatrix(MPI_Comm comm, HYPRE_BigInt glob_size,
/* Make sure that the first entry in each row is the diagonal one. */
HypreReadWrite();
hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(A));
#ifdef HYPRE_BIGINT
if (CanShallowCopy(diag->GetMemoryData(), GetHypreMemoryClass()))
{
CopyCSR_J(A->diag->num_nonzeros, mem_diag, diag->GetMemoryJ());
}
#endif
SyncBackCSR(diag, mem_diag); // update diag, if needed
hypre_MatvecCommPkgCreate(A);
}
@@ -842,12 +878,7 @@ HypreParMatrix::HypreParMatrix(MPI_Comm comm,
{
HypreReadWrite();
hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(A));
#ifdef HYPRE_BIGINT
if (CanShallowCopy(diag->GetMemoryData(), GetHypreMemoryClass()))
{
CopyCSR_J(A->diag->num_nonzeros, mem_diag, diag->GetMemoryJ());
}
#endif
SyncBackCSR(diag, mem_diag); // update diag, if needed
}
hypre_MatvecCommPkgCreate(A);
@@ -897,12 +928,7 @@ HypreParMatrix::HypreParMatrix(MPI_Comm comm,
{
HypreReadWrite();
hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(A));
#ifdef HYPRE_BIGINT
if (CanShallowCopy(diag->GetMemoryData(), GetHypreMemoryClass()))
{
CopyCSR_J(A->diag->num_nonzeros, mem_diag, diag->GetMemoryJ());
}
#endif
SyncBackCSR(diag, mem_diag); // update diag, if needed
}
hypre_MatvecCommPkgCreate(A);
@@ -1060,10 +1086,7 @@ HypreParMatrix::HypreParMatrix(MPI_Comm comm,
{
HypreReadWrite();
hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(A));
#ifdef HYPRE_BIGINT
// No need to sync the J array back to the Table diag.
// CopyCSR_J(A->diag->num_nonzeros, mem_diag, diag->GetJMemory());
#endif
SyncBackBoolCSR(diag, mem_diag); // update diag, if needed
}
hypre_MatvecCommPkgCreate(A);
miniapps/meshing/CMakeLists.txt
-3
View File
@@ -143,9 +143,6 @@ if (MFEM_USE_MPI)
)
# Meshing miniapps that return MFEM_SKIP_RETURN_VALUE in some cases:
set(SKIP_TESTS)
if (HYPRE_USING_CUDA OR HYPRE_USING_HIP)
list(APPEND SKIP_TESTS pmesh-fitting)
endif()
foreach(test ${PARALLEL_TESTS})
if (test IN_LIST SKIP_TESTS)
miniapps/meshing/pmesh-fitting.cpp
-5
View File
@@ -48,11 +48,6 @@ using namespace std;
int main (int argc, char *argv[])
{
#ifdef HYPRE_USING_GPU
cout << "\nThis miniapp is NOT supported with the GPU version of hypre.\n\n";
return MFEM_SKIP_RETURN_VALUE;
#endif
// 0. Initialize MPI and HYPRE.
Mpi::Init(argc, argv);
int myid = Mpi::WorldRank();
tests/unit/CMakeLists.txt
+1
View File
@@ -102,6 +102,7 @@ set(UNIT_TESTS_SRCS
fem/test_pa_kernels.cpp
fem/test_pgridfunc_save_serial.cpp
fem/test_project_bdr.cpp
fem/test_project_bdr_par.cpp
fem/test_quadf_coef.cpp
fem/test_quadinterpolator.cpp
fem/test_quadraturefunc.cpp
tests/unit/fem/test_project_bdr_par.cpp
+107
View File
@@ -0,0 +1,107 @@
// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
// LICENSE and NOTICE for details. LLNL-CODE-806117.
//
// This file is part of the MFEM library. For more information and source code
// availability visit https://mfem.org.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the BSD-3 license. We welcome feedback and contributions, see file
// CONTRIBUTING.md for details.
#include "mfem.hpp"
#include "unit_tests.hpp"
using namespace mfem;
#ifdef MFEM_USE_MPI
TEST_CASE("ProjectBdrCoefficient", "[Parallel]")
{
int num_procs = Mpi::WorldSize();
Mesh serial_mesh = Mesh::MakeCartesian3D(num_procs, num_procs, 1,
Element::HEXAHEDRON,
1.0, 1.0, 0.1, false);
// Assign alternating element attributes to each element to create a
// checkerboard pattern
for (int i = 0; i < serial_mesh.GetNE(); i++)
{
int attr = (i + (1 + num_procs % 2) * (i / num_procs)) % 2 + 1;
serial_mesh.SetAttribute(i, attr);
}
int bdr_max = serial_mesh.bdr_attributes.Max();
// Label all interior faces as boundary elements
Array<int> v(4);
for (int i=0; i < serial_mesh.GetNumFaces(); i++)
{
if (serial_mesh.FaceIsInterior(i))
{
serial_mesh.GetFaceVertices(i, v);
serial_mesh.AddBdrQuad(v, bdr_max + i + 1);
}
}
serial_mesh.FinalizeMesh();
serial_mesh.SetAttributes();
// Create an intentionally bad partitioning
Array<int> partitioning(num_procs * num_procs);
for (int i = 0; i < num_procs * num_procs; i++)
{
// The following creates a shifting pattern where neighboring elements
// are never owned by the same processor
partitioning[i] = (2 * num_procs - 1 - (i % num_procs) -
i / num_procs) % num_procs;
}
ParMesh par_mesh(MPI_COMM_WORLD, serial_mesh, partitioning);
H1_FECollection h1fec(2, par_mesh.Dimension());
ParFiniteElementSpace h1fes(&par_mesh, &h1fec);
ParGridFunction gf(&h1fes);
gf = 0.0;
int par_bdr_max = par_mesh.bdr_attributes.Max();
Array<int> all_bdr(par_bdr_max);
all_bdr = 1;
ConstantCoefficient coeff(123.456);
gf.ProjectBdrCoefficient(coeff, all_bdr);
// We projected a value to all interior and exterior boundary elements.
// An interior boundary element is only owned by one of the two sharing processors
// and we expect the GridFunction on each of the two processors to be the same value
// on that face. We test this by checking that each element sum is the same.
real_t local_sum_expected = 0.0;
real_t local_sum = 0.0;
for (int e = 0; e < par_mesh.GetNE(); e++)
{
Vector dof_vals;
gf.GetElementDofValues(e, dof_vals);
real_t e_sum = dof_vals.Sum();
if (e == 0) { local_sum_expected = e_sum * par_mesh.GetGlobalNE(); }
local_sum += e_sum;
}
real_t global_sum = 0.0, global_sum_expected = 0.0;
#ifdef MFEM_USE_DOUBLE
MPI_Allreduce(&local_sum, &global_sum, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
MPI_Allreduce(&local_sum_expected, &global_sum_expected, 1, MPI_DOUBLE, MPI_MAX,
MPI_COMM_WORLD);
#else
MPI_Allreduce(&local_sum, &global_sum, 1, MPI_FLOAT, MPI_SUM, MPI_COMM_WORLD);
MPI_Allreduce(&local_sum_expected, &global_sum_expected, 1, MPI_FLOAT, MPI_MAX,
MPI_COMM_WORLD);
#endif
REQUIRE(global_sum == global_sum_expected);
}
#endif
tests/unit/linalg/test_matrix_hypre.cpp
+139
View File
@@ -18,6 +18,145 @@ namespace mfem
#ifdef MFEM_USE_MPI
TEST_CASE("HypreParMatrixWrapConstructors-SyncChecks", "[Parallel], [CUDA]")
{
const int dim = 2;
const int n1d = 6;
const int p = 2;
Mesh smesh = Mesh::MakeCartesian2D(n1d, n1d, Element::QUADRILATERAL);
ParMesh mesh(MPI_COMM_WORLD, smesh);
smesh.Clear();
SECTION("SquareBlockDiagWrapConstructor")
{
H1_FECollection fec(p, dim);
ParFiniteElementSpace fespace(&mesh, &fec);
ParBilinearForm a(&fespace);
a.AddDomainIntegrator(new MassIntegrator);
a.Assemble();
a.Finalize();
SparseMatrix &spmat = a.SpMat();
const int height = spmat.Height();
const int nnz = spmat.NumNonZeroElems();
// Create a square block diagonal HypreParMatrix with blocks corresponding
// to the local sparse matrices, spmat. The constructed HypreParMatrix
// reuses the I, J and data arrays of spmat (with some exceptions).
// The constructor will also permute the entries of its J and data arrays
// to ensure that the diagonal entry is first in every row.
HypreParMatrix hpmat(mesh.GetComm(),
fespace.GlobalVSize(),
fespace.GetDofOffsets(),
&spmat);
// Verify that spmat's arrays are not out of sync:
REQUIRE(spmat.GetMemoryI().CompareHostAndDevice(height+1) == 0);
REQUIRE(spmat.GetMemoryJ().CompareHostAndDevice(nnz) == 0);
REQUIRE(spmat.GetMemoryData().CompareHostAndDevice(nnz) == 0);
}
SECTION("RectangularBlockDiagWrapConstructor")
{
H1_FECollection fec(p, dim);
ParFiniteElementSpace fespace(&mesh, &fec);
ParBilinearForm a(&fespace);
a.AddDomainIntegrator(new MassIntegrator);
a.Assemble();
a.Finalize();
SparseMatrix &spmat = a.SpMat();
const int height = spmat.Height();
const int nnz = spmat.NumNonZeroElems();
// Create a rectangular block diagonal HypreParMatrix with blocks
// corresponding to the local sparse matrices, spmat. The constructed
// HypreParMatrix reuses the I, J and data arrays of spmat (with some
// exceptions).
// When the row and column offsets are the same pointer, the constructor
// will also permute the entries of its J and data arrays to ensure that
// the diagonal entry is first in every row.
HypreParMatrix hpmat(mesh.GetComm(),
fespace.GlobalVSize(), // num rows
fespace.GlobalVSize(), // num cols
fespace.GetDofOffsets(), // row offsets
fespace.GetDofOffsets(), // col offsets
&spmat);
// Verify that spmat's arrays are not out of sync:
REQUIRE(spmat.GetMemoryI().CompareHostAndDevice(height+1) == 0);
REQUIRE(spmat.GetMemoryJ().CompareHostAndDevice(nnz) == 0);
REQUIRE(spmat.GetMemoryData().CompareHostAndDevice(nnz) == 0);
}
SECTION("RectangularDiagOffdWrapConstructor")
{
H1_FECollection fec(p, dim);
ParFiniteElementSpace fespace(&mesh, &fec);
ParBilinearForm a(&fespace);
a.AddDomainIntegrator(new MassIntegrator);
a.Assemble();
a.Finalize();
SparseMatrix &diag = a.SpMat();
const int height = diag.Height();
const int nnz = diag.NumNonZeroElems();
SparseMatrix offd(height, 0, 0); // height x 0 matrix
HYPRE_BigInt cmap = 0;
// Create a rectangular HypreParMatrix with diagonal blocks corresponding
// to the local sparse matrices, diag, and zero off-diagonal block, offd.
// The constructed HypreParMatrix reuses the I, J and data arrays of diag
// and offd (with some exceptions).
// When the row and column offsets are the same pointer, the constructor
// will also permute the entries of its block diagonal's J and data arrays
// to ensure that the diagonal entry is first in every row.
HypreParMatrix hpmat(mesh.GetComm(),
fespace.GlobalVSize(), // num rows
fespace.GlobalVSize(), // num cols
fespace.GetDofOffsets(), // row offsets
fespace.GetDofOffsets(), // col offsets
&diag,
&offd,
&cmap,
false);
// Verify that diag's arrays are not out of sync:
REQUIRE(diag.GetMemoryI().CompareHostAndDevice(height+1) == 0);
REQUIRE(diag.GetMemoryJ().CompareHostAndDevice(nnz) == 0);
REQUIRE(diag.GetMemoryData().CompareHostAndDevice(nnz) == 0);
}
SECTION("BooleanRectangularBlockDiagWrapConstructor")
{
H1_FECollection fec(p, dim);
ParFiniteElementSpace fespace(&mesh, &fec);
const Table &el_dof = fespace.GetElementToDofTable();
Table el_dof_t;
Transpose(el_dof, el_dof_t, fespace.GetNDofs());
Table dof_dof;
Mult(el_dof_t, el_dof, dof_dof);
const int height = dof_dof.Size();
const int nnz = dof_dof.Size_of_connections();
// Create a Boolean rectangular block diagonal HypreParMatrix with blocks
// corresponding to the local Table dof_dof. The constructed
// HypreParMatrix reuses the I and J arrays of dof_dof (with some
// exceptions).
// When the row and column offsets are the same pointer, the constructor
// will also permute the entries of its J and data arrays to ensure that
// the diagonal entry is first in every row.
HypreParMatrix hpm(mesh.GetComm(),
fespace.GlobalVSize(), // num rows
fespace.GlobalVSize(), // num cols
fespace.GetDofOffsets(), // row offsets
fespace.GetDofOffsets(), // col offsets
&dof_dof);
// Verify that dof_dof's arrays are not out of sync:
REQUIRE(dof_dof.GetIMemory().CompareHostAndDevice(height+1) == 0);
REQUIRE(dof_dof.GetJMemory().CompareHostAndDevice(nnz) == 0);
}
}
TEST_CASE("HypreParMatrixAbsMult", "[Parallel], [HypreParMatrixAbsMult]")
{
int rank;