From 7b3701254b5b15194869a5ee2159399003d91cbf Mon Sep 17 00:00:00 2001 From: li97 Date: Thu, 5 Jun 2025 17:19:59 -0700 Subject: [PATCH 01/11] serial dg works with new structure to support parallelization --- src/serac/numerics/functional/functional.hpp | 117 +++++++++++++----- .../functional/tests/functional_basic_dg.cpp | 4 +- 2 files changed, 86 insertions(+), 35 deletions(-) diff --git a/src/serac/numerics/functional/functional.hpp b/src/serac/numerics/functional/functional.hpp index f6d03cb5ae..4db55dca3e 100644 --- a/src/serac/numerics/functional/functional.hpp +++ b/src/serac/numerics/functional/functional.hpp @@ -247,23 +247,40 @@ class Functional { { SERAC_MARK_FUNCTION; - for (uint32_t i = 0; i < num_trial_spaces; i++) { - P_trial_[i] = trial_space_[i]->GetProlongationMatrix(); - - input_L_[i].SetSize(P_trial_[i]->Height(), mfem::Device::GetMemoryType()); - - // create the necessary number of empty mfem::Vectors, to be resized later - input_E_.push_back({}); - input_E_buffer_.push_back({}); - } - test_function_space_ = {test::family, test::order, test::components}; std::array trial_families = {trials::family...}; std::array trial_orders = {trials::order...}; std::array trial_components = {trials::components...}; + for (uint32_t i = 0; i < num_trial_spaces; i++) { trial_function_spaces_[i] = {trial_families[i], trial_orders[i], trial_components[i]}; + + P_trial_[i] = trial_space_[i]->GetProlongationMatrix(); + + // For L2 spaces, configure a ParGridFunction with external data to exchange face neighbor data + // This is for DG method where interior faces on the processor boundary need to access data on the neighrbor element + if (trial_function_spaces_[i].family == Family::L2) { + // Set the vector size to store tdof data + tdof_values_[i].SetSize(P_trial_[i]->Height(), mfem::Device::GetMemoryType()); + + // Link the ParGridFunction to external FE space and data by reference + trial_pargrid_functions_[i].MakeRef(const_cast(trial_space_[i]), tdof_values_[i].GetData()); + + // Exchange face nbr data to set the vector in ParGridFunction with the right size + trial_pargrid_functions_[i].ExchangeFaceNbrData(); + + // Set the local input vector size to store local + ghost data + input_L_[i].SetSize(tdof_values_[i].Size() + trial_pargrid_functions_[i].FaceNbrData().Size(), mfem::Device::GetMemoryType()); + } + else { + input_L_[i].SetSize(P_trial_[i]->Height(), mfem::Device::GetMemoryType()); + } + // std::cout << "local size = " << input_L_[i].Size() << std::endl; + + // create the necessary number of empty mfem::Vectors, to be resized later + input_E_.push_back({}); + input_E_buffer_.push_back({}); } // for (auto type : {Domain::Type::Elements, Domain::Type::BoundaryElements, Domain::Type::InteriorFaces}) { @@ -466,36 +483,64 @@ class Functional { // get the values for each local processor for (uint32_t i = 0; i < num_trial_spaces; i++) { -#if 0 +// #if 0 +// if (trial_function_spaces_[i].family == Family::L2) { +// +// // make a PGF on the fly +// // TODO: don't allocate/deallocate this data every invocation +// mfem::ParGridFunction X; +// X.MakeRef(trial_space_[i], input_L_[i].GetData()); +// +// // call exchange face nbr +// X.ExchangeFaceNbrData(); +// +// // copy input_L[i] and facenbrdata [i] into a common array like: +// // first part second part +// // [ --- L --- | --- FND --- ] +// mfem::Vector first_part; +// first_part.MakeRef(input_L_[i], 0); +// P_trial_[i]->Mult(*input_T[i], first_part); +// +// mfem::Vector second_part; +// second_part.MakeRef(input_L_[i], first_part.Size()); +// second_part = X.FaceNbrData(); +// +// } else { +// +// P_trial_[i]->Mult(*input_T[i], input_L_[i]); +// +// } +// #else +// P_trial_[i]->Mult(*input_T[i], input_L_[i]); +// #endif if (trial_function_spaces_[i].family == Family::L2) { - - // make a PGF on the fly - // TODO: don't allocate/deallocate this data every invocation - mfem::ParGridFunction X; - X.MakeRef(trial_space_[i], input_L_[i].GetData()); - - // call exchange face nbr - X.ExchangeFaceNbrData(); - // copy input_L[i] and facenbrdata [i] into a common array like: // first part second part // [ --- L --- | --- FND --- ] - mfem::Vector first_part; - first_part.MakeRef(input_L_[i], 0); - P_trial_[i]->Mult(*input_T[i], first_part); - - mfem::Vector second_part; - second_part.MakeRef(input_L_[i], first_part.Size()); - second_part = X.FaceNbrData(); - + P_trial_[i]->Mult(*input_T[i], tdof_values_[i]); + // MakeRef from mfem::Vector will actually delete the data in the original vector, so we assign the + // values directly by SetVector, which invokes additional operations. This can be optimized later. + input_L_[i].SetVector(tdof_values_[i], 0); + // // Debug + // std::cout << "TDOFs" << std::endl; + // for (int j = 0; j < tdof_values_[i].Size(); ++j) { + // std::cout << tdof_values_[i][j] << " "; + // } + // std::cout << std::endl; + + // Point ParGridFunction to the new data and call exchange face nbr data + trial_pargrid_functions_[i].MakeRef(const_cast(trial_space_[i]), tdof_values_[i].GetData()); + trial_pargrid_functions_[i].ExchangeFaceNbrData(); + input_L_[i].SetVector(trial_pargrid_functions_[i].FaceNbrData(), tdof_values_[i].Size()); } else { - P_trial_[i]->Mult(*input_T[i], input_L_[i]); - } -#else - P_trial_[i]->Mult(*input_T[i], input_L_[i]); -#endif + // // Debug + // std::cout << "Local vector" << std::endl; + // for (int j = 0; j < input_L_[i].Size(); ++j) { + // std::cout << input_L_[i][j] << " "; + // } + // std::cout << std::endl; } output_L_ = 0.0; @@ -823,6 +868,12 @@ class Functional { std::array trial_function_spaces_; FunctionSpace test_function_space_; + /// @brief For access of neighbor element info for processor boundary faces + std::array trial_pargrid_functions_; + + /// @brief Store values on true dofs of the local processor + std::array tdof_values_; + /** * @brief Operator that converts true (global) DOF values to local (current rank) DOF values * for the test space diff --git a/src/serac/numerics/functional/tests/functional_basic_dg.cpp b/src/serac/numerics/functional/tests/functional_basic_dg.cpp index 6705c8713d..927f9bcd04 100644 --- a/src/serac/numerics/functional/tests/functional_basic_dg.cpp +++ b/src/serac/numerics/functional/tests/functional_basic_dg.cpp @@ -69,8 +69,8 @@ void L2_test(std::string meshfile) double t = 0.0; - auto value = residual(t, U); - // check_gradient(residual, t, U); + // auto value = residual(t, U); + check_gradient(residual, t, U); } TEST(basic, L2_test_tris_and_quads_linear) { L2_test<2, 1>(SERAC_REPO_DIR "/data/meshes/patch2D_tris_and_quads.mesh"); } From b0e4bc13f7972260b150046e906ec56cce78bc70 Mon Sep 17 00:00:00 2001 From: li97 Date: Mon, 23 Jun 2025 15:09:00 -0700 Subject: [PATCH 02/11] Configured DG ghost data and indexing. Added test --- .../functional/element_restriction.cpp | 101 +++++---- src/serac/numerics/functional/functional.hpp | 155 ++++++++------ .../numerics/functional/geometric_factors.cpp | 8 +- .../numerics/functional/tests/CMakeLists.txt | 1 + .../functional/tests/dg_ghost_face_index.cpp | 101 +++++++++ .../tests/dg_restriction_operators.cpp | 34 ++-- .../functional/tests/domain_tests.cpp | 191 ++++++++++-------- .../tests/element_restriction_tests.cpp | 76 +++---- .../tests/geometric_factors_tests.cpp | 15 +- src/serac/numerics/functional/typedefs.hpp | 5 +- .../physics/state/finite_element_state.cpp | 4 +- 11 files changed, 438 insertions(+), 253 deletions(-) create mode 100644 src/serac/numerics/functional/tests/dg_ghost_face_index.cpp diff --git a/src/serac/numerics/functional/element_restriction.cpp b/src/serac/numerics/functional/element_restriction.cpp index 998a4d9e9b..3049b690d5 100644 --- a/src/serac/numerics/functional/element_restriction.cpp +++ b/src/serac/numerics/functional/element_restriction.cpp @@ -474,15 +474,6 @@ axom::Array GetFaceDofs(const serac::f std::vector > local_face_dofs = geom_local_face_dofs(p); std::vector > lex_perm = lexicographic_permutations(p); - // sam: this function contains several comments that relate to an investigation into adopting - // mfem's FaceNbrData pattern (which ended up being too invasive to implement initially). - // I leave some of the commented code here so Julian's recommendations are not lost. - // - // int components_per_node = fes->GetVDim(); - // bool by_vdim = fes->GetOrdering() == mfem::Ordering::byVDIM; - // fes->ExchangeFaceNbrData(); - // int LSize = fes->GetProlongationMatrix()->Height(); - uint64_t n = 0; for (int f : mfem_face_ids) { @@ -499,24 +490,6 @@ axom::Array GetFaceDofs(const serac::f mfem::Array elem_ids; face_to_elem->GetRow(f, elem_ids); -#if 0 - mpi::out << f << " elem ids: "; - for (auto elem : elem_ids) { - mpi::out << elem << " "; - } - mpi::out << std::endl; - - mfem::Array test_dofs; - fes->GetFaceDofs(f, test_dofs); - - mpi::out << " dofs (sz = " << test_dofs.Size() << "): "; - for (int z = 0; z < test_dofs.Size(); z++) { - mpi::out << test_dofs[z] << " "; - } - mpi::out << std::endl; - mpi::out << std::endl; -#endif - for (auto elem : elem_ids) { // 2a. get the list of faces (and their orientations) that belong to that element ... mfem::Array elem_side_ids, orientations; @@ -537,6 +510,7 @@ axom::Array GetFaceDofs(const serac::f } // 2b. ... and find `i` such that `elem_side_ids[i] == f` + // This generates a local_face_index that coincides with ones from GetFaceInformation int i; for (i = 0; i < elem_side_ids.Size(); i++) { if (elem_side_ids[i] == f) break; @@ -545,6 +519,11 @@ axom::Array GetFaceDofs(const serac::f // 3. get the dofs for the entire element mfem::Array elem_dof_ids; fes->GetElementDofs(elem, elem_dof_ids); + // // Debug hanyu + // mpi::out << " Shared face = " << info.IsShared() << " Face id = " << f << " Elem id = " << elem << std::endl + // << i << " " << info.element[0].local_face_id << " " << info.element[1].local_face_id << std::endl + // << orientations[i] << " " << info.element[0].orientation << " " << info.element[1].orientation << + // std::endl; mfem::Geometry::Type elem_geom = mesh->GetElementGeometry(elem); @@ -557,38 +536,75 @@ axom::Array GetFaceDofs(const serac::f int orientation = (mesh->Dimension() == 2 && info.IsBoundary()) ? 0 : orientations[i]; // 4. extract only the dofs that correspond to side `i` + // mpi::out << "face local dofs: "; for (auto k : face_perm(orientation)) { face_dofs.push_back(uint64_t(elem_dof_ids[local_face_dofs[uint32_t(elem_geom)](i, k)])); + // mpi::out << face_dofs.back().index() << " "; } + // mpi::out << std::endl; -#if 0 - // (5. optional) add remaining dofs that were omitted on shared faces + // 5. add remaining dofs that were omitted on shared faces + // FaceNbrData includes all dofs in the "volumetric" element adjacent to the shared faces + // We need to again find our face on the ghost element and extract only the face dofs + // This may not yet work for NC faces if (info.IsShared()) { + // FaceNbrData should be already exchanged in Functional when invoked by ParGridFunction + // which also invokes ExchangeFaceNbrData on subsequent ParFiniteElementSpace and ParMesh + int components_per_node = fes->GetVDim(); + bool by_vdim = fes->GetOrdering() == mfem::Ordering::byVDIM; + int LSize = fes->GetProlongationMatrix()->Height(); + + // Additional safeguard to make sure ByNODES ordering in vector L2 field raise an error + SLIC_ERROR_IF(!by_vdim && components_per_node > 1, "Unsupported: L2 vector field ordered by nodes"); mfem::Array shared_elem_vdof_ids; + mfem::Array shared_elem_dof_ids; // on my processor - // VVVVVVV + // VVVVVVV // dofs for the face [0 1 3 4 | 5 8 12 13] // ^^^^^^^^^ // on the other processor int other_element_id = info.element[1].index; - fes->GetFaceNbrElementVDofs(other_element_id, shared_elem_vdof_ids); // indices into vector from FaceNbrData + fes->GetFaceNbrElementVDofs(other_element_id, shared_elem_vdof_ids); // indices into vector from FaceNbrData - mpi::out << " this face is also shared so it has additional dofs: "; - int dofs_per_face = shared_elem_vdof_ids.Size() / components_per_node; - int stride = (by_vdim) ? components_per_node : 1; - - // sam: is this right? // byVDIM == x y z x y z x y z x y z // byNODES == x x x x y y y y z z z z - for (int k = 0; k < dofs_per_face; k++) { - face_dofs.push_back(uint64_t(fes->VDofToDof(shared_elem_vdof_ids[k * stride]) + LSize)); - mpi::out << face_dofs.back().index() << " "; + int dofs_per_element = shared_elem_vdof_ids.Size() / components_per_node; + int stride = (by_vdim) ? components_per_node : 1; + + shared_elem_dof_ids.SetSize(dofs_per_element); + for (int k = 0; k < dofs_per_element; ++k) { + shared_elem_dof_ids[k] = fes->VDofToDof(shared_elem_vdof_ids[k * stride]); // Change vdof to dof indices + } + + // Find the dofs on the shared face on the ghost element side + // Neighbor element local_face_index and orientation is already available from GetFaceInformation + // mpi::out << " this face is also shared so it has additional dofs: "; + mfem::Geometry::Type ghost_geom = fes->GetParMesh()->face_nbr_elements[other_element_id]->GetGeometryType(); + int j = info.element[1].local_face_id; + int ghost_orientation; + if (mesh->Dimension() == 2) { + // In 2D, orientations[i] sometimes is inverted as compared to the ones from GetFaceInformation + // In this case, we stay with the orientations constructed previously by GetElementEdges, which is strictly + // CCW So we set ghost orientation to be exactly the opposite of the original orientation (orientations[i]) + // side note: even if you use info.element[1].orientation, looks like it's still fine. + // The only difference is the order those face dof indices gets registered in std::vector face_dofs + ghost_orientation = (orientation == 0) ? 1 : 0; + } else { + // In 3D, I think it's consistently CCW. So we directly use the orientation from GetFaceInformation + // This orientation is already consistent with the permutation in {0, 1, ... , n} form + ghost_orientation = info.element[1].orientation; + } + + // extract only the dofs that correspond to side `j` + for (auto k : face_perm(ghost_orientation)) { + face_dofs.push_back(uint64_t(shared_elem_dof_ids[local_face_dofs[uint32_t(ghost_geom)](j, k)] + + LSize / components_per_node)); + // mpi::out << face_dofs.back().index() << " "; } - mpi::out << std::endl; + // mpi::out << std::endl; } -#endif } // H1 and Hcurl spaces are more straight-forward, since @@ -682,6 +698,9 @@ void ElementRestriction::Gather(const mfem::Vector& L_vector, mfem::Vector& E_ve for (uint64_t j = 0; j < nodes_per_elem; j++) { uint64_t E_id = (i * components + c) * nodes_per_elem + j; uint64_t L_id = GetVDof(dof_info(i, j), c).index(); + // // debug + // mpi::out << "elem #" << i << " " << "E_id = " << E_id << " " << dof_info(i, j).index() + // << " " << "L_id = " << L_id << std::endl; E_vector[int(E_id)] = L_vector[int(L_id)]; } } diff --git a/src/serac/numerics/functional/functional.hpp b/src/serac/numerics/functional/functional.hpp index 4db55dca3e..b038462640 100644 --- a/src/serac/numerics/functional/functional.hpp +++ b/src/serac/numerics/functional/functional.hpp @@ -166,6 +166,24 @@ generateParFiniteElementSpace(mfem::ParMesh* mesh) return std::pair(std::move(fes), std::move(fec)); } +/** + * @brief helper function to locally cast away const on FE space to update face neighbor data. + * this is ok because the original FE space is declared without const and we constrained + * the non-constness locally + */ +inline void updateFaceNbrData(const mfem::ParFiniteElementSpace* const_trial_space, mfem::ParGridFunction& trial_pgf, + mfem::Vector& trial_tdof_vals) +{ + mfem::ParFiniteElementSpace* nonconst_trial_space = const_cast(const_trial_space); + + // Link the ParGridFunction to external FE space and data by reference + trial_pgf.MakeRef(nonconst_trial_space, trial_tdof_vals.GetData()); + + // Exchange face nbr data to set the vector in ParGridFunction with the right size + // This should automatically invoke ExchangeFaceNbrData on the subsequent ParFESpace and ParMesh + trial_pgf.ExchangeFaceNbrData(); +} + /// @cond template class Functional; @@ -259,24 +277,20 @@ class Functional { P_trial_[i] = trial_space_[i]->GetProlongationMatrix(); // For L2 spaces, configure a ParGridFunction with external data to exchange face neighbor data - // This is for DG method where interior faces on the processor boundary need to access data on the neighrbor element + // This is for DG method where interior faces on the processor boundary need to access data on the neighrbor + // element if (trial_function_spaces_[i].family == Family::L2) { // Set the vector size to store tdof data - tdof_values_[i].SetSize(P_trial_[i]->Height(), mfem::Device::GetMemoryType()); - - // Link the ParGridFunction to external FE space and data by reference - trial_pargrid_functions_[i].MakeRef(const_cast(trial_space_[i]), tdof_values_[i].GetData()); + input_ldof_values_[i].SetSize(P_trial_[i]->Height(), mem_type); - // Exchange face nbr data to set the vector in ParGridFunction with the right size - trial_pargrid_functions_[i].ExchangeFaceNbrData(); + // Initialize face neighbor data vector + updateFaceNbrData(trial_space_[i], trial_pargrid_functions_[i], input_ldof_values_[i]); // Set the local input vector size to store local + ghost data - input_L_[i].SetSize(tdof_values_[i].Size() + trial_pargrid_functions_[i].FaceNbrData().Size(), mfem::Device::GetMemoryType()); - } - else { - input_L_[i].SetSize(P_trial_[i]->Height(), mfem::Device::GetMemoryType()); + input_L_[i].SetSize(input_ldof_values_[i].Size() + trial_pargrid_functions_[i].FaceNbrData().Size(), mem_type); + } else { + input_L_[i].SetSize(P_trial_[i]->Height(), mem_type); } - // std::cout << "local size = " << input_L_[i].Size() << std::endl; // create the necessary number of empty mfem::Vectors, to be resized later input_E_.push_back({}); @@ -374,10 +388,10 @@ class Functional { std::vector arg_vec = {args...}; for (uint32_t i : arg_vec) { domain.insert_restriction(trial_space_[i], trial_function_spaces_[i]); - check_interior_face_compatibility(domain.mesh_, trial_function_spaces_[i]); + // check_interior_face_compatibility(domain.mesh_, trial_function_spaces_[i]); } domain.insert_restriction(test_space_, test_function_space_); - check_interior_face_compatibility(domain.mesh_, test_function_space_); + // check_interior_face_compatibility(domain.mesh_, test_function_space_); using signature = test(decltype(serac::type(trial_spaces))...); integrals_.push_back( @@ -483,59 +497,75 @@ class Functional { // get the values for each local processor for (uint32_t i = 0; i < num_trial_spaces; i++) { -// #if 0 -// if (trial_function_spaces_[i].family == Family::L2) { -// -// // make a PGF on the fly -// // TODO: don't allocate/deallocate this data every invocation -// mfem::ParGridFunction X; -// X.MakeRef(trial_space_[i], input_L_[i].GetData()); -// -// // call exchange face nbr -// X.ExchangeFaceNbrData(); -// -// // copy input_L[i] and facenbrdata [i] into a common array like: -// // first part second part -// // [ --- L --- | --- FND --- ] -// mfem::Vector first_part; -// first_part.MakeRef(input_L_[i], 0); -// P_trial_[i]->Mult(*input_T[i], first_part); -// -// mfem::Vector second_part; -// second_part.MakeRef(input_L_[i], first_part.Size()); -// second_part = X.FaceNbrData(); -// -// } else { -// -// P_trial_[i]->Mult(*input_T[i], input_L_[i]); -// -// } -// #else -// P_trial_[i]->Mult(*input_T[i], input_L_[i]); -// #endif if (trial_function_spaces_[i].family == Family::L2) { // copy input_L[i] and facenbrdata [i] into a common array like: // first part second part // [ --- L --- | --- FND --- ] - P_trial_[i]->Mult(*input_T[i], tdof_values_[i]); + P_trial_[i]->Mult(*input_T[i], input_ldof_values_[i]); + // MakeRef from mfem::Vector will actually delete the data in the original vector, so we assign the // values directly by SetVector, which invokes additional operations. This can be optimized later. - input_L_[i].SetVector(tdof_values_[i], 0); - // // Debug - // std::cout << "TDOFs" << std::endl; - // for (int j = 0; j < tdof_values_[i].Size(); ++j) { - // std::cout << tdof_values_[i][j] << " "; - // } - // std::cout << std::endl; - - // Point ParGridFunction to the new data and call exchange face nbr data - trial_pargrid_functions_[i].MakeRef(const_cast(trial_space_[i]), tdof_values_[i].GetData()); - trial_pargrid_functions_[i].ExchangeFaceNbrData(); - input_L_[i].SetVector(trial_pargrid_functions_[i].FaceNbrData(), tdof_values_[i].Size()); + input_L_[i].SetVector(input_ldof_values_[i], 0); + + // Update again with new data + updateFaceNbrData(trial_space_[i], trial_pargrid_functions_[i], input_ldof_values_[i]); + + // Set the values of neighbor ghost dofs + // The entries in face_nbr_data vector is ALWAYS arranged byNODES per "volumetric" element. + // For example with two quadrilateral ghost elements we would have such face_nbr_data + // ----------- Ghost elem 1 ---------- ----------- Ghost elem 2 ---------- + // {X1 X1 X1 X1 Y1 Y1 Y1 Y1 Z1 Z1 Z1 Z1 X2 X2 X2 X2 Y2 Y2 Y2 Y2 Z2 Z2 Z2 Z2} + // + // This is because 1. mfem prepares face neighbor data element by element to communicate later (standard) + // 2. mfem uses GetElementVDofs to gather the data for communication which returns a set of indices ALWAYS + // ordered byNODES (???). If the ordering we set for ParFiniteElementSpace is byVDIM, we need the entries to be + // ----------- Ghost elem 1 ---------- ----------- Ghost elem 2 ---------- + // {X1 Y1 Z1 X1 Y1 Z1 X1 Y1 Z1 X1 Y1 Z1 X2 Y2 Z2 X2 Y2 Z2 X2 Y2 Z2 X2 Y2 Z2} + // so it's consistent with the local data. Therefore, we need to manually change this ordering + if (trial_space_[i]->GetVDim() == 1) { + // For scalar fields, this weird ordering doesn't cause any problems + input_L_[i].SetVector(trial_pargrid_functions_[i].FaceNbrData(), input_ldof_values_[i].Size()); + } else { + if (trial_space_[i]->GetOrdering() == mfem::Ordering::byVDIM) { + // For vector fields with byVDIM ordering, we manually set the entries in VDIM order. By the way this + // is really annoying because we have to get face neighbor vdof indices again in element restriction + int num_ghost_elem = trial_space_[i]->GetParMesh()->GetNFaceNeighborElements(); + int components_per_node = trial_space_[i]->GetVDim(); + const mfem::Vector& face_neighbor_data = trial_pargrid_functions_[i].FaceNbrData(); + + int offset = 0; + int local_size = input_ldof_values_[i].Size(); + + for (int n = 0; n < num_ghost_elem; ++n) { + mfem::Array old_shared_elem_vdof_ids; + trial_space_[i]->GetFaceNbrElementVDofs(n, old_shared_elem_vdof_ids); + int dofs_per_element = old_shared_elem_vdof_ids.Size() / components_per_node; + + // Put the entries into their new index + for (int m = 0; m < old_shared_elem_vdof_ids.Size(); ++m) { + // -------------- dofs before -------------- ---- component ----- + int new_id = (m % dofs_per_element) * components_per_node + m / dofs_per_element + offset; + input_L_[i][local_size + new_id] = face_neighbor_data(old_shared_elem_vdof_ids[m]); + } + // Increase offset for next ghost element + offset += old_shared_elem_vdof_ids.Size(); + } + } else { + // For vector fields with byNODES ordering, we need to prepare the local input vector in the form + // { true_comp_X ghost_comp_X true_comp_Y ghost_comp_Y true_comp_Z ghost_comp_Z } + // In this form, to ensure the correct mapping between dof and vdof, we need to change the ndofs in + // FiniteElementSpace to include ones from ghost element. Additionally, we need to prepare + // ghost_comp_X / ghost_comp_Y / ghost_comp_Z from element wise entries in face_nbr_data + // eg. ghost_comp_X includes X component entries of all ghost element dofs + // This requires significant renumbering of the dofs and therefore is not supported for now + SLIC_ERROR_ROOT("Unsupported: L2 vector field ordered by nodes"); + } + } } else { P_trial_[i]->Mult(*input_T[i], input_L_[i]); } // // Debug + // std::cout << "local size = " << input_L_[i].Size() << std::endl; // std::cout << "Local vector" << std::endl; // for (int j = 0; j < input_L_[i].Size(); ++j) { // std::cout << input_L_[i][j] << " "; @@ -555,6 +585,13 @@ class Functional { input_E_buffer_[i].SetSize(int(G_trial.ESize())); input_E_[i].Update(input_E_buffer_[i], G_trial.bOffsets()); G_trial.Gather(input_L_[i], input_E_[i]); + // // Debug + // std::cout << "element size = " << G_trial.ESize() << std::endl; + // std::cout << "Element vector" << std::endl; + // for (int j = 0; j < input_E_[i].Size(); ++j) { + // std::cout << input_E_[i][j] << " "; + // } + // std::cout << std::endl; } output_E_buffer_.SetSize(int(G_test.ESize())); @@ -871,8 +908,8 @@ class Functional { /// @brief For access of neighbor element info for processor boundary faces std::array trial_pargrid_functions_; - /// @brief Store values on true dofs of the local processor - std::array tdof_values_; + /// @brief Store values on L2 ldofs of the local processor + std::array input_ldof_values_; /** * @brief Operator that converts true (global) DOF values to local (current rank) DOF values diff --git a/src/serac/numerics/functional/geometric_factors.cpp b/src/serac/numerics/functional/geometric_factors.cpp index 972fef3d34..d41a7645b0 100644 --- a/src/serac/numerics/functional/geometric_factors.cpp +++ b/src/serac/numerics/functional/geometric_factors.cpp @@ -67,10 +67,10 @@ void compute_geometric_factors(mfem::Vector& positions_q, mfem::Vector& jacobian GeometricFactors::GeometricFactors(const Domain& domain, int q, mfem::Geometry::Type geom) { - // const mfem::ParGridFunction* nodes = static_cast< const mfem::ParGridFunction * >(domain.mesh_.GetNodes()); - // mfem::ParFiniteElementSpace * pfes = nodes->ParFESpace(); - const mfem::GridFunction* nodes = domain.mesh_.GetNodes(); - const mfem::FiniteElementSpace* fes = nodes->FESpace(); + const mfem::ParGridFunction* nodes = static_cast(domain.mesh_.GetNodes()); + mfem::ParFiniteElementSpace* fes = nodes->ParFESpace(); + // const mfem::GridFunction* nodes = domain.mesh_.GetNodes(); + // const mfem::FiniteElementSpace* fes = nodes->FESpace(); const std::vector& element_ids = domain.get_mfem_ids(geom); diff --git a/src/serac/numerics/functional/tests/CMakeLists.txt b/src/serac/numerics/functional/tests/CMakeLists.txt index a3f8b34eee..1357e65a1e 100644 --- a/src/serac/numerics/functional/tests/CMakeLists.txt +++ b/src/serac/numerics/functional/tests/CMakeLists.txt @@ -47,6 +47,7 @@ set(functional_parallel_test_sources functional_boundary_test.cpp functional_comparisons.cpp functional_comparison_L2.cpp + dg_ghost_face_index.cpp ) serac_add_tests(SOURCES ${functional_parallel_test_sources} diff --git a/src/serac/numerics/functional/tests/dg_ghost_face_index.cpp b/src/serac/numerics/functional/tests/dg_ghost_face_index.cpp new file mode 100644 index 0000000000..09ec22014a --- /dev/null +++ b/src/serac/numerics/functional/tests/dg_ghost_face_index.cpp @@ -0,0 +1,101 @@ +// Copyright (c) Lawrence Livermore National Security, LLC and +// other Serac Project Developers. See the top-level LICENSE file for +// details. +// +// SPDX-License-Identifier: (BSD-3-Clause) + +#include +#include + +#include "mfem.hpp" + +#include + +#include "axom/slic/core/SimpleLogger.hpp" +#include "serac/infrastructure/application_manager.hpp" +#include "serac/serac_config.hpp" +#include "serac/mesh_utils/mesh_utils_base.hpp" +#include "serac/numerics/stdfunction_operator.hpp" +#include "serac/numerics/functional/functional.hpp" +#include "serac/numerics/functional/tensor.hpp" + +using namespace serac; +using namespace serac::profiling; + +template +void L2_index_test(std::string meshfile) +{ + using test_space = L2; + using trial_space = L2; + + auto mesh = mesh::refineAndDistribute(buildMeshFromFile(meshfile), 0); + + auto fec = mfem::L2_FECollection(p, dim, mfem::BasisType::GaussLobatto); + mfem::ParFiniteElementSpace fespace(mesh.get(), &fec, dim, serac::ordering); + + mfem::ParGridFunction pgf(&fespace); + mfem::VectorFunctionCoefficient vcoef(dim, [](const mfem::Vector& X, mfem::Vector& F) { + int d = X.Size(); + F.SetSize(d); + for (int i = 0; i < d; ++i) { + F(i) = X(i); + } + }); + pgf.ProjectCoefficient(vcoef); + + mfem::Vector U(fespace.TrueVSize()); + U = pgf; + + // Construct the new functional object using the specified test and trial spaces + Functional residual(&fespace, {&fespace}); + + Domain interior_faces = InteriorFaces(*mesh); + + residual.AddInteriorFaceIntegral( + Dimension{}, DependsOn<0>{}, + [=](double /*t*/, auto X, auto velocity) { + // compute the surface normal + auto dX_dxi = get(X); + auto n = normalize(cross(dX_dxi)); + + // extract the velocity values from each side of the interface + // note: the orientation convention is such that the normal + // computed as above will point from from side 1->2 + auto [u_1, u_2] = velocity; + std::cout << "One side = " << u_1 << " The other side = " << u_2 << " Jump = " << u_1 - u_2 << std::endl; + + auto a = dot(u_2 - u_1, n); + + auto f_1 = u_1 * a; + auto f_2 = u_2 * a; + return serac::tuple{f_1, f_2}; + }, + interior_faces); + + double t = 0.0; + + auto value = residual(t, U); + EXPECT_NEAR(0., value.Norml2(), 1.e-12); +} +// TEST(index, L2_two_quad) { L2_index_test<2, 1>("./TwoQuad.mesh"); } +TEST(index, L2_test_tris_and_quads_linear) +{ + L2_index_test<2, 1>(SERAC_REPO_DIR "/data/meshes/patch2D_tris_and_quads.mesh"); +} +TEST(index, L2_test_tris_and_quads_quadratic) +{ + L2_index_test<2, 2>(SERAC_REPO_DIR "/data/meshes/patch2D_tris_and_quads.mesh"); +} + +TEST(index, L2_test_tets_linear) { L2_index_test<3, 1>(SERAC_REPO_DIR "/data/meshes/patch3D_tets.mesh"); } +TEST(index, L2_test_tets_quadratic) { L2_index_test<3, 2>(SERAC_REPO_DIR "/data/meshes/patch3D_tets.mesh"); } + +TEST(index, L2_test_hexes_linear) { L2_index_test<3, 1>(SERAC_REPO_DIR "/data/meshes/patch3D_hexes.mesh"); } +TEST(index, L2_test_hexes_quadratic) { L2_index_test<3, 2>(SERAC_REPO_DIR "/data/meshes/patch3D_hexes.mesh"); } + +int main(int argc, char* argv[]) +{ + ::testing::InitGoogleTest(&argc, argv); + serac::ApplicationManager applicationManager(argc, argv); + return RUN_ALL_TESTS(); +} diff --git a/src/serac/numerics/functional/tests/dg_restriction_operators.cpp b/src/serac/numerics/functional/tests/dg_restriction_operators.cpp index 49def5842c..518dd92cfc 100644 --- a/src/serac/numerics/functional/tests/dg_restriction_operators.cpp +++ b/src/serac/numerics/functional/tests/dg_restriction_operators.cpp @@ -7,6 +7,7 @@ #include #include "serac/mesh_utils/mesh_utils_base.hpp" +#include "serac/mesh_utils/mesh_utils.hpp" #include "serac/numerics/functional/domain.hpp" #include "serac/numerics/functional/element_restriction.hpp" #include "serac/numerics/functional/functional.hpp" @@ -286,13 +287,14 @@ void parametrized_test(int permutation) constexpr mfem::Geometry::Type face_geom = face_type(geom); constexpr int dim = dimension_of(geom); - mfem::Mesh mesh = generate_permuted_mesh(geom, permutation); + mfem::Mesh bmesh = generate_permuted_mesh(geom, permutation); + auto mesh = serac::mesh::refineAndDistribute(std::move(bmesh)); - Domain interior_faces = InteriorFaces(mesh); + Domain interior_faces = InteriorFaces(*mesh); // each one of these meshes should have two elements // and a single "face" that separates them - EXPECT_EQ(mesh.GetNE(), 2); + EXPECT_EQ(mesh->GetNE(), 2); std::vector face_ids; @@ -318,13 +320,13 @@ void parametrized_test(int permutation) auto Hcurl_fec = std::make_unique(polynomial_order, dim); auto L2_fec = std::make_unique(polynomial_order, dim, mfem::BasisType::GaussLobatto); - auto H1_fes = std::make_unique(&mesh, H1_fec.get()); - auto Hcurl_fes = std::make_unique(&mesh, Hcurl_fec.get()); - auto L2_fes = std::make_unique(&mesh, L2_fec.get()); + auto H1_fes = std::make_unique(mesh.get(), H1_fec.get()); + auto Hcurl_fes = std::make_unique(mesh.get(), Hcurl_fec.get()); + auto L2_fes = std::make_unique(mesh.get(), L2_fec.get()); - mfem::GridFunction H1_gf(H1_fes.get()); - mfem::GridFunction Hcurl_gf(Hcurl_fes.get()); - mfem::GridFunction L2_gf(L2_fes.get()); + mfem::ParGridFunction H1_gf(H1_fes.get()); + mfem::ParGridFunction Hcurl_gf(Hcurl_fes.get()); + mfem::ParGridFunction L2_gf(L2_fes.get()); mfem::FunctionCoefficient sfunc(scalar_func_mfem); mfem::VectorFunctionCoefficient vfunc(dim, vector_func_mfem); @@ -351,12 +353,12 @@ void parametrized_test(int permutation) using space = L2; - auto pmesh = mesh::refineAndDistribute(std::move(mesh), 0); - Domain pinterior_faces = InteriorFaces(*pmesh); + // auto pmesh = mesh::refineAndDistribute(std::move(mesh), 0); + // Domain pinterior_faces = InteriorFaces(*pmesh); - auto L2_pfes = mfem::ParFiniteElementSpace(pmesh.get(), L2_fec.get()); + // auto L2_pfes = mfem::ParFiniteElementSpace(pmesh.get(), L2_fec.get()); - Functional r({&L2_pfes}); + Functional r({L2_fes.get()}); r.AddInteriorFaceIntegral( Dimension{}, DependsOn<0>{}, @@ -369,10 +371,10 @@ void parametrized_test(int permutation) auto difference = rho1 - rho2; return difference * difference; }, - pinterior_faces); + interior_faces); - mfem::Vector U(L2_pfes.TrueVSize()); - mfem::ParGridFunction U_pgf(&L2_pfes); + mfem::Vector U(L2_fes->TrueVSize()); + mfem::ParGridFunction U_pgf(L2_fes.get()); U_pgf.ProjectCoefficient(sfunc); U_pgf.GetTrueDofs(U); diff --git a/src/serac/numerics/functional/tests/domain_tests.cpp b/src/serac/numerics/functional/tests/domain_tests.cpp index 35e9cb92eb..e8bf4dd82f 100644 --- a/src/serac/numerics/functional/tests/domain_tests.cpp +++ b/src/serac/numerics/functional/tests/domain_tests.cpp @@ -8,6 +8,7 @@ #include "serac/numerics/functional/domain.hpp" #include "serac/infrastructure/application_manager.hpp" +#include "serac/mesh_utils/mesh_utils.hpp" using namespace serac; @@ -31,16 +32,17 @@ mfem::Mesh import_mesh(std::string meshfile) TEST(domain, of_edges) { { - auto mesh = import_mesh("onehex.mesh"); - Domain d0 = Domain::ofEdges(mesh, std::function([](std::vector x) { - return (0.5 * (x[0][0] + x[1][0])) < 0.25; // x coordinate of edge midpoint - })); + // auto mesh = import_mesh("onehex.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(import_mesh("onehex.mesh"))); + Domain d0 = Domain::ofEdges(*mesh, std::function([](std::vector x) { + return (0.5 * (x[0][0] + x[1][0])) < 0.25; // x coordinate of edge midpoint + })); EXPECT_EQ(d0.edge_ids_.size(), 4); EXPECT_EQ(d0.dim_, 1); - Domain d1 = Domain::ofEdges(mesh, std::function([](std::vector x) { - return (0.5 * (x[0][1] + x[1][1])) < 0.25; // y coordinate of edge midpoint - })); + Domain d1 = Domain::ofEdges(*mesh, std::function([](std::vector x) { + return (0.5 * (x[0][1] + x[1][1])) < 0.25; // y coordinate of edge midpoint + })); EXPECT_EQ(d1.edge_ids_.size(), 4); EXPECT_EQ(d1.dim_, 1); @@ -58,16 +60,17 @@ TEST(domain, of_edges) } { - auto mesh = import_mesh("onetet.mesh"); - Domain d0 = Domain::ofEdges(mesh, std::function([](std::vector x) { - return (0.5 * (x[0][0] + x[1][0])) < 0.25; // x coordinate of edge midpoint - })); + // auto mesh = import_mesh("onetet.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(import_mesh("onetet.mesh"))); + Domain d0 = Domain::ofEdges(*mesh, std::function([](std::vector x) { + return (0.5 * (x[0][0] + x[1][0])) < 0.25; // x coordinate of edge midpoint + })); EXPECT_EQ(d0.edge_ids_.size(), 3); EXPECT_EQ(d0.dim_, 1); - Domain d1 = Domain::ofEdges(mesh, std::function([](std::vector x) { - return (0.5 * (x[0][1] + x[1][1])) < 0.25; // y coordinate of edge midpoint - })); + Domain d1 = Domain::ofEdges(*mesh, std::function([](std::vector x) { + return (0.5 * (x[0][1] + x[1][1])) < 0.25; // y coordinate of edge midpoint + })); EXPECT_EQ(d1.edge_ids_.size(), 3); EXPECT_EQ(d1.dim_, 1); @@ -86,17 +89,18 @@ TEST(domain, of_edges) { constexpr int dim = 2; - auto mesh = import_mesh("beam-quad.mesh"); - mesh.FinalizeQuadMesh(true); - Domain d0 = Domain::ofEdges(mesh, std::function([](std::vector x, int /* bdr_attr */) { - return (0.5 * (x[0][0] + x[1][0])) < 0.25; // x coordinate of edge midpoint - })); + auto bmesh = import_mesh("beam-quad.mesh"); + bmesh.FinalizeQuadMesh(true); + auto mesh = serac::mesh::refineAndDistribute(std::move(bmesh)); + Domain d0 = Domain::ofEdges(*mesh, std::function([](std::vector x, int /* bdr_attr */) { + return (0.5 * (x[0][0] + x[1][0])) < 0.25; // x coordinate of edge midpoint + })); EXPECT_EQ(d0.edge_ids_.size(), 1); EXPECT_EQ(d0.dim_, 1); - Domain d1 = Domain::ofEdges(mesh, std::function([](std::vector x, int /* bdr_attr */) { - return (0.5 * (x[0][1] + x[1][1])) < 0.25; // y coordinate of edge midpoint - })); + Domain d1 = Domain::ofEdges(*mesh, std::function([](std::vector x, int /* bdr_attr */) { + return (0.5 * (x[0][1] + x[1][1])) < 0.25; // y coordinate of edge midpoint + })); EXPECT_EQ(d1.edge_ids_.size(), 8); EXPECT_EQ(d1.dim_, 1); @@ -109,10 +113,10 @@ TEST(domain, of_edges) EXPECT_EQ(d3.dim_, 1); // check that by_attr compiles - Domain d4 = Domain::ofEdges(mesh, by_attr(3)); + Domain d4 = Domain::ofEdges(*mesh, by_attr(3)); EXPECT_EQ(d4.mfem_edge_ids_.size(), 16); - Domain d5 = Domain::ofBoundaryElements(mesh, [](std::vector, int) { return true; }); + Domain d5 = Domain::ofBoundaryElements(*mesh, [](std::vector, int) { return true; }); EXPECT_EQ(d5.edge_ids_.size(), 18); // 1x8 row of quads has 18 boundary edges } } @@ -121,16 +125,17 @@ TEST(domain, of_faces) { { constexpr int dim = 3; - auto mesh = import_mesh("onehex.mesh"); - Domain d0 = Domain::ofFaces(mesh, std::function([](std::vector vertices, int /*bdr_attr*/) { - return average(vertices)[0] < 0.25; // x coordinate of face center - })); + // auto mesh = import_mesh("onehex.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(import_mesh("onehex.mesh"))); + Domain d0 = Domain::ofFaces(*mesh, std::function([](std::vector vertices, int /*bdr_attr*/) { + return average(vertices)[0] < 0.25; // x coordinate of face center + })); EXPECT_EQ(d0.quad_ids_.size(), 1); EXPECT_EQ(d0.dim_, 2); - Domain d1 = Domain::ofFaces(mesh, std::function([](std::vector vertices, int /*bdr_attr*/) { - return average(vertices)[1] < 0.25; // y coordinate of face center - })); + Domain d1 = Domain::ofFaces(*mesh, std::function([](std::vector vertices, int /*bdr_attr*/) { + return average(vertices)[1] < 0.25; // y coordinate of face center + })); EXPECT_EQ(d1.quad_ids_.size(), 1); EXPECT_EQ(d1.dim_, 2); @@ -143,27 +148,28 @@ TEST(domain, of_faces) EXPECT_EQ(d3.dim_, 2); // check that by_attr compiles - Domain d4 = Domain::ofFaces(mesh, by_attr(3)); + Domain d4 = Domain::ofFaces(*mesh, by_attr(3)); - Domain d5 = Domain::ofBoundaryElements(mesh, [](std::vector, int) { return true; }); + Domain d5 = Domain::ofBoundaryElements(*mesh, [](std::vector, int) { return true; }); EXPECT_EQ(d5.quad_ids_.size(), 6); } { constexpr int dim = 3; - auto mesh = import_mesh("onetet.mesh"); - Domain d0 = Domain::ofFaces(mesh, std::function([](std::vector vertices, int /* bdr_attr */) { - // accept face if it contains a vertex whose x coordinate is less than 0.1 - for (auto v : vertices) { - if (v[0] < 0.1) return true; - } - return false; - })); + // auto mesh = import_mesh("onetet.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(import_mesh("onetet.mesh"))); + Domain d0 = Domain::ofFaces(*mesh, std::function([](std::vector vertices, int /* bdr_attr */) { + // accept face if it contains a vertex whose x coordinate is less than 0.1 + for (auto v : vertices) { + if (v[0] < 0.1) return true; + } + return false; + })); EXPECT_EQ(d0.tri_ids_.size(), 4); EXPECT_EQ(d0.dim_, 2); Domain d1 = Domain::ofFaces( - mesh, std::function([](std::vector x, int /* bdr_attr */) { return average(x)[1] < 0.1; })); + *mesh, std::function([](std::vector x, int /* bdr_attr */) { return average(x)[1] < 0.1; })); EXPECT_EQ(d1.tri_ids_.size(), 1); EXPECT_EQ(d1.dim_, 2); @@ -176,24 +182,25 @@ TEST(domain, of_faces) EXPECT_EQ(d3.dim_, 2); // check that by_attr compiles - Domain d4 = Domain::ofFaces(mesh, by_attr(3)); + Domain d4 = Domain::ofFaces(*mesh, by_attr(3)); - Domain d5 = Domain::ofBoundaryElements(mesh, [](std::vector, int) { return true; }); + Domain d5 = Domain::ofBoundaryElements(*mesh, [](std::vector, int) { return true; }); EXPECT_EQ(d5.tri_ids_.size(), 4); } { constexpr int dim = 2; - auto mesh = import_mesh("beam-quad.mesh"); - Domain d0 = Domain::ofFaces(mesh, std::function([](std::vector vertices, int /* attr */) { - return average(vertices)[0] < 2.25; // x coordinate of face center - })); + // auto mesh = import_mesh("beam-quad.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(import_mesh("beam-quad.mesh"))); + Domain d0 = Domain::ofFaces(*mesh, std::function([](std::vector vertices, int /* attr */) { + return average(vertices)[0] < 2.25; // x coordinate of face center + })); EXPECT_EQ(d0.quad_ids_.size(), 2); EXPECT_EQ(d0.dim_, 2); - Domain d1 = Domain::ofFaces(mesh, std::function([](std::vector vertices, int /* attr */) { - return average(vertices)[1] < 0.55; // y coordinate of face center - })); + Domain d1 = Domain::ofFaces(*mesh, std::function([](std::vector vertices, int /* attr */) { + return average(vertices)[1] < 0.55; // y coordinate of face center + })); EXPECT_EQ(d1.quad_ids_.size(), 8); EXPECT_EQ(d1.dim_, 2); @@ -206,7 +213,7 @@ TEST(domain, of_faces) EXPECT_EQ(d3.dim_, 2); // check that by_attr compiles - Domain d4 = Domain::ofFaces(mesh, by_attr(3)); + Domain d4 = Domain::ofFaces(*mesh, by_attr(3)); } } @@ -214,18 +221,19 @@ TEST(domain, of_elements) { { constexpr int dim = 3; - auto mesh = import_mesh("patch3D_tets_and_hexes.mesh"); - Domain d0 = Domain::ofElements(mesh, std::function([](std::vector vertices, int /*bdr_attr*/) { - return average(vertices)[0] < 0.7; // x coordinate of face center - })); + // auto mesh = import_mesh("patch3D_tets_and_hexes.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(import_mesh("patch3D_tets_and_hexes.mesh"))); + Domain d0 = Domain::ofElements(*mesh, std::function([](std::vector vertices, int /*bdr_attr*/) { + return average(vertices)[0] < 0.7; // x coordinate of face center + })); EXPECT_EQ(d0.tet_ids_.size(), 0); EXPECT_EQ(d0.hex_ids_.size(), 1); EXPECT_EQ(d0.dim_, 3); - Domain d1 = Domain::ofElements(mesh, std::function([](std::vector vertices, int /*bdr_attr*/) { - return average(vertices)[1] < 0.75; // y coordinate of face center - })); + Domain d1 = Domain::ofElements(*mesh, std::function([](std::vector vertices, int /*bdr_attr*/) { + return average(vertices)[1] < 0.75; // y coordinate of face center + })); EXPECT_EQ(d1.tet_ids_.size(), 6); EXPECT_EQ(d1.hex_ids_.size(), 1); EXPECT_EQ(d1.dim_, 3); @@ -241,20 +249,21 @@ TEST(domain, of_elements) EXPECT_EQ(d3.dim_, 3); // check that by_attr works - Domain d4 = Domain::ofElements(mesh, by_attr(3)); + Domain d4 = Domain::ofElements(*mesh, by_attr(3)); } { constexpr int dim = 2; - auto mesh = import_mesh("patch2D_tris_and_quads.mesh"); + // auto mesh = import_mesh("patch2D_tris_and_quads.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(import_mesh("patch2D_tris_and_quads.mesh"))); Domain d0 = Domain::ofElements( - mesh, std::function([](std::vector vertices, int /* attr */) { return average(vertices)[0] < 0.45; })); + *mesh, std::function([](std::vector vertices, int /* attr */) { return average(vertices)[0] < 0.45; })); EXPECT_EQ(d0.tri_ids_.size(), 1); EXPECT_EQ(d0.quad_ids_.size(), 1); EXPECT_EQ(d0.dim_, 2); Domain d1 = Domain::ofElements( - mesh, std::function([](std::vector vertices, int /* attr */) { return average(vertices)[1] < 0.45; })); + *mesh, std::function([](std::vector vertices, int /* attr */) { return average(vertices)[1] < 0.45; })); EXPECT_EQ(d1.tri_ids_.size(), 1); EXPECT_EQ(d1.quad_ids_.size(), 1); EXPECT_EQ(d1.dim_, 2); @@ -270,7 +279,7 @@ TEST(domain, of_elements) EXPECT_EQ(d3.dim_, 2); // check that by_attr compiles - Domain d4 = Domain::ofElements(mesh, by_attr(3)); + Domain d4 = Domain::ofElements(*mesh, by_attr(3)); } } @@ -278,16 +287,17 @@ TEST(domain, entireDomain2d) { constexpr int dim = 2; constexpr int p = 1; - auto mesh = import_mesh("patch2D_tris_and_quads.mesh"); + // auto mesh = import_mesh("patch2D_tris_and_quads.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(import_mesh("patch2D_tris_and_quads.mesh"))); - Domain d0 = EntireDomain(mesh); + Domain d0 = EntireDomain(*mesh); EXPECT_EQ(d0.dim_, 2); EXPECT_EQ(d0.tri_ids_.size(), 2); EXPECT_EQ(d0.quad_ids_.size(), 4); auto fec = mfem::H1_FECollection(p, dim); - auto fes = mfem::FiniteElementSpace(&mesh, &fec); + auto fes = mfem::ParFiniteElementSpace(mesh.get(), &fec); mfem::Array dof_indices = d0.dof_list(&fes); EXPECT_EQ(dof_indices.Size(), 8); @@ -297,16 +307,17 @@ TEST(domain, entireDomain3d) { constexpr int dim = 3; constexpr int p = 1; - auto mesh = import_mesh("patch3D_tets_and_hexes.mesh"); + // auto mesh = import_mesh("patch3D_tets_and_hexes.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(import_mesh("patch3D_tets_and_hexes.mesh"))); - Domain d0 = EntireDomain(mesh); + Domain d0 = EntireDomain(*mesh); EXPECT_EQ(d0.dim_, 3); EXPECT_EQ(d0.tet_ids_.size(), 12); EXPECT_EQ(d0.hex_ids_.size(), 7); auto fec = mfem::H1_FECollection(p, dim); - auto fes = mfem::FiniteElementSpace(&mesh, &fec); + auto fes = mfem::ParFiniteElementSpace(mesh.get(), &fec); mfem::Array dof_indices = d0.dof_list(&fes); EXPECT_EQ(dof_indices.Size(), 25); @@ -316,17 +327,18 @@ TEST(domain, of2dElementsFindsDofs) { constexpr int dim = 2; constexpr int p = 2; - auto mesh = import_mesh("patch2D_tris_and_quads.mesh"); + // auto mesh = import_mesh("patch2D_tris_and_quads.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(import_mesh("patch2D_tris_and_quads.mesh"))); auto fec = mfem::H1_FECollection(p, dim); - auto fes = mfem::FiniteElementSpace(&mesh, &fec); + auto fes = mfem::ParFiniteElementSpace(mesh.get(), &fec); auto find_element_0 = [](std::vector vertices, int /* attr */) { auto centroid = average(vertices); return (centroid[0] < 0.5) && (centroid[1] < 0.25); }; - Domain d0 = Domain::ofElements(mesh, find_element_0); + Domain d0 = Domain::ofElements(*mesh, find_element_0); mfem::Array dof_indices = d0.dof_list(&fes); @@ -339,7 +351,7 @@ TEST(domain, of2dElementsFindsDofs) tensor target{{0.533, 0.424}}; return norm(centroid - target) < 1e-2; }; - Domain d1 = Domain::ofElements(mesh, find_element_4); + Domain d1 = Domain::ofElements(*mesh, find_element_4); Domain elements_0_and_4 = d0 | d1; @@ -348,7 +360,7 @@ TEST(domain, of2dElementsFindsDofs) /////////////////////////////////////// - Domain d2 = EntireDomain(mesh) - elements_0_and_4; + Domain d2 = EntireDomain(*mesh) - elements_0_and_4; dof_indices = d2.dof_list(&fes); @@ -359,10 +371,11 @@ TEST(domain, of3dElementsFindsDofs) { constexpr int dim = 3; constexpr int p = 2; - auto mesh = import_mesh("patch3D_tets_and_hexes.mesh"); + // auto mesh = import_mesh("patch3D_tets_and_hexes.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(import_mesh("patch3D_tets_and_hexes.mesh"))); auto fec = mfem::H1_FECollection(p, dim); - auto fes = mfem::FiniteElementSpace(&mesh, &fec); + auto fes = mfem::ParFiniteElementSpace(mesh.get(), &fec); auto find_element_0 = [](std::vector vertices, int /* attr */) { auto centroid = average(vertices); @@ -370,7 +383,7 @@ TEST(domain, of3dElementsFindsDofs) return norm(centroid - target) < 1e-2; }; - Domain d0 = Domain::ofElements(mesh, find_element_0); + Domain d0 = Domain::ofElements(*mesh, find_element_0); mfem::Array dof_indices = d0.dof_list(&fes); @@ -384,7 +397,7 @@ TEST(domain, of3dElementsFindsDofs) vec3 target{{3.275, 1.2, 0.725}}; return norm(centroid - target) < 1e-2; }; - Domain d1 = Domain::ofElements(mesh, find_element_1); + Domain d1 = Domain::ofElements(*mesh, find_element_1); Domain elements_0_and_1 = d0 | d1; @@ -395,7 +408,7 @@ TEST(domain, of3dElementsFindsDofs) ///////////////////////////////////////// - Domain d2 = EntireDomain(mesh) - elements_0_and_1; + Domain d2 = EntireDomain(*mesh) - elements_0_and_1; dof_indices = d2.dof_list(&fes); @@ -406,17 +419,18 @@ TEST(domain, of2dBoundaryElementsFindsDofs) { constexpr int dim = 2; constexpr int p = 2; - auto mesh = import_mesh("patch2D_tris_and_quads.mesh"); + // auto mesh = import_mesh("patch2D_tris_and_quads.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(import_mesh("patch2D_tris_and_quads.mesh"))); auto find_right_boundary = [](std::vector vertices, int /* attr */) { return std::all_of(vertices.begin(), vertices.end(), [](vec2 X) { return X[0] > 1.0 - 1e-2; }); }; - Domain d0 = Domain::ofBoundaryElements(mesh, find_right_boundary); + Domain d0 = Domain::ofBoundaryElements(*mesh, find_right_boundary); EXPECT_EQ(d0.edge_ids_.size(), 1); auto fec = mfem::H1_FECollection(p, dim); - auto fes = mfem::FiniteElementSpace(&mesh, &fec); + auto fes = mfem::ParFiniteElementSpace(mesh.get(), &fec); mfem::Array dof_indices = d0.dof_list(&fes); @@ -426,7 +440,7 @@ TEST(domain, of2dBoundaryElementsFindsDofs) return std::all_of(vertices.begin(), vertices.end(), [](vec2 X) { return X[1] > 1.0 - 1e-2; }); }; - Domain d1 = Domain::ofBoundaryElements(mesh, find_top_boundary); + Domain d1 = Domain::ofBoundaryElements(*mesh, find_top_boundary); EXPECT_EQ(d1.edge_ids_.size(), 1); Domain d2 = d0 | d1; @@ -440,17 +454,18 @@ TEST(domain, of3dBoundaryElementsFindsDofs) { constexpr int dim = 3; constexpr int p = 2; - auto mesh = import_mesh("patch3D_tets.mesh"); + // auto mesh = import_mesh("patch3D_tets.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(import_mesh("patch3D_tets.mesh"))); auto find_xmax_boundary = [](std::vector vertices, int /* attr */) { return std::all_of(vertices.begin(), vertices.end(), [](vec3 X) { return X[0] > 1.0 - 1e-2; }); }; - Domain d0 = Domain::ofBoundaryElements(mesh, find_xmax_boundary); + Domain d0 = Domain::ofBoundaryElements(*mesh, find_xmax_boundary); EXPECT_EQ(d0.tri_ids_.size(), 2); auto fec = mfem::H1_FECollection(p, dim); - auto fes = mfem::FiniteElementSpace(&mesh, &fec); + auto fes = mfem::ParFiniteElementSpace(mesh.get(), &fec); mfem::Array dof_indices = d0.dof_list(&fes); @@ -460,7 +475,7 @@ TEST(domain, of3dBoundaryElementsFindsDofs) return std::all_of(vertices.begin(), vertices.end(), [](vec3 X) { return X[1] > 1.0 - 1e-2; }); }; - Domain d1 = Domain::ofBoundaryElements(mesh, find_ymax_boundary); + Domain d1 = Domain::ofBoundaryElements(*mesh, find_ymax_boundary); EXPECT_EQ(d1.tri_ids_.size(), 2); Domain d2 = d0 | d1; diff --git a/src/serac/numerics/functional/tests/element_restriction_tests.cpp b/src/serac/numerics/functional/tests/element_restriction_tests.cpp index 361ccfca14..02da2090c6 100644 --- a/src/serac/numerics/functional/tests/element_restriction_tests.cpp +++ b/src/serac/numerics/functional/tests/element_restriction_tests.cpp @@ -8,6 +8,7 @@ #include "serac/numerics/functional/domain.hpp" #include "serac/numerics/functional/element_restriction.hpp" +#include "serac/mesh_utils/mesh_utils.hpp" using namespace serac; @@ -26,18 +27,19 @@ TEST(patch_test_meshes, triangle_domains) { int p = 2; int dim = 2; - mfem::Mesh mesh(SERAC_REPO_DIR "/data/meshes/patch2D_tris.mesh"); + mfem::Mesh bmesh(SERAC_REPO_DIR "/data/meshes/patch2D_tris.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(bmesh)); auto H1_fec = std::make_unique(p, dim); - auto H1_fes = std::make_unique(&mesh, H1_fec.get()); + auto H1_fes = std::make_unique(mesh.get(), H1_fec.get()); auto Hcurl_fec = std::make_unique(p, dim); - auto Hcurl_fes = std::make_unique(&mesh, Hcurl_fec.get()); + auto Hcurl_fes = std::make_unique(mesh.get(), Hcurl_fec.get()); auto L2_fec = std::make_unique(p, dim, mfem::BasisType::GaussLobatto); - auto L2_fes = std::make_unique(&mesh, L2_fec.get()); + auto L2_fes = std::make_unique(mesh.get(), L2_fec.get()); - Domain whole = EntireDomain(mesh); + Domain whole = EntireDomain(*mesh); EXPECT_EQ(whole.mfem_edge_ids_.size(), 0); EXPECT_EQ(whole.mfem_tri_ids_.size(), 4); EXPECT_EQ(whole.mfem_quad_ids_.size(), 0); @@ -55,7 +57,7 @@ TEST(patch_test_meshes, triangle_domains) EXPECT_EQ(L2_BER.ESize(), 4 * 6); } - Domain boundary = EntireBoundary(mesh); + Domain boundary = EntireBoundary(*mesh); EXPECT_EQ(boundary.mfem_edge_ids_.size(), 4); EXPECT_EQ(boundary.mfem_tri_ids_.size(), 0); EXPECT_EQ(boundary.mfem_quad_ids_.size(), 0); @@ -73,7 +75,7 @@ TEST(patch_test_meshes, triangle_domains) EXPECT_EQ(L2_BER.ESize(), 4 * 3); } - Domain interior = InteriorFaces(mesh); + Domain interior = InteriorFaces(*mesh); EXPECT_EQ(interior.mfem_edge_ids_.size(), 4); EXPECT_EQ(interior.mfem_tri_ids_.size(), 0); EXPECT_EQ(interior.mfem_quad_ids_.size(), 0); @@ -96,17 +98,18 @@ TEST(patch_test_meshes, quadrilateral_domains) { int p = 2; int dim = 2; - mfem::Mesh mesh(SERAC_REPO_DIR "/data/meshes/patch2D_quads.mesh"); + mfem::Mesh bmesh(SERAC_REPO_DIR "/data/meshes/patch2D_quads.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(bmesh)); auto H1_fec = std::make_unique(p, dim); auto Hcurl_fec = std::make_unique(p, dim); auto L2_fec = std::make_unique(p, dim, mfem::BasisType::GaussLobatto); - auto H1_fes = std::make_unique(&mesh, H1_fec.get()); - auto Hcurl_fes = std::make_unique(&mesh, Hcurl_fec.get()); - auto L2_fes = std::make_unique(&mesh, L2_fec.get()); + auto H1_fes = std::make_unique(mesh.get(), H1_fec.get()); + auto Hcurl_fes = std::make_unique(mesh.get(), Hcurl_fec.get()); + auto L2_fes = std::make_unique(mesh.get(), L2_fec.get()); - Domain whole = EntireDomain(mesh); + Domain whole = EntireDomain(*mesh); EXPECT_EQ(whole.mfem_edge_ids_.size(), 0); EXPECT_EQ(whole.mfem_tri_ids_.size(), 0); EXPECT_EQ(whole.mfem_quad_ids_.size(), 5); @@ -124,7 +127,7 @@ TEST(patch_test_meshes, quadrilateral_domains) EXPECT_EQ(L2_BER.ESize(), 5 * 9); } - Domain boundary = EntireBoundary(mesh); + Domain boundary = EntireBoundary(*mesh); EXPECT_EQ(boundary.mfem_edge_ids_.size(), 4); EXPECT_EQ(boundary.mfem_tri_ids_.size(), 0); EXPECT_EQ(boundary.mfem_quad_ids_.size(), 0); @@ -142,7 +145,7 @@ TEST(patch_test_meshes, quadrilateral_domains) EXPECT_EQ(L2_BER.ESize(), 4 * 3); } - Domain interior = InteriorFaces(mesh); + Domain interior = InteriorFaces(*mesh); EXPECT_EQ(interior.mfem_edge_ids_.size(), 8); EXPECT_EQ(interior.mfem_tri_ids_.size(), 0); EXPECT_EQ(interior.mfem_quad_ids_.size(), 0); @@ -165,17 +168,18 @@ TEST(patch_test_meshes, triangle_and_quadrilateral_domains) { int p = 2; int dim = 2; - mfem::Mesh mesh(SERAC_REPO_DIR "/data/meshes/patch2D_tris_and_quads.mesh"); + mfem::Mesh bmesh(SERAC_REPO_DIR "/data/meshes/patch2D_tris_and_quads.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(bmesh)); auto H1_fec = std::make_unique(p, dim); auto Hcurl_fec = std::make_unique(p, dim); auto L2_fec = std::make_unique(p, dim, mfem::BasisType::GaussLobatto); - auto H1_fes = std::make_unique(&mesh, H1_fec.get()); - auto Hcurl_fes = std::make_unique(&mesh, Hcurl_fec.get()); - auto L2_fes = std::make_unique(&mesh, L2_fec.get()); + auto H1_fes = std::make_unique(mesh.get(), H1_fec.get()); + auto Hcurl_fes = std::make_unique(mesh.get(), Hcurl_fec.get()); + auto L2_fes = std::make_unique(mesh.get(), L2_fec.get()); - Domain whole = EntireDomain(mesh); + Domain whole = EntireDomain(*mesh); EXPECT_EQ(whole.mfem_edge_ids_.size(), 0); EXPECT_EQ(whole.mfem_tri_ids_.size(), 2); EXPECT_EQ(whole.mfem_quad_ids_.size(), 4); @@ -193,7 +197,7 @@ TEST(patch_test_meshes, triangle_and_quadrilateral_domains) EXPECT_EQ(L2_BER.ESize(), 2 * 6 + 4 * 9); } - Domain boundary = EntireBoundary(mesh); + Domain boundary = EntireBoundary(*mesh); EXPECT_EQ(boundary.mfem_edge_ids_.size(), 4); EXPECT_EQ(boundary.mfem_tri_ids_.size(), 0); EXPECT_EQ(boundary.mfem_quad_ids_.size(), 0); @@ -211,7 +215,7 @@ TEST(patch_test_meshes, triangle_and_quadrilateral_domains) EXPECT_EQ(L2_BER.ESize(), 4 * 3); } - Domain interior = InteriorFaces(mesh); + Domain interior = InteriorFaces(*mesh); EXPECT_EQ(interior.mfem_edge_ids_.size(), 9); EXPECT_EQ(interior.mfem_tri_ids_.size(), 0); EXPECT_EQ(interior.mfem_quad_ids_.size(), 0); @@ -234,17 +238,18 @@ TEST(patch_test_meshes, tetrahedron_domains) { int p = 2; int dim = 3; - mfem::Mesh mesh(SERAC_REPO_DIR "/data/meshes/patch3D_tets.mesh"); + mfem::Mesh bmesh(SERAC_REPO_DIR "/data/meshes/patch3D_tets.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(bmesh)); auto H1_fec = std::make_unique(p, dim); auto Hcurl_fec = std::make_unique(p, dim); auto L2_fec = std::make_unique(p, dim, mfem::BasisType::GaussLobatto); - auto H1_fes = std::make_unique(&mesh, H1_fec.get()); - auto Hcurl_fes = std::make_unique(&mesh, Hcurl_fec.get()); - auto L2_fes = std::make_unique(&mesh, L2_fec.get()); + auto H1_fes = std::make_unique(mesh.get(), H1_fec.get()); + auto Hcurl_fes = std::make_unique(mesh.get(), Hcurl_fec.get()); + auto L2_fes = std::make_unique(mesh.get(), L2_fec.get()); - Domain whole = EntireDomain(mesh); + Domain whole = EntireDomain(*mesh); EXPECT_EQ(whole.mfem_edge_ids_.size(), 0); EXPECT_EQ(whole.mfem_tri_ids_.size(), 0); EXPECT_EQ(whole.mfem_quad_ids_.size(), 0); @@ -262,7 +267,7 @@ TEST(patch_test_meshes, tetrahedron_domains) EXPECT_EQ(L2_BER.ESize(), 12 * 10); } - Domain boundary = EntireBoundary(mesh); + Domain boundary = EntireBoundary(*mesh); EXPECT_EQ(boundary.mfem_edge_ids_.size(), 0); EXPECT_EQ(boundary.mfem_tri_ids_.size(), 12); EXPECT_EQ(boundary.mfem_quad_ids_.size(), 0); @@ -280,7 +285,7 @@ TEST(patch_test_meshes, tetrahedron_domains) EXPECT_EQ(L2_BER.ESize(), 12 * 6); } - Domain interior = InteriorFaces(mesh); + Domain interior = InteriorFaces(*mesh); EXPECT_EQ(interior.mfem_edge_ids_.size(), 0); EXPECT_EQ(interior.mfem_tri_ids_.size(), 18); EXPECT_EQ(interior.mfem_quad_ids_.size(), 0); @@ -303,17 +308,18 @@ TEST(patch_test_meshes, hexahedron_domains) { int p = 2; int dim = 3; - mfem::Mesh mesh(SERAC_REPO_DIR "/data/meshes/patch3D_hexes.mesh"); + mfem::Mesh bmesh(SERAC_REPO_DIR "/data/meshes/patch3D_hexes.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(bmesh)); auto H1_fec = std::make_unique(p, dim); auto Hcurl_fec = std::make_unique(p, dim); auto L2_fec = std::make_unique(p, dim, mfem::BasisType::GaussLobatto); - auto H1_fes = std::make_unique(&mesh, H1_fec.get()); - auto Hcurl_fes = std::make_unique(&mesh, Hcurl_fec.get()); - auto L2_fes = std::make_unique(&mesh, L2_fec.get()); + auto H1_fes = std::make_unique(mesh.get(), H1_fec.get()); + auto Hcurl_fes = std::make_unique(mesh.get(), Hcurl_fec.get()); + auto L2_fes = std::make_unique(mesh.get(), L2_fec.get()); - Domain whole = EntireDomain(mesh); + Domain whole = EntireDomain(*mesh); EXPECT_EQ(whole.mfem_edge_ids_.size(), 0); EXPECT_EQ(whole.mfem_tri_ids_.size(), 0); EXPECT_EQ(whole.mfem_quad_ids_.size(), 0); @@ -331,7 +337,7 @@ TEST(patch_test_meshes, hexahedron_domains) EXPECT_EQ(L2_BER.ESize(), 7 * 27); } - Domain boundary = EntireBoundary(mesh); + Domain boundary = EntireBoundary(*mesh); EXPECT_EQ(boundary.mfem_edge_ids_.size(), 0); EXPECT_EQ(boundary.mfem_tri_ids_.size(), 0); EXPECT_EQ(boundary.mfem_quad_ids_.size(), 6); @@ -349,7 +355,7 @@ TEST(patch_test_meshes, hexahedron_domains) EXPECT_EQ(L2_BER.ESize(), 6 * 9); } - Domain interior = InteriorFaces(mesh); + Domain interior = InteriorFaces(*mesh); EXPECT_EQ(interior.mfem_edge_ids_.size(), 0); EXPECT_EQ(interior.mfem_tri_ids_.size(), 0); EXPECT_EQ(interior.mfem_quad_ids_.size(), 18); diff --git a/src/serac/numerics/functional/tests/geometric_factors_tests.cpp b/src/serac/numerics/functional/tests/geometric_factors_tests.cpp index 3a7a06f171..fdca0729b3 100644 --- a/src/serac/numerics/functional/tests/geometric_factors_tests.cpp +++ b/src/serac/numerics/functional/tests/geometric_factors_tests.cpp @@ -8,6 +8,7 @@ #include "serac/numerics/functional/geometric_factors.hpp" #include "serac/infrastructure/application_manager.hpp" +#include "serac/mesh_utils/mesh_utils.hpp" using namespace serac; @@ -30,11 +31,12 @@ mfem::Mesh import_mesh(std::string meshfile) TEST(geometric_factors, with_2D_domains) { - auto mesh = import_mesh("patch2D_tris_and_quads.mesh"); + // auto mesh = import_mesh("patch2D_tris_and_quads.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(import_mesh("patch2D_tris_and_quads.mesh"))); // `d` will consist of one tri and one quad Domain d = Domain::ofElements( - mesh, std::function([](std::vector vertices, int /* attr */) { return average(vertices)[0] < 0.45; })); + *mesh, std::function([](std::vector vertices, int /* attr */) { return average(vertices)[0] < 0.45; })); int q = 2; @@ -59,12 +61,13 @@ TEST(geometric_factors, with_2D_domains) TEST(geometric_factors, with_3D_domains) { - auto mesh = import_mesh("patch3D_tets_and_hexes.mesh"); + // auto mesh = import_mesh("patch3D_tets_and_hexes.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(import_mesh("patch3D_tets_and_hexes.mesh"))); // `d` will consist of 6 tets and 1 hex - Domain d = Domain::ofElements(mesh, std::function([](std::vector vertices, int /*bdr_attr*/) { - return average(vertices)[1] < 0.75; // y coordinate of face center - })); + Domain d = Domain::ofElements(*mesh, std::function([](std::vector vertices, int /*bdr_attr*/) { + return average(vertices)[1] < 0.75; // y coordinate of face center + })); int q = 2; diff --git a/src/serac/numerics/functional/typedefs.hpp b/src/serac/numerics/functional/typedefs.hpp index dad72993f0..62fb91adc4 100644 --- a/src/serac/numerics/functional/typedefs.hpp +++ b/src/serac/numerics/functional/typedefs.hpp @@ -12,9 +12,10 @@ namespace serac { // mfem::ParMesh / mfeme::ParFiniteElementSpace, but I'm not ready to pull the trigger on a big // interface change like that, so these typedefs mark the parts that would need to eventually change +// hanyu: I'm pulling the trigger for sam and the bullet might hit myself /// @cond -using mesh_t = mfem::Mesh; -using fes_t = mfem::FiniteElementSpace; +using mesh_t = mfem::ParMesh; +using fes_t = mfem::ParFiniteElementSpace; /// @endcond } // namespace serac diff --git a/src/serac/physics/state/finite_element_state.cpp b/src/serac/physics/state/finite_element_state.cpp index 0bb4d10d1c..29e81d7bda 100644 --- a/src/serac/physics/state/finite_element_state.cpp +++ b/src/serac/physics/state/finite_element_state.cpp @@ -75,7 +75,7 @@ void FiniteElementState::projectOnBoundary(mfem::VectorCoefficient& coef, const void FiniteElementState::project(mfem::Coefficient& coef, const Domain& domain) { - mfem::Array uniq_dof_ids = domain.dof_list(gridFunction().FESpace()); + mfem::Array uniq_dof_ids = domain.dof_list(gridFunction().ParFESpace()); mfem::ParGridFunction& grid_function = gridFunction(); grid_function.ProjectCoefficient(coef, uniq_dof_ids); setFromGridFunction(grid_function); @@ -83,7 +83,7 @@ void FiniteElementState::project(mfem::Coefficient& coef, const Domain& domain) void FiniteElementState::project(mfem::VectorCoefficient& coef, const Domain& domain) { - mfem::Array uniq_dof_ids = domain.dof_list(gridFunction().FESpace()); + mfem::Array uniq_dof_ids = domain.dof_list(gridFunction().ParFESpace()); mfem::ParGridFunction& grid_function = gridFunction(); grid_function.ProjectCoefficient(coef, uniq_dof_ids); setFromGridFunction(grid_function); From ccf832901029a0eea1b49cacf1c2cba7c7289473 Mon Sep 17 00:00:00 2001 From: li97 Date: Mon, 23 Jun 2025 19:14:13 -0700 Subject: [PATCH 03/11] Fixed three tests --- .../functional/tests/domain_tests.cpp | 52 +++++++++---------- .../tests/element_restriction_tests.cpp | 8 +++ .../tests/geometric_factors_tests.cpp | 8 +-- 3 files changed, 38 insertions(+), 30 deletions(-) diff --git a/src/serac/numerics/functional/tests/domain_tests.cpp b/src/serac/numerics/functional/tests/domain_tests.cpp index e8bf4dd82f..5ce35c8553 100644 --- a/src/serac/numerics/functional/tests/domain_tests.cpp +++ b/src/serac/numerics/functional/tests/domain_tests.cpp @@ -32,8 +32,8 @@ mfem::Mesh import_mesh(std::string meshfile) TEST(domain, of_edges) { { - // auto mesh = import_mesh("onehex.mesh"); - auto mesh = serac::mesh::refineAndDistribute(std::move(import_mesh("onehex.mesh"))); + auto bmesh = import_mesh("onehex.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(bmesh)); Domain d0 = Domain::ofEdges(*mesh, std::function([](std::vector x) { return (0.5 * (x[0][0] + x[1][0])) < 0.25; // x coordinate of edge midpoint })); @@ -60,8 +60,8 @@ TEST(domain, of_edges) } { - // auto mesh = import_mesh("onetet.mesh"); - auto mesh = serac::mesh::refineAndDistribute(std::move(import_mesh("onetet.mesh"))); + auto bmesh = import_mesh("onetet.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(bmesh)); Domain d0 = Domain::ofEdges(*mesh, std::function([](std::vector x) { return (0.5 * (x[0][0] + x[1][0])) < 0.25; // x coordinate of edge midpoint })); @@ -125,8 +125,8 @@ TEST(domain, of_faces) { { constexpr int dim = 3; - // auto mesh = import_mesh("onehex.mesh"); - auto mesh = serac::mesh::refineAndDistribute(std::move(import_mesh("onehex.mesh"))); + auto bmesh = import_mesh("onehex.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(bmesh)); Domain d0 = Domain::ofFaces(*mesh, std::function([](std::vector vertices, int /*bdr_attr*/) { return average(vertices)[0] < 0.25; // x coordinate of face center })); @@ -156,8 +156,8 @@ TEST(domain, of_faces) { constexpr int dim = 3; - // auto mesh = import_mesh("onetet.mesh"); - auto mesh = serac::mesh::refineAndDistribute(std::move(import_mesh("onetet.mesh"))); + auto bmesh = import_mesh("onetet.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(bmesh)); Domain d0 = Domain::ofFaces(*mesh, std::function([](std::vector vertices, int /* bdr_attr */) { // accept face if it contains a vertex whose x coordinate is less than 0.1 for (auto v : vertices) { @@ -190,8 +190,8 @@ TEST(domain, of_faces) { constexpr int dim = 2; - // auto mesh = import_mesh("beam-quad.mesh"); - auto mesh = serac::mesh::refineAndDistribute(std::move(import_mesh("beam-quad.mesh"))); + auto bmesh = import_mesh("beam-quad.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(bmesh)); Domain d0 = Domain::ofFaces(*mesh, std::function([](std::vector vertices, int /* attr */) { return average(vertices)[0] < 2.25; // x coordinate of face center })); @@ -221,8 +221,8 @@ TEST(domain, of_elements) { { constexpr int dim = 3; - // auto mesh = import_mesh("patch3D_tets_and_hexes.mesh"); - auto mesh = serac::mesh::refineAndDistribute(std::move(import_mesh("patch3D_tets_and_hexes.mesh"))); + auto bmesh = import_mesh("patch3D_tets_and_hexes.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(bmesh)); Domain d0 = Domain::ofElements(*mesh, std::function([](std::vector vertices, int /*bdr_attr*/) { return average(vertices)[0] < 0.7; // x coordinate of face center })); @@ -254,8 +254,8 @@ TEST(domain, of_elements) { constexpr int dim = 2; - // auto mesh = import_mesh("patch2D_tris_and_quads.mesh"); - auto mesh = serac::mesh::refineAndDistribute(std::move(import_mesh("patch2D_tris_and_quads.mesh"))); + auto bmesh = import_mesh("patch2D_tris_and_quads.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(bmesh)); Domain d0 = Domain::ofElements( *mesh, std::function([](std::vector vertices, int /* attr */) { return average(vertices)[0] < 0.45; })); EXPECT_EQ(d0.tri_ids_.size(), 1); @@ -287,8 +287,8 @@ TEST(domain, entireDomain2d) { constexpr int dim = 2; constexpr int p = 1; - // auto mesh = import_mesh("patch2D_tris_and_quads.mesh"); - auto mesh = serac::mesh::refineAndDistribute(std::move(import_mesh("patch2D_tris_and_quads.mesh"))); + auto bmesh = import_mesh("patch2D_tris_and_quads.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(bmesh)); Domain d0 = EntireDomain(*mesh); @@ -307,8 +307,8 @@ TEST(domain, entireDomain3d) { constexpr int dim = 3; constexpr int p = 1; - // auto mesh = import_mesh("patch3D_tets_and_hexes.mesh"); - auto mesh = serac::mesh::refineAndDistribute(std::move(import_mesh("patch3D_tets_and_hexes.mesh"))); + auto bmesh = import_mesh("patch3D_tets_and_hexes.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(bmesh)); Domain d0 = EntireDomain(*mesh); @@ -327,8 +327,8 @@ TEST(domain, of2dElementsFindsDofs) { constexpr int dim = 2; constexpr int p = 2; - // auto mesh = import_mesh("patch2D_tris_and_quads.mesh"); - auto mesh = serac::mesh::refineAndDistribute(std::move(import_mesh("patch2D_tris_and_quads.mesh"))); + auto bmesh = import_mesh("patch2D_tris_and_quads.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(bmesh)); auto fec = mfem::H1_FECollection(p, dim); auto fes = mfem::ParFiniteElementSpace(mesh.get(), &fec); @@ -371,8 +371,8 @@ TEST(domain, of3dElementsFindsDofs) { constexpr int dim = 3; constexpr int p = 2; - // auto mesh = import_mesh("patch3D_tets_and_hexes.mesh"); - auto mesh = serac::mesh::refineAndDistribute(std::move(import_mesh("patch3D_tets_and_hexes.mesh"))); + auto bmesh = import_mesh("patch3D_tets_and_hexes.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(bmesh)); auto fec = mfem::H1_FECollection(p, dim); auto fes = mfem::ParFiniteElementSpace(mesh.get(), &fec); @@ -419,8 +419,8 @@ TEST(domain, of2dBoundaryElementsFindsDofs) { constexpr int dim = 2; constexpr int p = 2; - // auto mesh = import_mesh("patch2D_tris_and_quads.mesh"); - auto mesh = serac::mesh::refineAndDistribute(std::move(import_mesh("patch2D_tris_and_quads.mesh"))); + auto bmesh = import_mesh("patch2D_tris_and_quads.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(bmesh)); auto find_right_boundary = [](std::vector vertices, int /* attr */) { return std::all_of(vertices.begin(), vertices.end(), [](vec2 X) { return X[0] > 1.0 - 1e-2; }); @@ -454,8 +454,8 @@ TEST(domain, of3dBoundaryElementsFindsDofs) { constexpr int dim = 3; constexpr int p = 2; - // auto mesh = import_mesh("patch3D_tets.mesh"); - auto mesh = serac::mesh::refineAndDistribute(std::move(import_mesh("patch3D_tets.mesh"))); + auto bmesh = import_mesh("patch3D_tets.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(bmesh)); auto find_xmax_boundary = [](std::vector vertices, int /* attr */) { return std::all_of(vertices.begin(), vertices.end(), [](vec3 X) { return X[0] > 1.0 - 1e-2; }); diff --git a/src/serac/numerics/functional/tests/element_restriction_tests.cpp b/src/serac/numerics/functional/tests/element_restriction_tests.cpp index 02da2090c6..2da708fb75 100644 --- a/src/serac/numerics/functional/tests/element_restriction_tests.cpp +++ b/src/serac/numerics/functional/tests/element_restriction_tests.cpp @@ -8,6 +8,7 @@ #include "serac/numerics/functional/domain.hpp" #include "serac/numerics/functional/element_restriction.hpp" +#include "serac/infrastructure/application_manager.hpp" #include "serac/mesh_utils/mesh_utils.hpp" using namespace serac; @@ -373,3 +374,10 @@ TEST(patch_test_meshes, hexahedron_domains) EXPECT_EQ(L2_BER.ESize(), 18 * (9 * 2)); } } + +int main(int argc, char* argv[]) +{ + ::testing::InitGoogleTest(&argc, argv); + serac::ApplicationManager applicationManager(argc, argv); + return RUN_ALL_TESTS(); +} diff --git a/src/serac/numerics/functional/tests/geometric_factors_tests.cpp b/src/serac/numerics/functional/tests/geometric_factors_tests.cpp index fdca0729b3..baea47021c 100644 --- a/src/serac/numerics/functional/tests/geometric_factors_tests.cpp +++ b/src/serac/numerics/functional/tests/geometric_factors_tests.cpp @@ -31,8 +31,8 @@ mfem::Mesh import_mesh(std::string meshfile) TEST(geometric_factors, with_2D_domains) { - // auto mesh = import_mesh("patch2D_tris_and_quads.mesh"); - auto mesh = serac::mesh::refineAndDistribute(std::move(import_mesh("patch2D_tris_and_quads.mesh"))); + auto bmesh = import_mesh("patch2D_tris_and_quads.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(bmesh)); // `d` will consist of one tri and one quad Domain d = Domain::ofElements( @@ -61,8 +61,8 @@ TEST(geometric_factors, with_2D_domains) TEST(geometric_factors, with_3D_domains) { - // auto mesh = import_mesh("patch3D_tets_and_hexes.mesh"); - auto mesh = serac::mesh::refineAndDistribute(std::move(import_mesh("patch3D_tets_and_hexes.mesh"))); + auto bmesh = import_mesh("patch3D_tets_and_hexes.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(bmesh)); // `d` will consist of 6 tets and 1 hex Domain d = Domain::ofElements(*mesh, std::function([](std::vector vertices, int /*bdr_attr*/) { From 6d088cf3adef2c318f76fda897a2725bda62085d Mon Sep 17 00:00:00 2001 From: li97 Date: Tue, 24 Jun 2025 09:30:42 -0700 Subject: [PATCH 04/11] Fixed output size --- .../functional/element_restriction.cpp | 6 +-- src/serac/numerics/functional/functional.hpp | 54 ++++++++++++------- .../functional/tests/dg_ghost_face_index.cpp | 18 +++---- 3 files changed, 46 insertions(+), 32 deletions(-) diff --git a/src/serac/numerics/functional/element_restriction.cpp b/src/serac/numerics/functional/element_restriction.cpp index 3049b690d5..b78e5e350b 100644 --- a/src/serac/numerics/functional/element_restriction.cpp +++ b/src/serac/numerics/functional/element_restriction.cpp @@ -698,9 +698,6 @@ void ElementRestriction::Gather(const mfem::Vector& L_vector, mfem::Vector& E_ve for (uint64_t j = 0; j < nodes_per_elem; j++) { uint64_t E_id = (i * components + c) * nodes_per_elem + j; uint64_t L_id = GetVDof(dof_info(i, j), c).index(); - // // debug - // mpi::out << "elem #" << i << " " << "E_id = " << E_id << " " << dof_info(i, j).index() - // << " " << "L_id = " << L_id << std::endl; E_vector[int(E_id)] = L_vector[int(L_id)]; } } @@ -714,6 +711,9 @@ void ElementRestriction::ScatterAdd(const mfem::Vector& E_vector, mfem::Vector& for (uint64_t j = 0; j < nodes_per_elem; j++) { uint64_t E_id = (i * components + c) * nodes_per_elem + j; uint64_t L_id = GetVDof(dof_info(i, j), c).index(); + // // debug + // mpi::out << "elem #" << i << " " << "E_id = " << E_id << " " << dof_info(i, j).index() + // << " " << "L_id = " << L_id << std::endl; L_vector[int(L_id)] += E_vector[int(E_id)]; } } diff --git a/src/serac/numerics/functional/functional.hpp b/src/serac/numerics/functional/functional.hpp index b038462640..bfe4b1909c 100644 --- a/src/serac/numerics/functional/functional.hpp +++ b/src/serac/numerics/functional/functional.hpp @@ -276,9 +276,9 @@ class Functional { P_trial_[i] = trial_space_[i]->GetProlongationMatrix(); - // For L2 spaces, configure a ParGridFunction with external data to exchange face neighbor data - // This is for DG method where interior faces on the processor boundary need to access data on the neighrbor - // element + // For L2 spaces, configure a ParGridFunction with external data to exchange face neighbor data. + // This is for DG method where interior faces on the processor boundary need to access data on + // the neighrbor element if (trial_function_spaces_[i].family == Family::L2) { // Set the vector size to store tdof data input_ldof_values_[i].SetSize(P_trial_[i]->Height(), mem_type); @@ -303,7 +303,17 @@ class Functional { P_test_ = test_space_->GetProlongationMatrix(); - output_L_.SetSize(P_test_->Height(), mem_type); + if (test_function_space_.family == Family::L2) { + // We only need these variables to set the output_L_ vector to the right size + mfem::ParGridFunction X_pgf; + mfem::Vector X_ldof_values; + X_ldof_values.SetSize(P_test_->Height(), mem_type); + + updateFaceNbrData(test_space_, X_pgf, X_ldof_values); + output_L_.SetSize(X_ldof_values.Size() + X_pgf.FaceNbrData().Size(), mem_type); + } else { + output_L_.SetSize(P_test_->Height(), mem_type); + } output_T_.SetSize(test_fes->GetTrueVSize(), mem_type); @@ -516,19 +526,19 @@ class Functional { // ----------- Ghost elem 1 ---------- ----------- Ghost elem 2 ---------- // {X1 X1 X1 X1 Y1 Y1 Y1 Y1 Z1 Z1 Z1 Z1 X2 X2 X2 X2 Y2 Y2 Y2 Y2 Z2 Z2 Z2 Z2} // - // This is because 1. mfem prepares face neighbor data element by element to communicate later (standard) - // 2. mfem uses GetElementVDofs to gather the data for communication which returns a set of indices ALWAYS + // This is because 1) mfem prepares face neighbor data element by element to communicate later (standard). + // 2) mfem uses GetElementVDofs to gather the data for communication which returns a set of indices ALWAYS // ordered byNODES (???). If the ordering we set for ParFiniteElementSpace is byVDIM, we need the entries to be // ----------- Ghost elem 1 ---------- ----------- Ghost elem 2 ---------- // {X1 Y1 Z1 X1 Y1 Z1 X1 Y1 Z1 X1 Y1 Z1 X2 Y2 Z2 X2 Y2 Z2 X2 Y2 Z2 X2 Y2 Z2} - // so it's consistent with the local data. Therefore, we need to manually change this ordering + // so it's consistent with the local data. Therefore, we need to manually change this ordering. if (trial_space_[i]->GetVDim() == 1) { - // For scalar fields, this weird ordering doesn't cause any problems + // For scalar fields, this weird ordering doesn't cause any problems. input_L_[i].SetVector(trial_pargrid_functions_[i].FaceNbrData(), input_ldof_values_[i].Size()); } else { if (trial_space_[i]->GetOrdering() == mfem::Ordering::byVDIM) { // For vector fields with byVDIM ordering, we manually set the entries in VDIM order. By the way this - // is really annoying because we have to get face neighbor vdof indices again in element restriction + // is really annoying because we have to get face neighbor vdof indices again in element restriction. int num_ghost_elem = trial_space_[i]->GetParMesh()->GetNFaceNeighborElements(); int components_per_node = trial_space_[i]->GetVDim(); const mfem::Vector& face_neighbor_data = trial_pargrid_functions_[i].FaceNbrData(); @@ -552,12 +562,12 @@ class Functional { } } else { // For vector fields with byNODES ordering, we need to prepare the local input vector in the form - // { true_comp_X ghost_comp_X true_comp_Y ghost_comp_Y true_comp_Z ghost_comp_Z } + // { true_comp_X ghost_comp_X true_comp_Y ghost_comp_Y true_comp_Z ghost_comp_Z }. // In this form, to ensure the correct mapping between dof and vdof, we need to change the ndofs in // FiniteElementSpace to include ones from ghost element. Additionally, we need to prepare - // ghost_comp_X / ghost_comp_Y / ghost_comp_Z from element wise entries in face_nbr_data - // eg. ghost_comp_X includes X component entries of all ghost element dofs - // This requires significant renumbering of the dofs and therefore is not supported for now + // ghost_comp_X / ghost_comp_Y / ghost_comp_Z from element wise entries in face_nbr_data, + // eg. ghost_comp_X includes X component entries of all ghost element dofs. + // This requires significant renumbering of the dofs and therefore is not supported for now. SLIC_ERROR_ROOT("Unsupported: L2 vector field ordered by nodes"); } } @@ -585,24 +595,28 @@ class Functional { input_E_buffer_[i].SetSize(int(G_trial.ESize())); input_E_[i].Update(input_E_buffer_[i], G_trial.bOffsets()); G_trial.Gather(input_L_[i], input_E_[i]); - // // Debug - // std::cout << "element size = " << G_trial.ESize() << std::endl; - // std::cout << "Element vector" << std::endl; - // for (int j = 0; j < input_E_[i].Size(); ++j) { - // std::cout << input_E_[i][j] << " "; - // } - // std::cout << std::endl; } output_E_buffer_.SetSize(int(G_test.ESize())); output_E_.Update(output_E_buffer_, G_test.bOffsets()); integral.Mult(t, input_E_, output_E_, wrt, update_qdata_); + // // Debug + // std::cout << "ESize = " << G_test.ESize() << std::endl; + // std::cout << "Element output vector" << std::endl; + // for (int j = 0; j < output_E_.Size(); ++j) { + // std::cout << output_E_[j] << " "; + // } + // std::cout << std::endl; // scatter-add to compute residuals on the local processor G_test.ScatterAdd(output_E_, output_L_); } // scatter-add to compute global residuals + // TODO: might need to extract a subvector from output_L_ representing local residuals excluding + // ghost dofs for L2 spaces. Right now this is not raising any error because P_test_ is not a + // ConformingProlongationOperator for L2 spaces, so it does not check if the dimensions of the + // vectors matches the height and width of the matrix. P_test_->MultTranspose(output_L_, output_T_); if constexpr (wrt != NO_DIFFERENTIATION) { diff --git a/src/serac/numerics/functional/tests/dg_ghost_face_index.cpp b/src/serac/numerics/functional/tests/dg_ghost_face_index.cpp index 09ec22014a..852eb816cc 100644 --- a/src/serac/numerics/functional/tests/dg_ghost_face_index.cpp +++ b/src/serac/numerics/functional/tests/dg_ghost_face_index.cpp @@ -28,12 +28,12 @@ void L2_index_test(std::string meshfile) using test_space = L2; using trial_space = L2; - auto mesh = mesh::refineAndDistribute(buildMeshFromFile(meshfile), 0); + auto mesh = mesh::refineAndDistribute(buildMeshFromFile(meshfile), 1); - auto fec = mfem::L2_FECollection(p, dim, mfem::BasisType::GaussLobatto); - mfem::ParFiniteElementSpace fespace(mesh.get(), &fec, dim, serac::ordering); + auto [test_fespace, test_fec] = serac::generateParFiniteElementSpace(mesh.get()); + auto [trial_fespace, trial_fec] = serac::generateParFiniteElementSpace(mesh.get()); - mfem::ParGridFunction pgf(&fespace); + mfem::ParGridFunction U_gf(trial_fespace.get()); mfem::VectorFunctionCoefficient vcoef(dim, [](const mfem::Vector& X, mfem::Vector& F) { int d = X.Size(); F.SetSize(d); @@ -41,13 +41,13 @@ void L2_index_test(std::string meshfile) F(i) = X(i); } }); - pgf.ProjectCoefficient(vcoef); + U_gf.ProjectCoefficient(vcoef); - mfem::Vector U(fespace.TrueVSize()); - U = pgf; + mfem::Vector U(trial_fespace->TrueVSize()); + U_gf.GetTrueDofs(U); // Construct the new functional object using the specified test and trial spaces - Functional residual(&fespace, {&fespace}); + Functional residual(test_fespace.get(), {trial_fespace.get()}); Domain interior_faces = InteriorFaces(*mesh); @@ -62,7 +62,7 @@ void L2_index_test(std::string meshfile) // note: the orientation convention is such that the normal // computed as above will point from from side 1->2 auto [u_1, u_2] = velocity; - std::cout << "One side = " << u_1 << " The other side = " << u_2 << " Jump = " << u_1 - u_2 << std::endl; + std::cout << "One side = " << u_1 << ", The other side = " << u_2 << ", Jump = " << u_1 - u_2 << std::endl; auto a = dot(u_2 - u_1, n); From abf5111000e865e48e06fee0053a17cbbb6a6879 Mon Sep 17 00:00:00 2001 From: li97 Date: Tue, 24 Jun 2025 21:17:50 -0700 Subject: [PATCH 05/11] Fixed true residual size --- src/serac/numerics/functional/functional.hpp | 19 ++++++++++++------- 1 file changed, 12 insertions(+), 7 deletions(-) diff --git a/src/serac/numerics/functional/functional.hpp b/src/serac/numerics/functional/functional.hpp index bfe4b1909c..c8883b3319 100644 --- a/src/serac/numerics/functional/functional.hpp +++ b/src/serac/numerics/functional/functional.hpp @@ -306,8 +306,7 @@ class Functional { if (test_function_space_.family == Family::L2) { // We only need these variables to set the output_L_ vector to the right size mfem::ParGridFunction X_pgf; - mfem::Vector X_ldof_values; - X_ldof_values.SetSize(P_test_->Height(), mem_type); + mfem::Vector X_ldof_values(P_test_->Height(), mem_type); updateFaceNbrData(test_space_, X_pgf, X_ldof_values); output_L_.SetSize(X_ldof_values.Size() + X_pgf.FaceNbrData().Size(), mem_type); @@ -613,11 +612,17 @@ class Functional { } // scatter-add to compute global residuals - // TODO: might need to extract a subvector from output_L_ representing local residuals excluding - // ghost dofs for L2 spaces. Right now this is not raising any error because P_test_ is not a - // ConformingProlongationOperator for L2 spaces, so it does not check if the dimensions of the - // vectors matches the height and width of the matrix. - P_test_->MultTranspose(output_L_, output_T_); + if (test_function_space_.family == Family::L2) { + // Extract a subvector from output_L_ for local residuals excluding ghost dofs for L2 spaces + // The output vector is in form [ --- L --- | --- FND --- ] + mfem::Vector output_ldof_values(P_test_->Height(), mem_type); + for (int j = 0; j < output_ldof_values.Size(); ++j) { + output_ldof_values[j] = output_L_[j]; + } + P_test_->MultTranspose(output_ldof_values, output_T_); + } else { + P_test_->MultTranspose(output_L_, output_T_); + } if constexpr (wrt != NO_DIFFERENTIATION) { // if the user has indicated they'd like to evaluate and differentiate w.r.t. From 42804174f653ba821530ff9d09210a401aa283e6 Mon Sep 17 00:00:00 2001 From: li97 Date: Wed, 25 Jun 2025 15:37:44 -0700 Subject: [PATCH 06/11] Refactor ghost data reordering. Fix ActionOfGradient --- src/serac/numerics/functional/functional.hpp | 90 ++++++++++++++------ 1 file changed, 65 insertions(+), 25 deletions(-) diff --git a/src/serac/numerics/functional/functional.hpp b/src/serac/numerics/functional/functional.hpp index c8883b3319..786f43fbb9 100644 --- a/src/serac/numerics/functional/functional.hpp +++ b/src/serac/numerics/functional/functional.hpp @@ -184,6 +184,33 @@ inline void updateFaceNbrData(const mfem::ParFiniteElementSpace* const_trial_spa trial_pgf.ExchangeFaceNbrData(); } +/** + * @brief helper functional to rearrange the ordering of FaceNbrData for L2 space to byVDIM + */ +inline void reorderFaceNbrData(const mfem::ParFiniteElementSpace* trial_space, const mfem::ParGridFunction& trial_pgf, + const int LSize, mfem::Vector& input_L) +{ + int num_ghost_elem = trial_space->GetParMesh()->GetNFaceNeighborElements(); + int components_per_node = trial_space->GetVDim(); + const mfem::Vector& face_neighbor_data = trial_pgf.FaceNbrData(); + + int offset = 0; + for (int n = 0; n < num_ghost_elem; ++n) { + mfem::Array old_shared_elem_vdof_ids; + trial_space->GetFaceNbrElementVDofs(n, old_shared_elem_vdof_ids); + int dofs_per_element = old_shared_elem_vdof_ids.Size() / components_per_node; + + // Put the entries into their new index + for (int m = 0; m < old_shared_elem_vdof_ids.Size(); ++m) { + // -------------- dofs before -------------- ---- component ----- + int new_id = (m % dofs_per_element) * components_per_node + m / dofs_per_element + offset; + input_L[LSize + new_id] = face_neighbor_data(old_shared_elem_vdof_ids[m]); + } + // Increase offset for next ghost element + offset += old_shared_elem_vdof_ids.Size(); + } +} + /// @cond template class Functional; @@ -459,7 +486,29 @@ class Functional { */ void ActionOfGradient(const mfem::Vector& input_T, mfem::Vector& output_T, uint32_t which) const { - P_trial_[which]->Mult(input_T, input_L_[which]); + // Please refer to 'operator()' below for detailed comments of this implementation + if (trial_function_spaces_[which].family == Family::L2) { + // copy input_L[which] and facenbrdata[which] into a common array like: + // first part second part + // [ --- L --- | --- FND --- ] + P_trial_[which]->Mult(input_T, input_ldof_values_[which]); + input_L_[which].SetVector(input_ldof_values_[which], 0); + updateFaceNbrData(trial_space_[which], trial_pargrid_functions_[which], input_ldof_values_[which]); + + // Ensure ghost data ordering is consistent with local data + if (trial_space_[which]->GetVDim() == 1) { + input_L_[which].SetVector(trial_pargrid_functions_[which].FaceNbrData(), input_ldof_values_[which].Size()); + } else { + if (trial_space_[which]->GetOrdering() == mfem::Ordering::byVDIM) { + int local_size = input_ldof_values_[which].Size(); + reorderFaceNbrData(trial_space_[which], trial_pargrid_functions_[which], local_size, input_L_[which]); + } else { + SLIC_ERROR_ROOT("Unsupported: L2 vector field ordered by nodes"); + } + } + } else { + P_trial_[which]->Mult(input_T, input_L_[which]); + } output_L_ = 0.0; @@ -484,7 +533,17 @@ class Functional { } // scatter-add to compute global residuals - P_test_->MultTranspose(output_L_, output_T); + if (test_function_space_.family == Family::L2) { + // Extract a subvector from output_L_ for local residuals excluding ghost dofs for L2 spaces + // The output vector is in form [ --- L --- | --- FND --- ] + mfem::Vector output_ldof_values(P_test_->Height(), mem_type); + for (int j = 0; j < output_ldof_values.Size(); ++j) { + output_ldof_values[j] = output_L_[j]; + } + P_test_->MultTranspose(output_ldof_values, output_T); + } else { + P_test_->MultTranspose(output_L_, output_T); + } } /** @@ -507,7 +566,7 @@ class Functional { // get the values for each local processor for (uint32_t i = 0; i < num_trial_spaces; i++) { if (trial_function_spaces_[i].family == Family::L2) { - // copy input_L[i] and facenbrdata [i] into a common array like: + // copy input_L[i] and facenbrdata[i] into a common array like: // first part second part // [ --- L --- | --- FND --- ] P_trial_[i]->Mult(*input_T[i], input_ldof_values_[i]); @@ -538,27 +597,8 @@ class Functional { if (trial_space_[i]->GetOrdering() == mfem::Ordering::byVDIM) { // For vector fields with byVDIM ordering, we manually set the entries in VDIM order. By the way this // is really annoying because we have to get face neighbor vdof indices again in element restriction. - int num_ghost_elem = trial_space_[i]->GetParMesh()->GetNFaceNeighborElements(); - int components_per_node = trial_space_[i]->GetVDim(); - const mfem::Vector& face_neighbor_data = trial_pargrid_functions_[i].FaceNbrData(); - - int offset = 0; int local_size = input_ldof_values_[i].Size(); - - for (int n = 0; n < num_ghost_elem; ++n) { - mfem::Array old_shared_elem_vdof_ids; - trial_space_[i]->GetFaceNbrElementVDofs(n, old_shared_elem_vdof_ids); - int dofs_per_element = old_shared_elem_vdof_ids.Size() / components_per_node; - - // Put the entries into their new index - for (int m = 0; m < old_shared_elem_vdof_ids.Size(); ++m) { - // -------------- dofs before -------------- ---- component ----- - int new_id = (m % dofs_per_element) * components_per_node + m / dofs_per_element + offset; - input_L_[i][local_size + new_id] = face_neighbor_data(old_shared_elem_vdof_ids[m]); - } - // Increase offset for next ghost element - offset += old_shared_elem_vdof_ids.Size(); - } + reorderFaceNbrData(trial_space_[i], trial_pargrid_functions_[i], local_size, input_L_[i]); } else { // For vector fields with byNODES ordering, we need to prepare the local input vector in the form // { true_comp_X ghost_comp_X true_comp_Y ghost_comp_Y true_comp_Z ghost_comp_Z }. @@ -925,10 +965,10 @@ class Functional { FunctionSpace test_function_space_; /// @brief For access of neighbor element info for processor boundary faces - std::array trial_pargrid_functions_; + mutable std::array trial_pargrid_functions_; /// @brief Store values on L2 ldofs of the local processor - std::array input_ldof_values_; + mutable std::array input_ldof_values_; /** * @brief Operator that converts true (global) DOF values to local (current rank) DOF values From f7a0b942bd215380e2f2874b81eac93689093210 Mon Sep 17 00:00:00 2001 From: li97 Date: Wed, 9 Jul 2025 16:03:36 -0700 Subject: [PATCH 07/11] fixed assemble, basic test passed in parallel --- .../functional/element_restriction.cpp | 8 - src/serac/numerics/functional/functional.hpp | 154 +++++++++++++++--- .../numerics/functional/tests/CMakeLists.txt | 2 +- .../functional/tests/dg_ghost_face_index.cpp | 2 +- .../tests/dg_restriction_operators.cpp | 5 - .../functional/tests/functional_basic_dg.cpp | 4 +- 6 files changed, 134 insertions(+), 41 deletions(-) diff --git a/src/serac/numerics/functional/element_restriction.cpp b/src/serac/numerics/functional/element_restriction.cpp index b78e5e350b..3fdfc14a59 100644 --- a/src/serac/numerics/functional/element_restriction.cpp +++ b/src/serac/numerics/functional/element_restriction.cpp @@ -519,11 +519,6 @@ axom::Array GetFaceDofs(const serac::f // 3. get the dofs for the entire element mfem::Array elem_dof_ids; fes->GetElementDofs(elem, elem_dof_ids); - // // Debug hanyu - // mpi::out << " Shared face = " << info.IsShared() << " Face id = " << f << " Elem id = " << elem << std::endl - // << i << " " << info.element[0].local_face_id << " " << info.element[1].local_face_id << std::endl - // << orientations[i] << " " << info.element[0].orientation << " " << info.element[1].orientation << - // std::endl; mfem::Geometry::Type elem_geom = mesh->GetElementGeometry(elem); @@ -711,9 +706,6 @@ void ElementRestriction::ScatterAdd(const mfem::Vector& E_vector, mfem::Vector& for (uint64_t j = 0; j < nodes_per_elem; j++) { uint64_t E_id = (i * components + c) * nodes_per_elem + j; uint64_t L_id = GetVDof(dof_info(i, j), c).index(); - // // debug - // mpi::out << "elem #" << i << " " << "E_id = " << E_id << " " << dof_info(i, j).index() - // << " " << "L_id = " << L_id << std::endl; L_vector[int(L_id)] += E_vector[int(E_id)]; } } diff --git a/src/serac/numerics/functional/functional.hpp b/src/serac/numerics/functional/functional.hpp index 786f43fbb9..a5e4f38f18 100644 --- a/src/serac/numerics/functional/functional.hpp +++ b/src/serac/numerics/functional/functional.hpp @@ -29,6 +29,9 @@ #include #include +// TODO REMOVE AFTER DEBUGGING +#include "serac/infrastructure/mpi_fstream.hpp" + namespace serac { /// @cond @@ -185,10 +188,11 @@ inline void updateFaceNbrData(const mfem::ParFiniteElementSpace* const_trial_spa } /** - * @brief helper functional to rearrange the ordering of FaceNbrData for L2 space to byVDIM + * @brief helper functional to reorder the ordering of FaceNbrData for L2 space to byVDIM and append this vector + * to the end of local dof vector, which will result in a vector in form [ --- L --- | --- FND --- ] */ -inline void reorderFaceNbrData(const mfem::ParFiniteElementSpace* trial_space, const mfem::ParGridFunction& trial_pgf, - const int LSize, mfem::Vector& input_L) +inline void appendFaceNbrData(const mfem::ParFiniteElementSpace* trial_space, const mfem::ParGridFunction& trial_pgf, + const int LSize, mfem::Vector& input_L) { int num_ghost_elem = trial_space->GetParMesh()->GetNFaceNeighborElements(); int components_per_node = trial_space->GetVDim(); @@ -211,6 +215,42 @@ inline void reorderFaceNbrData(const mfem::ParFiniteElementSpace* trial_space, c } } +/** + * @brief helper functional to reorder the face_nbr_glob_dof_map for L2 space to byVDIM + */ +inline void rearrangeFaceNbrDofGlobalIndex(const mfem::ParFiniteElementSpace* trial_space, + mfem::Array& face_nbr_glob_vdof_map) +{ + int components_per_node = trial_space->GetVDim(); + if (components_per_node == 1) { + face_nbr_glob_vdof_map = trial_space->face_nbr_glob_dof_map; + return; + } + + if (components_per_node > 1 && trial_space->GetOrdering() == mfem::Ordering::byNODES) { + SLIC_ERROR_ROOT("Unsupported: L2 vector field ordered by nodes"); + } + + face_nbr_glob_vdof_map.SetSize(trial_space->face_nbr_glob_dof_map.Size()); + int num_ghost_elem = trial_space->GetParMesh()->GetNFaceNeighborElements(); + const HYPRE_BigInt* face_nbr_glob_dof_index = trial_space->face_nbr_glob_dof_map; + + int offset = 0; + for (int n = 0; n < num_ghost_elem; ++n) { + mfem::Array old_shared_elem_vdof_ids; + trial_space->GetFaceNbrElementVDofs(n, old_shared_elem_vdof_ids); + int dofs_per_element = old_shared_elem_vdof_ids.Size() / components_per_node; + + // Put the entries into their new index + for (int m = 0; m < old_shared_elem_vdof_ids.Size(); ++m) { + int new_id = (m % dofs_per_element) * components_per_node + m / dofs_per_element + offset; + face_nbr_glob_vdof_map[new_id] = face_nbr_glob_dof_index[old_shared_elem_vdof_ids[m]]; + } + // Increase offset for next ghost element + offset += old_shared_elem_vdof_ids.Size(); + } +} + /// @cond template class Functional; @@ -315,6 +355,10 @@ class Functional { // Set the local input vector size to store local + ghost data input_L_[i].SetSize(input_ldof_values_[i].Size() + trial_pargrid_functions_[i].FaceNbrData().Size(), mem_type); + + // Rearrange the global index of face neighbor vdofs to the correct ordering. + // This only needs to be done once + rearrangeFaceNbrDofGlobalIndex(trial_space_[i], face_nbr_glob_vdof_maps_[i]); } else { input_L_[i].SetSize(P_trial_[i]->Height(), mem_type); } @@ -501,7 +545,7 @@ class Functional { } else { if (trial_space_[which]->GetOrdering() == mfem::Ordering::byVDIM) { int local_size = input_ldof_values_[which].Size(); - reorderFaceNbrData(trial_space_[which], trial_pargrid_functions_[which], local_size, input_L_[which]); + appendFaceNbrData(trial_space_[which], trial_pargrid_functions_[which], local_size, input_L_[which]); } else { SLIC_ERROR_ROOT("Unsupported: L2 vector field ordered by nodes"); } @@ -598,7 +642,7 @@ class Functional { // For vector fields with byVDIM ordering, we manually set the entries in VDIM order. By the way this // is really annoying because we have to get face neighbor vdof indices again in element restriction. int local_size = input_ldof_values_[i].Size(); - reorderFaceNbrData(trial_space_[i], trial_pargrid_functions_[i], local_size, input_L_[i]); + appendFaceNbrData(trial_space_[i], trial_pargrid_functions_[i], local_size, input_L_[i]); } else { // For vector fields with byNODES ordering, we need to prepare the local input vector in the form // { true_comp_X ghost_comp_X true_comp_Y ghost_comp_Y true_comp_Z ghost_comp_Z }. @@ -613,13 +657,6 @@ class Functional { } else { P_trial_[i]->Mult(*input_T[i], input_L_[i]); } - // // Debug - // std::cout << "local size = " << input_L_[i].Size() << std::endl; - // std::cout << "Local vector" << std::endl; - // for (int j = 0; j < input_L_[i].Size(); ++j) { - // std::cout << input_L_[i][j] << " "; - // } - // std::cout << std::endl; } output_L_ = 0.0; @@ -639,13 +676,6 @@ class Functional { output_E_buffer_.SetSize(int(G_test.ESize())); output_E_.Update(output_E_buffer_, G_test.bOffsets()); integral.Mult(t, input_E_, output_E_, wrt, update_qdata_); - // // Debug - // std::cout << "ESize = " << G_test.ESize() << std::endl; - // std::cout << "Element output vector" << std::endl; - // for (int j = 0; j < output_E_.Size(); ++j) { - // std::cout << output_E_[j] << " "; - // } - // std::cout << std::endl; // scatter-add to compute residuals on the local processor G_test.ScatterAdd(output_E_, output_L_); @@ -844,6 +874,14 @@ class Functional { if (row_ptr.empty()) { initialize_sparsity_pattern(); } + // // Debug + // mpi::out << "row offsets: " << (test_space_->GetDofOffsets())[0] << " " + // << (test_space_->GetDofOffsets())[1] << std::endl + // << "col offsets: " << (trial_space_->GetDofOffsets())[0] << " " + // << (trial_space_->GetDofOffsets())[1] << std::endl; + // HYPRE_BigInt row_offset = test_space_->GetMyDofOffset(); + // HYPRE_BigInt col_offset = trial_space_->GetMyDofOffset(); + // // // since we own the storage for row_ptr, col_ind, values, // we ask mfem to not deallocate those pointers in the SparseMatrix dtor @@ -904,6 +942,29 @@ class Functional { A_local.SearchRow(row, col) += K_e(e, i, j); } } + // // Debug + // for (uint32_t i = 0; i < rows_per_elem; i++) { + // int row = int(test_vdofs[i].index()); + // for (uint32_t j = 0; j < cols_per_elem; j++) { + // int col = int(trial_vdofs[j].index()); + // if (row < test_space_->GetVSize()) { + // mpi::out << "row = " << row + row_offset << " "; + // } else { + // mpi::out << "row = " << form_.face_nbr_glob_vdof_maps_[which_argument][row - + // test_space_->GetVSize()] << " "; + // } + + // if (col < trial_space_->GetVSize()) { + // mpi::out << "col = " << col + col_offset << " "; + // } else { + // mpi::out << "col = " << form_.face_nbr_glob_vdof_maps_[which_argument][col - + // trial_space_->GetVSize()] << " "; + // } + + // mpi::out << "A_ij = " << A_local.SearchRow(row, col) << std::endl; + // } + // } + // // } } } @@ -911,9 +972,51 @@ class Functional { auto* R = form_.test_space_->Dof_TrueDof_Matrix(); - auto* A_hypre = - new mfem::HypreParMatrix(test_space_->GetComm(), test_space_->GlobalVSize(), trial_space_->GlobalVSize(), - test_space_->GetDofOffsets(), trial_space_->GetDofOffsets(), &A_local); + // If either test_space_ or trial_space_ is L2, the local matrix (A_local) will includ ghost rows and columns like + // ------------ ------- + // | || | + // | diagonal || off-d | + // | block || block | + // | || | + // ------------ ------- + // ------------ ------- + // | ghost || XXX | + // | rows || XXX | + // ------------ ------- + // We only need + // ------------ ------- + // | || | + // | diagonal || off-d | + // | block || block | + // | || | + // ------------ ------- + // to compute local residual and the neighbor dof values will be communicated to multiply into off-digonal block. + // So we construct a HypreParMatrix that contains the off-diagonal block from the local sparse matrix + mfem::HypreParMatrix* A_hypre; + if (dynamic_cast(test_space_->FEColl()) || + dynamic_cast(trial_space_->FEColl())) { + int lrows = test_space_->GetVSize(); + int lcols = trial_space_->GetVSize(); + HYPRE_BigInt col_offset = trial_space_->GetMyDofOffset(); + mfem::Array glob_J(A_local.NumNonZeroElems()); + + int* J = A_local.GetJ(); + for (int j = 0; j < glob_J.Size(); ++j) { + if (J[j] < lcols) { + glob_J[j] = J[j] + col_offset; + } else { + glob_J[j] = form_.face_nbr_glob_vdof_maps_[which_argument][J[j] - lcols]; + } + } + A_hypre = new mfem::HypreParMatrix(test_space_->GetComm(), lrows, test_space_->GlobalVSize(), + trial_space_->GlobalVSize(), A_local.GetI(), glob_J, A_local.GetData(), + test_space_->GetDofOffsets(), trial_space_->GetDofOffsets()); + glob_J.DeleteAll(); + } else { + A_hypre = + new mfem::HypreParMatrix(test_space_->GetComm(), test_space_->GlobalVSize(), trial_space_->GlobalVSize(), + test_space_->GetDofOffsets(), trial_space_->GetDofOffsets(), &A_local); + } auto* P = trial_space_->Dof_TrueDof_Matrix(); @@ -964,10 +1067,13 @@ class Functional { std::array trial_function_spaces_; FunctionSpace test_function_space_; - /// @brief For access of neighbor element info for processor boundary faces + /// @brief Manage global index of face neighbor vdofs + std::array, num_trial_spaces> face_nbr_glob_vdof_maps_; + + /// @brief Cache for access of neighbor element info for processor boundary faces mutable std::array trial_pargrid_functions_; - /// @brief Store values on L2 ldofs of the local processor + /// @brief Cache to store values of L2 ldofs on the local processor mutable std::array input_ldof_values_; /** diff --git a/src/serac/numerics/functional/tests/CMakeLists.txt b/src/serac/numerics/functional/tests/CMakeLists.txt index 1357e65a1e..cf101153a4 100644 --- a/src/serac/numerics/functional/tests/CMakeLists.txt +++ b/src/serac/numerics/functional/tests/CMakeLists.txt @@ -20,7 +20,6 @@ set(functional_serial_test_sources simplex_basis_function_unit_tests.cpp element_restriction_tests.cpp dg_restriction_operators.cpp - functional_basic_dg.cpp bug_boundary_qoi.cpp domain_tests.cpp geometric_factors_tests.cpp @@ -41,6 +40,7 @@ set(functional_parallel_test_sources functional_with_domain.cpp functional_basic_h1_scalar.cpp functional_basic_h1_vector.cpp + functional_basic_dg.cpp functional_multiphysics.cpp functional_qoi.cpp functional_nonlinear.cpp diff --git a/src/serac/numerics/functional/tests/dg_ghost_face_index.cpp b/src/serac/numerics/functional/tests/dg_ghost_face_index.cpp index 852eb816cc..e60bec7258 100644 --- a/src/serac/numerics/functional/tests/dg_ghost_face_index.cpp +++ b/src/serac/numerics/functional/tests/dg_ghost_face_index.cpp @@ -77,7 +77,7 @@ void L2_index_test(std::string meshfile) auto value = residual(t, U); EXPECT_NEAR(0., value.Norml2(), 1.e-12); } -// TEST(index, L2_two_quad) { L2_index_test<2, 1>("./TwoQuad.mesh"); } + TEST(index, L2_test_tris_and_quads_linear) { L2_index_test<2, 1>(SERAC_REPO_DIR "/data/meshes/patch2D_tris_and_quads.mesh"); diff --git a/src/serac/numerics/functional/tests/dg_restriction_operators.cpp b/src/serac/numerics/functional/tests/dg_restriction_operators.cpp index 518dd92cfc..e3762202c1 100644 --- a/src/serac/numerics/functional/tests/dg_restriction_operators.cpp +++ b/src/serac/numerics/functional/tests/dg_restriction_operators.cpp @@ -353,11 +353,6 @@ void parametrized_test(int permutation) using space = L2; - // auto pmesh = mesh::refineAndDistribute(std::move(mesh), 0); - // Domain pinterior_faces = InteriorFaces(*pmesh); - - // auto L2_pfes = mfem::ParFiniteElementSpace(pmesh.get(), L2_fec.get()); - Functional r({L2_fes.get()}); r.AddInteriorFaceIntegral( diff --git a/src/serac/numerics/functional/tests/functional_basic_dg.cpp b/src/serac/numerics/functional/tests/functional_basic_dg.cpp index 927f9bcd04..de4ba2eedb 100644 --- a/src/serac/numerics/functional/tests/functional_basic_dg.cpp +++ b/src/serac/numerics/functional/tests/functional_basic_dg.cpp @@ -33,7 +33,7 @@ void L2_test(std::string meshfile) // int k = 0; // while (k == 0); - auto mesh = mesh::refineAndDistribute(buildMeshFromFile(meshfile), 0); + auto mesh = mesh::refineAndDistribute(buildMeshFromFile(meshfile), 1); auto fec = mfem::L2_FECollection(p, dim, mfem::BasisType::GaussLobatto); mfem::ParFiniteElementSpace fespace(mesh.get(), &fec, dim, serac::ordering); @@ -72,7 +72,7 @@ void L2_test(std::string meshfile) // auto value = residual(t, U); check_gradient(residual, t, U); } - +// TEST(basic, L2_two_quad) { L2_test<2, 1>("./TwoQuad.mesh"); } TEST(basic, L2_test_tris_and_quads_linear) { L2_test<2, 1>(SERAC_REPO_DIR "/data/meshes/patch2D_tris_and_quads.mesh"); } TEST(basic, L2_test_tris_and_quads_quadratic) { From b9c59586d2bbdb339f3c14e1cf95c9e15c7e60a4 Mon Sep 17 00:00:00 2001 From: li97 Date: Fri, 11 Jul 2025 16:03:53 -0700 Subject: [PATCH 08/11] dg parallelized --- src/serac/numerics/functional/domain.cpp | 53 ++++++++++ src/serac/numerics/functional/domain.hpp | 9 ++ .../functional/element_restriction.cpp | 9 -- src/serac/numerics/functional/functional.hpp | 55 ----------- .../numerics/functional/functional_qoi.inl | 98 ++++++++++++++++--- .../functional/tests/functional_basic_dg.cpp | 12 ++- 6 files changed, 158 insertions(+), 78 deletions(-) diff --git a/src/serac/numerics/functional/domain.cpp b/src/serac/numerics/functional/domain.cpp index dca0342319..ccb2d50b5c 100644 --- a/src/serac/numerics/functional/domain.cpp +++ b/src/serac/numerics/functional/domain.cpp @@ -535,6 +535,59 @@ void Domain::insert_restriction(const serac::fes_t* fes, FunctionSpace space) const BlockElementRestriction& Domain::get_restriction(FunctionSpace space) { return restriction_operators.at(space); }; +void Domain::compute_interior_face_qoi_weights() +{ + // Weights only need to be computed for Domain of InteriorFaces type + SLIC_ERROR_IF(type_ != Domain::Type::InteriorFaces, "This method is only for interior face domains"); + + // Weights only need to be computed once for the same domain + if (interior_face_qoi_weights.Size() == 0) { + if (dim_ == 1) { + interior_face_qoi_weights.SetSize(int(edge_ids_.size())); + + int i = 0; + for (int f : mfem_edge_ids_) { + mfem::Mesh::FaceInformation info = mesh_.GetFaceInformation(f); + + if (info.IsShared()) { + interior_face_qoi_weights[i] = 0.5; + } else { + interior_face_qoi_weights[i] = 1.0; + } + + ++i; + } + } else if (dim_ == 2) { + interior_face_qoi_weights.SetSize(int(tri_ids_.size() + quad_ids_.size())); + + int i = 0; + for (int f : mfem_tri_ids_) { + mfem::Mesh::FaceInformation info = mesh_.GetFaceInformation(f); + + if (info.IsShared()) { + interior_face_qoi_weights[i] = 0.5; + } else { + interior_face_qoi_weights[i] = 1.0; + } + + ++i; + } + + for (int f : mfem_quad_ids_) { + mfem::Mesh::FaceInformation info = mesh_.GetFaceInformation(f); + + if (info.IsShared()) { + interior_face_qoi_weights[i] = 0.5; + } else { + interior_face_qoi_weights[i] = 1.0; + } + + ++i; + } + } + } +} + /////////////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////////////// diff --git a/src/serac/numerics/functional/domain.hpp b/src/serac/numerics/functional/domain.hpp index cf5171b7ba..2e02469d53 100644 --- a/src/serac/numerics/functional/domain.hpp +++ b/src/serac/numerics/functional/domain.hpp @@ -244,6 +244,15 @@ struct Domain { /// thereby populate the element lists. void addElements(const std::vector& geom_id, const std::vector& elem_id, mfem::Geometry::Type element_geometry); + + /** + * @brief Compute face qoi weights so that interior faces shared by two processors + * are only integrated once. + */ + void compute_interior_face_qoi_weights(); + + /// @brief Interior face qoi weight vector + mfem::Vector interior_face_qoi_weights; }; /// @brief constructs a domain from all the elements in a mesh diff --git a/src/serac/numerics/functional/element_restriction.cpp b/src/serac/numerics/functional/element_restriction.cpp index 3fdfc14a59..8a5741933f 100644 --- a/src/serac/numerics/functional/element_restriction.cpp +++ b/src/serac/numerics/functional/element_restriction.cpp @@ -4,9 +4,6 @@ #include "serac/numerics/functional/geometry.hpp" -// TODO REMOVE AFTER DEBUGGING -#include "serac/infrastructure/mpi_fstream.hpp" - std::vector > lexicographic_permutations(int p) { // p == 0 is admissible for L2 spaces, but lexicographic permutations @@ -531,12 +528,9 @@ axom::Array GetFaceDofs(const serac::f int orientation = (mesh->Dimension() == 2 && info.IsBoundary()) ? 0 : orientations[i]; // 4. extract only the dofs that correspond to side `i` - // mpi::out << "face local dofs: "; for (auto k : face_perm(orientation)) { face_dofs.push_back(uint64_t(elem_dof_ids[local_face_dofs[uint32_t(elem_geom)](i, k)])); - // mpi::out << face_dofs.back().index() << " "; } - // mpi::out << std::endl; // 5. add remaining dofs that were omitted on shared faces // FaceNbrData includes all dofs in the "volumetric" element adjacent to the shared faces @@ -575,7 +569,6 @@ axom::Array GetFaceDofs(const serac::f // Find the dofs on the shared face on the ghost element side // Neighbor element local_face_index and orientation is already available from GetFaceInformation - // mpi::out << " this face is also shared so it has additional dofs: "; mfem::Geometry::Type ghost_geom = fes->GetParMesh()->face_nbr_elements[other_element_id]->GetGeometryType(); int j = info.element[1].local_face_id; int ghost_orientation; @@ -596,9 +589,7 @@ axom::Array GetFaceDofs(const serac::f for (auto k : face_perm(ghost_orientation)) { face_dofs.push_back(uint64_t(shared_elem_dof_ids[local_face_dofs[uint32_t(ghost_geom)](j, k)] + LSize / components_per_node)); - // mpi::out << face_dofs.back().index() << " "; } - // mpi::out << std::endl; } } diff --git a/src/serac/numerics/functional/functional.hpp b/src/serac/numerics/functional/functional.hpp index a5e4f38f18..7a8101c2ab 100644 --- a/src/serac/numerics/functional/functional.hpp +++ b/src/serac/numerics/functional/functional.hpp @@ -29,9 +29,6 @@ #include #include -// TODO REMOVE AFTER DEBUGGING -#include "serac/infrastructure/mpi_fstream.hpp" - namespace serac { /// @cond @@ -109,25 +106,6 @@ inline void check_for_unsupported_elements(const mfem::Mesh& mesh) } } -/** - * @brief function for verifying that DG spaces aren't used on interior face integrals over meshes that contain "shared" - * faces - * - * sam: I would like to support these "shared" faces, but apparently mfem handles them in a fundamentally different way - * than the "finite element operator decomposition" pattern used by everything else (see: "ExchangeFaceNbrData") - */ -inline void check_interior_face_compatibility(const mfem::Mesh& mesh, const FunctionSpace space) -{ - if (space.family == Family::L2) { - const mfem::ParMesh* pmesh = dynamic_cast(&mesh); - if (pmesh) { - SLIC_ERROR_IF( - pmesh->GetNSharedFaces() > 0, - "interior face integrals involving DG function spaces don't currently support meshes with shared faces"); - } - } -} - /** * @brief create an mfem::ParFiniteElementSpace from one of serac's * tag types: H1, Hcurl, L2 @@ -468,10 +446,8 @@ class Functional { std::vector arg_vec = {args...}; for (uint32_t i : arg_vec) { domain.insert_restriction(trial_space_[i], trial_function_spaces_[i]); - // check_interior_face_compatibility(domain.mesh_, trial_function_spaces_[i]); } domain.insert_restriction(test_space_, test_function_space_); - // check_interior_face_compatibility(domain.mesh_, test_function_space_); using signature = test(decltype(serac::type(trial_spaces))...); integrals_.push_back( @@ -874,14 +850,6 @@ class Functional { if (row_ptr.empty()) { initialize_sparsity_pattern(); } - // // Debug - // mpi::out << "row offsets: " << (test_space_->GetDofOffsets())[0] << " " - // << (test_space_->GetDofOffsets())[1] << std::endl - // << "col offsets: " << (trial_space_->GetDofOffsets())[0] << " " - // << (trial_space_->GetDofOffsets())[1] << std::endl; - // HYPRE_BigInt row_offset = test_space_->GetMyDofOffset(); - // HYPRE_BigInt col_offset = trial_space_->GetMyDofOffset(); - // // // since we own the storage for row_ptr, col_ind, values, // we ask mfem to not deallocate those pointers in the SparseMatrix dtor @@ -942,29 +910,6 @@ class Functional { A_local.SearchRow(row, col) += K_e(e, i, j); } } - // // Debug - // for (uint32_t i = 0; i < rows_per_elem; i++) { - // int row = int(test_vdofs[i].index()); - // for (uint32_t j = 0; j < cols_per_elem; j++) { - // int col = int(trial_vdofs[j].index()); - // if (row < test_space_->GetVSize()) { - // mpi::out << "row = " << row + row_offset << " "; - // } else { - // mpi::out << "row = " << form_.face_nbr_glob_vdof_maps_[which_argument][row - - // test_space_->GetVSize()] << " "; - // } - - // if (col < trial_space_->GetVSize()) { - // mpi::out << "col = " << col + col_offset << " "; - // } else { - // mpi::out << "col = " << form_.face_nbr_glob_vdof_maps_[which_argument][col - - // trial_space_->GetVSize()] << " "; - // } - - // mpi::out << "A_ij = " << A_local.SearchRow(row, col) << std::endl; - // } - // } - // // } } } diff --git a/src/serac/numerics/functional/functional_qoi.inl b/src/serac/numerics/functional/functional_qoi.inl index 76cfa97257..e53bcc3f06 100644 --- a/src/serac/numerics/functional/functional_qoi.inl +++ b/src/serac/numerics/functional/functional_qoi.inl @@ -88,29 +88,42 @@ class Functional { { SERAC_MARK_FUNCTION; + std::array trial_families = {trials::family...}; + std::array trial_orders = {trials::order...}; + std::array trial_components = {trials::components...}; + for (uint32_t i = 0; i < num_trial_spaces; i++) { + trial_function_spaces_[i] = {trial_families[i], trial_orders[i], trial_components[i]}; + P_trial_[i] = trial_space_[i]->GetProlongationMatrix(); - input_L_[i].SetSize(P_trial_[i]->Height(), mfem::Device::GetMemoryType()); + // For L2 spaces, configure a ParGridFunction with external data to exchange face neighbor data. + // This is for DG method where interior faces on the processor boundary need to access data on + // the neighrbor element + if (trial_function_spaces_[i].family == Family::L2) { + // Set the vector size to store tdof data + input_ldof_values_[i].SetSize(P_trial_[i]->Height(), mem_type); + + // Initialize face neighbor data vector + updateFaceNbrData(trial_space_[i], trial_pargrid_functions_[i], input_ldof_values_[i]); + + // Set the local input vector size to store local + ghost data + input_L_[i].SetSize(input_ldof_values_[i].Size() + trial_pargrid_functions_[i].FaceNbrData().Size(), mem_type); + } else { + input_L_[i].SetSize(P_trial_[i]->Height(), mem_type); + } // create the necessary number of empty mfem::Vectors, to be resized later input_E_.push_back({}); input_E_buffer_.push_back({}); } - std::array trial_families = {trials::family...}; - std::array trial_orders = {trials::order...}; - std::array trial_components = {trials::components...}; - for (uint32_t i = 0; i < num_trial_spaces; i++) { - trial_function_spaces_[i] = {trial_families[i], trial_orders[i], trial_components[i]}; - } - G_test_ = QoIElementRestriction(); P_test_ = QoIProlongation(trial_fes[0]->GetParMesh()->GetComm()); output_L_.SetSize(1, mem_type); - output_T_.SetSize(1, mfem::Device::GetMemoryType()); + output_T_.SetSize(1, mem_type); // gradient objects depend on some member variables in // Functional, so we initialize the gradient objects last @@ -202,9 +215,10 @@ class Functional { std::vector arg_vec = {args...}; for (uint32_t i : arg_vec) { domain.insert_restriction(trial_space_[i], trial_function_spaces_[i]); - check_interior_face_compatibility(domain.mesh_, trial_function_spaces_[i]); } + domain.compute_interior_face_qoi_weights(); + using signature = test(decltype(serac::type(trial_spaces))...); integrals_.push_back(MakeInteriorFaceIntegral(domain, integrand, arg_vec)); } @@ -262,7 +276,25 @@ class Functional { */ double ActionOfGradient(const mfem::Vector& input_T, uint32_t which) const { - P_trial_[which]->Mult(input_T, input_L_[which]); + // Please refer to 'ActionOfGradient' in 'functional.hpp' for detailed comments of this implementation + if (trial_function_spaces_[which].family == Family::L2) { + P_trial_[which]->Mult(input_T, input_ldof_values_[which]); + input_L_[which].SetVector(input_ldof_values_[which], 0); + updateFaceNbrData(trial_space_[which], trial_pargrid_functions_[which], input_ldof_values_[which]); + + if (trial_space_[which]->GetVDim() == 1) { + input_L_[which].SetVector(trial_pargrid_functions_[which].FaceNbrData(), input_ldof_values_[which].Size()); + } else { + if (trial_space_[which]->GetOrdering() == mfem::Ordering::byVDIM) { + int local_size = input_ldof_values_[which].Size(); + appendFaceNbrData(trial_space_[which], trial_pargrid_functions_[which], local_size, input_L_[which]); + } else { + SLIC_ERROR_ROOT("Unsupported: L2 vector field ordered by nodes"); + } + } + } else { + P_trial_[which]->Mult(input_T, input_L_[which]); + } output_L_ = 0.0; @@ -280,6 +312,11 @@ class Functional { integral.GradientMult(input_E_[which], output_E_, which); + // make sure shared interior faces are integrated only once + if (dom.type_ == Domain::Type::InteriorFaces) { + output_E_ *= dom.interior_face_qoi_weights; + } + // scatter-add to compute QoI value for the local processor G_test_.ScatterAdd(output_E_, output_L_); } @@ -306,9 +343,27 @@ class Functional { { const mfem::Vector* input_T[] = {&static_cast(args)...}; + // Please refer to 'operator()' in 'functional.hpp' for detailed comments of this implementation // get the values for each local processor for (uint32_t i = 0; i < num_trial_spaces; i++) { - P_trial_[i]->Mult(*input_T[i], input_L_[i]); + if (trial_function_spaces_[i].family == Family::L2) { + P_trial_[i]->Mult(*input_T[i], input_ldof_values_[i]); + input_L_[i].SetVector(input_ldof_values_[i], 0); + updateFaceNbrData(trial_space_[i], trial_pargrid_functions_[i], input_ldof_values_[i]); + + if (trial_space_[i]->GetVDim() == 1) { + input_L_[i].SetVector(trial_pargrid_functions_[i].FaceNbrData(), input_ldof_values_[i].Size()); + } else { + if (trial_space_[i]->GetOrdering() == mfem::Ordering::byVDIM) { + int local_size = input_ldof_values_[i].Size(); + appendFaceNbrData(trial_space_[i], trial_pargrid_functions_[i], local_size, input_L_[i]); + } else { + SLIC_ERROR_ROOT("Unsupported: L2 vector field ordered by nodes"); + } + } + } else { + P_trial_[i]->Mult(*input_T[i], input_L_[i]); + } } output_L_ = 0.0; @@ -329,6 +384,11 @@ class Functional { const bool update_qdata = false; integral.Mult(t, input_E_, output_E_, wrt, update_qdata); + // make sure shared interior faces are integrated only once + if (dom.type_ == Domain::Type::InteriorFaces) { + output_E_ *= dom.interior_face_qoi_weights; + } + // scatter-add to compute QoI value for the local processor G_test_.ScatterAdd(output_E_, output_L_); } @@ -427,6 +487,8 @@ class Functional { std::map> element_gradients[Domain::num_types]; + int lcols = form_.trial_space_[which_argument]->GetVSize(); + //////////////////////////////////////////////////////////////////////////////// for (auto& integral : form_.integrals_) { @@ -456,7 +518,11 @@ class Functional { for (uint32_t i = 0; i < cols_per_elem; i++) { int col = int(trial_vdofs[i].index()); - gradient_L_[col] += K_e(e, 0, i); + + // only add local DG dof gradients (other FE spaces satisfy col < lcols inherently) + if (col < lcols) { + gradient_L_[col] += K_e(e, 0, i); + } } } } @@ -492,6 +558,12 @@ class Functional { std::array trial_function_spaces_; + /// @brief Cache for access of neighbor element info for processor boundary faces + mutable std::array trial_pargrid_functions_; + + /// @brief Cache to store values of L2 ldofs on the local processor + mutable std::array input_ldof_values_; + /** * @brief Operator that converts true (global) DOF values to local (current rank) DOF values * for the test space diff --git a/src/serac/numerics/functional/tests/functional_basic_dg.cpp b/src/serac/numerics/functional/tests/functional_basic_dg.cpp index de4ba2eedb..6c307e65eb 100644 --- a/src/serac/numerics/functional/tests/functional_basic_dg.cpp +++ b/src/serac/numerics/functional/tests/functional_basic_dg.cpp @@ -72,7 +72,7 @@ void L2_test(std::string meshfile) // auto value = residual(t, U); check_gradient(residual, t, U); } -// TEST(basic, L2_two_quad) { L2_test<2, 1>("./TwoQuad.mesh"); } + TEST(basic, L2_test_tris_and_quads_linear) { L2_test<2, 1>(SERAC_REPO_DIR "/data/meshes/patch2D_tris_and_quads.mesh"); } TEST(basic, L2_test_tris_and_quads_quadratic) { @@ -195,11 +195,21 @@ TEST(basic, L2_mixed_scalar_test_tris_and_quads_linear) L2_scalar_valued_test<2, 1>(SERAC_REPO_DIR "/data/meshes/patch2D_tris_and_quads.mesh"); } +TEST(basic, L2_mixed_scalar_test_tris_and_quads_quadratic) +{ + L2_scalar_valued_test<2, 2>(SERAC_REPO_DIR "/data/meshes/patch2D_tris_and_quads.mesh"); +} + TEST(basic, L2_mixed_scalar_test_tets_and_hexes_linear) { L2_scalar_valued_test<3, 1>(SERAC_REPO_DIR "/data/meshes/patch3D_tets_and_hexes.mesh"); } +TEST(basic, L2_mixed_scalar_test_tets_and_hexes_quadratic) +{ + L2_scalar_valued_test<3, 2>(SERAC_REPO_DIR "/data/meshes/patch3D_tets_and_hexes.mesh"); +} + int main(int argc, char* argv[]) { ::testing::InitGoogleTest(&argc, argv); From eb535bb3474a40bbd8b26f4f0baae7543cb4211c Mon Sep 17 00:00:00 2001 From: li97 Date: Fri, 11 Jul 2025 23:33:54 -0700 Subject: [PATCH 09/11] update shared interior face list --- src/serac/numerics/functional/domain.cpp | 22 +++++-------------- src/serac/numerics/functional/domain.hpp | 12 +++++----- .../numerics/functional/functional_qoi.inl | 10 ++++++--- 3 files changed, 19 insertions(+), 25 deletions(-) diff --git a/src/serac/numerics/functional/domain.cpp b/src/serac/numerics/functional/domain.cpp index ccb2d50b5c..1b106b8a04 100644 --- a/src/serac/numerics/functional/domain.cpp +++ b/src/serac/numerics/functional/domain.cpp @@ -535,39 +535,31 @@ void Domain::insert_restriction(const serac::fes_t* fes, FunctionSpace space) const BlockElementRestriction& Domain::get_restriction(FunctionSpace space) { return restriction_operators.at(space); }; -void Domain::compute_interior_face_qoi_weights() +void Domain::insert_shared_interior_face_list() { // Weights only need to be computed for Domain of InteriorFaces type SLIC_ERROR_IF(type_ != Domain::Type::InteriorFaces, "This method is only for interior face domains"); - // Weights only need to be computed once for the same domain - if (interior_face_qoi_weights.Size() == 0) { + // make a list if we don't already have one + if (shared_interior_face_ids_.size() == 0) { if (dim_ == 1) { - interior_face_qoi_weights.SetSize(int(edge_ids_.size())); - int i = 0; for (int f : mfem_edge_ids_) { mfem::Mesh::FaceInformation info = mesh_.GetFaceInformation(f); if (info.IsShared()) { - interior_face_qoi_weights[i] = 0.5; - } else { - interior_face_qoi_weights[i] = 1.0; + shared_interior_face_ids_.push_back(i); } ++i; } } else if (dim_ == 2) { - interior_face_qoi_weights.SetSize(int(tri_ids_.size() + quad_ids_.size())); - int i = 0; for (int f : mfem_tri_ids_) { mfem::Mesh::FaceInformation info = mesh_.GetFaceInformation(f); if (info.IsShared()) { - interior_face_qoi_weights[i] = 0.5; - } else { - interior_face_qoi_weights[i] = 1.0; + shared_interior_face_ids_.push_back(i); } ++i; @@ -577,9 +569,7 @@ void Domain::compute_interior_face_qoi_weights() mfem::Mesh::FaceInformation info = mesh_.GetFaceInformation(f); if (info.IsShared()) { - interior_face_qoi_weights[i] = 0.5; - } else { - interior_face_qoi_weights[i] = 1.0; + shared_interior_face_ids_.push_back(i); } ++i; diff --git a/src/serac/numerics/functional/domain.hpp b/src/serac/numerics/functional/domain.hpp index 2e02469d53..53d76e7fa3 100644 --- a/src/serac/numerics/functional/domain.hpp +++ b/src/serac/numerics/functional/domain.hpp @@ -96,6 +96,9 @@ struct Domain { */ std::map restriction_operators; + /// @brief Shared interior face list + std::vector shared_interior_face_ids_; + /** * @brief empty Domain constructor, with connectivity info to be populated later */ @@ -246,13 +249,10 @@ struct Domain { mfem::Geometry::Type element_geometry); /** - * @brief Compute face qoi weights so that interior faces shared by two processors - * are only integrated once. + * @brief Find the list of interior faces shared by two processors to make sure + * these faces are only integrated once. */ - void compute_interior_face_qoi_weights(); - - /// @brief Interior face qoi weight vector - mfem::Vector interior_face_qoi_weights; + void insert_shared_interior_face_list(); }; /// @brief constructs a domain from all the elements in a mesh diff --git a/src/serac/numerics/functional/functional_qoi.inl b/src/serac/numerics/functional/functional_qoi.inl index e53bcc3f06..5675d2cc25 100644 --- a/src/serac/numerics/functional/functional_qoi.inl +++ b/src/serac/numerics/functional/functional_qoi.inl @@ -217,7 +217,7 @@ class Functional { domain.insert_restriction(trial_space_[i], trial_function_spaces_[i]); } - domain.compute_interior_face_qoi_weights(); + domain.insert_shared_interior_face_list(); using signature = test(decltype(serac::type(trial_spaces))...); integrals_.push_back(MakeInteriorFaceIntegral(domain, integrand, arg_vec)); @@ -314,7 +314,9 @@ class Functional { // make sure shared interior faces are integrated only once if (dom.type_ == Domain::Type::InteriorFaces) { - output_E_ *= dom.interior_face_qoi_weights; + for (int shared_id : dom.shared_interior_face_ids_) { + output_E_[shared_id] *= 0.5; + } } // scatter-add to compute QoI value for the local processor @@ -386,7 +388,9 @@ class Functional { // make sure shared interior faces are integrated only once if (dom.type_ == Domain::Type::InteriorFaces) { - output_E_ *= dom.interior_face_qoi_weights; + for (int shared_id : dom.shared_interior_face_ids_) { + output_E_[shared_id] *= 0.5; + } } // scatter-add to compute QoI value for the local processor From 3cfd82866c3e82a119a8a1921132d68f1d9b2cec Mon Sep 17 00:00:00 2001 From: li97 Date: Wed, 16 Jul 2025 13:29:30 -0700 Subject: [PATCH 10/11] Add suggested clean up --- src/serac/numerics/functional/tests/dg_ghost_face_index.cpp | 4 ++-- src/serac/numerics/functional/typedefs.hpp | 1 - 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/src/serac/numerics/functional/tests/dg_ghost_face_index.cpp b/src/serac/numerics/functional/tests/dg_ghost_face_index.cpp index e60bec7258..1eb13778bb 100644 --- a/src/serac/numerics/functional/tests/dg_ghost_face_index.cpp +++ b/src/serac/numerics/functional/tests/dg_ghost_face_index.cpp @@ -11,7 +11,6 @@ #include -#include "axom/slic/core/SimpleLogger.hpp" #include "serac/infrastructure/application_manager.hpp" #include "serac/serac_config.hpp" #include "serac/mesh_utils/mesh_utils_base.hpp" @@ -62,7 +61,8 @@ void L2_index_test(std::string meshfile) // note: the orientation convention is such that the normal // computed as above will point from from side 1->2 auto [u_1, u_2] = velocity; - std::cout << "One side = " << u_1 << ", The other side = " << u_2 << ", Jump = " << u_1 - u_2 << std::endl; + SLIC_INFO(axom::fmt::format("One size = {}, The other side = {}, Jump = {}", axom::fmt::streamed(u_1), + axom::fmt::streamed(u_2), axom::fmt::streamed(u_1 - u_2))); auto a = dot(u_2 - u_1, n); diff --git a/src/serac/numerics/functional/typedefs.hpp b/src/serac/numerics/functional/typedefs.hpp index 62fb91adc4..74eb0d4f9d 100644 --- a/src/serac/numerics/functional/typedefs.hpp +++ b/src/serac/numerics/functional/typedefs.hpp @@ -12,7 +12,6 @@ namespace serac { // mfem::ParMesh / mfeme::ParFiniteElementSpace, but I'm not ready to pull the trigger on a big // interface change like that, so these typedefs mark the parts that would need to eventually change -// hanyu: I'm pulling the trigger for sam and the bullet might hit myself /// @cond using mesh_t = mfem::ParMesh; using fes_t = mfem::ParFiniteElementSpace; From 4d8ade72f9b67ed9682b076fa0c64d23d1d29994 Mon Sep 17 00:00:00 2001 From: li97 Date: Mon, 4 Aug 2025 16:56:05 -0700 Subject: [PATCH 11/11] Add suggestions in PR --- src/serac/numerics/functional/domain.cpp | 6 ++++-- src/serac/numerics/functional/domain.hpp | 2 +- src/serac/numerics/functional/functional.hpp | 7 ++++--- src/serac/numerics/functional/functional_qoi.inl | 2 -- .../functional/tests/dg_ghost_face_index.cpp | 14 ++++++++++++++ 5 files changed, 23 insertions(+), 8 deletions(-) diff --git a/src/serac/numerics/functional/domain.cpp b/src/serac/numerics/functional/domain.cpp index 1b106b8a04..33f77ba516 100644 --- a/src/serac/numerics/functional/domain.cpp +++ b/src/serac/numerics/functional/domain.cpp @@ -538,10 +538,10 @@ const BlockElementRestriction& Domain::get_restriction(FunctionSpace space) { re void Domain::insert_shared_interior_face_list() { // Weights only need to be computed for Domain of InteriorFaces type - SLIC_ERROR_IF(type_ != Domain::Type::InteriorFaces, "This method is only for interior face domains"); + SLIC_ERROR_ROOT_IF(type_ != Domain::Type::InteriorFaces, "This method is only for interior face domains"); // make a list if we don't already have one - if (shared_interior_face_ids_.size() == 0) { + if (shared_interior_face_ids_.empty()) { if (dim_ == 1) { int i = 0; for (int f : mfem_edge_ids_) { @@ -647,6 +647,8 @@ Domain InteriorFaces(const mesh_t& mesh) } } + output.insert_shared_interior_face_list(); + return output; } diff --git a/src/serac/numerics/functional/domain.hpp b/src/serac/numerics/functional/domain.hpp index 53d76e7fa3..fb29cfec9e 100644 --- a/src/serac/numerics/functional/domain.hpp +++ b/src/serac/numerics/functional/domain.hpp @@ -96,7 +96,7 @@ struct Domain { */ std::map restriction_operators; - /// @brief Shared interior face list + /// @brief Ids of interior faces that lie on the boundary shared by two processors std::vector shared_interior_face_ids_; /** diff --git a/src/serac/numerics/functional/functional.hpp b/src/serac/numerics/functional/functional.hpp index 7a8101c2ab..c9fdd7e8ef 100644 --- a/src/serac/numerics/functional/functional.hpp +++ b/src/serac/numerics/functional/functional.hpp @@ -148,9 +148,10 @@ generateParFiniteElementSpace(mfem::ParMesh* mesh) } /** - * @brief helper function to locally cast away const on FE space to update face neighbor data. - * this is ok because the original FE space is declared without const and we constrained - * the non-constness locally + * @brief helper function to locally cast away const on FE space so we can update face neighbor + * data with ExchangeFaceNbrData. This is ok because : 1) the original trial FE space is declared + * without const; 2) we constrained the non-constness locally; 3) the locally owned data associated + * with the trial function space is NOT altered and ONLY ghost data is updated. */ inline void updateFaceNbrData(const mfem::ParFiniteElementSpace* const_trial_space, mfem::ParGridFunction& trial_pgf, mfem::Vector& trial_tdof_vals) diff --git a/src/serac/numerics/functional/functional_qoi.inl b/src/serac/numerics/functional/functional_qoi.inl index 5675d2cc25..6244ff136b 100644 --- a/src/serac/numerics/functional/functional_qoi.inl +++ b/src/serac/numerics/functional/functional_qoi.inl @@ -217,8 +217,6 @@ class Functional { domain.insert_restriction(trial_space_[i], trial_function_spaces_[i]); } - domain.insert_shared_interior_face_list(); - using signature = test(decltype(serac::type(trial_spaces))...); integrals_.push_back(MakeInteriorFaceIntegral(domain, integrand, arg_vec)); } diff --git a/src/serac/numerics/functional/tests/dg_ghost_face_index.cpp b/src/serac/numerics/functional/tests/dg_ghost_face_index.cpp index 1eb13778bb..7aaff324c4 100644 --- a/src/serac/numerics/functional/tests/dg_ghost_face_index.cpp +++ b/src/serac/numerics/functional/tests/dg_ghost_face_index.cpp @@ -21,6 +21,18 @@ using namespace serac; using namespace serac::profiling; +// This test initializes a DG field with nodal coordinates of the dofs, so that the +// discontinuous dof pairs across the interior faces have the value. For example +// on the interior face with the following dofs +// {1, 2} | {5, 6} +// | +// | +// {3, 4} | {7, 8} +// we have {1, 2} = {5, 6} and {3, 4} = {7, 8}. +// It then integrates the jump of dof values over all interior faces. +// If the ghost dof data is constructed correctly to align with locally owned data, +// then every entry in the residual vector should equal to zero. This is tested +// by the L2 norm of the residual equal to zero. template void L2_index_test(std::string meshfile) { @@ -32,6 +44,7 @@ void L2_index_test(std::string meshfile) auto [test_fespace, test_fec] = serac::generateParFiniteElementSpace(mesh.get()); auto [trial_fespace, trial_fec] = serac::generateParFiniteElementSpace(mesh.get()); + // Initialize the ParGridFunction by dof coordinates mfem::ParGridFunction U_gf(trial_fespace.get()); mfem::VectorFunctionCoefficient vcoef(dim, [](const mfem::Vector& X, mfem::Vector& F) { int d = X.Size(); @@ -50,6 +63,7 @@ void L2_index_test(std::string meshfile) Domain interior_faces = InteriorFaces(*mesh); + // Define the integral of jumps over all interior faces residual.AddInteriorFaceIntegral( Dimension{}, DependsOn<0>{}, [=](double /*t*/, auto X, auto velocity) {