diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 8000ea266..cd581f3bc 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -6,6 +6,7 @@ variables: #TODO: uncomment this when everyone has service user access across all the machines #LLNL_SERVICE_USER: asmith + GIT_STRATEGY: clone GIT_SUBMODULE_STRATEGY: recursive PROJECT_ALLOC_NAME: ${CI_PROJECT_NAME}_ci_${CI_PIPELINE_ID} BUILD_ROOT: ${CI_PROJECT_DIR} diff --git a/.gitlab/build_toss4.yml b/.gitlab/build_toss4.yml index 6a369bfab..c4f83eabc 100644 --- a/.gitlab/build_toss4.yml +++ b/.gitlab/build_toss4.yml @@ -61,11 +61,11 @@ toss4-gcc_10_3_1-src-no-tribol: ALLOC_TIME: "30" extends: .src_build_on_toss4 -toss4-gcc_10_3_1-src-no-sundials: +toss4-gcc_10_3_1-src-no-optional-solvers: variables: COMPILER: "gcc@10.3.1" HOST_CONFIG: "ruby-toss_4_x86_64_ib-${COMPILER}.cmake" - EXTRA_CMAKE_OPTIONS: "-DENABLE_BENCHMARKS=ON -USUNDIALS_DIR" + EXTRA_CMAKE_OPTIONS: "-DENABLE_BENCHMARKS=ON -USUNDIALS_DIR -UPETSC_DIR" ALLOC_NODES: "1" ALLOC_TIME: "20" extends: .src_build_on_toss4 diff --git a/cmake/blt b/cmake/blt index 66378271f..ca27776de 160000 --- a/cmake/blt +++ b/cmake/blt @@ -1 +1 @@ -Subproject commit 66378271f605de0b628ab93196fd41f5c7b59452 +Subproject commit ca27776de9626102725ebac14c37617a94e71a67 diff --git a/scripts/spack/specs.json b/scripts/spack/specs.json index 5a941d5d3..6f245df35 100644 --- a/scripts/spack/specs.json +++ b/scripts/spack/specs.json @@ -18,8 +18,5 @@ "clang@14.0.6+devtools+profiling" ], "blueos_3_ppc64le_ib_p9": - [ "clang@10.0.1+devtools+cuda+raja+umpire+profiling~strumpack cuda_arch=70" ], - - "darwin-x86_64": - [ "clang@9.0.0+devtools" ] + [ "clang@10.0.1+devtools+cuda+raja+umpire+profiling~strumpack cuda_arch=70" ] } diff --git a/src/docs/sphinx/dev_guide/profiling.rst b/src/docs/sphinx/dev_guide/profiling.rst index 7f5ea83a1..48df8a818 100644 --- a/src/docs/sphinx/dev_guide/profiling.rst +++ b/src/docs/sphinx/dev_guide/profiling.rst @@ -27,7 +27,7 @@ Introduction to SPOT -------------------- `SPOT `_ is a framework developed at -LLNL for vizualizing performance data. SPOT is an external tool and does not need to be +LLNL for visualizing performance data. SPOT is an external tool and does not need to be linked into Serac. TPL Build Instructions @@ -121,9 +121,21 @@ with benchmarking enabled (off by default). Then, run the build target ``run_ben cd make -j make run_benchmarks - find . -name "*.cali" -print0 | xargs -0 mv -t . pwd This will run all of Serac's benchmarks multiple times with varying MPI task counts, and generate a Caliper file for -each benchmark run. The ``find`` command afterwards ensures all Caliper files are moved to the same directory. Now, you -can visualize the results with SPOT, entering the path printed from ``pwd``. +each benchmark run at ``PROJECT_BINARY_DIR``. Now, you can visualize the results with SPOT, entering the path printed +from ``pwd``. + +Visualizing Benchmarks using SPOT +--------------------------------- + +If you have access to LC, you can go to the following website and enter a directory in CZ/ RZ that contains Caliper +files: + +- `SPOT CZ `_ +- `SPOT RZ `_ + +.. note:: + There is a bug in SPOT where if you remove Caliper files from a directory, they still show up on SPOT - if you've + visualized them previously. The current workaround is by removing the ``llnl.gov`` site cache manually. diff --git a/src/docs/sphinx/dev_guide/testing.rst b/src/docs/sphinx/dev_guide/testing.rst index 1962d9603..7718c1c0a 100644 --- a/src/docs/sphinx/dev_guide/testing.rst +++ b/src/docs/sphinx/dev_guide/testing.rst @@ -11,7 +11,7 @@ Testing Serac has two levels of tests, unit and integration. Unit tests are used to test individual components of code, such as a class or function. While integration tests -are for testing the code as a whole. For example, testing the `serac` driver with +are for testing the code as a whole. For example, testing the ``serac`` driver with an input file against blessed answers. Unit Tests @@ -41,7 +41,7 @@ Requirements: >>> socket.gethostname().rstrip('1234567890') >>> exit() - Currently, there are configuration json files for Toss3 and BlueOS which can be + Currently, there are configuration json files for TOSS4 and BlueOS which can be used as reference. #. **Build the code.** @@ -54,7 +54,7 @@ Requirements: # BlueOS $ lalloc 2 ./ats.sh - # Toss3 + # TOSS4 $ salloc -N2 ./ats.sh # Personal Machine (currently runs subset of tests) @@ -72,6 +72,22 @@ Requirements: ATS also outputs both a ``.log`` and ``.log.err`` for each test and checker that is run. +#. **Rebaseline tests (as needed).** + If tolerances to tests need to be updated, first ensure you've generated new tolerances by running the integration + tests like mentioned above. Then, use the ``-b`` option:: + + # Single baseline + $ ./ats.sh -b dyn_solve_serial + + # Comma-separated baselines + $ ./ats.sh -b dyn_solve_serial,dyn_solve_parallel + + # All baselines + $ ./ats.sh -b all + + This will update the json files located in the `serac_tests `_ submodule. To + avoid Caliper files from additionally being generated, configure with ``-DENABLE_BENCHMARKS=OFF``. + Installing ATS -------------- diff --git a/src/docs/sphinx/quickstart.rst b/src/docs/sphinx/quickstart.rst index 0b3ffd25a..f9fe3b4fd 100644 --- a/src/docs/sphinx/quickstart.rst +++ b/src/docs/sphinx/quickstart.rst @@ -78,7 +78,7 @@ For example on **Ubuntu 20.04**: python3 scripts/uberenv/uberenv.py --project-json=scripts/spack/devtools.json --spack-env-file=scripts/spack/configs/linux_ubuntu_20/spack.yaml --prefix=../path/to/install Unlike Serac's library dependencies, our developer tools can be built with any compiler because -they are not linked into the serac executable. We recommend GCC 8 because we have tested that they all +they are not linked into the serac executable. We recommend GCC 10 because we have tested that they all build with that compiler. Building Serac's Dependencies via Spack/uberenv @@ -97,7 +97,7 @@ This has been encapsulated using `Uberenv `_. U doing the following: * Pulls a blessed version of Spack locally -* If you are on a known operating system (like TOSS3), we have defined Spack configuration files +* If you are on a known operating system (like TOSS4), we have defined Spack configuration files to keep Spack from building the world * Installs our Spack packages into the local Spack * Simplifies whole dependency build into one command @@ -120,18 +120,18 @@ Serac. .. note:: On LC machines, it is good practice to do the build step in parallel on a compute node. - Here is an example command: ``salloc -ppdebug -N1-1 python3 scripts/uberenv/uberenv.py`` + Here is an example command: ``salloc -ppdebug -N1 python3 scripts/uberenv/uberenv.py`` Unless otherwise specified Spack will default to a compiler. This is generally not a good idea when developing large codes. To specify which compiler to use add the compiler specification to the ``--spec`` Uberenv -command line option. On TOSS3, we recommend and have tested ``--spec=%clang@10.0.0``. More compiler specs +command line option. On TOSS4, we recommend and have tested ``--spec=%clang@14.0.6``. More compiler specs can be found in the Spack compiler files in our repository: ``scripts/spack/configs//spack.yaml``. We currently regularly test the following Spack configuration files: -* Linux Ubuntu 20.04 (via Windows WSL 2) -* TOSS 3 (On Ruby at LC) +* Linux Ubuntu 22.04 (via Azure) +* TOSS4 (On Ruby at LC) * BlueOS (On Lassen at LC) To install Serac on a new platform, it is a good idea to start with a known Spack environments file, or ``spack.yaml`` file, @@ -144,7 +144,7 @@ in the `Spack docs `_ If you do not have a ``spack.yaml`` already, you can leave off that command line option from ``uberenv`` and Spack will generate a new one for you. Uberenv will copy it where you ran your uberenv command for future use. .. note:: - A newer vesion of cmake (>=3.20) and llvm (>=14) may be required. + A newer version of cmake (>=3.20) and llvm (>=14) may be required. Some helpful uberenv options include : @@ -152,11 +152,11 @@ Some helpful uberenv options include : * ``--spec=" build_type=Debug"`` (build the MFEM and Hypre libraries with debug symbols) * ``--spec=+profiling`` (build the Adiak and Caliper libraries) * ``--spec=+devtools`` (also build the devtools with one command) -* ``--spec=%clang@10.0.0`` (build with a specific compiler as defined in the ``spack.yaml`` file) +* ``--spec=%clang@14.0.6`` (build with a specific compiler as defined in the ``spack.yaml`` file) * ``--spack-env-file=`` (use specific Spack environment configuration file) * ``--prefix=`` (required, build and install the dependencies in a particular location) - this *must be outside* of your local Serac repository -The modifiers to the Spack specification ``spec`` can be chained together, e.g. ``--spec='%clang@10.0.0+devtools build_type=Debug'``. +The modifiers to the Spack specification ``spec`` can be chained together, e.g. ``--spec='%clang@14.0.6+devtools build_type=Debug'``. If you already have a Spack instance from another project that you would like to reuse, you can do so by changing the uberenv command as follows: @@ -169,7 +169,7 @@ Building Serac's Dependencies by Hand ------------------------------------- To build Serac's dependencies by hand, use of a ``host-config`` CMake configuration file is -stongly encouraged. A good place to start is by copying an existing host config in the +strongly encouraged. A good place to start is by copying an existing host config in the ``host-config`` directory and modifying it according to your system setup. .. _build-label: diff --git a/src/serac/infrastructure/accelerator.cpp b/src/serac/infrastructure/accelerator.cpp index 818993efc..cdce5eba4 100644 --- a/src/serac/infrastructure/accelerator.cpp +++ b/src/serac/infrastructure/accelerator.cpp @@ -25,7 +25,7 @@ void initializeDevice() { SLIC_ERROR_ROOT_IF(device, "serac::accelerator::initializeDevice cannot be called more than once"); device = std::make_unique(); -#ifdef MFEM_USE_CUDA +#if defined(MFEM_USE_CUDA) && defined(SERAC_USE_CUDA_KERNEL_EVALUATION) device->Configure("cuda"); #endif } diff --git a/src/serac/infrastructure/tests/error_handling.cpp b/src/serac/infrastructure/tests/error_handling.cpp index c221ad027..7b63a82ae 100644 --- a/src/serac/infrastructure/tests/error_handling.cpp +++ b/src/serac/infrastructure/tests/error_handling.cpp @@ -57,39 +57,6 @@ TEST(ErrorHandling, BcOneComponentVectorCoefDofs) EXPECT_THROW(BoundaryCondition(coef, 0, *space, dofs), SlicErrorException); } -TEST(ErrorHandling, BcRetrieveScalarCoef) -{ - auto mesh = mesh::refineAndDistribute(buildDiskMesh(10), 0, 0); - - auto [space, coll] = serac::generateParFiniteElementSpace>(mesh.get()); - - auto coef = std::make_shared(1.0); - BoundaryCondition bc(coef, -1, *space, std::set{}); - EXPECT_NO_THROW(bc.scalarCoefficient()); - EXPECT_THROW(bc.vectorCoefficient(), SlicErrorException); - - const auto& const_bc = bc; - EXPECT_NO_THROW(const_bc.scalarCoefficient()); - EXPECT_THROW(const_bc.vectorCoefficient(), SlicErrorException); -} - -TEST(ErrorHandling, BcRetrieveVecCoef) -{ - auto mesh = mesh::refineAndDistribute(buildDiskMesh(10), 0, 0); - - auto [space, coll] = serac::generateParFiniteElementSpace>(mesh.get()); - - mfem::Vector vec; - auto coef = std::make_shared(vec); - BoundaryCondition bc(coef, {}, *space, std::set{}); - EXPECT_NO_THROW(bc.vectorCoefficient()); - EXPECT_THROW(bc.scalarCoefficient(), SlicErrorException); - - const auto& const_bc = bc; - EXPECT_NO_THROW(const_bc.vectorCoefficient()); - EXPECT_THROW(const_bc.scalarCoefficient(), SlicErrorException); -} - TEST(ErrorHandling, InvalidCmdlineArg) { // The command is actually --input-file diff --git a/src/serac/numerics/functional/shape_aware_functional.hpp b/src/serac/numerics/functional/shape_aware_functional.hpp index e1c91aacf..2e2cfc172 100644 --- a/src/serac/numerics/functional/shape_aware_functional.hpp +++ b/src/serac/numerics/functional/shape_aware_functional.hpp @@ -13,6 +13,7 @@ #pragma once +#include "serac/infrastructure/accelerator.hpp" #include "serac/numerics/functional/functional.hpp" namespace serac { @@ -41,7 +42,7 @@ struct ShapeCorrection { * * @param p The current shape displacement field at the underlying quadrature point */ - ShapeCorrection(Dimension, shape_type p) + SERAC_HOST_DEVICE ShapeCorrection(Dimension, shape_type p) : J_(Identity() + get(p)), detJ_(det(J_)), inv_J_(inv(J_)), inv_JT_(transpose(inv_J_)) { } @@ -267,7 +268,7 @@ SERAC_HOST_DEVICE auto apply_shape_aware_qf_helper_with_state(const lambda& qf, } // namespace detail /// @cond -template +template class ShapeAwareFunctional; /// @endcond @@ -356,6 +357,86 @@ class ShapeAwareFunctional { functional_ = std::make_unique>(prepended_spaces); } + /** + * @brief Functor representing a shape-aware integrand. Used instead of an extended generic + * lambda for compatibility with NVCC. + */ + template + struct ShapeAwareIntegrandWrapper { + /// @brief Constructor for functor + ShapeAwareIntegrandWrapper(Integrand integrand) : integrand_(integrand) {} + /// @brief Integrand + Integrand integrand_; + /** + * @brief integrand call + * + * @tparam PositionType position type + * @tparam ShapeValueType shape value type + * @tparam QFuncArgs type of variadic pack to forward to qfunc + * @param[in] time time + * @param[in] x position + * @param[in] shape_val shape + * @param[in] qfunc_args qfunc parameter pack + * @return Shape aware qf value + */ + template + SERAC_HOST_DEVICE auto operator()(double time, PositionType x, ShapeValueType shape_val, + QFuncArgs... qfunc_args) const + { + auto qfunc_tuple = make_tuple(qfunc_args...); + auto reduced_trial_space_tuple = make_tuple(get(trial_spaces)...); + + detail::ShapeCorrection shape_correction(Dimension{}, shape_val); + // TODO(CUDA): When this is compiled to device code, the below make_integer_sequence will + // need to change to a camp integer sequence. + auto unmodified_qf_return = detail::apply_shape_aware_qf_helper( + integrand_, time, x, shape_val, reduced_trial_space_tuple, qfunc_tuple, shape_correction, + std::make_integer_sequence{}); + return shape_correction.modify_shape_aware_qf_return(test_space, unmodified_qf_return); + } + }; + + /** + * @brief Functor representing a shape-aware integrand with state. Used instead of an extended generic + * lambda for compatibility with NVCC. + */ + template + struct ShapeAwareIntegrandWrapperWithState { + /// @brief Constructor for functor + ShapeAwareIntegrandWrapperWithState(Integrand integrand) : integrand_(integrand) {} + /// @brief integrand + Integrand integrand_; + /** + * @brief integrand call + * + * @tparam PositionType position type + * @tparam StateType state type + * @tparam ShapeValueType shape value type + * @tparam QFuncArgs type of variadic pack to forward to qfunc + * @param[in] time time + * @param[in] x position + * @param[in] state reference to state + * @param[in] shape_val shape + * @param[in] qfunc_args qfunc parameter pack + * @return shape aware integrand value + */ + template + SERAC_HOST_DEVICE auto operator()(double time, PositionType x, StateType& state, ShapeValueType shape_val, + QFuncArgs... qfunc_args) const + { + auto qfunc_tuple = make_tuple(qfunc_args...); + auto reduced_trial_space_tuple = make_tuple(get(trial_spaces)...); + + detail::ShapeCorrection shape_correction(Dimension{}, shape_val); + // TODO(CUDA): When this is compiled to device code, the below make_integer_sequence will + // need to change to a camp integer sequence. + auto unmodified_qf_return = detail::apply_shape_aware_qf_helper_with_state( + integrand_, time, x, state, shape_val, reduced_trial_space_tuple, qfunc_tuple, shape_correction, + std::make_integer_sequence{}); + return shape_correction.modify_shape_aware_qf_return(test_space, unmodified_qf_return); + } + }; + /** * @brief Adds a domain integral term to the weak formulation of the PDE * @@ -377,40 +458,48 @@ class ShapeAwareFunctional { std::shared_ptr> qdata = NoQData) { if constexpr (std::is_same_v) { - functional_->AddDomainIntegral( - Dimension{}, DependsOn<0, (args + 1)...>{}, - [integrand](double time, auto x, auto shape_val, auto... qfunc_args) { - auto qfunc_tuple = make_tuple(qfunc_args...); - auto reduced_trial_space_tuple = make_tuple(get(trial_spaces)...); - - detail::ShapeCorrection shape_correction(Dimension{}, shape_val); - - auto unmodified_qf_return = detail::apply_shape_aware_qf_helper( - integrand, time, x, shape_val, reduced_trial_space_tuple, qfunc_tuple, shape_correction, - std::make_integer_sequence{}); - return shape_correction.modify_shape_aware_qf_return(test_space, unmodified_qf_return); - }, - domain, qdata); + functional_->AddDomainIntegral(Dimension{}, DependsOn<0, (args + 1)...>{}, + ShapeAwareIntegrandWrapper(integrand), domain, qdata); } else { - functional_->AddDomainIntegral( - Dimension{}, DependsOn<0, (args + 1)...>{}, - [integrand](double time, auto x, auto& state, auto shape_val, auto... qfunc_args) { - auto qfunc_tuple = make_tuple(qfunc_args...); - auto reduced_trial_space_tuple = make_tuple(get(trial_spaces)...); - - detail::ShapeCorrection shape_correction(Dimension{}, shape_val); - - auto unmodified_qf_return = detail::apply_shape_aware_qf_helper_with_state( - integrand, time, x, state, shape_val, reduced_trial_space_tuple, qfunc_tuple, shape_correction, - std::make_integer_sequence{}); - return shape_correction.modify_shape_aware_qf_return(test_space, unmodified_qf_return); - }, - domain, qdata); + functional_->AddDomainIntegral(Dimension{}, DependsOn<0, (args + 1)...>{}, + ShapeAwareIntegrandWrapperWithState(integrand), domain, + qdata); } } /** - * @brief Adds a boundary integral term to the weak formulation of the PDE + * @brief Functor representing a shape-aware integrand. Used instead of an extended generic + * lambda for compatibility with NVCC. + */ + template + struct ShapeAwareBoundaryIntegrandWrapper { + /// @brief Constructor for functor + ShapeAwareBoundaryIntegrandWrapper(Integrand integrand) : integrand_(integrand) {} + /// @brief integrand + Integrand integrand_; + /** + * @brief integrand call + * + * @tparam PositionType position type + * @tparam ShapeValueType shape value type + * @tparam QFuncArgs type of variadic pack to forward to qfunc + * @param[in] time time + * @param[in] x position + * @param[in] shape_val shape + * @param[in] qfunc_args qfunc parameter pack + * @return qf function corrected for boundary area + */ + template + SERAC_HOST_DEVICE auto operator()(double time, PositionType x, ShapeValueType shape_val, + QFuncArgs... qfunc_args) const + { + auto unmodified_qf_return = integrand_(time, x + shape_val, qfunc_args...); + return unmodified_qf_return * detail::compute_boundary_area_correction(x, shape_val); + } + }; + + /** + * @brief Adds a boundary integral term to the weak formulationx of the PDE * * @tparam dim The dimension of the element (2 for quad, 3 for hex, etc) * @tparam args The type of the trial function input arguments @@ -426,14 +515,8 @@ class ShapeAwareFunctional { template void AddBoundaryIntegral(Dimension, DependsOn, const lambda& integrand, domain_type& domain) { - functional_->AddBoundaryIntegral( - Dimension{}, DependsOn<0, (args + 1)...>{}, - [integrand](double time, auto x, auto shape_val, auto... qfunc_args) { - auto unmodified_qf_return = integrand(time, x + shape_val, qfunc_args...); - - return unmodified_qf_return * detail::compute_boundary_area_correction(x, shape_val); - }, - domain); + functional_->AddBoundaryIntegral(Dimension{}, DependsOn<0, (args + 1)...>{}, + ShapeAwareBoundaryIntegrandWrapper(integrand), domain); } /** diff --git a/src/serac/numerics/functional/tests/bug_boundary_qoi.cpp b/src/serac/numerics/functional/tests/bug_boundary_qoi.cpp index d0e8ce78c..97e6a6159 100644 --- a/src/serac/numerics/functional/tests/bug_boundary_qoi.cpp +++ b/src/serac/numerics/functional/tests/bug_boundary_qoi.cpp @@ -24,6 +24,14 @@ using namespace serac::profiling; double t = 0.0; +struct IdentityFunctor { + template + SERAC_HOST_DEVICE auto operator()(Arg1, Arg2) const + { + return 1.0; + } +}; + int num_procs, my_rank; TEST(BoundaryIntegralQOI, AttrBug) @@ -41,8 +49,7 @@ TEST(BoundaryIntegralQOI, AttrBug) auto [shape_fes, shape_fec] = serac::generateParFiniteElementSpace(pmesh.get()); serac::ShapeAwareFunctional totalSurfArea(shape_fes.get(), {}); - totalSurfArea.AddBoundaryIntegral( - serac::Dimension<2 - 1>{}, serac::DependsOn<>{}, [](auto, auto) { return 1.0; }, *pmesh); + totalSurfArea.AddBoundaryIntegral(serac::Dimension<2 - 1>{}, serac::DependsOn<>{}, IdentityFunctor{}, *pmesh); serac::FiniteElementState shape(*shape_fes); double totalSurfaceArea = totalSurfArea(0.0, shape); @@ -50,8 +57,7 @@ TEST(BoundaryIntegralQOI, AttrBug) serac::Domain attr1 = serac::Domain::ofBoundaryElements(*pmesh, serac::by_attr<2>(1)); serac::ShapeAwareFunctional attr1SurfArea(shape_fes.get(), {}); - attr1SurfArea.AddBoundaryIntegral( - serac::Dimension<2 - 1>{}, serac::DependsOn<>{}, [](auto, auto) { return 1.0; }, attr1); + attr1SurfArea.AddBoundaryIntegral(serac::Dimension<2 - 1>{}, serac::DependsOn<>{}, IdentityFunctor{}, attr1); double attr1SurfaceArea = attr1SurfArea(0.0, shape); EXPECT_NEAR(attr1SurfaceArea, 1.0, 1.0e-14); diff --git a/src/serac/numerics/functional/tests/check_gradient.hpp b/src/serac/numerics/functional/tests/check_gradient.hpp index 3292be88d..fb48070f3 100644 --- a/src/serac/numerics/functional/tests/check_gradient.hpp +++ b/src/serac/numerics/functional/tests/check_gradient.hpp @@ -9,24 +9,34 @@ #include "mfem.hpp" +#include "serac/infrastructure/accelerator.hpp" #include "serac/serac_config.hpp" #include "serac/numerics/functional/functional.hpp" -template -void check_gradient(serac::Functional& f, double t, const mfem::Vector& U, double epsilon = 1.0e-4) +template +void check_gradient(serac::Functional& f, double t, mfem::Vector& U, double epsilon = 1.0e-4) { int seed = 42; mfem::Vector dU(U.Size()); + // Set memory backend to device if using GPU execution. + if constexpr (exec == serac::ExecutionSpace::GPU) { + dU.UseDevice(true); + } dU.Randomize(seed); auto [value, dfdU] = f(t, serac::differentiate_wrt(U)); std::unique_ptr dfdU_matrix = assemble(dfdU); // jacobian vector products - mfem::Vector df_jvp1 = dfdU(dU); // matrix-free + mfem::Vector df_jvp1; + df_jvp1.UseDevice(true); + df_jvp1 = dfdU(dU); // matrix-free mfem::Vector df_jvp2(df_jvp1.Size()); + if constexpr (exec == serac::ExecutionSpace::GPU) { + df_jvp2.UseDevice(true); + } dfdU_matrix->Mult(dU, df_jvp2); // sparse matvec if (df_jvp1.Norml2() != 0) { @@ -214,8 +224,8 @@ void check_gradient(serac::Functional& f, double t, const mfem::Vector& U, co // qoi overloads // /////////////////// -template -void check_gradient(serac::Functional& f, double t, const mfem::HypreParVector& U) +template +void check_gradient(serac::Functional& f, double t, const mfem::HypreParVector& U) { int seed = 42; diff --git a/src/serac/numerics/functional/tests/functional_basic_h1_scalar.cpp b/src/serac/numerics/functional/tests/functional_basic_h1_scalar.cpp index 09d63a8bf..7a4fc6118 100644 --- a/src/serac/numerics/functional/tests/functional_basic_h1_scalar.cpp +++ b/src/serac/numerics/functional/tests/functional_basic_h1_scalar.cpp @@ -6,9 +6,8 @@ #include #include - +#include "serac/infrastructure/accelerator.hpp" #include "mfem.hpp" - #include #include "axom/slic/core/SimpleLogger.hpp" @@ -24,6 +23,12 @@ using namespace serac; using namespace serac::profiling; +#ifdef SERAC_USE_CUDA_KERNEL_EVALUATION +constexpr auto exec_space = serac::ExecutionSpace::GPU; +#else +constexpr auto exec_space = serac::ExecutionSpace::CPU; +#endif + template struct TestThermalModelOne { template @@ -71,13 +76,14 @@ void thermal_test_impl(std::unique_ptr& mesh) U_gf.GetTrueDofs(U); // Construct the new functional object using the known test and trial spaces - Functional residual(test_fespace.get(), {trial_fespace.get()}); + Functional residual(test_fespace.get(), {trial_fespace.get()}); residual.AddDomainIntegral(Dimension{}, DependsOn<0>{}, TestThermalModelOne{}, *mesh); residual.AddBoundaryIntegral(Dimension{}, DependsOn<0>{}, TestThermalModelTwo{}, *mesh); double t = 0.0; + check_gradient(residual, t, U); } @@ -95,19 +101,21 @@ void thermal_test(std::string meshfile) } } +TEST(basic, thermal_hexes) { thermal_test<1, 1>("/data/meshes/patch3D_hexes.mesh"); } +#ifndef SERAC_USE_CUDA_KERNEL_EVALUATION TEST(basic, thermal_tris) { thermal_test<1, 1>("/data/meshes/patch2D_tris.mesh"); } +TEST(basic, thermal_tets) { thermal_test<1, 1>("/data/meshes/patch3D_tets.mesh"); } TEST(basic, thermal_quads) { thermal_test<1, 1>("/data/meshes/patch2D_quads.mesh"); } TEST(basic, thermal_tris_and_quads) { thermal_test<1, 1>("/data/meshes/patch2D_tris_and_quads.mesh"); } -TEST(basic, thermal_tets) { thermal_test<1, 1>("/data/meshes/patch3D_tets.mesh"); } -TEST(basic, thermal_hexes) { thermal_test<1, 1>("/data/meshes/patch3D_hexes.mesh"); } TEST(basic, thermal_tets_and_hexes) { thermal_test<1, 1>("/data/meshes/patch3D_tets_and_hexes.mesh"); } - TEST(mixed, thermal_tris_and_quads) { thermal_test<2, 1>("/data/meshes/patch2D_tris_and_quads.mesh"); } TEST(mixed, thermal_tets_and_hexes) { thermal_test<2, 1>("/data/meshes/patch3D_tets_and_hexes.mesh"); } +#endif int main(int argc, char* argv[]) { + serac::accelerator::initializeDevice(); ::testing::InitGoogleTest(&argc, argv); int num_procs, myid; @@ -120,6 +128,7 @@ int main(int argc, char* argv[]) int result = RUN_ALL_TESTS(); MPI_Finalize(); + serac::accelerator::terminateDevice(); return result; } diff --git a/src/serac/numerics/functional/tests/functional_basic_h1_vector.cpp b/src/serac/numerics/functional/tests/functional_basic_h1_vector.cpp index 1d2287751..a69aefec6 100644 --- a/src/serac/numerics/functional/tests/functional_basic_h1_vector.cpp +++ b/src/serac/numerics/functional/tests/functional_basic_h1_vector.cpp @@ -24,6 +24,12 @@ using namespace serac; using namespace serac::profiling; +#ifdef SERAC_USE_CUDA_KERNEL_EVALUATION +constexpr auto exec_space = serac::ExecutionSpace::GPU; +#else +constexpr auto exec_space = serac::ExecutionSpace::CPU; +#endif + template struct MixedModelOne { template @@ -89,9 +95,9 @@ void weird_mixed_test(std::unique_ptr& mesh) auto [test_fes, test_col] = generateParFiniteElementSpace(mesh.get()); mfem::Vector U(trial_fes->TrueVSize()); - U.Randomize(); - Functional residual(test_fes.get(), {trial_fes.get()}); + Functional residual(test_fes.get(), {trial_fes.get()}); + U.Randomize(); // note: this is not really an elasticity problem, it's testing source and flux // terms that have the appropriate shapes to ensure that all the differentiation @@ -117,7 +123,7 @@ void elasticity_test(std::unique_ptr& mesh) mfem::Vector U(trial_fes->TrueVSize()); U.Randomize(); - Functional residual(test_fes.get(), {trial_fes.get()}); + Functional residual(test_fes.get(), {trial_fes.get()}); // note: this is not really an elasticity problem, it's testing source and flux // terms that have the appropriate shapes to ensure that all the differentiation @@ -150,19 +156,22 @@ void test_suite(std::string meshfile) } } +TEST(VectorValuedH1, test_suite_hexes) { test_suite("/data/meshes/patch3D_hexes.mesh"); } + +#ifndef SERAC_USE_CUDA_KERNEL_EVALUATION TEST(VectorValuedH1, test_suite_tris) { test_suite("/data/meshes/patch2D_tris.mesh"); } TEST(VectorValuedH1, test_suite_quads) { test_suite("/data/meshes/patch2D_quads.mesh"); } TEST(VectorValuedH1, test_suite_tris_and_quads) { test_suite("/data/meshes/patch2D_tris_and_quads.mesh"); } - TEST(VectorValuedH1, test_suite_tets) { test_suite("/data/meshes/patch3D_tets.mesh"); } -TEST(VectorValuedH1, test_suite_hexes) { test_suite("/data/meshes/patch3D_hexes.mesh"); } TEST(VectorValuedH1, test_suite_tets_and_hexes) { test_suite("/data/meshes/patch3D_tets_and_hexes.mesh"); } +#endif int main(int argc, char* argv[]) { int num_procs, myid; ::testing::InitGoogleTest(&argc, argv); + serac::accelerator::initializeDevice(); MPI_Init(&argc, &argv); MPI_Comm_size(MPI_COMM_WORLD, &num_procs); @@ -173,6 +182,6 @@ int main(int argc, char* argv[]) int result = RUN_ALL_TESTS(); MPI_Finalize(); - + serac::accelerator::terminateDevice(); return result; } diff --git a/src/serac/numerics/functional/tests/functional_comparison_L2.cpp b/src/serac/numerics/functional/tests/functional_comparison_L2.cpp index cf1906092..3678fa213 100644 --- a/src/serac/numerics/functional/tests/functional_comparison_L2.cpp +++ b/src/serac/numerics/functional/tests/functional_comparison_L2.cpp @@ -25,6 +25,12 @@ constexpr bool verbose = true; std::unique_ptr mesh2D; std::unique_ptr mesh3D; +#ifdef SERAC_USE_CUDA_KERNEL_EVALUATION +constexpr auto exec_space = serac::ExecutionSpace::GPU; +#else +constexpr auto exec_space = serac::ExecutionSpace::CPU; +#endif + template struct hcurl_qfunction { template @@ -100,7 +106,7 @@ void functional_test(mfem::ParMesh& mesh, L2

test, L2

trial, Dimension residual(fespace.get(), {fespace.get()}); + Functional residual(fespace.get(), {fespace.get()}); // Add the total domain residual term to the weak form residual.AddDomainIntegral(Dimension{}, DependsOn<0>{}, test_qfunction{}, mesh); diff --git a/src/serac/numerics/functional/tests/functional_comparisons.cpp b/src/serac/numerics/functional/tests/functional_comparisons.cpp index 17bb76a67..99f4ae9a2 100644 --- a/src/serac/numerics/functional/tests/functional_comparisons.cpp +++ b/src/serac/numerics/functional/tests/functional_comparisons.cpp @@ -10,6 +10,7 @@ #include "axom/slic/core/SimpleLogger.hpp" #include "mfem.hpp" +#include "serac/infrastructure/accelerator.hpp" #include "serac/infrastructure/input.hpp" #include "serac/serac_config.hpp" #include "serac/mesh/mesh_utils_base.hpp" @@ -17,6 +18,7 @@ #include "serac/numerics/functional/functional.hpp" #include "serac/numerics/functional/tensor.hpp" #include +#include "serac/infrastructure/debug_print.hpp" using namespace serac; @@ -30,6 +32,12 @@ std::unique_ptr mesh3D; static constexpr double a = 1.7; static constexpr double b = 2.1; +#ifdef SERAC_USE_CUDA_KERNEL_EVALUATION +constexpr auto exec_space = serac::ExecutionSpace::GPU; +#else +constexpr auto exec_space = serac::ExecutionSpace::CPU; +#endif + template struct hcurl_qfunction { template @@ -42,6 +50,32 @@ struct hcurl_qfunction { } }; +template +struct StressFunctor { + template + SERAC_HOST_DEVICE auto operator()(double /*t*/, Position, Displacement displacement) const + { + static constexpr auto I = DenseIdentity(); + auto [u, du_dx] = displacement; + auto body_force = a * u + I[0]; + auto strain = 0.5 * (du_dx + transpose(du_dx)); + auto stress = b * tr(strain) * I + 2.0 * b * strain; + return serac::tuple{body_force, stress}; + } +}; + +struct SourceFluxFunctor { + template + SERAC_HOST_DEVICE auto operator()(double /*t*/, Position position, Temperature temperature) const + { + auto [X, dX_dxi] = position; + auto [u, du_dX] = temperature; + auto source = a * u - (100 * X[0] * X[1]); + auto flux = b * du_dX; + return serac::tuple{source, flux}; + } +}; + // this test sets up a toy "thermal" problem where the residual includes contributions // from a temperature-dependent source term and a temperature-gradient-dependent flux // @@ -102,21 +136,15 @@ void functional_test(mfem::ParMesh& mesh, H1

test, H1

trial, Dimension residual(fespace.get(), {fespace.get()}); + Functional residual(fespace.get(), {fespace.get()}); // Add the total domain residual term to the functional - residual.AddDomainIntegral( - Dimension{}, DependsOn<0>{}, - [=](double /*t*/, auto position, auto temperature) { - // get the value and the gradient from the input tuple - auto [X, dX_dxi] = position; - auto [u, du_dX] = temperature; - auto source = a * u - (100 * X[0] * X[1]); - auto flux = b * du_dX; - return serac::tuple{source, flux}; - }, - mesh); + residual.AddDomainIntegral(Dimension{}, DependsOn<0>{}, SourceFluxFunctor{}, mesh); // Compute the residual using standard MFEM methods // mfem::Vector r1 = (*J_mfem) * U - (*F); @@ -229,20 +257,13 @@ void functional_test(mfem::ParMesh& mesh, H1 test, H1 trial, Dim mfem::Vector U(fespace->TrueVSize()); u_global.GetTrueDofs(U); + using test_space = decltype(test); + using trial_space = decltype(trial); [[maybe_unused]] static constexpr auto I = DenseIdentity(); - Functional residual(fespace.get(), {fespace.get()}); + Functional residual(fespace.get(), {fespace.get()}); - residual.AddDomainIntegral( - Dimension{}, DependsOn<0>{}, - [=](double /*t*/, auto /*x*/, auto displacement) { - auto [u, du_dx] = displacement; - auto body_force = a * u + I[0]; - auto strain = 0.5 * (du_dx + transpose(du_dx)); - auto stress = b * tr(strain) * I + 2.0 * b * strain; - return serac::tuple{body_force, stress}; - }, - mesh); + residual.AddDomainIntegral(Dimension{}, DependsOn<0>{}, StressFunctor{}, mesh); // mfem::Vector r1 = (*J_mfem) * U - (*F); mfem::Vector r1(U.Size()); @@ -351,7 +372,10 @@ void functional_test(mfem::ParMesh& mesh, Hcurl

test, Hcurl

trial, Dimensi mfem::Vector U(fespace->TrueVSize()); u_global.GetTrueDofs(U); - Functional residual(fespace.get(), {fespace.get()}); + using test_space = decltype(test); + using trial_space = decltype(trial); + + Functional residual(fespace.get(), {fespace.get()}); residual.AddDomainIntegral(Dimension{}, DependsOn<0>{}, hcurl_qfunction{}, mesh); @@ -403,14 +427,23 @@ void functional_test(mfem::ParMesh& mesh, Hcurl

test, Hcurl

trial, Dimensi EXPECT_NEAR(0., diff2.Norml2() / g1.Norml2(), 1.e-13); } +#ifndef SERAC_USE_CUDA_KERNEL_EVALUATION TEST(Thermal, 2DLinear) { functional_test(*mesh2D, H1<1>{}, H1<1>{}, Dimension<2>{}); } TEST(Thermal, 2DQuadratic) { functional_test(*mesh2D, H1<2>{}, H1<2>{}, Dimension<2>{}); } TEST(Thermal, 2DCubic) { functional_test(*mesh2D, H1<3>{}, H1<3>{}, Dimension<2>{}); } +TEST(Elasticity, 2DLinear) { functional_test(*mesh2D, H1<1, 2>{}, H1<1, 2>{}, Dimension<2>{}); } +TEST(Elasticity, 2DQuadratic) { functional_test(*mesh2D, H1<2, 2>{}, H1<2, 2>{}, Dimension<2>{}); } +TEST(Elasticity, 2DCubic) { functional_test(*mesh2D, H1<3, 2>{}, H1<3, 2>{}, Dimension<2>{}); } +#endif TEST(Thermal, 3DLinear) { functional_test(*mesh3D, H1<1>{}, H1<1>{}, Dimension<3>{}); } TEST(Thermal, 3DQuadratic) { functional_test(*mesh3D, H1<2>{}, H1<2>{}, Dimension<3>{}); } TEST(Thermal, 3DCubic) { functional_test(*mesh3D, H1<3>{}, H1<3>{}, Dimension<3>{}); } +TEST(Elasticity, 3DLinear) { functional_test(*mesh3D, H1<1, 3>{}, H1<1, 3>{}, Dimension<3>{}); } +TEST(Elasticity, 3DQuadratic) { functional_test(*mesh3D, H1<2, 3>{}, H1<2, 3>{}, Dimension<3>{}); } +TEST(Elasticity, 3DCubic) { functional_test(*mesh3D, H1<3, 3>{}, H1<3, 3>{}, Dimension<3>{}); } + // TODO: reenable these once hcurl implements of simplex elements is finished // TEST(Hcurl, 2DLinear) { functional_test(*mesh2D, Hcurl<1>{}, Hcurl<1>{}, Dimension<2>{}); } // TEST(Hcurl, 2DQuadratic) { functional_test(*mesh2D, Hcurl<2>{}, Hcurl<2>{}, Dimension<2>{}); } @@ -420,16 +453,9 @@ TEST(Thermal, 3DCubic) { functional_test(*mesh3D, H1<3>{}, H1<3>{}, Dimension<3> // TEST(Hcurl, 3DQuadratic) { functional_test(*mesh3D, Hcurl<2>{}, Hcurl<2>{}, Dimension<3>{}); } // TEST(Hcurl, 3DCubic) { functional_test(*mesh3D, Hcurl<3>{}, Hcurl<3>{}, Dimension<3>{}); } -TEST(Elasticity, 2DLinear) { functional_test(*mesh2D, H1<1, 2>{}, H1<1, 2>{}, Dimension<2>{}); } -TEST(Elasticity, 2DQuadratic) { functional_test(*mesh2D, H1<2, 2>{}, H1<2, 2>{}, Dimension<2>{}); } -TEST(Elasticity, 2DCubic) { functional_test(*mesh2D, H1<3, 2>{}, H1<3, 2>{}, Dimension<2>{}); } - -TEST(Elasticity, 3DLinear) { functional_test(*mesh3D, H1<1, 3>{}, H1<1, 3>{}, Dimension<3>{}); } -TEST(Elasticity, 3DQuadratic) { functional_test(*mesh3D, H1<2, 3>{}, H1<2, 3>{}, Dimension<3>{}); } -TEST(Elasticity, 3DCubic) { functional_test(*mesh3D, H1<3, 3>{}, H1<3, 3>{}, Dimension<3>{}); } - int main(int argc, char* argv[]) { + serac::accelerator::initializeDevice(); ::testing::InitGoogleTest(&argc, argv); MPI_Init(&argc, &argv); MPI_Comm_size(MPI_COMM_WORLD, &num_procs); @@ -468,5 +494,6 @@ int main(int argc, char* argv[]) int result = RUN_ALL_TESTS(); MPI_Finalize(); + serac::accelerator::terminateDevice(); return result; } diff --git a/src/serac/numerics/functional/tests/functional_qoi.cpp b/src/serac/numerics/functional/tests/functional_qoi.cpp index ddace8569..50808c267 100644 --- a/src/serac/numerics/functional/tests/functional_qoi.cpp +++ b/src/serac/numerics/functional/tests/functional_qoi.cpp @@ -11,6 +11,7 @@ #include #include "mfem.hpp" +#include "serac/infrastructure/accelerator.hpp" #include "serac/serac_config.hpp" #include "serac/mesh/mesh_utils_base.hpp" #include "serac/numerics/functional/functional.hpp" diff --git a/src/serac/physics/benchmarks/CMakeLists.txt b/src/serac/physics/benchmarks/CMakeLists.txt index 7eb946021..137a89978 100644 --- a/src/serac/physics/benchmarks/CMakeLists.txt +++ b/src/serac/physics/benchmarks/CMakeLists.txt @@ -26,6 +26,7 @@ foreach(physics_benchmark ${physics_benchmark_targets}) blt_add_benchmark(NAME ${physics_benchmark}_${task_count}_task_count COMMAND ${physics_benchmark} NUM_MPI_TASKS ${task_count} + WORKING_DIRECTORY ${PROJECT_BINARY_DIR} ) endforeach() endforeach() diff --git a/src/serac/physics/boundary_conditions/boundary_condition.cpp b/src/serac/physics/boundary_conditions/boundary_condition.cpp index a25936baf..4a216c78c 100644 --- a/src/serac/physics/boundary_conditions/boundary_condition.cpp +++ b/src/serac/physics/boundary_conditions/boundary_condition.cpp @@ -53,21 +53,6 @@ void BoundaryCondition::setTrueDofList(const mfem::Array& true_dofs) mfem::FiniteElementSpace::MarkerToList(local_dof_marker, local_dofs_); } -void BoundaryCondition::setLocalDofList(const mfem::Array& local_dofs) -{ - local_dofs_ = local_dofs; - - // Create a marker arrays for the true and local dofs - mfem::Array true_dof_marker(space_.GetTrueVSize()); - mfem::Array local_dof_marker(space_.GetVSize()); - - mfem::FiniteElementSpace::ListToMarker(local_dofs_, space_.GetVSize(), local_dof_marker); - - space_.GetRestrictionMatrix()->BooleanMult(local_dof_marker, true_dof_marker); - - mfem::FiniteElementSpace::MarkerToList(true_dof_marker, true_dofs_); -} - void BoundaryCondition::setDofListsFromAttributeMarkers() { auto& mutable_space = const_cast(space_); @@ -133,48 +118,4 @@ void BoundaryCondition::apply(mfem::HypreParMatrix& k_mat, mfem::Vector& rhs, mf mfem::EliminateBC(k_mat, *eliminated_entries, true_dofs_, state, rhs); } -const mfem::Coefficient& BoundaryCondition::scalarCoefficient() const -{ - auto scalar_coef = get_if>(&coef_); - if (scalar_coef) { - return **scalar_coef; - } else { - SLIC_ERROR_ROOT("Asking for a scalar coefficient on a BoundaryCondition that contains a vector coefficient."); - exit(1); - } -} - -mfem::Coefficient& BoundaryCondition::scalarCoefficient() -{ - auto scalar_coef = get_if>(&coef_); - if (scalar_coef) { - return **scalar_coef; - } else { - SLIC_ERROR_ROOT("Asking for a scalar coefficient on a BoundaryCondition that contains a vector coefficient."); - exit(1); - } -} - -const mfem::VectorCoefficient& BoundaryCondition::vectorCoefficient() const -{ - auto vec_coef = get_if>(&coef_); - if (vec_coef) { - return **vec_coef; - } else { - SLIC_ERROR_ROOT("Asking for a vector coefficient on a BoundaryCondition that contains a scalar coefficient."); - exit(1); - } -} - -mfem::VectorCoefficient& BoundaryCondition::vectorCoefficient() -{ - auto vec_coef = get_if>(&coef_); - if (vec_coef) { - return **vec_coef; - } else { - SLIC_ERROR_ROOT("Asking for a vector coefficient on a BoundaryCondition that contains a scalar coefficient."); - exit(1); - } -} - } // namespace serac diff --git a/src/serac/physics/boundary_conditions/boundary_condition.hpp b/src/serac/physics/boundary_conditions/boundary_condition.hpp index 2ecd19453..a09cea690 100644 --- a/src/serac/physics/boundary_conditions/boundary_condition.hpp +++ b/src/serac/physics/boundary_conditions/boundary_condition.hpp @@ -93,46 +93,6 @@ class BoundaryCondition { */ mfem::Array& markers() { return attr_markers_; } - /** - * @brief Accessor for the underlying vector coefficient - * - * This method performs an internal check to verify the underlying GeneralCoefficient - * is in fact a vector. - * - * @return A non-owning reference to the underlying vector coefficient - */ - const mfem::VectorCoefficient& vectorCoefficient() const; - - /** - * @brief Accessor for the underlying vector coefficient - * - * This method performs an internal check to verify the underlying GeneralCoefficient - * is in fact a vector. - * - * @return A non-owning reference to the underlying vector coefficient - */ - mfem::VectorCoefficient& vectorCoefficient(); - - /** - * @brief Accessor for the underlying scalar coefficient - * - * This method performs an internal check to verify the underlying GeneralCoefficient - * is in fact a scalar. - * - * @return A non-owning reference to the underlying scalar coefficient - */ - const mfem::Coefficient& scalarCoefficient() const; - - /** - * @brief Accessor for the underlying scalar coefficient - * - * This method performs an internal check to verify the underlying GeneralCoefficient - * is in fact a scalar. - * - * @return A non-owning reference to the underlying scalar coefficient - */ - mfem::Coefficient& scalarCoefficient(); - /** * @brief Returns the DOF indices for an essential boundary condition * @return A non-owning reference to the array of indices @@ -200,16 +160,6 @@ class BoundaryCondition { */ void setTrueDofList(const mfem::Array& true_dofs); - /** - * @brief "Manually" set the DOF indices without specifying the field to which they apply - * @param[in] local_dofs The local (finite element/grid function) indices of the DOFs constrained by the boundary - * condition - * - * @note This will set both the true and local internal dof index arrays. - * @note True and local dofs are described in the MFEM documentation - */ - void setLocalDofList(const mfem::Array& local_dofs); - /** * @brief A coefficient containing either a mfem::Coefficient or an mfem::VectorCoefficient */ diff --git a/src/serac/physics/boundary_conditions/boundary_condition_manager.cpp b/src/serac/physics/boundary_conditions/boundary_condition_manager.cpp index e23e84204..512c0baa1 100644 --- a/src/serac/physics/boundary_conditions/boundary_condition_manager.cpp +++ b/src/serac/physics/boundary_conditions/boundary_condition_manager.cpp @@ -13,13 +13,6 @@ namespace serac { -void BoundaryConditionManager::addNatural(const std::set& nat_bdr, serac::GeneralCoefficient nat_bdr_coef, - mfem::ParFiniteElementSpace& space, const std::optional component) -{ - nat_bdr_.emplace_back(nat_bdr_coef, component, space, nat_bdr); - all_dofs_valid_ = false; -} - void BoundaryConditionManager::addEssential(const std::set& ess_bdr, serac::GeneralCoefficient ess_bdr_coef, mfem::ParFiniteElementSpace& space, const std::optional component) { diff --git a/src/serac/physics/boundary_conditions/boundary_condition_manager.hpp b/src/serac/physics/boundary_conditions/boundary_condition_manager.hpp index a8e5c5348..9f128d8da 100644 --- a/src/serac/physics/boundary_conditions/boundary_condition_manager.hpp +++ b/src/serac/physics/boundary_conditions/boundary_condition_manager.hpp @@ -20,136 +20,6 @@ namespace serac { -/** - * @brief A "view" for filtering a container - * @note Will be made obsolete by C++20 - * @see std::ranges::views::filter - */ -template -class FilterView { -public: - /** - * @brief An iterator over a filtered view - */ - class FilterViewIterator { - public: - /** - * @brief Constructs a new iterator object - * @param[in] curr The element in the container that should be initially "pointed to" - * @param[in] end The element "one past the end" of the container - * @param[in] pred The predicate to filter with - */ - FilterViewIterator(Iter curr, Iter end, const Pred& pred) : curr_(curr), end_(end), pred_(pred) {} - - /** - * @brief Advances the pointed-to container element to the next element that - * satisfies the predicate - */ - FilterViewIterator& operator++() - { - // Move forward once to advance the element, then continue - // advancing until a predicate-satisfying element is found - ++curr_; - while ((curr_ != end_) && (!pred_(*curr_))) { - ++curr_; - } - return *this; - } - - /** - * @brief Dereferences the iterator - * @return A non-owning reference to the pointed-to element - */ - const auto& operator*() const { return *curr_; } - - /** - * @overload - */ - auto& operator*() { return *curr_; } - - /** - * @brief Comparison operation, checks for iterator inequality - */ - bool operator!=(const FilterViewIterator& other) const { return curr_ != other.curr_; } - - private: - /** - * @brief The currently pointed to element - */ - Iter curr_; - - /** - * @brief One past the last element of the container, used for bounds checking - */ - Iter end_; - - /** - * @brief A reference for the predicate to filter with - */ - const Pred& pred_; - }; - - /** - * @brief Constructs a new lazily-evaluated filtering view over a container - * @param[in] begin The begin() iterator to the container - * @param[in] end The end() iterator to the container - * @param[in] pred The predicate for the filter - */ - FilterView(Iter begin, Iter end, Pred&& pred) : begin_(begin), end_(end), pred_(std::move(pred)) - { - // Skip to the first element that meets the predicate, making sure not to deref the "end" iterator - while ((begin_ != end_) && (!pred_(*begin_))) { - ++begin_; - } - } - - /** - * @brief Returns the first filtered element, i.e., the first element in the - * underlying container that satisfies the predicate - */ - FilterViewIterator begin() { return FilterViewIterator(begin_, end_, pred_); } - - /** - * @overload - */ - const FilterViewIterator begin() const { return FilterViewIterator(begin_, end_, pred_); } - - /** - * @brief Returns one past the end of the container, primarily for bounds-checking - */ - FilterViewIterator end() { return FilterViewIterator(end_, end_, pred_); } - - /** - * @overload - */ - const FilterViewIterator end() const { return FilterViewIterator(end_, end_, pred_); } - -private: - /** - * @brief begin() iterator to the underlying container - */ - Iter begin_; - - /** - * @brief end() iterator to the underlying container - */ - Iter end_; - - /** - * @brief Predicate to filter with - */ - Pred pred_; -}; - -/** - * @brief Deduction guide - iterator and lambda types must be deduced, so this mitigates a "builder" function - * - * @tparam Iter Iterator for the view - * @tparam Pred Predicate for the view - */ -template -FilterView(Iter, Iter, Pred&&) -> FilterView; - /** * @brief A container for the boundary condition information relating to a specific physics module */ @@ -162,17 +32,6 @@ class BoundaryConditionManager { */ explicit BoundaryConditionManager(const mfem::ParMesh& mesh) : num_attrs_(mesh.bdr_attributes.Max()) {} - /** - * @brief Set the natural boundary conditions from a list of boundary markers and a coefficient - * - * @param[in] nat_bdr The set of mesh attributes denoting a natural boundary - * @param[in] nat_bdr_coef The coefficient defining the natural boundary function - * @param[in] space The finite element space to which the BC should be applied - * @param[in] component The component to set (null implies all components are set) - */ - void addNatural(const std::set& nat_bdr, serac::GeneralCoefficient nat_bdr_coef, - mfem::ParFiniteElementSpace& space, const std::optional component = {}); - /** * @brief Set the essential boundary conditions from a list of boundary markers and a coefficient * @@ -198,26 +57,6 @@ class BoundaryConditionManager { void addEssential(const mfem::Array& true_dofs, std::shared_ptr ess_bdr_coef, mfem::ParFiniteElementSpace& space); - /** - * @brief Set a generic boundary condition from a list of boundary markers and a coefficient - * - * @tparam The type of the tag to use - * @param[in] bdr_attr The set of mesh attributes denoting a natural boundary - * @param[in] bdr_coef The coefficient defining the natural boundary function - * @param[in] tag The tag for the generic boundary condition, for identification purposes - * @param[in] space The finite element space to which the BC should be applied - * @param[in] component The component to set (null implies all components are set) - * @pre Template type "Tag" must be an enumeration - */ - template - void addGeneric(const std::set& bdr_attr, serac::GeneralCoefficient bdr_coef, const Tag tag, - mfem::ParFiniteElementSpace& space, const std::optional component = {}) - { - other_bdr_.emplace_back(bdr_coef, component, space, bdr_attr); - other_bdr_.back().setTag(tag); - all_dofs_valid_ = false; - } - /** * @brief Returns all the true degrees of freedom associated with all the essential BCs * @return A const reference to the list of true DOF indices, without duplicates and sorted @@ -258,40 +97,11 @@ class BoundaryConditionManager { * @brief Accessor for the essential BC objects */ std::vector& essentials() { return ess_bdr_; } - /** - * @brief Accessor for the natural BC objects - */ - std::vector& naturals() { return nat_bdr_; } - /** - * @brief Accessor for the generic BC objects - */ - std::vector& generics() { return other_bdr_; } /** * @brief Accessor for the essential BC objects */ const std::vector& essentials() const { return ess_bdr_; } - /** - * @brief Accessor for the natural BC objects - */ - const std::vector& naturals() const { return nat_bdr_; } - /** - * @brief Accessor for the generic BC objects - */ - const std::vector& generics() const { return other_bdr_; } - - /** - * @brief View over all "other"/generic boundary conditions with a specific tag - * @tparam Tag The template type for the tag - * @param tag The tag to filter with - * @pre Tag must be an enumeration type - */ - template - auto genericsWithTag(const Tag tag) - { - static_assert(std::is_enum_v, "Only enumerations can be used to tag a boundary condition."); - return FilterView(other_bdr_.begin(), other_bdr_.end(), [tag](const auto& bc) { return bc.tagEquals(tag); }); - } private: /** @@ -309,16 +119,6 @@ class BoundaryConditionManager { */ std::vector ess_bdr_; - /** - * @brief The vector of natural boundary conditions - */ - std::vector nat_bdr_; - - /** - * @brief The vector of generic (not Dirichlet or Neumann) boundary conditions - */ - std::vector other_bdr_; - /** * @brief The set of boundary attributes associated with * already-registered BCs diff --git a/src/serac/physics/boundary_conditions/tests/boundary_cond.cpp b/src/serac/physics/boundary_conditions/tests/boundary_cond.cpp index 869899770..71ea40c53 100644 --- a/src/serac/physics/boundary_conditions/tests/boundary_cond.cpp +++ b/src/serac/physics/boundary_conditions/tests/boundary_cond.cpp @@ -86,54 +86,6 @@ TEST(BoundaryCond, DirectTrueDofs) } } -enum TestTag -{ - Tag1 = 0, - Tag2 = 1 -}; - -enum OtherTag -{ - Fake1 = 0, - Fake2 = 1 -}; - -TEST(BoundaryCond, FilterGenerics) -{ - MPI_Barrier(MPI_COMM_WORLD); - constexpr int N = 15; - auto mesh = mfem::Mesh::MakeCartesian2D(N, N, mfem::Element::TRIANGLE); - mfem::ParMesh par_mesh(MPI_COMM_WORLD, mesh); - - auto [space, coll] = serac::generateParFiniteElementSpace>(&par_mesh); - - BoundaryConditionManager bcs(par_mesh); - auto coef = std::make_shared(1); - for (int i = 0; i < N; i++) { - bcs.addGeneric({}, coef, TestTag::Tag1, *space, 1); - bcs.addGeneric({}, coef, TestTag::Tag2, *space, 1); - } - - int bcs_with_tag1 = 0; - for (const auto& bc : bcs.genericsWithTag(TestTag::Tag1)) { - EXPECT_TRUE(bc.tagEquals(TestTag::Tag1)); - // Also check that a different enum with the same underlying value will fail - EXPECT_FALSE(bc.tagEquals(OtherTag::Fake1)); - bcs_with_tag1++; - } - EXPECT_EQ(bcs_with_tag1, N); - - int bcs_with_tag2 = 0; - for (const auto& bc : bcs.genericsWithTag(TestTag::Tag2)) { - EXPECT_TRUE(bc.tagEquals(TestTag::Tag2)); - EXPECT_FALSE(bc.tagEquals(OtherTag::Fake2)); - bcs_with_tag2++; - } - EXPECT_EQ(bcs_with_tag2, N); - - MPI_Barrier(MPI_COMM_WORLD); -} - TEST(BoundaryCondHelper, ElementAttributeDofListScalar) { MPI_Barrier(MPI_COMM_WORLD); diff --git a/src/serac/physics/state/CMakeLists.txt b/src/serac/physics/state/CMakeLists.txt index 9115d9b1b..b935113a3 100644 --- a/src/serac/physics/state/CMakeLists.txt +++ b/src/serac/physics/state/CMakeLists.txt @@ -18,6 +18,7 @@ set(state_sources ) set(state_depends serac_infrastructure serac_functional) +blt_list_append(TO state_depends ELEMENTS cuda IF ENABLE_CUDA) blt_add_library( NAME serac_state diff --git a/src/serac/physics/tests/parameterized_thermomechanics_example.cpp b/src/serac/physics/tests/parameterized_thermomechanics_example.cpp index 6886d7ddd..37b086681 100644 --- a/src/serac/physics/tests/parameterized_thermomechanics_example.cpp +++ b/src/serac/physics/tests/parameterized_thermomechanics_example.cpp @@ -37,7 +37,7 @@ struct ParameterizedThermoelasticMaterial { double theta_ref; ///< datum temperature for thermal expansion template - auto operator()(State& /*state*/, const tensor& grad_u, T2 temperature, + auto operator()(const State& /*state*/, const tensor& grad_u, T2 temperature, T3 coefficient_of_thermal_expansion) const { auto theta = get(temperature); diff --git a/src/serac/physics/tests/solid_statics_patch.cpp b/src/serac/physics/tests/solid_statics_patch.cpp index 0066c8b1f..05b4540c4 100644 --- a/src/serac/physics/tests/solid_statics_patch.cpp +++ b/src/serac/physics/tests/solid_statics_patch.cpp @@ -101,11 +101,11 @@ class AffineSolution { /** * @brief Exact displacement solution of the form: - * + * * p == 1: u(X) := A.X + b (affine) - * p == 2: u(X) := A.diag(X).X + b + * p == 2: u(X) := A.diag(X).X + b * p == 3: u(X) := A.diag(X).diag(X).X + b - * + * * @tparam dim number of spatial dimensions */ template @@ -153,7 +153,7 @@ class ManufacturedSolution { template < typename T > auto gradient(const tensor & X) const { - return make_tensor([&](int i, int j){ + return make_tensor([&](int i, int j){ using std::pow; return A(i, j) * p * pow(X[j], p-1); }); @@ -190,7 +190,7 @@ class ManufacturedSolution { typename material_type::State state{}; auto sigma = material(state, H); auto P = solid_mechanics::CauchyToPiola(sigma, H); - return dot(P, n0); + return dot(P, n0); }; sf.setTraction(traction, EntireBoundary(sf.mesh())); @@ -287,7 +287,7 @@ double solution_error(PatchBoundaryCondition bc) constexpr int dim = dimension_of(element_type::geometry); // BT: shouldn't this assertion be in the physics module? - // Putting it here prevents tests from having a nonsensical spatial dimension value, + // Putting it here prevents tests from having a nonsensical spatial dimension value, // but the physics module should be catching this error to protect users. static_assert(dim == 2 || dim == 3, "Dimension must be 2 or 3 for solid test"); @@ -295,10 +295,10 @@ double solution_error(PatchBoundaryCondition bc) // In the future, it would be better to generalize it // so that solution_polynomial_order = p, but my initial // attempt had some issues. I'll revisit it at a later date - // + // // relevant issue: https://github.com/LLNL/serac/issues/926 constexpr int solution_polynomial_order = 1; - auto exact_displacement = ManufacturedSolution(solution_polynomial_order); + auto exact_displacement = ManufacturedSolution(solution_polynomial_order); std::string meshdir = std::string(SERAC_REPO_DIR) + "/data/meshes/"; std::string filename; @@ -308,9 +308,9 @@ double solution_error(PatchBoundaryCondition bc) case mfem::Geometry::TETRAHEDRON: filename = meshdir + "patch3D_tets.mesh"; break; case mfem::Geometry::CUBE: filename = meshdir + "patch3D_hexes.mesh"; break; default: SLIC_ERROR_ROOT("unsupported element type for patch test"); break; - } + } auto mesh = mesh::refineAndDistribute(buildMeshFromFile(filename)); - + std::string mesh_tag{"mesh_tag"}; auto& pmesh = serac::StateManager::setMesh(std::move(mesh), mesh_tag); @@ -383,11 +383,11 @@ double pressure_error() case mfem::Geometry::TETRAHEDRON: filename = meshdir + "patch3D_tets.mesh"; break; case mfem::Geometry::CUBE: filename = meshdir + "patch3D_hexes.mesh"; break; default: SLIC_ERROR_ROOT("unsupported element type for patch test"); break; - } + } auto mesh = mesh::refineAndDistribute(buildMeshFromFile(filename)); - + std::string mesh_tag{"mesh"}; - + auto& pmesh = serac::StateManager::setMesh(std::move(mesh), mesh_tag); // Construct a solid mechanics solver @@ -405,13 +405,13 @@ double pressure_error() solid.setMaterial(mat); typename solid_mechanics::NeoHookean::State state; - auto H = make_tensor([](int i, int j) { + auto H = make_tensor([](int i, int j) { if ( i == 0 && j == 0) { return -0.1; } return 0.0; }); - + tensor sigma = mat(state, H); auto P = solid_mechanics::CauchyToPiola(sigma, H); double pressure = -1.0 * P(0,0); @@ -429,7 +429,7 @@ double pressure_error() solid.setDisplacementBCs({1,3}, [](auto&){return 0.0; }, 1); } else if (element_type::geometry == mfem::Geometry::SQUARE) { solid.setDisplacementBCs({1}, exact_uniaxial_strain); - solid.setDisplacementBCs({2,4}, [](auto&){return 0.0; }, 1); + solid.setDisplacementBCs({2,4}, [](auto&){return 0.0; }, 1); } } else if (dim == 3) { solid.setDisplacementBCs({1}, exact_uniaxial_strain); diff --git a/tests b/tests index 3dd50b00c..ba34e958e 160000 --- a/tests +++ b/tests @@ -1 +1 @@ -Subproject commit 3dd50b00cc8f5b5dad308c75a9fb7102272aa481 +Subproject commit ba34e958ec41f4a8d5d7197056cdb4acb38dd734