diff --git a/src/serac/numerics/functional/quadrature_data.hpp b/src/serac/numerics/functional/quadrature_data.hpp index b3ace3c4ec..1bc07f9c89 100644 --- a/src/serac/numerics/functional/quadrature_data.hpp +++ b/src/serac/numerics/functional/quadrature_data.hpp @@ -99,6 +99,13 @@ class ArrayView { } // namespace axom namespace serac { +namespace detail { +constexpr std::array qdata_geometries = {mfem::Geometry::SEGMENT, mfem::Geometry::TRIANGLE, + mfem::Geometry::SQUARE, mfem::Geometry::TETRAHEDRON, + mfem::Geometry::CUBE}; +constexpr std::array qdata_geometry_names = {"Segment", "Triangle", "Square", "Tetrahedron", + "Cube"}; +} // namespace detail /** * @brief A class for storing and access user-defined types at quadrature points @@ -122,10 +129,7 @@ struct QuadratureData { */ QuadratureData(geom_array_t elements, geom_array_t qpts_per_element, T value = T{}) { - constexpr std::array geometries = {mfem::Geometry::SEGMENT, mfem::Geometry::TRIANGLE, mfem::Geometry::SQUARE, - mfem::Geometry::TETRAHEDRON, mfem::Geometry::CUBE}; - - for (auto geom : geometries) { + for (auto geom : detail::qdata_geometries) { if (elements[uint32_t(geom)] > 0) { data[geom] = axom::Array(elements[uint32_t(geom)], qpts_per_element[uint32_t(geom)]); data[geom].fill(value); @@ -140,7 +144,7 @@ struct QuadratureData { axom::ArrayView operator[](mfem::Geometry::Type geom) { return axom::ArrayView(data.at(geom)); } /// @brief a 3D array indexed by (which geometry, which element, which quadrature point) - std::map > data; + std::map> data; }; /// @cond @@ -168,7 +172,7 @@ struct QuadratureData { /// @endcond /// these values exist to serve as default arguments for materials without material state -extern std::shared_ptr > NoQData; -extern std::shared_ptr > EmptyQData; +extern std::shared_ptr> NoQData; +extern std::shared_ptr> EmptyQData; } // namespace serac diff --git a/src/serac/physics/solid_mechanics.hpp b/src/serac/physics/solid_mechanics.hpp index eec39eb6f6..c29f265ad4 100644 --- a/src/serac/physics/solid_mechanics.hpp +++ b/src/serac/physics/solid_mechanics.hpp @@ -417,7 +417,7 @@ class SolidMechanics, std::integer_se qdata_type createQuadratureDataBuffer(T initial_state, const std::optional& optional_domain = std::nullopt) { Domain domain = (optional_domain) ? *optional_domain : EntireDomain(mesh_); - return StateManager::newQuadratureDataBuffer(domain, order, dim, initial_state); + return StateManager::newQuadratureDataBuffer(mesh_tag_, domain, order, dim, initial_state); } /** diff --git a/src/serac/physics/state/CMakeLists.txt b/src/serac/physics/state/CMakeLists.txt index 863df13db3..90eb49d1ca 100644 --- a/src/serac/physics/state/CMakeLists.txt +++ b/src/serac/physics/state/CMakeLists.txt @@ -32,3 +32,8 @@ install(TARGETS serac_state EXPORT serac-targets DESTINATION lib ) + +if(SERAC_ENABLE_TESTS) + add_subdirectory(tests) +endif() + diff --git a/src/serac/physics/state/state_manager.cpp b/src/serac/physics/state/state_manager.cpp index 0d40a2673c..754c214a5a 100644 --- a/src/serac/physics/state/state_manager.cpp +++ b/src/serac/physics/state/state_manager.cpp @@ -20,10 +20,10 @@ std::string StateManag std::unordered_map StateManager::named_states_; std::unordered_map StateManager::named_duals_; -double StateManager::newDataCollection(const std::string& name, const std::optional cycle_to_load) +double StateManager::newDataCollection(const std::string& mesh_tag, const std::optional cycle_to_load) { SLIC_ERROR_ROOT_IF(!ds_, "Cannot construct a DataCollection without a DataStore"); - std::string coll_name = name + "_datacoll"; + std::string coll_name = getCollectionName(mesh_tag); auto global_grp = ds_->getRoot()->createGroup(coll_name + "_global"); auto bp_index_grp = global_grp->createGroup("blueprint_index/" + coll_name); @@ -31,7 +31,7 @@ double StateManager::newDataCollection(const std::string& name, const std::optio // Needs to be configured to own the mesh data so all mesh data is saved to datastore/output file constexpr bool owns_mesh_data = true; - auto [iter, _] = datacolls_.emplace(std::piecewise_construct, std::forward_as_tuple(name), + auto [iter, _] = datacolls_.emplace(std::piecewise_construct, std::forward_as_tuple(mesh_tag), std::forward_as_tuple(coll_name, bp_index_grp, domain_grp, owns_mesh_data)); auto& datacoll = iter->second; datacoll.SetComm(MPI_COMM_WORLD); @@ -56,7 +56,7 @@ double StateManager::newDataCollection(const std::string& name, const std::optio // indicates that the mesh is periodic and the new nodal grid function must also // be discontinuous. bool is_discontinuous = false; - auto nodes = mesh(name).GetNodes(); + auto nodes = mesh(mesh_tag).GetNodes(); if (nodes) { is_discontinuous = nodes->FESpace()->FEColl()->GetContType() == mfem::FiniteElementCollection::DISCONTINUOUS; SLIC_WARNING_ROOT_IF( @@ -71,17 +71,17 @@ double StateManager::newDataCollection(const std::string& name, const std::optio // 2. Uses the existing continuity of the mesh finite element space (periodic meshes are discontinuous) // 3. Uses the spatial dimension as the mesh dimension (i.e. it is not a lower dimension manifold) // 4. Uses the ordering set by serac::ordering - mesh(name).SetCurvature(1, is_discontinuous, -1, serac::ordering); + mesh(mesh_tag).SetCurvature(1, is_discontinuous, -1, serac::ordering); // Sidre will destruct the nodal grid function instead of the mesh - mesh(name).SetNodesOwner(false); + mesh(mesh_tag).SetNodesOwner(false); // Generate the face neighbor information in the mesh. This is needed by the face restriction // operators used by Functional - mesh(name).ExchangeFaceNbrData(); + mesh(mesh_tag).ExchangeFaceNbrData(); // Construct and store the shape displacement fields and sensitivities associated with this mesh - constructShapeFields(name); + constructShapeFields(mesh_tag); } else { datacoll.SetCycle(0); // Iteration counter @@ -94,10 +94,10 @@ double StateManager::newDataCollection(const std::string& name, const std::optio void StateManager::loadCheckpointedStates(int cycle_to_load, std::vector states_to_load) { SERAC_MARK_FUNCTION; - mfem::ParMesh* meshPtr = &(*states_to_load.begin())->mesh(); - std::string mesh_name = collectionID(meshPtr); + mfem::ParMesh* meshPtr = &(*states_to_load.begin())->mesh(); + std::string mesh_tag = collectionID(meshPtr); - std::string coll_name = mesh_name + "_datacoll"; + std::string coll_name = getCollectionName(mesh_tag); axom::sidre::MFEMSidreDataCollection previous_datacoll(coll_name); @@ -107,7 +107,7 @@ void StateManager::loadCheckpointedStates(int cycle_to_load, std::vectormesh(); - SLIC_ERROR_ROOT_IF(collectionID(meshPtr) != mesh_name, + SLIC_ERROR_ROOT_IF(collectionID(meshPtr) != mesh_tag, "Loading FiniteElementStates from two different meshes at one time is not allowed."); mfem::ParGridFunction* datacoll_owned_grid_function = previous_datacoll.GetParField(state->name()); diff --git a/src/serac/physics/state/state_manager.hpp b/src/serac/physics/state/state_manager.hpp index 34207a8b76..9c99b1789e 100644 --- a/src/serac/physics/state/state_manager.hpp +++ b/src/serac/physics/state/state_manager.hpp @@ -95,10 +95,110 @@ class StateManager { */ static void storeState(FiniteElementState& state); + /** + * @brief Store a pre-constructed Quadrature Data in the state manager + * + * @tparam T the type to be created at each quadrature point + * @param mesh_tag The tag for the stored mesh used to locate the datacollection + * @param qdata The quadrature data to store + */ + template + static void storeQuadratureData(const std::string& mesh_tag, std::shared_ptr> qdata) + { + SLIC_ERROR_ROOT_IF(!ds_, "Serac's data store was not initialized - call StateManager::initialize first"); + SLIC_ERROR_ROOT_IF(!hasMesh(mesh_tag), + axom::fmt::format("Serac's state manager does not have a mesh with given tag '{}'", mesh_tag)); + + constexpr const char* qds_group_name = "quadraturedatas"; + + // Get Sidre location for quadrature data inside data collection + auto& datacoll = datacolls_.at(mesh_tag); + axom::sidre::Group* bp_group = datacoll.GetBPGroup(); // mesh_datacoll + // For each geometry type, use i to get both type and name from matching arrays + for (std::size_t i = 0; i < detail::qdata_geometries.size(); ++i) { + auto geom_type = detail::qdata_geometries[i]; + + // Check if geometry type has any data + if ((*qdata).data.find(geom_type) != (*qdata).data.end()) { + auto geom_name = detail::qdata_geometry_names[i]; + + // Get axom::Array of states in map + auto states = (*qdata)[geom_type]; + + // Get various size information + auto num_states = static_cast(states.size()); + SLIC_ERROR_ROOT_IF(num_states == 0, "Number of States should be more than 0 at this point."); + auto state_size = static_cast(sizeof(*(states.begin()))); + auto total_size = num_states * state_size; + // Sidre treats information as an array of uint8s + auto num_uint8s = total_size / static_cast(sizeof(std::uint8_t)); + + if (!is_restart_) { + axom::sidre::Group* qdatas_group = bp_group->createGroup(qds_group_name); + + // Create Sidre group, store basic information, and point Sidre at the array external to Sidre + // Note: Sidre will not own this data. + axom::sidre::Group* geom_group = qdatas_group->createGroup(std::string(geom_name)); + geom_group->createViewScalar("num_states", num_states); + geom_group->createViewScalar("state_size", state_size); + geom_group->createViewScalar("total_size", total_size); + + // Tell Sidre where the external array is, how large it is (calculated above), and what is in it (faking + // uint8) + axom::sidre::View* states_view = geom_group->createView("states"); + states_view->setExternalDataPtr(axom::sidre::UINT8_ID, num_uint8s, states.data()); + } else { + // Get Sidre group of where the states were stored. + // Note: this data is not owned by Sidre and the array should have been created at this point but + // the previous data has not been loaded yet into the array. + SLIC_ERROR_ROOT_IF(!bp_group->hasGroup(qds_group_name), + axom::fmt::format("Loaded Sidre Datastore did not have group for Quadrature Datas")); + axom::sidre::Group* qdatas_group = bp_group->getGroup(qds_group_name); + SLIC_ERROR_ROOT_IF( + !qdatas_group->hasGroup(std::string(geom_name)), + axom::fmt::format("Loaded Sidre Datastore did not have group for Quadrature Data geometry type '{}'", + std::string(geom_name))); + axom::sidre::Group* geom_group = qdatas_group->getGroup(std::string(geom_name)); + + // Verify size correctness + auto verify_size = [](axom::sidre::Group* group, int value, const std::string& view_name, + const std::string& err_msg) { + SLIC_ERROR_IF( + !group->hasView(view_name), + axom::fmt::format("Loaded Sidre Datastore does not have value '{}' for Quadrature Data.", view_name)); + auto prev_value = group->getView(view_name)->getData(); + SLIC_ERROR_IF(value != prev_value, axom::fmt::format(err_msg, value, prev_value)); + }; + verify_size(geom_group, num_states, "num_states", + "Current number of Quadrature Data States '{}' does not match value in restart '{}'."); + verify_size(geom_group, state_size, "state_size", + "Current size of Quadrature Data State '{}' does not match value in restart '{}'."); + verify_size(geom_group, total_size, "total_size", + "Current total size of Quadrature Data States '{}' does not match value in restart '{}'."); + + // Tell Sidre where the external array is + SLIC_ERROR_ROOT_IF(!geom_group->hasView("states"), + "Loaded Quadrature Data geometry Sidre view did not have 'states'"); + axom::sidre::View* states_view = geom_group->getView("states"); + states_view->setExternalDataPtr(states.data()); + } + } + } + + if (is_restart_) { + // NOTE: This call will reload all external buffers from file stored in the DataStore + // TODO: This should be changed to load only the current material quadrature data after + // MFEMSidreDatacollection::LoadExternalData and SPIO is enhanced to allow loading the + // external data piecemeal + datacoll.LoadExternalData(); + } + } + /** * @brief Create a shared ptr to a quadrature data buffer for the given material type * * @tparam T the type to be created at each quadrature point + * @param mesh_tag The tag for the stored mesh used to locate the datacollection * @param domain The spatial domain over which to allocate the quadrature data * @param order The order of the discretization of the primal fields * @param dim The spatial dimension of the mesh @@ -106,8 +206,8 @@ class StateManager { * @return shared pointer to quadrature data buffer */ template - static std::shared_ptr> newQuadratureDataBuffer(const Domain& domain, int order, int dim, - T initial_state) + static std::shared_ptr> newQuadratureDataBuffer(const std::string& mesh_tag, const Domain& domain, + int order, int dim, T initial_state) { int Q = order + 1; @@ -125,7 +225,9 @@ class StateManager { qpts_per_elem[size_t(geom)] = uint32_t(num_quadrature_points(geom, Q)); } - return std::make_shared>(elems, qpts_per_elem, initial_state); + auto qdata = std::make_shared>(elems, qpts_per_elem, initial_state); + storeQuadratureData(mesh_tag, qdata); + return qdata; } /** @@ -336,11 +438,18 @@ class StateManager { private: /** * @brief Creates a new datacollection based on a registered mesh - * @param[in] name The name of the new datacollection + * @param[in] mesh_tag The mesh name used to name the new datacollection * @param[in] cycle_to_load What cycle to load the DataCollection from, if applicable * @return The time from specified restart cycle. Otherwise zero. */ - static double newDataCollection(const std::string& name, const std::optional cycle_to_load = {}); + static double newDataCollection(const std::string& mesh_tag, const std::optional cycle_to_load = {}); + + /** + * @brief Returns the name of the data collection's name for a given mesh_tag + * @param[in] mesh_tag The mesh name used to name the new datacollection + * @return The name of the data collection for the given mesh_tag + */ + static std::string getCollectionName(const std::string& mesh_tag) { return mesh_tag + "_datacoll"; } /** * @brief Construct the shape displacement field for the requested mesh diff --git a/src/serac/physics/state/tests/CMakeLists.txt b/src/serac/physics/state/tests/CMakeLists.txt new file mode 100644 index 0000000000..d6fd8cd798 --- /dev/null +++ b/src/serac/physics/state/tests/CMakeLists.txt @@ -0,0 +1,23 @@ +# Copyright (c) 2019-2024, Lawrence Livermore National Security, LLC and +# other Serac Project Developers. See the top-level LICENSE file for +# details. +# +# SPDX-License-Identifier: (BSD-3-Clause) + +set(state_manager_tests_depends serac_state serac_physics_materials serac_functional serac_mesh gtest) + +set(state_manager_serial_test_sources + state_manager.cpp + ) + +serac_add_tests(SOURCES ${state_manager_serial_test_sources} + DEPENDS_ON ${state_manager_tests_depends} + NUM_MPI_TASKS 1) + +# set(state_manager_parallel_test_sources +# state_manager.cpp +# ) + +# serac_add_tests(SOURCES ${state_manager_parallel_test_sources} +# DEPENDS_ON ${state_manager_tests_depends} +# NUM_MPI_TASKS 2) diff --git a/src/serac/physics/state/tests/state_manager.cpp b/src/serac/physics/state/tests/state_manager.cpp new file mode 100644 index 0000000000..919c41afa4 --- /dev/null +++ b/src/serac/physics/state/tests/state_manager.cpp @@ -0,0 +1,154 @@ +// Copyright (c) 2019-2024, 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 "axom/slic/core/SimpleLogger.hpp" +#include "gtest/gtest.h" + +#include "serac/mesh/mesh_utils.hpp" +#include "serac/numerics/functional/domain.hpp" +#include "serac/numerics/functional/quadrature_data.hpp" +#include "serac/physics/state/state_manager.hpp" +#include "serac/physics/materials/solid_material.hpp" +#include "serac/serac_config.hpp" + +namespace serac { + +namespace detail { + +template +void apply_function_to_quadrature_data_states(const double starting_value, std::shared_ptr qdata, + std::function& apply_function) +{ + for (std::size_t i = 0; i < detail::qdata_geometries.size(); ++i) { + auto geom_type = detail::qdata_geometries[i]; + + // Check if geometry type has any data + if ((*qdata).data.find(geom_type) != (*qdata).data.end()) { + // Get axom::Array of states in map + auto states = (*qdata)[geom_type]; + double curr_value = starting_value; + for (auto& state : states) { + apply_function(curr_value, state); + } + } + } +} + +template +bool compare_tensors(const tensor& a, const tensor& b) +{ + for (int i = 0; i < M; ++i) { + for (int j = 0; j < N; ++j) { + if (a(i, j) != b(i, j)) { + return false; + } + } + } + return true; +} + +} // namespace detail + +TEST(state_manager, basic) +{ + MPI_Barrier(MPI_COMM_WORLD); + + // TODO: are these right? + constexpr int dim = 3; + constexpr int order = 2; + + // Info about this test's Quadrature data state + /* + struct State { + tensor Fpinv = DenseIdentity<3>(); ///< inverse of plastic distortion tensor + double accumulated_plastic_strain; ///< uniaxial equivalent plastic strain + }; + */ + using State = serac::solid_mechanics::J2::State; + + //--------------------------------- Helper functions for this test + // Lamda to check the state against a starting value which is incremented after each check + std::function check_state = [](double& curr_value, State& state) { + tensor expected_tensor = make_tensor([&](int i, int j) { return i + curr_value * j; }); + EXPECT_TRUE(detail::compare_tensors(state.Fpinv, expected_tensor)); + EXPECT_DOUBLE_EQ(state.accumulated_plastic_strain, curr_value); + curr_value++; + }; + + // Lamda to fill the state against a starting value which is incremented after each check + std::function fill_state = [](double& curr_value, State& state) { + state.Fpinv = make_tensor([&](int i, int j) { return i + curr_value * j; }); + state.accumulated_plastic_strain = curr_value; + curr_value++; + }; + //--------------------------------- + + // Create DataStore + std::string name = "basic"; + axom::sidre::DataStore datastore; + StateManager::initialize(datastore, name + "_data"); + + // Construct the appropriate dimension mesh and give it to the StateManager + std::string filename = SERAC_REPO_DIR "/data/meshes/ball.mesh"; + std::string mesh_tag = "ball_mesh"; + auto mesh = mesh::refineAndDistribute(buildMeshFromFile(filename), 1, 0); + StateManager::setMesh(std::move(mesh), mesh_tag); + + // Create and store the initial state of the quadrature data in sidre + SLIC_INFO("Creating Quadrature Data with initial state"); + Domain domain = EntireDomain(StateManager::mesh(mesh_tag)); + State initial_state{}; + std::shared_ptr> qdata = + StateManager::newQuadratureDataBuffer(mesh_tag, domain, order, dim, initial_state); + + // Change data + SLIC_INFO("Populating Quadrature Data"); + constexpr double good_starting_value = 1.0; + detail::apply_function_to_quadrature_data_states(good_starting_value, qdata, fill_state); + SLIC_INFO("Verifying populated Quadrature Data"); + detail::apply_function_to_quadrature_data_states(good_starting_value, qdata, check_state); + + // Save to disk and simulate a restart + const int cycle = 1; + const double time_saved = 1.5; + SLIC_INFO(axom::fmt::format("Saving mesh restart '{0}' at cycle '{1}'", mesh_tag, cycle)); + StateManager::save(time_saved, cycle, mesh_tag); + + // Reset StateManager then load from disk + SLIC_INFO("Clearing current and loading previously saved State Manager"); + StateManager::reset(); + axom::sidre::DataStore new_datastore; + StateManager::initialize(new_datastore, name + "_data"); + StateManager::load(cycle, mesh_tag); + + // Load data from disk + SLIC_INFO("Loading previously saved Quadrature Data"); + Domain new_domain = EntireDomain(StateManager::mesh(mesh_tag)); + std::shared_ptr> new_qdata = + StateManager::newQuadratureDataBuffer(mesh_tag, new_domain, order, dim, initial_state); + + // Verify data has reloaded to restart data + SLIC_INFO("Verifying loaded Quadrature Data"); + detail::apply_function_to_quadrature_data_states(good_starting_value, new_qdata, check_state); +} + +} // namespace serac + +int main(int argc, char* argv[]) +{ + ::testing::InitGoogleTest(&argc, argv); + MPI_Init(&argc, &argv); + + axom::slic::SimpleLogger logger; + + int result = RUN_ALL_TESTS(); + MPI_Finalize(); + + return result; +} diff --git a/src/serac/physics/tests/solid_robin_condition.cpp b/src/serac/physics/tests/solid_robin_condition.cpp index 7b90bc2529..d335e6280d 100644 --- a/src/serac/physics/tests/solid_robin_condition.cpp +++ b/src/serac/physics/tests/solid_robin_condition.cpp @@ -44,7 +44,6 @@ void functional_solid_test_robin_condition() auto mesh = mesh::refineAndDistribute(buildMeshFromFile(filename), serial_refinement, parallel_refinement); serac::StateManager::setMesh(std::move(mesh), mesh_tag); - // _solver_params_start serac::NonlinearSolverOptions nonlinear_options{.nonlin_solver = NonlinearSolver::Newton, .relative_tol = 1.0e-12, .absolute_tol = 1.0e-12, @@ -53,7 +52,6 @@ void functional_solid_test_robin_condition() SolidMechanics solid_solver(nonlinear_options, solid_mechanics::default_linear_options, solid_mechanics::default_quasistatic_options, "solid_mechanics", mesh_tag); - // _solver_params_end solid_mechanics::LinearIsotropic mat{ 1.0, // mass density diff --git a/src/serac/physics/tests/thermal_nonlinear_solve.cpp b/src/serac/physics/tests/thermal_nonlinear_solve.cpp index 32cc4a1c5f..4cea8ad4c0 100644 --- a/src/serac/physics/tests/thermal_nonlinear_solve.cpp +++ b/src/serac/physics/tests/thermal_nonlinear_solve.cpp @@ -43,7 +43,6 @@ void functional_thermal_test_nonlinear() auto mesh = mesh::refineAndDistribute(buildMeshFromFile(filename), serial_refinement, parallel_refinement); serac::StateManager::setMesh(std::move(mesh), mesh_tag); - // _solver_params_start serac::NonlinearSolverOptions nonlinear_options{.nonlin_solver = NonlinearSolver::NewtonLineSearch, .relative_tol = 1.0e-12, .absolute_tol = 1.0e-12, @@ -52,7 +51,6 @@ void functional_thermal_test_nonlinear() HeatTransfer thermal_solver(nonlinear_options, heat_transfer::default_linear_options, heat_transfer::default_static_options, "heat_transfer", mesh_tag); - // _solver_params_end heat_transfer::LinearIsotropicConductor mat{ 1.0, // mass density diff --git a/src/serac/physics/tests/thermal_robin_condition.cpp b/src/serac/physics/tests/thermal_robin_condition.cpp index 6a2d7b6146..5fd9871292 100644 --- a/src/serac/physics/tests/thermal_robin_condition.cpp +++ b/src/serac/physics/tests/thermal_robin_condition.cpp @@ -43,7 +43,6 @@ void functional_thermal_test_robin_condition() auto mesh = mesh::refineAndDistribute(buildMeshFromFile(filename), serial_refinement, parallel_refinement); serac::StateManager::setMesh(std::move(mesh), mesh_tag); - // _solver_params_start serac::NonlinearSolverOptions nonlinear_options{.nonlin_solver = NonlinearSolver::Newton, .relative_tol = 1.0e-12, .absolute_tol = 1.0e-12, @@ -52,7 +51,6 @@ void functional_thermal_test_robin_condition() HeatTransfer thermal_solver(nonlinear_options, heat_transfer::default_linear_options, heat_transfer::default_static_options, "heat_transfer", mesh_tag); - // _solver_params_end heat_transfer::LinearIsotropicConductor mat{ 1.0, // mass density