From 260ee5b2c59208dcfd9ab89e92cece16e87ae0c1 Mon Sep 17 00:00:00 2001 From: EB Chin Date: Tue, 3 Jun 2025 12:21:49 -0700 Subject: [PATCH 01/71] initial changes --- cmake/thirdparty/SetupSeracThirdParty.cmake | 140 ++++++++-------- ...x86_64_ib_cray-rocmcc@6.2.1_hip_dfem.cmake | 149 ++++++++++++++++++ ...hippet-toss_4_x86_64_ib-clang@19.1.7.cmake | 143 +++++++++++++++++ mfem | 2 +- .../spack/configs/toss_4_x86_64_ib/spack.yaml | 17 ++ scripts/spack/specs.json | 3 +- src/serac/physics/contact/contact_data.hpp | 4 + 7 files changed, 391 insertions(+), 67 deletions(-) create mode 100644 host-configs/rzvernal-toss_4_x86_64_ib_cray-rocmcc@6.2.1_hip_dfem.cmake create mode 100644 host-configs/rzwhippet-toss_4_x86_64_ib-clang@19.1.7.cmake diff --git a/cmake/thirdparty/SetupSeracThirdParty.cmake b/cmake/thirdparty/SetupSeracThirdParty.cmake index d9fba1287..0f5e80d1a 100644 --- a/cmake/thirdparty/SetupSeracThirdParty.cmake +++ b/cmake/thirdparty/SetupSeracThirdParty.cmake @@ -195,6 +195,29 @@ if (NOT SERAC_THIRD_PARTY_LIBRARIES_FOUND) endif() message(STATUS "ARPACK support is ${ARPACK_FOUND}") + #------------------------------------------------------------------------------ + # Enzyme + #------------------------------------------------------------------------------ + if (ENZYME_DIR) + serac_assert_is_directory(DIR_VARIABLE ENZYME_DIR) + set(Enzyme_ROOT ${ENZYME_DIR} CACHE PATH "") + find_dependency(Enzyme REQUIRED) + + # serac_assert_find_succeeded(PROJECT_NAME Enzyme + # TARGET LLDEnzymeFlags + # DIR_VARIABLE ENZYME_DIR) + + serac_assert_find_succeeded(PROJECT_NAME Enzyme + TARGET ClangEnzymeFlags + DIR_VARIABLE ENZYME_DIR) + + message(STATUS "Enzyme support is ON") + set(ENZYME_FOUND TRUE) + else() + message(STATUS "Enzyme support is OFF") + set(ENZYME_FOUND FALSE) + endif() + #------------------------------------------------------------------------------ # MFEM #------------------------------------------------------------------------------ @@ -214,7 +237,7 @@ if (NOT SERAC_THIRD_PARTY_LIBRARIES_FOUND) message(STATUS "Using MFEM submodule") #### Store Data that MFEM clears - set(tpls_to_save ADIAK AMGX AXOM CALIPER CAMP CONDUIT HDF5 + set(tpls_to_save ADIAK AMGX AXOM CALIPER CAMP CONDUIT ENZYME HDF5 HYPRE LUA METIS MFEM NETCDF PARMETIS PETSC RAJA SLEPC SUPERLU_DIST STRUMPACK SUNDIALS TRIBOL UMPIRE) @@ -292,6 +315,10 @@ if (NOT SERAC_THIRD_PARTY_LIBRARIES_FOUND) endif() set(MFEM_USE_UMPIRE OFF CACHE BOOL "") set(MFEM_USE_ZLIB ON CACHE BOOL "") + if(ENZYME_DIR) + serac_assert_is_directory(DIR_VARIABLE ENZYME_DIR) + set(MFEM_USE_ENZYME ON CACHE BOOL "") + endif() #### MFEM Configuration Options @@ -469,78 +496,61 @@ if (NOT SERAC_THIRD_PARTY_LIBRARIES_FOUND) set(ENABLE_FORTRAN ON CACHE BOOL "" FORCE) endif() - #------------------------------------------------------------------------------ - # Enzyme (used by Tribol) - #------------------------------------------------------------------------------ - if (ENZYME_DIR) - serac_assert_is_directory(DIR_VARIABLE ENZYME_DIR) - set(Enzyme_ROOT ${ENZYME_DIR} CACHE PATH "") - find_dependency(Enzyme REQUIRED) - - serac_assert_find_succeeded(PROJECT_NAME Enzyme - TARGET LLDEnzymeFlags - DIR_VARIABLE ENZYME_DIR) - - message(STATUS "Enzyme support is ON") - set(ENZYME_FOUND TRUE) - else() - message(STATUS "Enzyme support is OFF") - set(ENZYME_FOUND FALSE) - endif() - #------------------------------------------------------------------------------ # Tribol #------------------------------------------------------------------------------ - if (NOT SERAC_ENABLE_CODEVELOP) - if(TRIBOL_DIR) - serac_assert_is_directory(DIR_VARIABLE TRIBOL_DIR) + if (NOT SERAC_DISABLE_TRIBOL) + if (NOT SERAC_ENABLE_CODEVELOP) + if(TRIBOL_DIR) + serac_assert_is_directory(DIR_VARIABLE TRIBOL_DIR) - find_dependency(tribol REQUIRED PATHS "${TRIBOL_DIR}/lib/cmake") + find_dependency(tribol REQUIRED PATHS "${TRIBOL_DIR}/lib/cmake") - serac_assert_find_succeeded(PROJECT_NAME Tribol - TARGET tribol - DIR_VARIABLE TRIBOL_DIR) - blt_convert_to_system_includes(TARGET tribol) - set(TRIBOL_FOUND ON) + serac_assert_find_succeeded(PROJECT_NAME Tribol + TARGET tribol + DIR_VARIABLE TRIBOL_DIR) + blt_convert_to_system_includes(TARGET tribol) + set(TRIBOL_FOUND ON) + else() + set(TRIBOL_FOUND OFF) + endif() + + message(STATUS "Tribol support is " ${TRIBOL_FOUND}) else() - set(TRIBOL_FOUND OFF) - endif() - - message(STATUS "Tribol support is " ${TRIBOL_FOUND}) - else() - set(ENABLE_FORTRAN OFF CACHE BOOL "" FORCE) - # Otherwise we use the submodule - message(STATUS "Using Tribol submodule") - set(BUILD_REDECOMP ${SERAC_ENABLE_MPI} CACHE BOOL "") - set(TRIBOL_USE_MPI ${SERAC_ENABLE_MPI} CACHE BOOL "") - set(TRIBOL_ENABLE_TESTS OFF CACHE BOOL "") - set(TRIBOL_ENABLE_EXAMPLES OFF CACHE BOOL "") - set(TRIBOL_ENABLE_DOCS OFF CACHE BOOL "") + set(ENABLE_FORTRAN OFF CACHE BOOL "" FORCE) + # Otherwise we use the submodule + message(STATUS "Using Tribol submodule") + set(BUILD_REDECOMP ${SERAC_ENABLE_MPI} CACHE BOOL "") + set(TRIBOL_USE_MPI ${SERAC_ENABLE_MPI} CACHE BOOL "") + set(TRIBOL_ENABLE_TESTS OFF CACHE BOOL "") + set(TRIBOL_ENABLE_EXAMPLES OFF CACHE BOOL "") + set(TRIBOL_ENABLE_DOCS OFF CACHE BOOL "") + + if(${PROJECT_NAME} STREQUAL "smith") + set(tribol_repo_dir "${PROJECT_SOURCE_DIR}/serac/tribol") + else() + set(tribol_repo_dir "${PROJECT_SOURCE_DIR}/tribol") + endif() - if(${PROJECT_NAME} STREQUAL "smith") - set(tribol_repo_dir "${PROJECT_SOURCE_DIR}/serac/tribol") - else() - set(tribol_repo_dir "${PROJECT_SOURCE_DIR}/tribol") + add_subdirectory(${tribol_repo_dir} ${CMAKE_BINARY_DIR}/tribol) + + target_include_directories(redecomp PUBLIC + $ + ) + target_include_directories(tribol PUBLIC + $ + $ + $ + ) + target_include_directories(tribol_shared PUBLIC + $ + $ + $ + ) + + set(TRIBOL_FOUND TRUE CACHE BOOL "" FORCE) + set(ENABLE_FORTRAN ON CACHE BOOL "" FORCE) endif() - - add_subdirectory(${tribol_repo_dir} ${CMAKE_BINARY_DIR}/tribol) - - target_include_directories(redecomp PUBLIC - $ - ) - target_include_directories(tribol PUBLIC - $ - $ - $ - ) - target_include_directories(tribol_shared PUBLIC - $ - $ - $ - ) - - set(TRIBOL_FOUND TRUE CACHE BOOL "" FORCE) - set(ENABLE_FORTRAN ON CACHE BOOL "" FORCE) endif() #--------------------------------------------------------------------------- diff --git a/host-configs/rzvernal-toss_4_x86_64_ib_cray-rocmcc@6.2.1_hip_dfem.cmake b/host-configs/rzvernal-toss_4_x86_64_ib_cray-rocmcc@6.2.1_hip_dfem.cmake new file mode 100644 index 000000000..f7132a3f6 --- /dev/null +++ b/host-configs/rzvernal-toss_4_x86_64_ib_cray-rocmcc@6.2.1_hip_dfem.cmake @@ -0,0 +1,149 @@ +#------------------------------------------------------------------------------ +# !!!! This is a generated file, edit at own risk !!!! +#------------------------------------------------------------------------------ +# CMake executable path: /usr/tce/bin/cmake +#------------------------------------------------------------------------------ + +set(CMAKE_PREFIX_PATH "/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_43/rocmcc-6.2.1/caliper-2.12.1-h2elmogt2mtorvunf67jrcxqdcrzjjem;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_43/rocmcc-6.2.1/tribol-0.1.0.18-neize4eyd3wxmeuui44az2x5ecqfe2pd;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_43/rocmcc-6.2.1/adiak-0.4.1-djvwhzv7uucbw342y2lxn6mjv4onnzeh;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_43/rocmcc-6.2.1/libunwind-1.8.1-667lulihqw77lixyhoyezlqadqmw3ufq;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_43/rocmcc-6.2.1/axom-0.10.1.1-4udyvzqnaoic74uv3zcumyguv5mbqqcg;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_43/rocmcc-6.2.1/conduit-0.9.2-wo23sn3yekayvxjphe2eu4vvmvkqrovj;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_43/rocmcc-6.2.1/mfem-4.8.0.1-s6vr55m7kq5gylxswlaqy5usltbrslzd;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_43/rocmcc-6.2.1/raja-2024.07.0-m322mdffprdfso4vg7iox3n5goelhtyp;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_43/rocmcc-6.2.1/umpire-2024.07.0-h4aizqln4js5uo6o5dnkkjadcq5apw37;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_43/rocmcc-6.2.1/netcdf-c-4.7.4-m3moxxpcdu3id35a6ib662pmcpe6hlxt;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_43/rocmcc-6.2.1/strumpack-8.0.0-w5kwnpidy7gzm4q2kplzw7ip2hebkojg;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_43/rocmcc-6.2.1/sundials-6.7.0-ihjcn5jozjo6hqxvueljdahsvpg4fpyr;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_43/rocmcc-6.2.1/superlu-dist-8.1.2-ulxns2zfpqyuezhrgxjh7do45ezer2m5;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_43/rocmcc-6.2.1/camp-2024.07.0-223y52kngyoyilxr7gz3i4pe5btnm5uq;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_43/rocmcc-6.2.1/fmt-11.0.2-haismsyfch472jneh7tj5hoz5faug2wp;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_43/rocmcc-6.2.1/hdf5-1.8.23-zhqxrmbx5hqezzazejozk5l7wuz3ijx6;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_43/rocmcc-6.2.1/netlib-scalapack-2.2.2-zinikcerv4cwfrzt4kbhby6xm76fjmvq;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_43/rocmcc-6.2.1/hypre-2.26.0-sffmaz74r5zmnzir6l73deslw36yypsh;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_43/rocmcc-6.2.1/parmetis-4.0.3-7o2cbbpbawncpn24yftcj24gbopz4aev;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_43/rocmcc-6.2.1/metis-5.1.0-q4s4zr5oj2vsyc4eqjy22wrwndgb4fnh;/usr/tce;/usr/tce/packages/cray-mpich-tce/cray-mpich-8.1.29-rocmcc-6.2.1;/opt/rocm-6.2.1;/opt/rocm-6.2.1;/opt/rocm-6.2.1;/opt/rocm-6.2.1/llvm;/opt/rocm-6.2.1;/opt/rocm-6.2.1;/opt/rocm-6.2.1;/opt/rocm-6.2.1" CACHE STRING "") + +set(CMAKE_INSTALL_RPATH_USE_LINK_PATH "ON" CACHE STRING "") + +set(CMAKE_BUILD_RPATH "/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_43/rocmcc-6.2.1/serac-develop-og4fmn5dkubewiybduaji237suiuzssb/lib;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_43/rocmcc-6.2.1/serac-develop-og4fmn5dkubewiybduaji237suiuzssb/lib64;;/opt/rh/gcc-toolset-12/root/usr/lib/gcc/x86_64-redhat-linux/12" CACHE STRING "") + +set(CMAKE_INSTALL_RPATH "/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_43/rocmcc-6.2.1/serac-develop-og4fmn5dkubewiybduaji237suiuzssb/lib;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_43/rocmcc-6.2.1/serac-develop-og4fmn5dkubewiybduaji237suiuzssb/lib64;;/opt/rh/gcc-toolset-12/root/usr/lib/gcc/x86_64-redhat-linux/12" CACHE STRING "") + +set(CMAKE_BUILD_TYPE "Release" CACHE STRING "") + +#------------------------------------------------------------------------------ +# Compilers +#------------------------------------------------------------------------------ +# Compiler Spec: rocmcc@=6.2.1 +#------------------------------------------------------------------------------ +if(DEFINED ENV{SPACK_CC}) + + set(CMAKE_C_COMPILER "/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_43/spack/lib/spack/env/rocmcc/amdclang" CACHE PATH "") + + set(CMAKE_CXX_COMPILER "/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_43/spack/lib/spack/env/rocmcc/amdclang++" CACHE PATH "") + + set(CMAKE_Fortran_COMPILER "/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_43/spack/lib/spack/env/rocmcc/amdflang" CACHE PATH "") + +else() + + set(CMAKE_C_COMPILER "/opt/rocm-6.2.1/llvm/bin/amdclang" CACHE PATH "") + + set(CMAKE_CXX_COMPILER "/opt/rocm-6.2.1/llvm/bin/amdclang++" CACHE PATH "") + + set(CMAKE_Fortran_COMPILER "/opt/rocm-6.2.1/llvm/bin/amdflang" CACHE PATH "") + +endif() + +#------------------------------------------------------------------------------ +# MPI +#------------------------------------------------------------------------------ + +set(MPI_C_COMPILER "/usr/tce/packages/cray-mpich-tce/cray-mpich-8.1.29-rocmcc-6.2.1/bin/mpicc" CACHE PATH "") + +set(MPI_CXX_COMPILER "/usr/tce/packages/cray-mpich-tce/cray-mpich-8.1.29-rocmcc-6.2.1/bin/mpicxx" CACHE PATH "") + +set(MPI_Fortran_COMPILER "/usr/tce/packages/cray-mpich-tce/cray-mpich-8.1.29-rocmcc-6.2.1/bin/mpif90" CACHE PATH "") + +set(MPIEXEC_NUMPROC_FLAG "-n" CACHE STRING "") + +set(ENABLE_MPI ON CACHE BOOL "") + +set(MPIEXEC_EXECUTABLE "/usr/global/tools/flux_wrappers/bin/srun" CACHE PATH "") + +#------------------------------------------------------------------------------ +# Hardware +#------------------------------------------------------------------------------ + +#------------------------------------------------ +# ROCm +#------------------------------------------------ + +set(HIP_ROOT_DIR "/opt/rocm-6.2.1" CACHE PATH "") + +set(CMAKE_HIP_COMPILER "/opt/rocm-6.2.1/llvm/bin/clang++" CACHE FILEPATH "") + +set(CMAKE_HIP_ARCHITECTURES "gfx90a" CACHE STRING "") + +set(AMDGPU_TARGETS "gfx90a" CACHE STRING "") + +set(GPU_TARGETS "gfx90a" CACHE STRING "") + +set(ENABLE_OPENMP OFF CACHE BOOL "") + +set(ENABLE_HIP ON CACHE BOOL "") + +set(ROCM_PATH "/opt/rocm-6.2.1" CACHE PATH "") + +set(CMAKE_EXE_LINKER_FLAGS "-L/opt/rocm-6.2.1/lib -Wl,-rpath,/opt/rocm-6.2.1/lib -L/opt/rocm-6.2.1/llvm/lib -Wl,-rpath,/opt/rocm-6.2.1/llvm/lib -lxpmem -L/opt/cray/pe/mpich/8.1.29/gtl/lib -Wl,-rpath,/opt/cray/pe/mpich/8.1.29/gtl/lib -lmpi_gtl_hsa -Wl,--disable-new-dtags -lflang -lflangrti -lamdhip64 -lhsakmt -lhsa-runtime64 -lamd_comgr -lpgmath -lhipblas" CACHE STRING "") + +#------------------------------------------------------------------------------ +# TPLs +#------------------------------------------------------------------------------ + +set(TPL_ROOT "/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_43/rocmcc-6.2.1" CACHE PATH "") + +set(AXOM_DIR "${TPL_ROOT}/axom-0.10.1.1-4udyvzqnaoic74uv3zcumyguv5mbqqcg" CACHE PATH "") + +set(CAMP_DIR "${TPL_ROOT}/camp-2024.07.0-223y52kngyoyilxr7gz3i4pe5btnm5uq" CACHE PATH "") + +set(CONDUIT_DIR "${TPL_ROOT}/conduit-0.9.2-wo23sn3yekayvxjphe2eu4vvmvkqrovj" CACHE PATH "") + +set(LUA_DIR "/usr" CACHE PATH "") + +set(MFEM_DIR "${TPL_ROOT}/mfem-4.8.0.1-s6vr55m7kq5gylxswlaqy5usltbrslzd" CACHE PATH "") + +set(HDF5_DIR "${TPL_ROOT}/hdf5-1.8.23-zhqxrmbx5hqezzazejozk5l7wuz3ijx6" CACHE PATH "") + +set(HYPRE_DIR "${TPL_ROOT}/hypre-2.26.0-sffmaz74r5zmnzir6l73deslw36yypsh" CACHE PATH "") + +set(METIS_DIR "${TPL_ROOT}/metis-5.1.0-q4s4zr5oj2vsyc4eqjy22wrwndgb4fnh" CACHE PATH "") + +set(PARMETIS_DIR "${TPL_ROOT}/parmetis-4.0.3-7o2cbbpbawncpn24yftcj24gbopz4aev" CACHE PATH "") + +set(NETCDF_DIR "${TPL_ROOT}/netcdf-c-4.7.4-m3moxxpcdu3id35a6ib662pmcpe6hlxt" CACHE PATH "") + +set(SUPERLUDIST_DIR "${TPL_ROOT}/superlu-dist-8.1.2-ulxns2zfpqyuezhrgxjh7do45ezer2m5" CACHE PATH "") + +set(ADIAK_DIR "${TPL_ROOT}/adiak-0.4.1-djvwhzv7uucbw342y2lxn6mjv4onnzeh" CACHE PATH "") + +# AMGX not built + +set(CALIPER_DIR "${TPL_ROOT}/caliper-2.12.1-h2elmogt2mtorvunf67jrcxqdcrzjjem" CACHE PATH "") + +# PETSC not built + +set(RAJA_DIR "${TPL_ROOT}/raja-2024.07.0-m322mdffprdfso4vg7iox3n5goelhtyp" CACHE PATH "") + +# SLEPC not built + +set(STRUMPACK_DIR "${TPL_ROOT}/strumpack-8.0.0-w5kwnpidy7gzm4q2kplzw7ip2hebkojg" CACHE PATH "") + +set(SUNDIALS_DIR "${TPL_ROOT}/sundials-6.7.0-ihjcn5jozjo6hqxvueljdahsvpg4fpyr" CACHE PATH "") + +set(UMPIRE_DIR "${TPL_ROOT}/umpire-2024.07.0-h4aizqln4js5uo6o5dnkkjadcq5apw37" CACHE PATH "") + +set(TRIBOL_DIR "${TPL_ROOT}/tribol-0.1.0.18-neize4eyd3wxmeuui44az2x5ecqfe2pd" CACHE PATH "") + +set(ENZYME_DIR "/usr/WS2/smithdev/toss_4_x86_64_ib_cray/enzyme/rocm-6.2.1/0.0.180" CACHE PATH "") + +#------------------------------------------------------------------------------ +# Devtools +#------------------------------------------------------------------------------ + +# Code checks disabled due to disabled devtools + +set(SERAC_ENABLE_CODE_CHECKS OFF CACHE BOOL "") + +set(ENABLE_CLANGFORMAT OFF CACHE BOOL "") + +set(ENABLE_CLANGTIDY OFF CACHE BOOL "") + +set(ENABLE_DOCS OFF CACHE BOOL "") + +set(SERAC_ENABLE_CODEVELOP ON CACHE BOOL "") + +set(SERAC_DISABLE_TRIBOL ON CACHE BOOL "") + diff --git a/host-configs/rzwhippet-toss_4_x86_64_ib-clang@19.1.7.cmake b/host-configs/rzwhippet-toss_4_x86_64_ib-clang@19.1.7.cmake new file mode 100644 index 000000000..747b9ce1a --- /dev/null +++ b/host-configs/rzwhippet-toss_4_x86_64_ib-clang@19.1.7.cmake @@ -0,0 +1,143 @@ +#------------------------------------------------------------------------------ +# !!!! This is a generated file, edit at own risk !!!! +#------------------------------------------------------------------------------ +# CMake executable path: /usr/tce/bin/cmake +#------------------------------------------------------------------------------ + +set(CMAKE_PREFIX_PATH "/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_05_27_15_56_53/clang-19.1.7/caliper-2.12.1-a4chnmcpoktn7vlyrjxbw2qskkkytkj3;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_05_27_15_56_53/clang-19.1.7/tribol-0.1.0.18-my4bcnwnewxcnewg3xgbzvym64dfizuj;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_05_27_15_56_53/clang-19.1.7/adiak-0.4.1-5wjrv4xwbm3y7yh7745ijc45sbcwqwqi;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_05_27_15_56_53/clang-19.1.7/libunwind-1.8.1-tgmphz65mn62bautt2ksl45ckecw5e2c;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_05_27_15_56_53/clang-19.1.7/axom-0.10.1.1-tkmjd26d72eqk4tts7hwlwognvwuyj5s;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_05_27_15_56_53/clang-19.1.7/conduit-0.9.2-t2zobt76hkwap7d2ppddp3traorxpffl;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_05_27_15_56_53/clang-19.1.7/mfem-4.8.0.1-gbv5tbeosslcqh2zhkcrgkwtqxob257j;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_05_27_15_56_53/clang-19.1.7/raja-2024.07.0-t5ay2skmkw3ybywrq33tcyl5xbn5g3gv;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_05_27_15_56_53/clang-19.1.7/umpire-2024.07.0-47yl3s6ob3ru4g6qnma4tiv6ob5fr6p6;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_05_27_15_56_53/clang-19.1.7/netcdf-c-4.7.4-62ri63jaeyajwrrvwk75dbwvrnj6fony;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_05_27_15_56_53/clang-19.1.7/slepc-3.21.2-rji2oyfq6x4lyovnrkfqnleoc3yrqrm6;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_05_27_15_56_53/clang-19.1.7/sundials-6.7.0-wwas67vc5xfgayl2uimersm3xszlqv62;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_05_27_15_56_53/clang-19.1.7/camp-2024.07.0-afpyaxeo4ulr2lkljpy2zr7rjeq64cox;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_05_27_15_56_53/clang-19.1.7/fmt-11.0.2-xpfumigjok66664dhruf6tfaa7gwcn3n;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_05_27_15_56_53/clang-19.1.7/hdf5-1.8.23-pu3kif4ltqdw2kcwvlswfze6pmmsjs6q;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_05_27_15_56_53/clang-19.1.7/arpack-ng-3.9.0-ahgmmge3jagsx3n72wbxtcznbubcclca;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_05_27_15_56_53/clang-19.1.7/petsc-3.21.6-elp7iw6xe64daspgy5ctubfo3jj4nd43;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_05_27_15_56_53/clang-19.1.7/hypre-2.26.0-ta6n6j6ziyuw2tfejt27nkdulvjq27hu;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_05_27_15_56_53/clang-19.1.7/strumpack-8.0.0-zmfwdjw54im5o6m3gqssnkvuagdzljwe;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_05_27_15_56_53/clang-19.1.7/superlu-dist-8.1.2-k53l7hvjwbqnepeb7etlujfgzrzfliua;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_05_27_15_56_53/clang-19.1.7/netlib-scalapack-2.2.2-fuk7tt7il4e6tekxxctc625g4xk4hqeu;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_05_27_15_56_53/clang-19.1.7/parmetis-4.0.3-morhecieft72ri6gi2r6hsatzrgdhe5n;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_05_27_15_56_53/clang-19.1.7/metis-5.1.0-ebbueyypx5bn5mekybxavligkykhm5fy;/usr/workspace/smithdev/devtools/toss_4_x86_64_ib/latest/python-3.11.7;/usr/workspace/smithdev/devtools/toss_4_x86_64_ib/latest/python-3.11.7;/usr/tce;/usr/workspace/smithdev/devtools/toss_4_x86_64_ib/latest/cppcheck-2.9;/usr/workspace/smithdev/devtools/toss_4_x86_64_ib/latest/doxygen-1.9.8;/usr/tce/packages/clang/clang-14.0.6;/usr/tce/packages/mvapich2/mvapich2-2.3.7-gcc-10.3.1;/usr/tce/packages/python/python-3.9.12" CACHE STRING "") + +set(CMAKE_INSTALL_RPATH_USE_LINK_PATH "ON" CACHE STRING "") + +set(CMAKE_BUILD_RPATH "/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_05_27_15_56_53/clang-19.1.7/serac-develop-t2mywg4ueanbqllbapg5wvrivg3umjoj/lib;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_05_27_15_56_53/clang-19.1.7/serac-develop-t2mywg4ueanbqllbapg5wvrivg3umjoj/lib64;;/usr/tce/packages/gcc/gcc-10.3.1/lib/gcc/x86_64-redhat-linux/10" CACHE STRING "") + +set(CMAKE_INSTALL_RPATH "/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_05_27_15_56_53/clang-19.1.7/serac-develop-t2mywg4ueanbqllbapg5wvrivg3umjoj/lib;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_05_27_15_56_53/clang-19.1.7/serac-develop-t2mywg4ueanbqllbapg5wvrivg3umjoj/lib64;;/usr/tce/packages/gcc/gcc-10.3.1/lib/gcc/x86_64-redhat-linux/10" CACHE STRING "") + +set(CMAKE_BUILD_TYPE "Release" CACHE STRING "") + +#------------------------------------------------------------------------------ +# Compilers +#------------------------------------------------------------------------------ +# Compiler Spec: clang@=19.1.7 +#------------------------------------------------------------------------------ +if(DEFINED ENV{SPACK_CC}) + + set(CMAKE_C_COMPILER "/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_05_27_15_56_53/spack/lib/spack/env/clang/clang" CACHE PATH "") + + set(CMAKE_CXX_COMPILER "/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_05_27_15_56_53/spack/lib/spack/env/clang/clang++" CACHE PATH "") + + set(CMAKE_Fortran_COMPILER "/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_05_27_15_56_53/spack/lib/spack/env/clang/gfortran" CACHE PATH "") + +else() + + set(CMAKE_C_COMPILER "/usr/WS2/smithdev/toss_4_x86_64_ib/llvm/19.1.7/bin/clang" CACHE PATH "") + + set(CMAKE_CXX_COMPILER "/usr/WS2/smithdev/toss_4_x86_64_ib/llvm/19.1.7/bin/clang++" CACHE PATH "") + + set(CMAKE_Fortran_COMPILER "/usr/tce/packages/gcc/gcc-10.3.1/bin/gfortran" CACHE PATH "") + +endif() + +set(CMAKE_C_FLAGS "--gcc-toolchain=/usr/tce/packages/gcc/gcc-10.3.1" CACHE STRING "") + +set(CMAKE_CXX_FLAGS "--gcc-toolchain=/usr/tce/packages/gcc/gcc-10.3.1" CACHE STRING "") + +set(CMAKE_Fortran_FLAGS "-fPIE -fPIC" CACHE STRING "") + +#------------------------------------------------------------------------------ +# MPI +#------------------------------------------------------------------------------ + +set(MPI_C_COMPILER "/usr/tce/packages/mvapich2/mvapich2-2.3.7-gcc-10.3.1/bin/mpicc" CACHE PATH "") + +set(MPI_CXX_COMPILER "/usr/tce/packages/mvapich2/mvapich2-2.3.7-gcc-10.3.1/bin/mpicxx" CACHE PATH "") + +set(MPI_Fortran_COMPILER "/usr/tce/packages/mvapich2/mvapich2-2.3.7-gcc-10.3.1/bin/mpif90" CACHE PATH "") + +set(MPIEXEC_NUMPROC_FLAG "-n" CACHE STRING "") + +set(ENABLE_MPI ON CACHE BOOL "") + +set(MPIEXEC_EXECUTABLE "/usr/bin/srun" CACHE PATH "") + +#------------------------------------------------------------------------------ +# Hardware +#------------------------------------------------------------------------------ + +set(ENABLE_OPENMP OFF CACHE BOOL "") + +#------------------------------------------------------------------------------ +# TPLs +#------------------------------------------------------------------------------ + +set(TPL_ROOT "/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_05_27_15_56_53/clang-19.1.7" CACHE PATH "") + +set(AXOM_DIR "${TPL_ROOT}/axom-0.10.1.1-tkmjd26d72eqk4tts7hwlwognvwuyj5s" CACHE PATH "") + +set(CAMP_DIR "${TPL_ROOT}/camp-2024.07.0-afpyaxeo4ulr2lkljpy2zr7rjeq64cox" CACHE PATH "") + +set(CONDUIT_DIR "${TPL_ROOT}/conduit-0.9.2-t2zobt76hkwap7d2ppddp3traorxpffl" CACHE PATH "") + +set(LUA_DIR "/usr" CACHE PATH "") + +set(MFEM_DIR "${TPL_ROOT}/mfem-4.8.0.1-gbv5tbeosslcqh2zhkcrgkwtqxob257j" CACHE PATH "") + +set(HDF5_DIR "${TPL_ROOT}/hdf5-1.8.23-pu3kif4ltqdw2kcwvlswfze6pmmsjs6q" CACHE PATH "") + +set(HYPRE_DIR "${TPL_ROOT}/hypre-2.26.0-ta6n6j6ziyuw2tfejt27nkdulvjq27hu" CACHE PATH "") + +set(METIS_DIR "${TPL_ROOT}/metis-5.1.0-ebbueyypx5bn5mekybxavligkykhm5fy" CACHE PATH "") + +set(PARMETIS_DIR "${TPL_ROOT}/parmetis-4.0.3-morhecieft72ri6gi2r6hsatzrgdhe5n" CACHE PATH "") + +set(NETCDF_DIR "${TPL_ROOT}/netcdf-c-4.7.4-62ri63jaeyajwrrvwk75dbwvrnj6fony" CACHE PATH "") + +set(SUPERLUDIST_DIR "${TPL_ROOT}/superlu-dist-8.1.2-k53l7hvjwbqnepeb7etlujfgzrzfliua" CACHE PATH "") + +set(ARPACK_DIR "${TPL_ROOT}/arpack-ng-3.9.0-ahgmmge3jagsx3n72wbxtcznbubcclca" CACHE PATH "") + +set(ADIAK_DIR "${TPL_ROOT}/adiak-0.4.1-5wjrv4xwbm3y7yh7745ijc45sbcwqwqi" CACHE PATH "") + +# AMGX not built + +set(CALIPER_DIR "${TPL_ROOT}/caliper-2.12.1-a4chnmcpoktn7vlyrjxbw2qskkkytkj3" CACHE PATH "") + +set(PETSC_DIR "${TPL_ROOT}/petsc-3.21.6-elp7iw6xe64daspgy5ctubfo3jj4nd43" CACHE PATH "") + +set(RAJA_DIR "${TPL_ROOT}/raja-2024.07.0-t5ay2skmkw3ybywrq33tcyl5xbn5g3gv" CACHE PATH "") + +set(SLEPC_DIR "${TPL_ROOT}/slepc-3.21.2-rji2oyfq6x4lyovnrkfqnleoc3yrqrm6" CACHE PATH "") + +set(STRUMPACK_DIR "${TPL_ROOT}/strumpack-8.0.0-zmfwdjw54im5o6m3gqssnkvuagdzljwe" CACHE PATH "") + +set(SUNDIALS_DIR "${TPL_ROOT}/sundials-6.7.0-wwas67vc5xfgayl2uimersm3xszlqv62" CACHE PATH "") + +set(UMPIRE_DIR "${TPL_ROOT}/umpire-2024.07.0-47yl3s6ob3ru4g6qnma4tiv6ob5fr6p6" CACHE PATH "") + +set(TRIBOL_DIR "${TPL_ROOT}/tribol-0.1.0.18-my4bcnwnewxcnewg3xgbzvym64dfizuj" CACHE PATH "") + +set(ENZYME_DIR "/usr/WS2/smithdev/toss_4_x86_64_ib/enzyme/llvm_19.1.7/0.0.180" CACHE PATH "") + +#------------------------------------------------------------------------------ +# Devtools +#------------------------------------------------------------------------------ + +set(DEVTOOLS_ROOT "/usr/WS2/smithdev/devtools/toss_4_x86_64_ib/2024_05_30_15_09_55/._view/rkqkran3ydsuprr2wip5pdnz5wh7xwnr" CACHE PATH "") + +set(ATS_EXECUTABLE "${DEVTOOLS_ROOT}/python-3.11.7/bin/ats" CACHE PATH "") + +set(CLANGFORMAT_EXECUTABLE "/usr/tce/packages/clang/clang-14.0.6/bin/clang-format" CACHE PATH "") + +set(CLANGTIDY_EXECUTABLE "/usr/tce/packages/clang/clang-14.0.6/bin/clang-tidy" CACHE PATH "") + +set(ENABLE_DOCS ON CACHE BOOL "") + +set(SPHINX_EXECUTABLE "${DEVTOOLS_ROOT}/python-3.11.7/bin/sphinx-build" CACHE PATH "") + +set(CPPCHECK_EXECUTABLE "${DEVTOOLS_ROOT}/cppcheck-2.9/bin/cppcheck" CACHE PATH "") + +set(DOXYGEN_EXECUTABLE "${DEVTOOLS_ROOT}/doxygen-1.9.8/bin/doxygen" CACHE PATH "") + +set(SERAC_ENABLE_CODEVELOP ON CACHE BOOL "") + +set(SERAC_DISABLE_TRIBOL ON CACHE BOOL "") + diff --git a/mfem b/mfem index d9c1c34fd..3e3a5e5d1 160000 --- a/mfem +++ b/mfem @@ -1 +1 @@ -Subproject commit d9c1c34fdfaf3f7a9f56dfc82f7c083082a36fca +Subproject commit 3e3a5e5d14003d7856f779292fe032067c10fc1b diff --git a/scripts/spack/configs/toss_4_x86_64_ib/spack.yaml b/scripts/spack/configs/toss_4_x86_64_ib/spack.yaml index dedf5afc0..5b7a7f061 100644 --- a/scripts/spack/configs/toss_4_x86_64_ib/spack.yaml +++ b/scripts/spack/configs/toss_4_x86_64_ib/spack.yaml @@ -58,6 +58,23 @@ spack: fc: /usr/tce/packages/gcc/gcc-10.3.1/bin/gfortran spec: clang@16.0.6 target: x86_64 + - compiler: + environment: {} + extra_rpaths: [] + flags: + cflags: --gcc-toolchain=/usr/tce/packages/gcc/gcc-10.3.1 + cxxflags: --gcc-toolchain=/usr/tce/packages/gcc/gcc-10.3.1 + fflags: -fPIE -fPIC + modules: + - clang/19.1.7 + operating_system: rhel8 + paths: + cc: /usr/WS2/smithdev/toss_4_x86_64_ib/llvm/19.1.7/bin/clang + cxx: /usr/WS2/smithdev/toss_4_x86_64_ib/llvm/19.1.7/bin/clang++ + f77: /usr/tce/packages/gcc/gcc-10.3.1/bin/gfortran + fc: /usr/tce/packages/gcc/gcc-10.3.1/bin/gfortran + spec: clang@19.1.7 + target: x86_64 packages: all: diff --git a/scripts/spack/specs.json b/scripts/spack/specs.json index c128493d4..c5a1fe6a0 100644 --- a/scripts/spack/specs.json +++ b/scripts/spack/specs.json @@ -16,7 +16,8 @@ "toss_4_x86_64_ib": [ "gcc@10.3.1+devtools+profiling", "clang@14.0.6+devtools+profiling", - "clang@16.0.6+devtools+profiling~openmp" ], + "clang@16.0.6+devtools+profiling~openmp", + "clang@19.1.7+devtools+profiling~openmp" ], "toss_4_x86_64_ib_cray": [ "rocmcc@6.2.1~openmp+profiling+rocm+raja+umpire+strumpack~petsc amdgpu_target=gfx90a ^hip@6.2.1 ^rocprim@6.2.1 ^hsa-rocr-dev@6.2.1 ^llvm-amdgpu@6.2.1 ^hdf5 cflags=-Wno-int-conversion"], diff --git a/src/serac/physics/contact/contact_data.hpp b/src/serac/physics/contact/contact_data.hpp index 01335cc92..2bd963d1d 100644 --- a/src/serac/physics/contact/contact_data.hpp +++ b/src/serac/physics/contact/contact_data.hpp @@ -247,6 +247,7 @@ class ContactData { */ int num_pressure_dofs_; +#ifdef SERAC_USE_TRIBOL /** * @brief Tracks whether the Jacobian and pressure offsets need to be updated * @@ -254,6 +255,7 @@ class ContactData { * */ mutable bool offsets_up_to_date_; +#endif /** * @brief Offsets giving size of each block Jacobian contribution @@ -284,9 +286,11 @@ class ContactData { */ mutable mfem::Array global_pressure_dof_offsets_; +#ifdef SERAC_USE_TRIBOL int cycle_{0}; double time_{0.0}; double dt_{1.0}; +#endif }; } // namespace serac From 0ad28904f342ad642c380561f398de1a97bb0043 Mon Sep 17 00:00:00 2001 From: "Eric B. Chin" Date: Tue, 3 Jun 2025 17:06:16 -0700 Subject: [PATCH 02/71] tioga host config --- cmake/thirdparty/SetupSeracThirdParty.cmake | 6 +- ...x86_64_ib_cray-rocmcc@6.2.1_hip_dfem.cmake | 0 ...hippet-toss_4_x86_64_ib-clang@19.1.7.cmake | 0 ...x86_64_ib_cray-rocmcc@6.2.1_hip_dfem.cmake | 149 ++++++++++++++++++ 4 files changed, 152 insertions(+), 3 deletions(-) rename host-configs/{ => dfem}/rzvernal-toss_4_x86_64_ib_cray-rocmcc@6.2.1_hip_dfem.cmake (100%) rename host-configs/{ => dfem}/rzwhippet-toss_4_x86_64_ib-clang@19.1.7.cmake (100%) create mode 100644 host-configs/dfem/tioga-toss_4_x86_64_ib_cray-rocmcc@6.2.1_hip_dfem.cmake diff --git a/cmake/thirdparty/SetupSeracThirdParty.cmake b/cmake/thirdparty/SetupSeracThirdParty.cmake index 0f5e80d1a..8592d1ce6 100644 --- a/cmake/thirdparty/SetupSeracThirdParty.cmake +++ b/cmake/thirdparty/SetupSeracThirdParty.cmake @@ -203,9 +203,9 @@ if (NOT SERAC_THIRD_PARTY_LIBRARIES_FOUND) set(Enzyme_ROOT ${ENZYME_DIR} CACHE PATH "") find_dependency(Enzyme REQUIRED) - # serac_assert_find_succeeded(PROJECT_NAME Enzyme - # TARGET LLDEnzymeFlags - # DIR_VARIABLE ENZYME_DIR) + serac_assert_find_succeeded(PROJECT_NAME Enzyme + TARGET LLDEnzymeFlags + DIR_VARIABLE ENZYME_DIR) serac_assert_find_succeeded(PROJECT_NAME Enzyme TARGET ClangEnzymeFlags diff --git a/host-configs/rzvernal-toss_4_x86_64_ib_cray-rocmcc@6.2.1_hip_dfem.cmake b/host-configs/dfem/rzvernal-toss_4_x86_64_ib_cray-rocmcc@6.2.1_hip_dfem.cmake similarity index 100% rename from host-configs/rzvernal-toss_4_x86_64_ib_cray-rocmcc@6.2.1_hip_dfem.cmake rename to host-configs/dfem/rzvernal-toss_4_x86_64_ib_cray-rocmcc@6.2.1_hip_dfem.cmake diff --git a/host-configs/rzwhippet-toss_4_x86_64_ib-clang@19.1.7.cmake b/host-configs/dfem/rzwhippet-toss_4_x86_64_ib-clang@19.1.7.cmake similarity index 100% rename from host-configs/rzwhippet-toss_4_x86_64_ib-clang@19.1.7.cmake rename to host-configs/dfem/rzwhippet-toss_4_x86_64_ib-clang@19.1.7.cmake diff --git a/host-configs/dfem/tioga-toss_4_x86_64_ib_cray-rocmcc@6.2.1_hip_dfem.cmake b/host-configs/dfem/tioga-toss_4_x86_64_ib_cray-rocmcc@6.2.1_hip_dfem.cmake new file mode 100644 index 000000000..4a75768ea --- /dev/null +++ b/host-configs/dfem/tioga-toss_4_x86_64_ib_cray-rocmcc@6.2.1_hip_dfem.cmake @@ -0,0 +1,149 @@ +#------------------------------------------------------------------------------ +# !!!! This is a generated file, edit at own risk !!!! +#------------------------------------------------------------------------------ +# CMake executable path: /usr/tce/bin/cmake +#------------------------------------------------------------------------------ + +set(CMAKE_PREFIX_PATH "/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_36/rocmcc-6.2.1/caliper-2.12.1-h2elmogt2mtorvunf67jrcxqdcrzjjem;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_36/rocmcc-6.2.1/tribol-0.1.0.18-neize4eyd3wxmeuui44az2x5ecqfe2pd;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_36/rocmcc-6.2.1/adiak-0.4.1-djvwhzv7uucbw342y2lxn6mjv4onnzeh;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_36/rocmcc-6.2.1/libunwind-1.8.1-667lulihqw77lixyhoyezlqadqmw3ufq;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_36/rocmcc-6.2.1/axom-0.10.1.1-4udyvzqnaoic74uv3zcumyguv5mbqqcg;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_36/rocmcc-6.2.1/conduit-0.9.2-wo23sn3yekayvxjphe2eu4vvmvkqrovj;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_36/rocmcc-6.2.1/mfem-4.8.0.1-s6vr55m7kq5gylxswlaqy5usltbrslzd;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_36/rocmcc-6.2.1/raja-2024.07.0-m322mdffprdfso4vg7iox3n5goelhtyp;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_36/rocmcc-6.2.1/umpire-2024.07.0-h4aizqln4js5uo6o5dnkkjadcq5apw37;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_36/rocmcc-6.2.1/netcdf-c-4.7.4-m3moxxpcdu3id35a6ib662pmcpe6hlxt;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_36/rocmcc-6.2.1/strumpack-8.0.0-w5kwnpidy7gzm4q2kplzw7ip2hebkojg;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_36/rocmcc-6.2.1/sundials-6.7.0-ihjcn5jozjo6hqxvueljdahsvpg4fpyr;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_36/rocmcc-6.2.1/superlu-dist-8.1.2-ulxns2zfpqyuezhrgxjh7do45ezer2m5;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_36/rocmcc-6.2.1/camp-2024.07.0-223y52kngyoyilxr7gz3i4pe5btnm5uq;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_36/rocmcc-6.2.1/fmt-11.0.2-haismsyfch472jneh7tj5hoz5faug2wp;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_36/rocmcc-6.2.1/hdf5-1.8.23-zhqxrmbx5hqezzazejozk5l7wuz3ijx6;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_36/rocmcc-6.2.1/netlib-scalapack-2.2.2-zinikcerv4cwfrzt4kbhby6xm76fjmvq;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_36/rocmcc-6.2.1/hypre-2.26.0-sffmaz74r5zmnzir6l73deslw36yypsh;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_36/rocmcc-6.2.1/parmetis-4.0.3-7o2cbbpbawncpn24yftcj24gbopz4aev;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_36/rocmcc-6.2.1/metis-5.1.0-q4s4zr5oj2vsyc4eqjy22wrwndgb4fnh;/usr/tce;/usr/tce/packages/cray-mpich-tce/cray-mpich-8.1.29-rocmcc-6.2.1;/opt/rocm-6.2.1;/opt/rocm-6.2.1;/opt/rocm-6.2.1;/opt/rocm-6.2.1/llvm;/opt/rocm-6.2.1;/opt/rocm-6.2.1;/opt/rocm-6.2.1;/opt/rocm-6.2.1" CACHE STRING "") + +set(CMAKE_INSTALL_RPATH_USE_LINK_PATH "ON" CACHE STRING "") + +set(CMAKE_BUILD_RPATH "/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_36/rocmcc-6.2.1/serac-develop-og4fmn5dkubewiybduaji237suiuzssb/lib;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_36/rocmcc-6.2.1/serac-develop-og4fmn5dkubewiybduaji237suiuzssb/lib64;;/opt/rh/gcc-toolset-12/root/usr/lib/gcc/x86_64-redhat-linux/12" CACHE STRING "") + +set(CMAKE_INSTALL_RPATH "/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_36/rocmcc-6.2.1/serac-develop-og4fmn5dkubewiybduaji237suiuzssb/lib;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_36/rocmcc-6.2.1/serac-develop-og4fmn5dkubewiybduaji237suiuzssb/lib64;;/opt/rh/gcc-toolset-12/root/usr/lib/gcc/x86_64-redhat-linux/12" CACHE STRING "") + +set(CMAKE_BUILD_TYPE "Release" CACHE STRING "") + +#------------------------------------------------------------------------------ +# Compilers +#------------------------------------------------------------------------------ +# Compiler Spec: rocmcc@=6.2.1 +#------------------------------------------------------------------------------ +if(DEFINED ENV{SPACK_CC}) + + set(CMAKE_C_COMPILER "/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_36/spack/lib/spack/env/rocmcc/amdclang" CACHE PATH "") + + set(CMAKE_CXX_COMPILER "/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_36/spack/lib/spack/env/rocmcc/amdclang++" CACHE PATH "") + + set(CMAKE_Fortran_COMPILER "/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_36/spack/lib/spack/env/rocmcc/amdflang" CACHE PATH "") + +else() + + set(CMAKE_C_COMPILER "/opt/rocm-6.2.1/llvm/bin/amdclang" CACHE PATH "") + + set(CMAKE_CXX_COMPILER "/opt/rocm-6.2.1/llvm/bin/amdclang++" CACHE PATH "") + + set(CMAKE_Fortran_COMPILER "/opt/rocm-6.2.1/llvm/bin/amdflang" CACHE PATH "") + +endif() + +#------------------------------------------------------------------------------ +# MPI +#------------------------------------------------------------------------------ + +set(MPI_C_COMPILER "/usr/tce/packages/cray-mpich-tce/cray-mpich-8.1.29-rocmcc-6.2.1/bin/mpicc" CACHE PATH "") + +set(MPI_CXX_COMPILER "/usr/tce/packages/cray-mpich-tce/cray-mpich-8.1.29-rocmcc-6.2.1/bin/mpicxx" CACHE PATH "") + +set(MPI_Fortran_COMPILER "/usr/tce/packages/cray-mpich-tce/cray-mpich-8.1.29-rocmcc-6.2.1/bin/mpif90" CACHE PATH "") + +set(MPIEXEC_NUMPROC_FLAG "-n" CACHE STRING "") + +set(ENABLE_MPI ON CACHE BOOL "") + +set(MPIEXEC_EXECUTABLE "/usr/global/tools/flux_wrappers/bin/srun" CACHE PATH "") + +#------------------------------------------------------------------------------ +# Hardware +#------------------------------------------------------------------------------ + +#------------------------------------------------ +# ROCm +#------------------------------------------------ + +set(HIP_ROOT_DIR "/opt/rocm-6.2.1" CACHE PATH "") + +set(CMAKE_HIP_COMPILER "/opt/rocm-6.2.1/llvm/bin/clang++" CACHE FILEPATH "") + +set(CMAKE_HIP_ARCHITECTURES "gfx90a" CACHE STRING "") + +set(AMDGPU_TARGETS "gfx90a" CACHE STRING "") + +set(GPU_TARGETS "gfx90a" CACHE STRING "") + +set(ENABLE_OPENMP OFF CACHE BOOL "") + +set(ENABLE_HIP ON CACHE BOOL "") + +set(ROCM_PATH "/opt/rocm-6.2.1" CACHE PATH "") + +set(CMAKE_EXE_LINKER_FLAGS "-L/opt/rocm-6.2.1/lib -Wl,-rpath,/opt/rocm-6.2.1/lib -L/opt/rocm-6.2.1/llvm/lib -Wl,-rpath,/opt/rocm-6.2.1/llvm/lib -lxpmem -L/opt/cray/pe/mpich/8.1.29/gtl/lib -Wl,-rpath,/opt/cray/pe/mpich/8.1.29/gtl/lib -lmpi_gtl_hsa -Wl,--disable-new-dtags -lflang -lflangrti -lamdhip64 -lhsakmt -lhsa-runtime64 -lamd_comgr -lpgmath -lhipblas" CACHE STRING "") + +#------------------------------------------------------------------------------ +# TPLs +#------------------------------------------------------------------------------ + +set(TPL_ROOT "/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib_cray/2025_04_16_13_46_36/rocmcc-6.2.1" CACHE PATH "") + +set(AXOM_DIR "${TPL_ROOT}/axom-0.10.1.1-4udyvzqnaoic74uv3zcumyguv5mbqqcg" CACHE PATH "") + +set(CAMP_DIR "${TPL_ROOT}/camp-2024.07.0-223y52kngyoyilxr7gz3i4pe5btnm5uq" CACHE PATH "") + +set(CONDUIT_DIR "${TPL_ROOT}/conduit-0.9.2-wo23sn3yekayvxjphe2eu4vvmvkqrovj" CACHE PATH "") + +set(LUA_DIR "/usr" CACHE PATH "") + +set(MFEM_DIR "${TPL_ROOT}/mfem-4.8.0.1-s6vr55m7kq5gylxswlaqy5usltbrslzd" CACHE PATH "") + +set(HDF5_DIR "${TPL_ROOT}/hdf5-1.8.23-zhqxrmbx5hqezzazejozk5l7wuz3ijx6" CACHE PATH "") + +set(HYPRE_DIR "${TPL_ROOT}/hypre-2.26.0-sffmaz74r5zmnzir6l73deslw36yypsh" CACHE PATH "") + +set(METIS_DIR "${TPL_ROOT}/metis-5.1.0-q4s4zr5oj2vsyc4eqjy22wrwndgb4fnh" CACHE PATH "") + +set(PARMETIS_DIR "${TPL_ROOT}/parmetis-4.0.3-7o2cbbpbawncpn24yftcj24gbopz4aev" CACHE PATH "") + +set(NETCDF_DIR "${TPL_ROOT}/netcdf-c-4.7.4-m3moxxpcdu3id35a6ib662pmcpe6hlxt" CACHE PATH "") + +set(SUPERLUDIST_DIR "${TPL_ROOT}/superlu-dist-8.1.2-ulxns2zfpqyuezhrgxjh7do45ezer2m5" CACHE PATH "") + +set(ADIAK_DIR "${TPL_ROOT}/adiak-0.4.1-djvwhzv7uucbw342y2lxn6mjv4onnzeh" CACHE PATH "") + +# AMGX not built + +set(CALIPER_DIR "${TPL_ROOT}/caliper-2.12.1-h2elmogt2mtorvunf67jrcxqdcrzjjem" CACHE PATH "") + +# PETSC not built + +set(RAJA_DIR "${TPL_ROOT}/raja-2024.07.0-m322mdffprdfso4vg7iox3n5goelhtyp" CACHE PATH "") + +# SLEPC not built + +set(STRUMPACK_DIR "${TPL_ROOT}/strumpack-8.0.0-w5kwnpidy7gzm4q2kplzw7ip2hebkojg" CACHE PATH "") + +set(SUNDIALS_DIR "${TPL_ROOT}/sundials-6.7.0-ihjcn5jozjo6hqxvueljdahsvpg4fpyr" CACHE PATH "") + +set(UMPIRE_DIR "${TPL_ROOT}/umpire-2024.07.0-h4aizqln4js5uo6o5dnkkjadcq5apw37" CACHE PATH "") + +set(TRIBOL_DIR "${TPL_ROOT}/tribol-0.1.0.18-neize4eyd3wxmeuui44az2x5ecqfe2pd" CACHE PATH "") + +set(ENZYME_DIR "/usr/WS2/smithdev/toss_4_x86_64_ib_cray/enzyme/rocm-6.2.1/0.0.180" CACHE PATH "") + +#------------------------------------------------------------------------------ +# Devtools +#------------------------------------------------------------------------------ + +# Code checks disabled due to disabled devtools + +set(SERAC_ENABLE_CODE_CHECKS OFF CACHE BOOL "") + +set(ENABLE_CLANGFORMAT OFF CACHE BOOL "") + +set(ENABLE_CLANGTIDY OFF CACHE BOOL "") + +set(ENABLE_DOCS OFF CACHE BOOL "") + +set(SERAC_ENABLE_CODEVELOP ON CACHE BOOL "") + +set(SERAC_DISABLE_TRIBOL ON CACHE BOOL "") + From 94ea8368150c396cc4688a9814822f110d28bfa0 Mon Sep 17 00:00:00 2001 From: "Eric B. Chin" Date: Tue, 3 Jun 2025 22:09:04 -0700 Subject: [PATCH 03/71] ruby dfem host config --- .../ruby-toss_4_x86_64_ib-clang@19.1.7.cmake | 138 ++++++++++++++++++ 1 file changed, 138 insertions(+) create mode 100644 host-configs/dfem/ruby-toss_4_x86_64_ib-clang@19.1.7.cmake diff --git a/host-configs/dfem/ruby-toss_4_x86_64_ib-clang@19.1.7.cmake b/host-configs/dfem/ruby-toss_4_x86_64_ib-clang@19.1.7.cmake new file mode 100644 index 000000000..fa0cd4152 --- /dev/null +++ b/host-configs/dfem/ruby-toss_4_x86_64_ib-clang@19.1.7.cmake @@ -0,0 +1,138 @@ +#------------------------------------------------------------------------------ +# !!!! This is a generated file, edit at own risk !!!! +#------------------------------------------------------------------------------ +# CMake executable path: /usr/tce/bin/cmake +#------------------------------------------------------------------------------ + +set(CMAKE_PREFIX_PATH "/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_06_03_16_52_06/clang-19.1.7/caliper-2.12.1-5o3zfbeosigjwrp5ffizjs6sf7pyoqjv;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_06_03_16_52_06/clang-19.1.7/tribol-0.1.0.18-my4bcnwnewxcnewg3xgbzvym64dfizuj;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_06_03_16_52_06/clang-19.1.7/adiak-0.4.1-5wjrv4xwbm3y7yh7745ijc45sbcwqwqi;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_06_03_16_52_06/clang-19.1.7/libunwind-1.8.1-tgmphz65mn62bautt2ksl45ckecw5e2c;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_06_03_16_52_06/clang-19.1.7/axom-0.10.1.1-tkmjd26d72eqk4tts7hwlwognvwuyj5s;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_06_03_16_52_06/clang-19.1.7/conduit-0.9.2-t2zobt76hkwap7d2ppddp3traorxpffl;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_06_03_16_52_06/clang-19.1.7/mfem-4.8.0.1-gbv5tbeosslcqh2zhkcrgkwtqxob257j;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_06_03_16_52_06/clang-19.1.7/raja-2024.07.0-t5ay2skmkw3ybywrq33tcyl5xbn5g3gv;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_06_03_16_52_06/clang-19.1.7/umpire-2024.07.0-47yl3s6ob3ru4g6qnma4tiv6ob5fr6p6;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_06_03_16_52_06/clang-19.1.7/netcdf-c-4.7.4-62ri63jaeyajwrrvwk75dbwvrnj6fony;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_06_03_16_52_06/clang-19.1.7/slepc-3.21.2-rji2oyfq6x4lyovnrkfqnleoc3yrqrm6;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_06_03_16_52_06/clang-19.1.7/sundials-6.7.0-wwas67vc5xfgayl2uimersm3xszlqv62;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_06_03_16_52_06/clang-19.1.7/camp-2024.07.0-afpyaxeo4ulr2lkljpy2zr7rjeq64cox;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_06_03_16_52_06/clang-19.1.7/fmt-11.0.2-xpfumigjok66664dhruf6tfaa7gwcn3n;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_06_03_16_52_06/clang-19.1.7/hdf5-1.8.23-pu3kif4ltqdw2kcwvlswfze6pmmsjs6q;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_06_03_16_52_06/clang-19.1.7/arpack-ng-3.9.0-ahgmmge3jagsx3n72wbxtcznbubcclca;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_06_03_16_52_06/clang-19.1.7/petsc-3.21.6-elp7iw6xe64daspgy5ctubfo3jj4nd43;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_06_03_16_52_06/clang-19.1.7/hypre-2.26.0-ta6n6j6ziyuw2tfejt27nkdulvjq27hu;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_06_03_16_52_06/clang-19.1.7/strumpack-8.0.0-zmfwdjw54im5o6m3gqssnkvuagdzljwe;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_06_03_16_52_06/clang-19.1.7/superlu-dist-8.1.2-k53l7hvjwbqnepeb7etlujfgzrzfliua;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_06_03_16_52_06/clang-19.1.7/netlib-scalapack-2.2.2-fuk7tt7il4e6tekxxctc625g4xk4hqeu;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_06_03_16_52_06/clang-19.1.7/parmetis-4.0.3-morhecieft72ri6gi2r6hsatzrgdhe5n;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_06_03_16_52_06/clang-19.1.7/metis-5.1.0-ebbueyypx5bn5mekybxavligkykhm5fy;/usr/workspace/smithdev/devtools/toss_4_x86_64_ib/latest/python-3.11.7;/usr/workspace/smithdev/devtools/toss_4_x86_64_ib/latest/python-3.11.7;/usr/tce;/usr/workspace/smithdev/devtools/toss_4_x86_64_ib/latest/cppcheck-2.9;/usr/workspace/smithdev/devtools/toss_4_x86_64_ib/latest/doxygen-1.9.8;/usr/tce/packages/clang/clang-14.0.6;/usr/tce/packages/mvapich2/mvapich2-2.3.7-gcc-10.3.1;/usr/tce/packages/python/python-3.9.12" CACHE STRING "") + +set(CMAKE_INSTALL_RPATH_USE_LINK_PATH "ON" CACHE STRING "") + +set(CMAKE_BUILD_RPATH "/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_06_03_16_52_06/clang-19.1.7/serac-develop-rxm4riztuqljjlnrhyomn4cqusqy5g6z/lib;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_06_03_16_52_06/clang-19.1.7/serac-develop-rxm4riztuqljjlnrhyomn4cqusqy5g6z/lib64;;/usr/tce/packages/gcc/gcc-10.3.1/lib/gcc/x86_64-redhat-linux/10" CACHE STRING "") + +set(CMAKE_INSTALL_RPATH "/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_06_03_16_52_06/clang-19.1.7/serac-develop-rxm4riztuqljjlnrhyomn4cqusqy5g6z/lib;/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_06_03_16_52_06/clang-19.1.7/serac-develop-rxm4riztuqljjlnrhyomn4cqusqy5g6z/lib64;;/usr/tce/packages/gcc/gcc-10.3.1/lib/gcc/x86_64-redhat-linux/10" CACHE STRING "") + +set(CMAKE_BUILD_TYPE "Release" CACHE STRING "") + +#------------------------------------------------------------------------------ +# Compilers +#------------------------------------------------------------------------------ +# Compiler Spec: clang@=19.1.7 +#------------------------------------------------------------------------------ +if(DEFINED ENV{SPACK_CC}) + + set(CMAKE_C_COMPILER "/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_06_03_16_52_06/spack/lib/spack/env/clang/clang" CACHE PATH "") + + set(CMAKE_CXX_COMPILER "/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_06_03_16_52_06/spack/lib/spack/env/clang/clang++" CACHE PATH "") + + set(CMAKE_Fortran_COMPILER "/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_06_03_16_52_06/spack/lib/spack/env/clang/gfortran" CACHE PATH "") + +else() + + set(CMAKE_C_COMPILER "/usr/WS2/smithdev/toss_4_x86_64_ib/llvm/19.1.7/bin/clang" CACHE PATH "") + + set(CMAKE_CXX_COMPILER "/usr/WS2/smithdev/toss_4_x86_64_ib/llvm/19.1.7/bin/clang++" CACHE PATH "") + + set(CMAKE_Fortran_COMPILER "/usr/tce/packages/gcc/gcc-10.3.1/bin/gfortran" CACHE PATH "") + +endif() + +set(CMAKE_C_FLAGS "--gcc-toolchain=/usr/tce/packages/gcc/gcc-10.3.1" CACHE STRING "") + +set(CMAKE_CXX_FLAGS "--gcc-toolchain=/usr/tce/packages/gcc/gcc-10.3.1" CACHE STRING "") + +set(CMAKE_Fortran_FLAGS "-fPIE -fPIC" CACHE STRING "") + +#------------------------------------------------------------------------------ +# MPI +#------------------------------------------------------------------------------ + +set(MPI_C_COMPILER "/usr/tce/packages/mvapich2/mvapich2-2.3.7-gcc-10.3.1/bin/mpicc" CACHE PATH "") + +set(MPI_CXX_COMPILER "/usr/tce/packages/mvapich2/mvapich2-2.3.7-gcc-10.3.1/bin/mpicxx" CACHE PATH "") + +set(MPI_Fortran_COMPILER "/usr/tce/packages/mvapich2/mvapich2-2.3.7-gcc-10.3.1/bin/mpif90" CACHE PATH "") + +set(MPIEXEC_NUMPROC_FLAG "-n" CACHE STRING "") + +set(ENABLE_MPI ON CACHE BOOL "") + +set(MPIEXEC_EXECUTABLE "/usr/bin/srun" CACHE PATH "") + +#------------------------------------------------------------------------------ +# Hardware +#------------------------------------------------------------------------------ + +set(ENABLE_OPENMP OFF CACHE BOOL "") + +#------------------------------------------------------------------------------ +# TPLs +#------------------------------------------------------------------------------ + +set(TPL_ROOT "/usr/WS2/smithdev/libs/serac/toss_4_x86_64_ib/2025_06_03_16_52_06/clang-19.1.7" CACHE PATH "") + +set(AXOM_DIR "${TPL_ROOT}/axom-0.10.1.1-tkmjd26d72eqk4tts7hwlwognvwuyj5s" CACHE PATH "") + +set(CAMP_DIR "${TPL_ROOT}/camp-2024.07.0-afpyaxeo4ulr2lkljpy2zr7rjeq64cox" CACHE PATH "") + +set(CONDUIT_DIR "${TPL_ROOT}/conduit-0.9.2-t2zobt76hkwap7d2ppddp3traorxpffl" CACHE PATH "") + +set(LUA_DIR "/usr" CACHE PATH "") + +set(MFEM_DIR "${TPL_ROOT}/mfem-4.8.0.1-gbv5tbeosslcqh2zhkcrgkwtqxob257j" CACHE PATH "") + +set(HDF5_DIR "${TPL_ROOT}/hdf5-1.8.23-pu3kif4ltqdw2kcwvlswfze6pmmsjs6q" CACHE PATH "") + +set(HYPRE_DIR "${TPL_ROOT}/hypre-2.26.0-ta6n6j6ziyuw2tfejt27nkdulvjq27hu" CACHE PATH "") + +set(METIS_DIR "${TPL_ROOT}/metis-5.1.0-ebbueyypx5bn5mekybxavligkykhm5fy" CACHE PATH "") + +set(PARMETIS_DIR "${TPL_ROOT}/parmetis-4.0.3-morhecieft72ri6gi2r6hsatzrgdhe5n" CACHE PATH "") + +set(NETCDF_DIR "${TPL_ROOT}/netcdf-c-4.7.4-62ri63jaeyajwrrvwk75dbwvrnj6fony" CACHE PATH "") + +set(SUPERLUDIST_DIR "${TPL_ROOT}/superlu-dist-8.1.2-k53l7hvjwbqnepeb7etlujfgzrzfliua" CACHE PATH "") + +set(ARPACK_DIR "${TPL_ROOT}/arpack-ng-3.9.0-ahgmmge3jagsx3n72wbxtcznbubcclca" CACHE PATH "") + +set(ADIAK_DIR "${TPL_ROOT}/adiak-0.4.1-5wjrv4xwbm3y7yh7745ijc45sbcwqwqi" CACHE PATH "") + +# AMGX not built + +set(CALIPER_DIR "${TPL_ROOT}/caliper-2.12.1-5o3zfbeosigjwrp5ffizjs6sf7pyoqjv" CACHE PATH "") + +set(PETSC_DIR "${TPL_ROOT}/petsc-3.21.6-elp7iw6xe64daspgy5ctubfo3jj4nd43" CACHE PATH "") + +set(RAJA_DIR "${TPL_ROOT}/raja-2024.07.0-t5ay2skmkw3ybywrq33tcyl5xbn5g3gv" CACHE PATH "") + +set(SLEPC_DIR "${TPL_ROOT}/slepc-3.21.2-rji2oyfq6x4lyovnrkfqnleoc3yrqrm6" CACHE PATH "") + +set(STRUMPACK_DIR "${TPL_ROOT}/strumpack-8.0.0-zmfwdjw54im5o6m3gqssnkvuagdzljwe" CACHE PATH "") + +set(SUNDIALS_DIR "${TPL_ROOT}/sundials-6.7.0-wwas67vc5xfgayl2uimersm3xszlqv62" CACHE PATH "") + +set(UMPIRE_DIR "${TPL_ROOT}/umpire-2024.07.0-47yl3s6ob3ru4g6qnma4tiv6ob5fr6p6" CACHE PATH "") + +set(TRIBOL_DIR "${TPL_ROOT}/tribol-0.1.0.18-my4bcnwnewxcnewg3xgbzvym64dfizuj" CACHE PATH "") + +#------------------------------------------------------------------------------ +# Devtools +#------------------------------------------------------------------------------ + +set(DEVTOOLS_ROOT "/usr/WS2/smithdev/devtools/toss_4_x86_64_ib/2024_05_30_15_02_20/._view/psk2dcrijss6s4i6qmxplzthrzm3y7nh" CACHE PATH "") + +set(ATS_EXECUTABLE "${DEVTOOLS_ROOT}/python-3.11.7/bin/ats" CACHE PATH "") + +set(CLANGFORMAT_EXECUTABLE "/usr/tce/packages/clang/clang-14.0.6/bin/clang-format" CACHE PATH "") + +set(CLANGTIDY_EXECUTABLE "/usr/tce/packages/clang/clang-14.0.6/bin/clang-tidy" CACHE PATH "") + +set(ENABLE_DOCS ON CACHE BOOL "") + +set(SPHINX_EXECUTABLE "${DEVTOOLS_ROOT}/python-3.11.7/bin/sphinx-build" CACHE PATH "") + +set(CPPCHECK_EXECUTABLE "${DEVTOOLS_ROOT}/cppcheck-2.9/bin/cppcheck" CACHE PATH "") + +set(DOXYGEN_EXECUTABLE "${DEVTOOLS_ROOT}/doxygen-1.9.8/bin/doxygen" CACHE PATH "") + + From bffad920fce2e134113b4efdca893003acddda52 Mon Sep 17 00:00:00 2001 From: "E. B. Chin" Date: Wed, 18 Jun 2025 21:55:18 -0700 Subject: [PATCH 04/71] more general dfem qfunction material --- mfem | 2 +- src/serac/numerics/functional/tensor.hpp | 8 +- src/serac/physics/dfem_residual.hpp | 123 +++++++ src/serac/physics/solid_dfem_residual.hpp | 103 ++++++ .../physics/tests/test_dfem_vs_functional.cpp | 340 ++++++++++++++++++ 5 files changed, 571 insertions(+), 5 deletions(-) create mode 100644 src/serac/physics/dfem_residual.hpp create mode 100644 src/serac/physics/solid_dfem_residual.hpp create mode 100644 src/serac/physics/tests/test_dfem_vs_functional.cpp diff --git a/mfem b/mfem index 3e3a5e5d1..e72ed572c 160000 --- a/mfem +++ b/mfem @@ -1 +1 @@ -Subproject commit 3e3a5e5d14003d7856f779292fe032067c10fc1b +Subproject commit e72ed572c4023ba6a3f0b3a28452890894d0580a diff --git a/src/serac/numerics/functional/tensor.hpp b/src/serac/numerics/functional/tensor.hpp index 83e5a6b41..7ce1082c9 100644 --- a/src/serac/numerics/functional/tensor.hpp +++ b/src/serac/numerics/functional/tensor.hpp @@ -1318,8 +1318,8 @@ SERAC_HOST_DEVICE constexpr auto det(const tensor& A) * @param A Input matrix * @return det(A + I) - 1, where I is the identity matrix */ -template -SERAC_HOST_DEVICE constexpr auto detApIm1(const tensor& A) +template class TensorType> +SERAC_HOST_DEVICE constexpr auto detApIm1(const TensorType& A) { // From the Cayley-Hamilton theorem, we get that for any N by N matrix A, // det(A - I) - 1 = I1(A) + I2(A) + ... + IN(A), @@ -1331,8 +1331,8 @@ SERAC_HOST_DEVICE constexpr auto detApIm1(const tensor& A) } /// @overload -template -SERAC_HOST_DEVICE constexpr auto detApIm1(const tensor& A) +template class TensorType> +SERAC_HOST_DEVICE constexpr auto detApIm1(const TensorType& A) { // For notes on the implementation, see the 2x2 version. diff --git a/src/serac/physics/dfem_residual.hpp b/src/serac/physics/dfem_residual.hpp new file mode 100644 index 000000000..302323e2e --- /dev/null +++ b/src/serac/physics/dfem_residual.hpp @@ -0,0 +1,123 @@ +// 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) + +/** + * @file dfem_residual.hpp + */ + +#pragma once + +#include "mfem.hpp" + +#include "serac/physics/residual.hpp" +#include "serac/physics/mesh.hpp" +#include "serac/physics/state/finite_element_state.hpp" +#include "serac/physics/state/finite_element_dual.hpp" + +namespace serac { + +/** + * @brief The nonlinear residual class + * + * This uses dFEM to compute fairly arbitrary residuals and tangent + * stiffness matrices based on body and boundary integrals. + * + */ +class DfemResidual : public Residual { + public: + using SpacesT = std::vector; ///< typedef + + /** + * @brief Construct a new DfemResidual object + * + * @param physics_name A name for the physics module instance + * @param mesh The serac mesh + * @param diff_op A differentiable operator that computes the residual and jacobian + */ + DfemResidual(std::string physics_name, std::shared_ptr mesh, mfem::future::DifferentiableOperator&& diff_op) + : Residual(physics_name), mesh_(mesh), diff_op_(std::move(diff_op)) + { + } + + /** + * @brief Add a body integral contribution to the residual + * + * @tparam BodyIntegralType The type of the body integral + * // DependsOn can be indices into fields which the body integral may depend on + * @param body_name The name of the registered domain over which the body force is applied. If nothing is supplied + * the entire domain is + * @param body_integral A function describing the body force applied + * used. + * @pre body_integral must be a object that can be called with the following arguments: + * 1. `double t` the time + * 2. `tensor x` the spatial coordinates for the quadrature point + * 3. `tuple{value, derivative}`, a variadic list of tuples (each with a values and derivative), + * one tuple for each of the trial spaces specified in the `DependsOn<...>` argument. + * @note The actual types of these arguments passed will be `double`, `tensor` or tuples thereof + * when doing direct evaluation. When differentiating with respect to one of the inputs, its stored + * values will change to `dual` numbers rather than `double`. (e.g. `tensor` becomes `tensor, + * 3>`) + * + */ + template + void addBodyIntegral(QFunctionType qfunction, InputType inputs, OutputType outputs, + const mfem::IntegrationRule& integration_rule, mfem::Array domain_attributes, + std::integer_sequence derivative_ids) + { + diff_op_.AddDomainIntegrator(qfunction, inputs, outputs, integration_rule, domain_attributes, derivative_ids); + } + + /// @overload + mfem::Vector residual(double, double, const std::vector& fields, int block_row = 0) const override + { + SLIC_ERROR_IF(block_row != 0, "Invalid block row and column requested in fieldJacobian for FunctionalResidual"); + diff_op_.SetParameters({mesh_nodes_}); + mfem::Vector ret(fields[0]->space().GetVSize()); + diff_op_.Mult(fields[0]->gridFunction(), ret); + return ret; + } + + /// @overload + std::unique_ptr jacobian(double, double, const std::vector&, + const std::vector&, int) const override + { + return std::make_unique(); + } + + /// @overload + void jvp(double, double, const std::vector& fields, const std::vector& v_fields, + const std::vector& jvp_reactions) const override + { + SLIC_ERROR_IF(v_fields.size() != fields.size(), + "Invalid number of field sensitivities relative to the number of fields"); + SLIC_ERROR_IF(jvp_reactions.size() != 1, "DfemResidual nonlinear systems only supports 1 output residual"); + + auto grad_op = diff_op_.GetDerivative(0, {fields[0]}, {mesh_nodes_}); + grad_op->AddMult(*v_fields[0], *jvp_reactions[0]); + } + + /// @overload + void vjp(double, double, const std::vector& fields, const std::vector& v_fields, + const std::vector& vjp_sensitivities) const override + { + SLIC_ERROR_IF(vjp_sensitivities.size() != fields.size(), + "Invalid number of field sensitivities relative to the number of fields"); + SLIC_ERROR_IF(v_fields.size() != 1, "FunctionalResidual nonlinear systems only supports 1 output residual"); + + auto grad_op = diff_op_.GetDerivative(0, {fields[0]}, {mesh_nodes_}); + grad_op->AddMultTranspose(*v_fields[0], *vjp_sensitivities[0]); + } + + protected: + /// @brief primary mesh + std::shared_ptr mesh_; + + mfem::ParGridFunction* mesh_nodes_; ///< grid function for the mesh nodes + + mutable mfem::future::DifferentiableOperator diff_op_; ///< differentiable operator for the residual +}; + +} // namespace serac diff --git a/src/serac/physics/solid_dfem_residual.hpp b/src/serac/physics/solid_dfem_residual.hpp new file mode 100644 index 000000000..d707ff52b --- /dev/null +++ b/src/serac/physics/solid_dfem_residual.hpp @@ -0,0 +1,103 @@ +// 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) + +/** + * @file solid_dfem_residual.hpp + * + * @brief Implements the residual interface for solid mechanics physics. + * Derives from dfem_residual. + */ + +#pragma once + +#include "serac/physics/dfem_residual.hpp" + +namespace serac { + +template +struct MomentumQFunction { + SERAC_HOST_DEVICE inline auto operator()( + mfem::real_t dt, const mfem::future::tensor& internal_state, + const mfem::future::tensor& du_dxi, + const mfem::future::tensor& dv_dxi, const mfem::future::tensor&, + const mfem::future::tensor& dX_dxi, mfem::real_t weight) const + { + auto dxi_dX = mfem::future::inv(dX_dxi); + auto du_dX = mfem::future::dot(du_dxi, dxi_dX); + auto dv_dX = mfem::future::dot(dv_dxi, dxi_dX); + auto P = mfem::future::get<0>(material.pkStress(dt, internal_state, du_dX, dv_dxi)); + auto JxW = mfem::future::det(dX_dxi) * weight * mfem::future::transpose(dxi_dX); + return mfem::future::tuple{P * JxW}; + } + + Material material; ///< the material model to use for computing the stress +}; + +/** + * @brief The nonlinear residual class + * + * This uses Functional to compute the solid mechanics residuals and tangent + * stiffness matrices. + * + * @tparam order The order of the discretization of the displacement and velocity fields + * @tparam dim The spatial dimension of the mesh + */ +template +class SolidDfemResidual : public DfemResidual { + public: + /// @brief disp, velo, accel + static constexpr int NUM_STATE_VARS = 3; + + /** + * @brief Construct a new SolidResidual object + * + * @param physics_name A name for the physics module instance + * @param mesh The serac Mesh + * @param diff_op A differentiable operator that computes the residual and jacobian + */ + SolidDfemResidual(std::string physics_name, std::shared_ptr mesh, + mfem::future::DifferentiableOperator&& diff_op) + : DfemResidual(physics_name, mesh, std::move(diff_op)) + { + } + + /** + * @brief Set the material stress response and mass properties for the physics module + * + * @tparam MaterialType The solid material type + * @tparam StateType the type that contains the internal variables for MaterialType + * @param body_name string name for a registered body Domain on the mesh + * @param material A material that provides a function to evaluate stress + * @pre material must be a object that can be called with the following arguments: + * 1. `MaterialType::State & state` an mutable reference to the internal variables for this quadrature point + * 2. `tensor du_dx` the displacement gradient at this quadrature point + * 3. `tuple{value, derivative}`, a tuple of values and derivatives for each parameter field + * specified in the `DependsOn<...>` argument. + * + * @note The actual types of these arguments passed will be `double`, `tensor` or tuples thereof + * when doing direct evaluation. When differentiating with respect to one of the inputs, its stored + * values will change to `dual` numbers rather than `double`. (e.g. `tensor` becomes `tensor, + * 3>`) + * + * @param qdata the buffer of material internal variables at each quadrature point + * + * @pre MaterialType must have a public method `density` which can take FiniteElementState parameter inputs + * @pre MaterialType must have a public method 'pkStress' which returns the first Piola-Kirchhoff stress + * + */ + template + void setMaterial(std::string, const MaterialType& material, const mfem::IntegrationRule& displacement_ir) + { + mfem::future::tuple inputs{mfem::future::Gradient<0>{}, mfem::future::Gradient<1>{}, mfem::future::Gradient<2>{}, + mfem::future::Gradient<3>{}, mfem::future::Weight{}}; + mfem::future::tuple outputs{mfem::future::Gradient<4>{}}; + mfem::Array solid_domain_attributes(DfemResidual::mesh_->mfemParMesh().attributes.Max()); + auto momentum_qf = MomentumQFunction{.material = material}; + DfemResidual::diff_op_.AddDomainIntegrator(momentum_qf, inputs, outputs, displacement_ir, solid_domain_attributes); + } +}; + +} // namespace serac diff --git a/src/serac/physics/tests/test_dfem_vs_functional.cpp b/src/serac/physics/tests/test_dfem_vs_functional.cpp new file mode 100644 index 000000000..4accf3544 --- /dev/null +++ b/src/serac/physics/tests/test_dfem_vs_functional.cpp @@ -0,0 +1,340 @@ +// 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 "mfem.hpp" +#include "serac/infrastructure/application_manager.hpp" +#include "serac/mesh_utils/mesh_utils.hpp" +#include "serac/physics/materials/solid_material.hpp" +#include "serac/physics/mesh.hpp" +#include "serac/physics/state/state_manager.hpp" +#include "serac/physics/solid_residual.hpp" + +auto element_shape = mfem::Element::QUADRILATERAL; + +template +auto getPointers(std::vector& values) +{ + std::vector pointers; + for (auto& t : values) { + pointers.push_back(&t); + } + return pointers; +} + +template +auto getPointers(std::vector& states, std::vector& params) +{ + assert(params.size() > 0); + std::vector pointers; + for (auto& t : states) { + pointers.push_back(&t); + } + for (auto& t : params) { + pointers.push_back(&t); + } + return pointers; +} + +template +auto getPointers(T& v) +{ + return std::vector{&v}; +} + +void pseudoRand(serac::FiniteElementState& dual) +{ + int sz = dual.Size(); + for (int i = 0; i < sz; ++i) { + dual(i) = -1.2 + 2.02 * (double(i) / sz); + } +} + +void pseudoRand(serac::FiniteElementDual& dual) +{ + int sz = dual.Size(); + for (int i = 0; i < sz; ++i) { + dual(i) = -1.2 + 2.02 * (double(i) / sz); + } +} + +namespace serac { + +struct NeoHookeanWithFieldWithRateFunctional { + using State = Empty; ///< this material has no internal variables + + template + SERAC_HOST_DEVICE auto pkStress(double /*dt*/, State& /* state */, const tensor& du_dX, + const tensor& /*dv_dX*/) const + { + using std::log1p; + constexpr auto I = Identity(); + auto lambda = K - (2.0 / 3.0) * G; + auto B_minus_I = dot(du_dX, transpose(du_dX)) + transpose(du_dX) + du_dX; + + auto logJ = log1p(detApIm1(du_dX)); + // Kirchoff stress, in form that avoids cancellation error when F is near I + auto TK = lambda * logJ * I + G * B_minus_I; + + // Pull back to Piola + auto F = du_dX + I; + return dot(TK, inv(transpose(F))); + } + + SERAC_HOST_DEVICE auto density() const { return Rho; } + + double K; ///< bulk modulus + double G; ///< shear modulus + double Rho; ///< density +}; + +struct NeoHookeanWithFieldWithRateDfem { + static constexpr int state_size = 0; + + template + SERAC_HOST_DEVICE auto pkStress(double /*dt*/, const mfem::future::tensor& /* state */, + const mfem::future::tensor& du_dX, + const mfem::future::tensor& /*dv_dX*/) const + { + using std::log1p; + constexpr auto I = mfem::future::IdentityMatrix(); + auto lambda = K - (2.0 / 3.0) * G; + auto B_minus_I = mfem::future::dot(du_dX, mfem::future::transpose(du_dX)) + mfem::future::transpose(du_dX) + du_dX; + + auto logJ = log1p(detApIm1(du_dX)); + // Kirchoff stress, in form that avoids cancellation error when F is near I + auto TK = lambda * logJ * I + G * B_minus_I; + + // Pull back to Piola + auto F = du_dX + I; + return mfem::future::dot(TK, mfem::future::inv(mfem::future::transpose(F))); + } + + SERAC_HOST_DEVICE auto density() const { return Rho; } + + double K; ///< bulk modulus + double G; ///< shear modulus + double Rho; ///< density +}; + +} // namespace serac + +struct ResidualFixture : public testing::Test { + static constexpr int dim = 2; + static constexpr int disp_order = 1; + + using VectorSpace = serac::H1; + using DensitySpace = serac::L2; + + // using SolidMaterial = serac::solid_mechanics::NeoHookeanWithFieldDensity; + using FunctionalMaterial = serac::NeoHookeanWithFieldWithRateFunctional; + using DfemMaterial = serac::NeoHookeanWithFieldWithRateDfem; + + enum StateOrder + { + SHAPE, + DISP, + VELO, + ACCEL + }; + + void SetUp() + { + MPI_Barrier(MPI_COMM_WORLD); + serac::StateManager::initialize(datastore, "solid_dynamics"); + + double length = 0.5; + double width = 2.0; + mesh = std::make_shared(mfem::Mesh::MakeCartesian2D(6, 20, element_shape, true, length, width), + "this_mesh_name", 0, 0); + + serac::FiniteElementState disp = serac::StateManager::newState(VectorSpace{}, "displacement", mesh->tag()); + serac::FiniteElementState velo = serac::StateManager::newState(VectorSpace{}, "velocity", mesh->tag()); + serac::FiniteElementState accel = serac::StateManager::newState(VectorSpace{}, "acceleration", mesh->tag()); + serac::FiniteElementState shape_disp = + serac::StateManager::newState(VectorSpace{}, "shape_displacement", mesh->tag()); + serac::FiniteElementState density = serac::StateManager::newState(DensitySpace{}, "density", mesh->tag()); + + states = {shape_disp, disp, velo, accel}; + params = {density}; + + for (auto s : states) { + dual_states.push_back(serac::FiniteElementDual(s.space(), s.name() + "_dual")); + } + for (auto p : params) { + dual_params.push_back(serac::FiniteElementDual(p.space(), p.name() + "_dual")); + } + + v_rhs_states = states; + v_rhs_params = params; + + std::string physics_name = "solid"; + + using SolidResidualT = serac::SolidResidual>; + auto solid_mechanics_residual = std::make_shared(physics_name, mesh, states[SHAPE].space(), + states[DISP].space(), getSpaces(params)); + // SolidMaterial mat; + // mat.K = 1.0; + // mat.G = 0.5; + SolidRateMaterial rate_mat; + rate_mat.K = 1.0; + rate_mat.G = 0.5; + rate_mat.Rho = 1.5; + + // solid_mechanics_residual->setMaterial(serac::DependsOn<0>{}, mesh->entireBodyName(), mat); + solid_mechanics_residual->setRateMaterial(serac::DependsOn<>{}, mesh->entireBodyName(), rate_mat); + + // apply traction boundary conditions + std::string surface_name = "side"; + mesh->addDomainOfBoundaryElements(surface_name, serac::by_attr(1)); + solid_mechanics_residual->addBoundaryIntegral(surface_name, [](auto /*t*/, auto /*x*/, auto n) { return 1.0 * n; }); + solid_mechanics_residual->addPressure(surface_name, [](auto /*t*/, auto /*x*/) { return 0.6; }); + + // initialize fields for testing + for (auto& s : v_rhs_states) { + pseudoRand(s); + } + for (auto& p : v_rhs_params) { + pseudoRand(p); + } + + dual_states[0] = 1.0; // used to test that vjp acts via +=, add initial value to shape displacement dual + + states[DISP].setFromFieldFunction([](serac::tensor x) { + auto u = 0.1 * x; + return u; + }); + states[VELO].setFromFieldFunction([](serac::tensor x) { + auto u = -0.1 * x; + return u; + }); + states[ACCEL].setFromFieldFunction([](serac::tensor x) { + auto u = -0.01 * x; + return u; + }); + states[SHAPE] = 0.0; + params[0] = 1.2; + + // residual is abstract Residual class to ensure usage only through BasePhysics interface + residual = solid_mechanics_residual; + } + + const double time = 0.0; + const double dt = 1.0; + + std::string velo_name = "solid_velocity"; + + axom::sidre::DataStore datastore; + std::shared_ptr mesh; + std::shared_ptr residual; + + std::vector states; + std::vector params; + + std::vector dual_states; + std::vector dual_params; + + std::vector v_rhs_states; + std::vector v_rhs_params; +}; + +TEST_F(ResidualFixture, VjpConsistency) +{ + // initialize the displacement and acceleration to a non-trivial field + auto all_states = getPointers(states, params); + + serac::FiniteElementDual res_vector(states[DISP].space(), "residual"); + res_vector = residual->residual(time, dt, all_states); + ASSERT_NE(0.0, res_vector.Norml2()); + + auto jacobian_weights = [&](size_t i) { + std::vector tangents(all_states.size()); + tangents[i] = 1.0; + return tangents; + }; + + // test vjp + serac::FiniteElementState v(res_vector.space(), "v"); + pseudoRand(v); + auto all_vjps = getPointers(dual_states, dual_params); + + residual->vjp(time, dt, all_states, getPointers(v), all_vjps); + + for (size_t i = 0; i < all_states.size(); ++i) { + serac::FiniteElementState vjp = *all_states[i]; + vjp = 0.0; + auto J = residual->jacobian(time, dt, all_states, jacobian_weights(i)); + J->MultTranspose(v, vjp); + if (i == 0) vjp += 1.0; // make sure vjp uses += + EXPECT_NEAR(vjp.Norml2(), all_vjps[i]->Norml2(), 1e-12); + } +} + +TEST_F(ResidualFixture, JvpConsistency) +{ + // initialize the displacement and acceleration to a non-trivial field + auto all_states = getPointers(states, params); + + serac::FiniteElementDual res_vector(states[DISP].space(), "residual"); + res_vector = residual->residual(time, dt, all_states); + ASSERT_NE(0.0, res_vector.Norml2()); + + auto jacobianWeights = [&](size_t i) { + std::vector tangents(all_states.size()); + tangents[i] = 1.0; + return tangents; + }; + + auto selectStates = [&](size_t i) { + auto pts = getPointers(v_rhs_states, v_rhs_params); + for (size_t j = 0; j < pts.size(); ++j) { + if (i != j) { + pts[j] = nullptr; + } + } + return pts; + }; + + serac::FiniteElementDual jvp_slow(states[DISP].space(), "jvp_slow"); + serac::FiniteElementDual jvp(states[DISP].space(), "jvp"); + jvp = 4.0; // set to some value to test jvp resets these values + std::vector jvps = getPointers(jvp); + + auto all_v_rhs_states = getPointers(v_rhs_states, v_rhs_params); + + for (size_t i = 0; i < all_states.size(); ++i) { + auto J = residual->jacobian(time, dt, all_states, jacobianWeights(i)); + J->Mult(*all_v_rhs_states[i], jvp_slow); + residual->jvp(time, dt, all_states, selectStates(i), jvps); + EXPECT_NEAR(jvp_slow.Norml2(), jvp.Norml2(), 1e-12); + } + + // test jacobians in weighted combinations + { + all_v_rhs_states[SHAPE] = nullptr; + all_v_rhs_states[VELO] = nullptr; + all_v_rhs_states[4] = nullptr; + + double acceleration_factor = 0.2; + std::vector jacobian_weights = {0.0, 1.0, 0.0, acceleration_factor, 0.0}; + + auto J = residual->jacobian(time, dt, all_states, jacobian_weights); + J->Mult(*all_v_rhs_states[DISP], jvp_slow); + + *all_v_rhs_states[ACCEL] = *all_v_rhs_states[DISP]; + *all_v_rhs_states[ACCEL] *= acceleration_factor; + + residual->jvp(time, dt, all_states, all_v_rhs_states, jvps); + EXPECT_NEAR(jvp_slow.Norml2(), jvp.Norml2(), 1e-12); + } +} + +int main(int argc, char* argv[]) +{ + ::testing::InitGoogleTest(&argc, argv); + serac::ApplicationManager applicationManager(argc, argv); + return RUN_ALL_TESTS(); +} From 34f86a9a63a71a5a4f6f4a9ae2b55e6002eb3d6f Mon Sep 17 00:00:00 2001 From: "E. B. Chin" Date: Fri, 20 Jun 2025 22:58:56 -0700 Subject: [PATCH 05/71] add more features --- src/serac/physics/solid_dfem_residual.hpp | 124 +++++++++++++++--- .../physics/tests/test_dfem_vs_functional.cpp | 14 +- 2 files changed, 116 insertions(+), 22 deletions(-) diff --git a/src/serac/physics/solid_dfem_residual.hpp b/src/serac/physics/solid_dfem_residual.hpp index d707ff52b..76d038c5c 100644 --- a/src/serac/physics/solid_dfem_residual.hpp +++ b/src/serac/physics/solid_dfem_residual.hpp @@ -17,18 +17,21 @@ namespace serac { -template -struct MomentumQFunction { +template +struct StressDivQFunction { SERAC_HOST_DEVICE inline auto operator()( - mfem::real_t dt, const mfem::future::tensor& internal_state, - const mfem::future::tensor& du_dxi, - const mfem::future::tensor& dv_dxi, const mfem::future::tensor&, - const mfem::future::tensor& dX_dxi, mfem::real_t weight) const + // mfem::real_t dt, // TODO: figure out how to pass this in + const mfem::future::tensor& du_dxi, + const mfem::future::tensor& dv_dxi, + const mfem::future::tensor&, + const mfem::future::tensor& dX_dxi, mfem::real_t weight, + Parameters... params) const { auto dxi_dX = mfem::future::inv(dX_dxi); auto du_dX = mfem::future::dot(du_dxi, dxi_dX); auto dv_dX = mfem::future::dot(dv_dxi, dxi_dX); - auto P = mfem::future::get<0>(material.pkStress(dt, internal_state, du_dX, dv_dxi)); + double dt = 1.0; // TODO: figure out how to pass this in to the qfunction + auto P = mfem::future::get<0>(material.pkStress(dt, du_dX, dv_dxi, params...)); auto JxW = mfem::future::det(dX_dxi) * weight * mfem::future::transpose(dxi_dX); return mfem::future::tuple{P * JxW}; } @@ -49,7 +52,17 @@ template class SolidDfemResidual : public DfemResidual { public: /// @brief disp, velo, accel - static constexpr int NUM_STATE_VARS = 3; + static constexpr int NUM_STATE_VARS = 4; + + enum FieldIDs + { + DISP, ///< displacement + VELO, ///< velocity + ACCEL, ///< acceleration + COORD, ///< coordinates + PARAM, ///< parameter field, if any + STATE ///< internal state field, if any + }; /** * @brief Construct a new SolidResidual object @@ -58,9 +71,14 @@ class SolidDfemResidual : public DfemResidual { * @param mesh The serac Mesh * @param diff_op A differentiable operator that computes the residual and jacobian */ - SolidDfemResidual(std::string physics_name, std::shared_ptr mesh, - mfem::future::DifferentiableOperator&& diff_op) - : DfemResidual(physics_name, mesh, std::move(diff_op)) + SolidDfemResidual(std::string physics_name, std::shared_ptr mesh, const mfem::ParFiniteElementSpace& test_space, + std::vector parameter_fe_spaces = {}) + : DfemResidual( + physics_name, mesh, + mfem::future::DifferentiableOperator(solutionFieldDescriptors(mesh->mfemParMesh(), test_space), + parameterFieldDescriptors(parameter_fe_spaces), mesh->mfemParMesh())), + solution_fields_(solutionFieldDescriptors(mesh->mfemParMesh(), test_space)), + parameter_fields_(parameterFieldDescriptors(parameter_fe_spaces)) { } @@ -89,15 +107,87 @@ class SolidDfemResidual : public DfemResidual { * */ template - void setMaterial(std::string, const MaterialType& material, const mfem::IntegrationRule& displacement_ir) + void setMaterial(std::string /*body_name*/, const MaterialType& material, + const mfem::IntegrationRule& displacement_ir) { - mfem::future::tuple inputs{mfem::future::Gradient<0>{}, mfem::future::Gradient<1>{}, mfem::future::Gradient<2>{}, - mfem::future::Gradient<3>{}, mfem::future::Weight{}}; - mfem::future::tuple outputs{mfem::future::Gradient<4>{}}; + mfem::future::tuple outputs{mfem::future::Gradient{}}; + // TODO: find out the right attributes from the body name mfem::Array solid_domain_attributes(DfemResidual::mesh_->mfemParMesh().attributes.Max()); - auto momentum_qf = MomentumQFunction{.material = material}; - DfemResidual::diff_op_.AddDomainIntegrator(momentum_qf, inputs, outputs, displacement_ir, solid_domain_attributes); + if constexpr (MaterialType::has_state && MaterialType::has_parameters) { + // build new diff_op_ with state and parameter fields + mfem::future::tuple inputs{// TODO: figure out how to pass in dt + mfem::future::Gradient{}, mfem::future::Gradient{}, + mfem::future::Gradient{}, mfem::future::Gradient{}, + mfem::future::Weight{}, mfem::future::Value{}, + mfem::future::Value{}}; + auto stress_div_qf = + StressDivQFunction{ + .material = material}; + DfemResidual::diff_op_.AddDomainIntegrator(stress_div_qf, inputs, outputs, displacement_ir, + solid_domain_attributes); + } else if constexpr (MaterialType::has_parameters) { + // build new diff_op_ with parameter fields + mfem::future::tuple inputs{// TODO: figure out how to pass in dt + mfem::future::Gradient{}, mfem::future::Gradient{}, + mfem::future::Gradient{}, mfem::future::Gradient{}, + mfem::future::Weight{}, mfem::future::Value{}}; + auto stress_div_qf = StressDivQFunction{.material = material}; + DfemResidual::diff_op_.AddDomainIntegrator(stress_div_qf, inputs, outputs, displacement_ir, + solid_domain_attributes); + } else if constexpr (MaterialType::has_state) { + // build new diff_op_ with state fields + mfem::future::tuple inputs{// TODO: figure out how to pass in dt + mfem::future::Gradient{}, mfem::future::Gradient{}, + mfem::future::Gradient{}, mfem::future::Gradient{}, + mfem::future::Weight{}, mfem::future::Value{}}; + auto stress_div_qf = StressDivQFunction{.material = material}; + DfemResidual::diff_op_.AddDomainIntegrator(stress_div_qf, inputs, outputs, displacement_ir, + solid_domain_attributes); + } else { + // existing diff_op_ should work fine + mfem::future::tuple inputs{// TODO: figure out how to pass in dt + mfem::future::Gradient{}, mfem::future::Gradient{}, + mfem::future::Gradient{}, mfem::future::Gradient{}, + mfem::future::Weight{}}; + auto stress_div_qf = StressDivQFunction{.material = material}; + DfemResidual::diff_op_.AddDomainIntegrator(stress_div_qf, inputs, outputs, displacement_ir, + solid_domain_attributes); + } + } + + protected: + /** + * @brief Construct the field descriptors for the differentiable operator + * + * @param parameter_fe_spaces The parameter finite element spaces + * @return A vector of field descriptors for the differentiable operator + */ + static std::vector parameterFieldDescriptors( + std::vector parameter_fe_spaces) + { + std::vector field_descriptors; + field_descriptors.reserve(parameter_fe_spaces.size()); + for (const auto& space : parameter_fe_spaces) { + field_descriptors.emplace_back(NUM_STATE_VARS + field_descriptors.size(), space); + } + return field_descriptors; + } + /** + * @brief Construct the field descriptors for the solution fields + * + */ + static std::vector solutionFieldDescriptors( + const mfem::ParMesh& mesh, const mfem::ParFiniteElementSpace& test_space) + { + return { + {DISP, &test_space}, + {VELO, &test_space}, + {ACCEL, &test_space}, + {COORD, static_cast(mesh.GetNodalFESpace())}, + }; } + std::vector solution_fields_; + std::vector parameter_fields_; }; } // namespace serac diff --git a/src/serac/physics/tests/test_dfem_vs_functional.cpp b/src/serac/physics/tests/test_dfem_vs_functional.cpp index 4accf3544..54f6c073d 100644 --- a/src/serac/physics/tests/test_dfem_vs_functional.cpp +++ b/src/serac/physics/tests/test_dfem_vs_functional.cpp @@ -93,11 +93,15 @@ struct NeoHookeanWithFieldWithRateFunctional { struct NeoHookeanWithFieldWithRateDfem { static constexpr int state_size = 0; + static constexpr bool has_parameters = true; + static constexpr bool has_state = false; - template - SERAC_HOST_DEVICE auto pkStress(double /*dt*/, const mfem::future::tensor& /* state */, - const mfem::future::tensor& du_dX, - const mfem::future::tensor& /*dv_dX*/) const + using param_type = double; + std::unique_ptr fe_collection = std::make_unique(0); + + template typename TensorT> + SERAC_HOST_DEVICE auto pkStress(double /*dt*/, const TensorT& du_dX, + const TensorT& /*dv_dX*/, double K) const { using std::log1p; constexpr auto I = mfem::future::IdentityMatrix(); @@ -115,7 +119,7 @@ struct NeoHookeanWithFieldWithRateDfem { SERAC_HOST_DEVICE auto density() const { return Rho; } - double K; ///< bulk modulus + // double K; ///< bulk modulus double G; ///< shear modulus double Rho; ///< density }; From 734a6de29bea787437fba44b041e26ecf95c5401 Mon Sep 17 00:00:00 2001 From: "E. B. Chin" Date: Mon, 23 Jun 2025 17:33:47 -0700 Subject: [PATCH 06/71] start testing a diff_op for each term --- src/serac/physics/dfem_residual2.hpp | 153 ++++++++++++++ src/serac/physics/solid_dfem_residual2.hpp | 189 ++++++++++++++++++ .../physics/tests/test_dfem_vs_functional.cpp | 4 +- 3 files changed, 344 insertions(+), 2 deletions(-) create mode 100644 src/serac/physics/dfem_residual2.hpp create mode 100644 src/serac/physics/solid_dfem_residual2.hpp diff --git a/src/serac/physics/dfem_residual2.hpp b/src/serac/physics/dfem_residual2.hpp new file mode 100644 index 000000000..476a0eb3d --- /dev/null +++ b/src/serac/physics/dfem_residual2.hpp @@ -0,0 +1,153 @@ +// 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) + +/** + * @file dfem_residual2.hpp + */ + +#pragma once + +#include "mfem.hpp" + +#include "serac/physics/residual.hpp" +#include "serac/physics/mesh.hpp" +#include "serac/physics/state/finite_element_state.hpp" +#include "serac/physics/state/finite_element_dual.hpp" + +namespace serac { + +/** + * @brief The nonlinear residual class + * + * This uses dFEM to compute fairly arbitrary residuals and tangent + * stiffness matrices based on body and boundary integrals. + * + */ +class DfemResidual2 : public Residual { + public: + using SpacesT = std::vector; ///< typedef + + /** + * @brief Construct a new DfemResidual object + * + * @param physics_name A name for the physics module instance + * @param mesh The serac mesh + * @param diff_op A differentiable operator that computes the residual and jacobian + */ + DfemResidual2(std::string physics_name, std::shared_ptr mesh) : Residual(physics_name), mesh_(mesh) {} + + /** + * @brief Add a body integral contribution to the residual + * + * @tparam BodyIntegralType The type of the body integral + * // DependsOn can be indices into fields which the body integral may depend on + * @param body_name The name of the registered domain over which the body force is applied. If nothing is supplied + * the entire domain is + * @param body_integral A function describing the body force applied + * used. + * @pre body_integral must be a object that can be called with the following arguments: + * 1. `double t` the time + * 2. `tensor x` the spatial coordinates for the quadrature point + * 3. `tuple{value, derivative}`, a variadic list of tuples (each with a values and derivative), + * one tuple for each of the trial spaces specified in the `DependsOn<...>` argument. + * @note The actual types of these arguments passed will be `double`, `tensor` or tuples thereof + * when doing direct evaluation. When differentiating with respect to one of the inputs, its stored + * values will change to `dual` numbers rather than `double`. (e.g. `tensor` becomes `tensor, + * 3>`) + * + */ + template + void addBodyIntegral(QFunctionType qfunction, const std::vector& solution_fields, + const std::vector& param_fields, InputType qfunction_inputs, + OutputType qfunction_outputs, const mfem::IntegrationRule& integration_rule, + mfem::Array domain_attributes, std::integer_sequence derivative_ids) + { + // get field IDs + std::vector field_ids; + field_ids.reserve(solution_fields.size() + param_fields.size()); + for (const auto& field : solution_fields) { + field_ids.push_back(field.id); + } + for (const auto& field : param_fields) { + field_ids.push_back(field.id); + } + // just keep the first input field as a solution field, make the rest parameter fields to simplify the Mult() call + std::vector other_fields; + other_fields.reserve(param_fields.size() + solution_fields.size() - 1); + for (size_t i = 1; i < solution_fields.size(); ++i) { + other_fields.push_back(solution_fields[i]); + } + for (const auto& field : param_fields) { + other_fields.push_back(field); + } + diff_ops_with_field_ids_.emplace_back( + std::make_tuple(mfem::future::DifferentiableOperator({solution_fields[0]}, other_fields, mesh_->mfemParMesh()), + std::move(field_ids))); + std::get<0>(diff_ops_with_field_ids_.back()) + .AddDomainIntegrator(qfunction, qfunction_inputs, qfunction_outputs, integration_rule, domain_attributes, + derivative_ids); + } + + /// @overload + mfem::Vector residual(double, double, const std::vector& fields, int block_row = 0) const override + { + SLIC_ERROR_IF(block_row != 0, "Invalid block row and column requested in fieldJacobian for DfemResidual2"); + mfem::Vector resid(fields[0]->space().GetVSize()); + resid = 0.0; + for (auto& diff_op_with_field_ids : diff_ops_with_field_ids_) { + std::vector other_fields; + auto& field_ids = std::get<1>(diff_op_with_field_ids); + other_fields.reserve(field_ids.size() - 1); + for (size_t i = 1; i < field_ids.size(); ++i) { + other_fields.push_back(&fields[field_ids[i]]->gridFunction()); + } + std::get<0>(diff_op_with_field_ids).SetParameters(other_fields); + mfem::Vector term_resid(fields[0]->space().GetTrueVSize()); + std::get<0>(diff_op_with_field_ids).Mult(*fields[0], term_resid); + resid += term_resid; + } + return resid; + } + + /// @overload + std::unique_ptr jacobian(double, double, const std::vector&, + const std::vector&, int) const override + { + return std::make_unique(); + } + + /// @overload + void jvp(double, double, const std::vector& fields, const std::vector& v_fields, + const std::vector& jvp_reactions) const override + { + SLIC_ERROR_IF(v_fields.size() != fields.size(), + "Invalid number of field sensitivities relative to the number of fields"); + SLIC_ERROR_IF(jvp_reactions.size() != 1, "DfemResidual nonlinear systems only supports 1 output residual"); + + auto grad_op = diff_op_.GetDerivative(0, {fields[0]}, {mesh_nodes_}); + grad_op->AddMult(*v_fields[0], *jvp_reactions[0]); + } + + /// @overload + void vjp(double, double, const std::vector& fields, const std::vector& v_fields, + const std::vector& vjp_sensitivities) const override + { + SLIC_ERROR_IF(vjp_sensitivities.size() != fields.size(), + "Invalid number of field sensitivities relative to the number of fields"); + SLIC_ERROR_IF(v_fields.size() != 1, "FunctionalResidual nonlinear systems only supports 1 output residual"); + + auto grad_op = diff_op_.GetDerivative(0, {fields[0]}, {mesh_nodes_}); + grad_op->AddMultTranspose(*v_fields[0], *vjp_sensitivities[0]); + } + + protected: + /// @brief primary mesh + std::shared_ptr mesh_; + + mutable std::vector>> diff_ops_with_field_ids_; +}; + +} // namespace serac diff --git a/src/serac/physics/solid_dfem_residual2.hpp b/src/serac/physics/solid_dfem_residual2.hpp new file mode 100644 index 000000000..079392620 --- /dev/null +++ b/src/serac/physics/solid_dfem_residual2.hpp @@ -0,0 +1,189 @@ +// 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) + +/** + * @file solid_dfem_residual.hpp + * + * @brief Implements the residual interface for solid mechanics physics. + * Derives from dfem_residual. + */ + +#pragma once + +#include "serac/physics/dfem_residual2.hpp" + +namespace serac { + +template +struct StressDivQFunction { + SERAC_HOST_DEVICE inline auto operator()( + // mfem::real_t dt, // TODO: figure out how to pass this in + const mfem::future::tensor& du_dxi, + const mfem::future::tensor& dv_dxi, + const mfem::future::tensor&, + const mfem::future::tensor& dX_dxi, mfem::real_t weight, + Parameters... params) const + { + auto dxi_dX = mfem::future::inv(dX_dxi); + auto du_dX = mfem::future::dot(du_dxi, dxi_dX); + auto dv_dX = mfem::future::dot(dv_dxi, dxi_dX); + double dt = 1.0; // TODO: figure out how to pass this in to the qfunction + auto P = mfem::future::get<0>(material.pkStress(dt, du_dX, dv_dxi, params...)); + auto JxW = mfem::future::det(dX_dxi) * weight * mfem::future::transpose(dxi_dX); + return mfem::future::tuple{P * JxW}; + } + + Material material; ///< the material model to use for computing the stress +}; + +/** + * @brief The nonlinear residual class + * + * This uses Functional to compute the solid mechanics residuals and tangent + * stiffness matrices. + * + * @tparam order The order of the discretization of the displacement and velocity fields + * @tparam dim The spatial dimension of the mesh + */ +template +class SolidDfemResidual2 : public DfemResidual2 { + public: + /// @brief disp, velo, accel + static constexpr int NUM_STATE_VARS = 4; + + enum FieldIDs + { + DISP, ///< displacement + VELO, ///< velocity + ACCEL, ///< acceleration + COORD, ///< coordinates + PARAM, ///< parameter field, if any + STATE ///< internal state field, if any + }; + + /** + * @brief Construct a new SolidResidual object + * + * @param physics_name A name for the physics module instance + * @param mesh The serac Mesh + */ + SolidDfemResidual2(std::string physics_name, std::shared_ptr mesh, + const mfem::ParFiniteElementSpace& test_space) + : DfemResidual2(physics_name, mesh), test_space_(test_space) + { + } + + /** + * @brief Set the material stress response and mass properties for the physics module + * + * @tparam MaterialType The solid material type + * @tparam StateType the type that contains the internal variables for MaterialType + * @param body_name string name for a registered body Domain on the mesh + * @param material A material that provides a function to evaluate stress + * @pre material must be a object that can be called with the following arguments: + * 1. `MaterialType::State & state` an mutable reference to the internal variables for this quadrature point + * 2. `tensor du_dx` the displacement gradient at this quadrature point + * 3. `tuple{value, derivative}`, a tuple of values and derivatives for each parameter field + * specified in the `DependsOn<...>` argument. + * + * @note The actual types of these arguments passed will be `double`, `tensor` or tuples thereof + * when doing direct evaluation. When differentiating with respect to one of the inputs, its stored + * values will change to `dual` numbers rather than `double`. (e.g. `tensor` becomes `tensor, + * 3>`) + * + * @param qdata the buffer of material internal variables at each quadrature point + * + * @pre MaterialType must have a public method `density` which can take FiniteElementState parameter inputs + * @pre MaterialType must have a public method 'pkStress' which returns the first Piola-Kirchhoff stress + * + */ + template + void setMaterial(std::string /*body_name*/, const MaterialType& material, + const mfem::IntegrationRule& displacement_ir) + { + mfem::future::tuple outputs{mfem::future::Gradient{}}; + // TODO: find out the right attributes from the body name + mfem::Array solid_domain_attributes(DfemResidual2::mesh_->mfemParMesh().attributes.Max()); + // if constexpr (MaterialType::has_state && MaterialType::has_parameters) { + // // build new diff_op_ with state and parameter fields + // mfem::future::tuple inputs{// TODO: figure out how to pass in dt + // mfem::future::Gradient{}, mfem::future::Gradient{}, + // mfem::future::Gradient{}, mfem::future::Gradient{}, + // mfem::future::Weight{}, mfem::future::Value{}, + // mfem::future::Value{}}; + // auto stress_div_qf = + // StressDivQFunction{ + // .material = material}; + // DfemResidual2::addBodyIntegral(stress_div_qf, ) + // DfemResidual::diff_op_.AddDomainIntegrator(stress_div_qf, inputs, outputs, displacement_ir, + // solid_domain_attributes); + // } else + if constexpr (MaterialType::num_parameters > 0) { + mfem::future::tuple inputs{// TODO: figure out how to pass in dt + mfem::future::Gradient{}, mfem::future::Gradient{}, + mfem::future::Gradient{}, mfem::future::Gradient{}, + mfem::future::Weight{}, mfem::future::Value{}}; + auto stress_div_qf = StressDivQFunction{.material = material}; + DfemResidual::diff_op_.AddDomainIntegrator(stress_div_qf, inputs, outputs, displacement_ir, + solid_domain_attributes); + } + // else if constexpr (MaterialType::has_state) { + // // build new diff_op_ with state fields + // mfem::future::tuple inputs{// TODO: figure out how to pass in dt + // mfem::future::Gradient{}, mfem::future::Gradient{}, + // mfem::future::Gradient{}, mfem::future::Gradient{}, + // mfem::future::Weight{}, mfem::future::Value{}}; + // auto stress_div_qf = StressDivQFunction{.material = material}; + // DfemResidual::diff_op_.AddDomainIntegrator(stress_div_qf, inputs, outputs, displacement_ir, + // solid_domain_attributes); + // } + else { + // existing diff_op_ should work fine + mfem::future::tuple inputs{// TODO: figure out how to pass in dt + mfem::future::Gradient{}, mfem::future::Gradient{}, + mfem::future::Gradient{}, mfem::future::Gradient{}, + mfem::future::Weight{}}; + auto stress_div_qf = StressDivQFunction{.material = material}; + DfemResidual::diff_op_.AddDomainIntegrator(stress_div_qf, inputs, outputs, displacement_ir, + solid_domain_attributes); + } + } + + protected: + /** + * @brief Construct the field descriptors for the differentiable operator + * + * @param parameter_fe_spaces The parameter finite element spaces + * @return A vector of field descriptors for the differentiable operator + */ + static std::vector parameterFieldDescriptors( + std::vector parameter_fe_spaces) + { + std::vector field_descriptors; + field_descriptors.reserve(parameter_fe_spaces.size()); + for (const auto& space : parameter_fe_spaces) { + field_descriptors.emplace_back(NUM_STATE_VARS + field_descriptors.size(), space); + } + return field_descriptors; + } + /** + * @brief Construct the field descriptors for the solution fields + * + */ + static std::vector solutionFieldDescriptors( + const mfem::ParMesh& mesh, const mfem::ParFiniteElementSpace& test_space) + { + return { + {DISP, &test_space}, + {VELO, &test_space}, + {ACCEL, &test_space}, + {COORD, static_cast(mesh.GetNodalFESpace())}, + }; + } + const mfem::ParFiniteElementSpace& test_space_; +}; + +} // namespace serac diff --git a/src/serac/physics/tests/test_dfem_vs_functional.cpp b/src/serac/physics/tests/test_dfem_vs_functional.cpp index 54f6c073d..b9b74ae9c 100644 --- a/src/serac/physics/tests/test_dfem_vs_functional.cpp +++ b/src/serac/physics/tests/test_dfem_vs_functional.cpp @@ -93,8 +93,8 @@ struct NeoHookeanWithFieldWithRateFunctional { struct NeoHookeanWithFieldWithRateDfem { static constexpr int state_size = 0; - static constexpr bool has_parameters = true; - static constexpr bool has_state = false; + static constexpr int num_parameters = 1; + static constexpr int num_states = 0; using param_type = double; std::unique_ptr fe_collection = std::make_unique(0); From ff9074896ee0e6f968e9dec452f656e8754b4465 Mon Sep 17 00:00:00 2001 From: "E. B. Chin" Date: Tue, 24 Jun 2025 16:31:49 -0700 Subject: [PATCH 07/71] set params at construction --- src/serac/physics/solid_dfem_residual2.hpp | 81 +++++++--------------- 1 file changed, 25 insertions(+), 56 deletions(-) diff --git a/src/serac/physics/solid_dfem_residual2.hpp b/src/serac/physics/solid_dfem_residual2.hpp index 079392620..c5b1f471b 100644 --- a/src/serac/physics/solid_dfem_residual2.hpp +++ b/src/serac/physics/solid_dfem_residual2.hpp @@ -59,9 +59,7 @@ class SolidDfemResidual2 : public DfemResidual2 { DISP, ///< displacement VELO, ///< velocity ACCEL, ///< acceleration - COORD, ///< coordinates - PARAM, ///< parameter field, if any - STATE ///< internal state field, if any + COORD ///< coordinates }; /** @@ -71,9 +69,14 @@ class SolidDfemResidual2 : public DfemResidual2 { * @param mesh The serac Mesh */ SolidDfemResidual2(std::string physics_name, std::shared_ptr mesh, - const mfem::ParFiniteElementSpace& test_space) + const mfem::ParFiniteElementSpace& test_space, + const std::vector& parameter_fe_spaces = {}) : DfemResidual2(physics_name, mesh), test_space_(test_space) { + parameter_fields_.reserve(parameter_fe_spaces.size()); + for (auto space : parameter_fe_spaces) { + parameter_fields_.emplace_back(NUM_STATE_VARS + parameter_fields_.size(), space); + } } /** @@ -100,56 +103,22 @@ class SolidDfemResidual2 : public DfemResidual2 { * @pre MaterialType must have a public method 'pkStress' which returns the first Piola-Kirchhoff stress * */ - template - void setMaterial(std::string /*body_name*/, const MaterialType& material, + template + void setMaterial(DependsOn, std::string /*body_name*/, const MaterialType& material, const mfem::IntegrationRule& displacement_ir) { - mfem::future::tuple outputs{mfem::future::Gradient{}}; + mfem::future::tuple qfunction_inputs{ + // TODO: figure out how to pass in dt + mfem::future::Gradient{}, mfem::future::Gradient{}, + mfem::future::Gradient{}, mfem::future::Gradient{}, + mfem::future::Weight{}, mfem::future::Value{}...}; + mfem::future::tuple qfunction_outputs{mfem::future::Gradient{}}; // TODO: find out the right attributes from the body name mfem::Array solid_domain_attributes(DfemResidual2::mesh_->mfemParMesh().attributes.Max()); - // if constexpr (MaterialType::has_state && MaterialType::has_parameters) { - // // build new diff_op_ with state and parameter fields - // mfem::future::tuple inputs{// TODO: figure out how to pass in dt - // mfem::future::Gradient{}, mfem::future::Gradient{}, - // mfem::future::Gradient{}, mfem::future::Gradient{}, - // mfem::future::Weight{}, mfem::future::Value{}, - // mfem::future::Value{}}; - // auto stress_div_qf = - // StressDivQFunction{ - // .material = material}; - // DfemResidual2::addBodyIntegral(stress_div_qf, ) - // DfemResidual::diff_op_.AddDomainIntegrator(stress_div_qf, inputs, outputs, displacement_ir, - // solid_domain_attributes); - // } else - if constexpr (MaterialType::num_parameters > 0) { - mfem::future::tuple inputs{// TODO: figure out how to pass in dt - mfem::future::Gradient{}, mfem::future::Gradient{}, - mfem::future::Gradient{}, mfem::future::Gradient{}, - mfem::future::Weight{}, mfem::future::Value{}}; - auto stress_div_qf = StressDivQFunction{.material = material}; - DfemResidual::diff_op_.AddDomainIntegrator(stress_div_qf, inputs, outputs, displacement_ir, - solid_domain_attributes); - } - // else if constexpr (MaterialType::has_state) { - // // build new diff_op_ with state fields - // mfem::future::tuple inputs{// TODO: figure out how to pass in dt - // mfem::future::Gradient{}, mfem::future::Gradient{}, - // mfem::future::Gradient{}, mfem::future::Gradient{}, - // mfem::future::Weight{}, mfem::future::Value{}}; - // auto stress_div_qf = StressDivQFunction{.material = material}; - // DfemResidual::diff_op_.AddDomainIntegrator(stress_div_qf, inputs, outputs, displacement_ir, - // solid_domain_attributes); - // } - else { - // existing diff_op_ should work fine - mfem::future::tuple inputs{// TODO: figure out how to pass in dt - mfem::future::Gradient{}, mfem::future::Gradient{}, - mfem::future::Gradient{}, mfem::future::Gradient{}, - mfem::future::Weight{}}; - auto stress_div_qf = StressDivQFunction{.material = material}; - DfemResidual::diff_op_.AddDomainIntegrator(stress_div_qf, inputs, outputs, displacement_ir, - solid_domain_attributes); - } + auto stress_div_qf = StressDivQFunction{.material = material}; + DfemResidual2::addBodyIntegral(stress_div_qf, solutionFieldDescriptors(), parameter_fields_, qfunction_inputs, + qfunction_outputs, displacement_ir, solid_domain_attributes, + std::integer_sequence{}); } protected: @@ -173,17 +142,17 @@ class SolidDfemResidual2 : public DfemResidual2 { * @brief Construct the field descriptors for the solution fields * */ - static std::vector solutionFieldDescriptors( - const mfem::ParMesh& mesh, const mfem::ParFiniteElementSpace& test_space) + std::vector solutionFieldDescriptors() { return { - {DISP, &test_space}, - {VELO, &test_space}, - {ACCEL, &test_space}, - {COORD, static_cast(mesh.GetNodalFESpace())}, + {DISP, &test_space_}, + {VELO, &test_space_}, + {ACCEL, &test_space_}, + {COORD, static_cast(DfemResidual2::mesh_->mfemParMesh().GetNodalFESpace())}, }; } const mfem::ParFiniteElementSpace& test_space_; + std::vector parameter_fields_; }; } // namespace serac From a22d0fa899ca5f51c84c69ccaba0eea2462c77b7 Mon Sep 17 00:00:00 2001 From: "E. B. Chin" Date: Tue, 24 Jun 2025 23:05:53 -0700 Subject: [PATCH 08/71] cleanup --- src/serac/physics/dfem_residual.hpp | 138 ++++++++++++++---- src/serac/physics/dfem_residual2.hpp | 153 -------------------- src/serac/physics/solid_dfem_residual.hpp | 123 ++++++---------- src/serac/physics/solid_dfem_residual2.hpp | 158 --------------------- 4 files changed, 149 insertions(+), 423 deletions(-) delete mode 100644 src/serac/physics/dfem_residual2.hpp delete mode 100644 src/serac/physics/solid_dfem_residual2.hpp diff --git a/src/serac/physics/dfem_residual.hpp b/src/serac/physics/dfem_residual.hpp index 302323e2e..21d08c7cb 100644 --- a/src/serac/physics/dfem_residual.hpp +++ b/src/serac/physics/dfem_residual.hpp @@ -5,7 +5,7 @@ // SPDX-License-Identifier: (BSD-3-Clause) /** - * @file dfem_residual.hpp + * @file dfem_residual2.hpp */ #pragma once @@ -37,10 +37,7 @@ class DfemResidual : public Residual { * @param mesh The serac mesh * @param diff_op A differentiable operator that computes the residual and jacobian */ - DfemResidual(std::string physics_name, std::shared_ptr mesh, mfem::future::DifferentiableOperator&& diff_op) - : Residual(physics_name), mesh_(mesh), diff_op_(std::move(diff_op)) - { - } + DfemResidual(std::string physics_name, std::shared_ptr mesh) : Residual(physics_name), mesh_(mesh) {} /** * @brief Add a body integral contribution to the residual @@ -63,51 +60,131 @@ class DfemResidual : public Residual { * */ template - void addBodyIntegral(QFunctionType qfunction, InputType inputs, OutputType outputs, - const mfem::IntegrationRule& integration_rule, mfem::Array domain_attributes, - std::integer_sequence derivative_ids) + void addBodyIntegral(QFunctionType qfunction, const std::vector& solution_fields, + const std::vector& param_fields, InputType qfunction_inputs, + OutputType qfunction_outputs, const mfem::IntegrationRule& integration_rule, + mfem::Array domain_attributes, std::integer_sequence derivative_ids) { - diff_op_.AddDomainIntegrator(qfunction, inputs, outputs, integration_rule, domain_attributes, derivative_ids); + // get field IDs + std::vector field_ids; + field_ids.reserve(solution_fields.size() + param_fields.size()); + for (const auto& field : solution_fields) { + field_ids.push_back(field.id); + } + for (const auto& field : param_fields) { + field_ids.push_back(field.id); + } + // just keep the first input field as a solution field, make the rest parameter fields to simplify the Mult() call + std::vector other_fields; + other_fields.reserve(param_fields.size() + solution_fields.size() - 1); + for (size_t i = 1; i < solution_fields.size(); ++i) { + other_fields.push_back(solution_fields[i]); + } + for (const auto& field : param_fields) { + other_fields.push_back(field); + } + residual_terms_with_field_ids_.emplace_back( + std::make_tuple(mfem::future::DifferentiableOperator({solution_fields[0]}, other_fields, mesh_->mfemParMesh()), + std::move(field_ids))); + std::get<0>(residual_terms_with_field_ids_.back()) + .AddDomainIntegrator(qfunction, qfunction_inputs, qfunction_outputs, integration_rule, domain_attributes, + derivative_ids); } /// @overload - mfem::Vector residual(double, double, const std::vector& fields, int block_row = 0) const override + mfem::Vector residual(double, double, const std::vector& fields, int block_row = 0) const override { - SLIC_ERROR_IF(block_row != 0, "Invalid block row and column requested in fieldJacobian for FunctionalResidual"); - diff_op_.SetParameters({mesh_nodes_}); - mfem::Vector ret(fields[0]->space().GetVSize()); - diff_op_.Mult(fields[0]->gridFunction(), ret); - return ret; + SLIC_ERROR_ROOT_IF(block_row != 0, "Invalid block row and column requested in fieldJacobian for DfemResidual"); + mfem::Vector resid(fields[0]->space().GetVSize()); + resid = 0.0; + for (auto& residual_term_with_field_ids : residual_terms_with_field_ids_) { + std::vector other_fields; + auto& field_ids = std::get<1>(residual_term_with_field_ids); + other_fields.reserve(field_ids.size() - 1); + for (size_t i = 1; i < field_ids.size(); ++i) { + other_fields.push_back(&fields[field_ids[i]]->gridFunction()); + } + std::get<0>(residual_term_with_field_ids).SetParameters(other_fields); + mfem::Vector term_resid(fields[0]->space().GetTrueVSize()); + std::get<0>(residual_term_with_field_ids).Mult(*fields[0], term_resid); + resid += term_resid; + } + return resid; } /// @overload - std::unique_ptr jacobian(double, double, const std::vector&, + std::unique_ptr jacobian(double, double, const std::vector&, const std::vector&, int) const override { + SLIC_ERROR_ROOT("DfemResidual does not support matrix assembly"); return std::make_unique(); } /// @overload - void jvp(double, double, const std::vector& fields, const std::vector& v_fields, + void jvp(double, double, const std::vector& fields, const std::vector& v_fields, const std::vector& jvp_reactions) const override { - SLIC_ERROR_IF(v_fields.size() != fields.size(), - "Invalid number of field sensitivities relative to the number of fields"); - SLIC_ERROR_IF(jvp_reactions.size() != 1, "DfemResidual nonlinear systems only supports 1 output residual"); - - auto grad_op = diff_op_.GetDerivative(0, {fields[0]}, {mesh_nodes_}); - grad_op->AddMult(*v_fields[0], *jvp_reactions[0]); + SLIC_ERROR_ROOT_IF(v_fields.size() != fields.size(), + "Invalid number of field sensitivities relative to the number of fields"); + SLIC_ERROR_ROOT_IF(jvp_reactions.size() != 1, "DfemResidual nonlinear systems only supports 1 output residual"); + + *jvp_reactions[0] = 0.0; + + for (size_t input_col = 0; input_col < fields.size(); ++input_col) { + if (v_fields[input_col] == nullptr) { + continue; + } + for (auto& residual_term_with_field_ids : residual_terms_with_field_ids_) { + const auto& field_ids = std::get<1>(residual_term_with_field_ids); + auto field_id_it = std::find(field_ids.begin(), field_ids.end(), input_col); + // this residual term does not depend on this input field + if (field_id_it == field_ids.end()) { + continue; + } + // if (field_ids[0] != input_col) { + // continue; + // } + // auto& residual_term_ = std::get<0>(residual_term_with_field_ids); + // std::vector other_fields; + // other_fields.reserve(field_ids.size() - 1); + // for (size_t i = 1; i < field_ids.size(); ++i) { + // other_fields.push_back(&fields[field_ids[i]]->gridFunction()); + // } + // residual_term_.SetParameters(other_fields); + // mfem::Vector K(fields[0]->space().GetTrueVSize()); + // residual_term_.AddMult(*v_fields[input_col], K); + // (*jvp_reactions[0]) += K; + // } + // // find the residual term that corresponds to this input field + // auto it = std::find_if(residual_terms_with_field_ids_.begin(), residual_terms_with_field_ids_.end(), + // [&](const auto& term_and_ids) { + // const auto& field_ids = std::get<1>(term_and_ids); + // return field_ids[0] == input_col; + // }); + // if (it == residual_terms_with_field_ids_.end()) { + // continue; + // } + // const auto& residual_term_ = std::get<0>(*it); + // std::vector other_fields; + // const auto& field_ids = std::get<1>(*it); + // other_fields.reserve(field_ids.size() - 1); + // for (size_t i = 1; i < field_ids.size(); ++i) { + // other_fields.push_back(&fields[field_ids[i]]->gridFunction()); + // } + // residual_term_.SetParameters(other_fields); + } + } } /// @overload - void vjp(double, double, const std::vector& fields, const std::vector& v_fields, + void vjp(double, double, const std::vector& fields, const std::vector& v_fields, const std::vector& vjp_sensitivities) const override { - SLIC_ERROR_IF(vjp_sensitivities.size() != fields.size(), - "Invalid number of field sensitivities relative to the number of fields"); - SLIC_ERROR_IF(v_fields.size() != 1, "FunctionalResidual nonlinear systems only supports 1 output residual"); + SLIC_ERROR_ROOT_IF(vjp_sensitivities.size() != fields.size(), + "Invalid number of field sensitivities relative to the number of fields"); + SLIC_ERROR_ROOT_IF(v_fields.size() != 1, "FunctionalResidual nonlinear systems only supports 1 output residual"); - auto grad_op = diff_op_.GetDerivative(0, {fields[0]}, {mesh_nodes_}); + auto grad_op = residual_term_.GetDerivative(0, {fields[0]}, {mesh_nodes_}); grad_op->AddMultTranspose(*v_fields[0], *vjp_sensitivities[0]); } @@ -115,9 +192,8 @@ class DfemResidual : public Residual { /// @brief primary mesh std::shared_ptr mesh_; - mfem::ParGridFunction* mesh_nodes_; ///< grid function for the mesh nodes - - mutable mfem::future::DifferentiableOperator diff_op_; ///< differentiable operator for the residual + mutable std::vector>> + residual_terms_with_field_ids_; }; } // namespace serac diff --git a/src/serac/physics/dfem_residual2.hpp b/src/serac/physics/dfem_residual2.hpp deleted file mode 100644 index 476a0eb3d..000000000 --- a/src/serac/physics/dfem_residual2.hpp +++ /dev/null @@ -1,153 +0,0 @@ -// 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) - -/** - * @file dfem_residual2.hpp - */ - -#pragma once - -#include "mfem.hpp" - -#include "serac/physics/residual.hpp" -#include "serac/physics/mesh.hpp" -#include "serac/physics/state/finite_element_state.hpp" -#include "serac/physics/state/finite_element_dual.hpp" - -namespace serac { - -/** - * @brief The nonlinear residual class - * - * This uses dFEM to compute fairly arbitrary residuals and tangent - * stiffness matrices based on body and boundary integrals. - * - */ -class DfemResidual2 : public Residual { - public: - using SpacesT = std::vector; ///< typedef - - /** - * @brief Construct a new DfemResidual object - * - * @param physics_name A name for the physics module instance - * @param mesh The serac mesh - * @param diff_op A differentiable operator that computes the residual and jacobian - */ - DfemResidual2(std::string physics_name, std::shared_ptr mesh) : Residual(physics_name), mesh_(mesh) {} - - /** - * @brief Add a body integral contribution to the residual - * - * @tparam BodyIntegralType The type of the body integral - * // DependsOn can be indices into fields which the body integral may depend on - * @param body_name The name of the registered domain over which the body force is applied. If nothing is supplied - * the entire domain is - * @param body_integral A function describing the body force applied - * used. - * @pre body_integral must be a object that can be called with the following arguments: - * 1. `double t` the time - * 2. `tensor x` the spatial coordinates for the quadrature point - * 3. `tuple{value, derivative}`, a variadic list of tuples (each with a values and derivative), - * one tuple for each of the trial spaces specified in the `DependsOn<...>` argument. - * @note The actual types of these arguments passed will be `double`, `tensor` or tuples thereof - * when doing direct evaluation. When differentiating with respect to one of the inputs, its stored - * values will change to `dual` numbers rather than `double`. (e.g. `tensor` becomes `tensor, - * 3>`) - * - */ - template - void addBodyIntegral(QFunctionType qfunction, const std::vector& solution_fields, - const std::vector& param_fields, InputType qfunction_inputs, - OutputType qfunction_outputs, const mfem::IntegrationRule& integration_rule, - mfem::Array domain_attributes, std::integer_sequence derivative_ids) - { - // get field IDs - std::vector field_ids; - field_ids.reserve(solution_fields.size() + param_fields.size()); - for (const auto& field : solution_fields) { - field_ids.push_back(field.id); - } - for (const auto& field : param_fields) { - field_ids.push_back(field.id); - } - // just keep the first input field as a solution field, make the rest parameter fields to simplify the Mult() call - std::vector other_fields; - other_fields.reserve(param_fields.size() + solution_fields.size() - 1); - for (size_t i = 1; i < solution_fields.size(); ++i) { - other_fields.push_back(solution_fields[i]); - } - for (const auto& field : param_fields) { - other_fields.push_back(field); - } - diff_ops_with_field_ids_.emplace_back( - std::make_tuple(mfem::future::DifferentiableOperator({solution_fields[0]}, other_fields, mesh_->mfemParMesh()), - std::move(field_ids))); - std::get<0>(diff_ops_with_field_ids_.back()) - .AddDomainIntegrator(qfunction, qfunction_inputs, qfunction_outputs, integration_rule, domain_attributes, - derivative_ids); - } - - /// @overload - mfem::Vector residual(double, double, const std::vector& fields, int block_row = 0) const override - { - SLIC_ERROR_IF(block_row != 0, "Invalid block row and column requested in fieldJacobian for DfemResidual2"); - mfem::Vector resid(fields[0]->space().GetVSize()); - resid = 0.0; - for (auto& diff_op_with_field_ids : diff_ops_with_field_ids_) { - std::vector other_fields; - auto& field_ids = std::get<1>(diff_op_with_field_ids); - other_fields.reserve(field_ids.size() - 1); - for (size_t i = 1; i < field_ids.size(); ++i) { - other_fields.push_back(&fields[field_ids[i]]->gridFunction()); - } - std::get<0>(diff_op_with_field_ids).SetParameters(other_fields); - mfem::Vector term_resid(fields[0]->space().GetTrueVSize()); - std::get<0>(diff_op_with_field_ids).Mult(*fields[0], term_resid); - resid += term_resid; - } - return resid; - } - - /// @overload - std::unique_ptr jacobian(double, double, const std::vector&, - const std::vector&, int) const override - { - return std::make_unique(); - } - - /// @overload - void jvp(double, double, const std::vector& fields, const std::vector& v_fields, - const std::vector& jvp_reactions) const override - { - SLIC_ERROR_IF(v_fields.size() != fields.size(), - "Invalid number of field sensitivities relative to the number of fields"); - SLIC_ERROR_IF(jvp_reactions.size() != 1, "DfemResidual nonlinear systems only supports 1 output residual"); - - auto grad_op = diff_op_.GetDerivative(0, {fields[0]}, {mesh_nodes_}); - grad_op->AddMult(*v_fields[0], *jvp_reactions[0]); - } - - /// @overload - void vjp(double, double, const std::vector& fields, const std::vector& v_fields, - const std::vector& vjp_sensitivities) const override - { - SLIC_ERROR_IF(vjp_sensitivities.size() != fields.size(), - "Invalid number of field sensitivities relative to the number of fields"); - SLIC_ERROR_IF(v_fields.size() != 1, "FunctionalResidual nonlinear systems only supports 1 output residual"); - - auto grad_op = diff_op_.GetDerivative(0, {fields[0]}, {mesh_nodes_}); - grad_op->AddMultTranspose(*v_fields[0], *vjp_sensitivities[0]); - } - - protected: - /// @brief primary mesh - std::shared_ptr mesh_; - - mutable std::vector>> diff_ops_with_field_ids_; -}; - -} // namespace serac diff --git a/src/serac/physics/solid_dfem_residual.hpp b/src/serac/physics/solid_dfem_residual.hpp index 76d038c5c..af9d7e467 100644 --- a/src/serac/physics/solid_dfem_residual.hpp +++ b/src/serac/physics/solid_dfem_residual.hpp @@ -49,7 +49,7 @@ struct StressDivQFunction { * @tparam dim The spatial dimension of the mesh */ template -class SolidDfemResidual : public DfemResidual { +class SolidDfemResidual2 : public DfemResidual { public: /// @brief disp, velo, accel static constexpr int NUM_STATE_VARS = 4; @@ -59,9 +59,7 @@ class SolidDfemResidual : public DfemResidual { DISP, ///< displacement VELO, ///< velocity ACCEL, ///< acceleration - COORD, ///< coordinates - PARAM, ///< parameter field, if any - STATE ///< internal state field, if any + COORD ///< coordinates }; /** @@ -69,17 +67,16 @@ class SolidDfemResidual : public DfemResidual { * * @param physics_name A name for the physics module instance * @param mesh The serac Mesh - * @param diff_op A differentiable operator that computes the residual and jacobian */ - SolidDfemResidual(std::string physics_name, std::shared_ptr mesh, const mfem::ParFiniteElementSpace& test_space, - std::vector parameter_fe_spaces = {}) - : DfemResidual( - physics_name, mesh, - mfem::future::DifferentiableOperator(solutionFieldDescriptors(mesh->mfemParMesh(), test_space), - parameterFieldDescriptors(parameter_fe_spaces), mesh->mfemParMesh())), - solution_fields_(solutionFieldDescriptors(mesh->mfemParMesh(), test_space)), - parameter_fields_(parameterFieldDescriptors(parameter_fe_spaces)) + SolidDfemResidual2(std::string physics_name, std::shared_ptr mesh, + const mfem::ParFiniteElementSpace& test_space, + const std::vector& parameter_fe_spaces = {}) + : DfemResidual(physics_name, mesh), test_space_(test_space) { + parameter_fields_.reserve(parameter_fe_spaces.size()); + for (auto space : parameter_fe_spaces) { + parameter_fields_.emplace_back(NUM_STATE_VARS + parameter_fields_.size(), space); + } } /** @@ -106,87 +103,51 @@ class SolidDfemResidual : public DfemResidual { * @pre MaterialType must have a public method 'pkStress' which returns the first Piola-Kirchhoff stress * */ - template - void setMaterial(std::string /*body_name*/, const MaterialType& material, - const mfem::IntegrationRule& displacement_ir) + template + void setMaterial(DependsOn, const mfem::Array& domain_attributes, + const MaterialType& material, const mfem::IntegrationRule& displacement_ir) { - mfem::future::tuple outputs{mfem::future::Gradient{}}; - // TODO: find out the right attributes from the body name - mfem::Array solid_domain_attributes(DfemResidual::mesh_->mfemParMesh().attributes.Max()); - if constexpr (MaterialType::has_state && MaterialType::has_parameters) { - // build new diff_op_ with state and parameter fields - mfem::future::tuple inputs{// TODO: figure out how to pass in dt - mfem::future::Gradient{}, mfem::future::Gradient{}, - mfem::future::Gradient{}, mfem::future::Gradient{}, - mfem::future::Weight{}, mfem::future::Value{}, - mfem::future::Value{}}; - auto stress_div_qf = - StressDivQFunction{ - .material = material}; - DfemResidual::diff_op_.AddDomainIntegrator(stress_div_qf, inputs, outputs, displacement_ir, - solid_domain_attributes); - } else if constexpr (MaterialType::has_parameters) { - // build new diff_op_ with parameter fields - mfem::future::tuple inputs{// TODO: figure out how to pass in dt - mfem::future::Gradient{}, mfem::future::Gradient{}, - mfem::future::Gradient{}, mfem::future::Gradient{}, - mfem::future::Weight{}, mfem::future::Value{}}; - auto stress_div_qf = StressDivQFunction{.material = material}; - DfemResidual::diff_op_.AddDomainIntegrator(stress_div_qf, inputs, outputs, displacement_ir, - solid_domain_attributes); - } else if constexpr (MaterialType::has_state) { - // build new diff_op_ with state fields - mfem::future::tuple inputs{// TODO: figure out how to pass in dt - mfem::future::Gradient{}, mfem::future::Gradient{}, - mfem::future::Gradient{}, mfem::future::Gradient{}, - mfem::future::Weight{}, mfem::future::Value{}}; - auto stress_div_qf = StressDivQFunction{.material = material}; - DfemResidual::diff_op_.AddDomainIntegrator(stress_div_qf, inputs, outputs, displacement_ir, - solid_domain_attributes); - } else { - // existing diff_op_ should work fine - mfem::future::tuple inputs{// TODO: figure out how to pass in dt - mfem::future::Gradient{}, mfem::future::Gradient{}, - mfem::future::Gradient{}, mfem::future::Gradient{}, - mfem::future::Weight{}}; - auto stress_div_qf = StressDivQFunction{.material = material}; - DfemResidual::diff_op_.AddDomainIntegrator(stress_div_qf, inputs, outputs, displacement_ir, - solid_domain_attributes); - } + std::vector active_parameter_fields({parameter_fields_[active_parameters]...}); + // NOTE: this is used to size all of the arguments. for scalars, the qfunction argument can just be a double, but + // for tensors, it has to be a mfem::future::tensor + mfem::future::tuple qfunction_inputs{ + // TODO: figure out how to pass in dt + mfem::future::Gradient{}, mfem::future::Gradient{}, + mfem::future::Gradient{}, mfem::future::Gradient{}, + mfem::future::Weight{}, mfem::future::Value{}...}; + mfem::future::tuple qfunction_outputs{mfem::future::Gradient{}}; + // TODO: use serac::Domain instead of mfem attributes + // NOTE: MaterialParameters is used to set the qfunction parameters. TODO: potentially embed this into the material + // itself through e.g. a tuple? + auto stress_div_qfunction = StressDivQFunction{.material = material}; + DfemResidual::addBodyIntegral(stress_div_qfunction, solutionFieldDescriptors(), active_parameter_fields, + qfunction_inputs, qfunction_outputs, displacement_ir, domain_attributes, + std::integer_sequence{}); } - protected: - /** - * @brief Construct the field descriptors for the differentiable operator - * - * @param parameter_fe_spaces The parameter finite element spaces - * @return A vector of field descriptors for the differentiable operator - */ - static std::vector parameterFieldDescriptors( - std::vector parameter_fe_spaces) + /// @overload + template + void setMaterial(const mfem::Array& domain_attributes, const MaterialType& material, + const mfem::IntegrationRule& displacement_ir) { - std::vector field_descriptors; - field_descriptors.reserve(parameter_fe_spaces.size()); - for (const auto& space : parameter_fe_spaces) { - field_descriptors.emplace_back(NUM_STATE_VARS + field_descriptors.size(), space); - } - return field_descriptors; + setMaterial(DependsOn<>{}, domain_attributes, material, displacement_ir); } + + protected: /** * @brief Construct the field descriptors for the solution fields * */ - static std::vector solutionFieldDescriptors( - const mfem::ParMesh& mesh, const mfem::ParFiniteElementSpace& test_space) + std::vector solutionFieldDescriptors() { return { - {DISP, &test_space}, - {VELO, &test_space}, - {ACCEL, &test_space}, - {COORD, static_cast(mesh.GetNodalFESpace())}, + {DISP, &test_space_}, + {VELO, &test_space_}, + {ACCEL, &test_space_}, + {COORD, static_cast(DfemResidual::mesh_->mfemParMesh().GetNodalFESpace())}, }; } - std::vector solution_fields_; + const mfem::ParFiniteElementSpace& test_space_; std::vector parameter_fields_; }; diff --git a/src/serac/physics/solid_dfem_residual2.hpp b/src/serac/physics/solid_dfem_residual2.hpp deleted file mode 100644 index c5b1f471b..000000000 --- a/src/serac/physics/solid_dfem_residual2.hpp +++ /dev/null @@ -1,158 +0,0 @@ -// 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) - -/** - * @file solid_dfem_residual.hpp - * - * @brief Implements the residual interface for solid mechanics physics. - * Derives from dfem_residual. - */ - -#pragma once - -#include "serac/physics/dfem_residual2.hpp" - -namespace serac { - -template -struct StressDivQFunction { - SERAC_HOST_DEVICE inline auto operator()( - // mfem::real_t dt, // TODO: figure out how to pass this in - const mfem::future::tensor& du_dxi, - const mfem::future::tensor& dv_dxi, - const mfem::future::tensor&, - const mfem::future::tensor& dX_dxi, mfem::real_t weight, - Parameters... params) const - { - auto dxi_dX = mfem::future::inv(dX_dxi); - auto du_dX = mfem::future::dot(du_dxi, dxi_dX); - auto dv_dX = mfem::future::dot(dv_dxi, dxi_dX); - double dt = 1.0; // TODO: figure out how to pass this in to the qfunction - auto P = mfem::future::get<0>(material.pkStress(dt, du_dX, dv_dxi, params...)); - auto JxW = mfem::future::det(dX_dxi) * weight * mfem::future::transpose(dxi_dX); - return mfem::future::tuple{P * JxW}; - } - - Material material; ///< the material model to use for computing the stress -}; - -/** - * @brief The nonlinear residual class - * - * This uses Functional to compute the solid mechanics residuals and tangent - * stiffness matrices. - * - * @tparam order The order of the discretization of the displacement and velocity fields - * @tparam dim The spatial dimension of the mesh - */ -template -class SolidDfemResidual2 : public DfemResidual2 { - public: - /// @brief disp, velo, accel - static constexpr int NUM_STATE_VARS = 4; - - enum FieldIDs - { - DISP, ///< displacement - VELO, ///< velocity - ACCEL, ///< acceleration - COORD ///< coordinates - }; - - /** - * @brief Construct a new SolidResidual object - * - * @param physics_name A name for the physics module instance - * @param mesh The serac Mesh - */ - SolidDfemResidual2(std::string physics_name, std::shared_ptr mesh, - const mfem::ParFiniteElementSpace& test_space, - const std::vector& parameter_fe_spaces = {}) - : DfemResidual2(physics_name, mesh), test_space_(test_space) - { - parameter_fields_.reserve(parameter_fe_spaces.size()); - for (auto space : parameter_fe_spaces) { - parameter_fields_.emplace_back(NUM_STATE_VARS + parameter_fields_.size(), space); - } - } - - /** - * @brief Set the material stress response and mass properties for the physics module - * - * @tparam MaterialType The solid material type - * @tparam StateType the type that contains the internal variables for MaterialType - * @param body_name string name for a registered body Domain on the mesh - * @param material A material that provides a function to evaluate stress - * @pre material must be a object that can be called with the following arguments: - * 1. `MaterialType::State & state` an mutable reference to the internal variables for this quadrature point - * 2. `tensor du_dx` the displacement gradient at this quadrature point - * 3. `tuple{value, derivative}`, a tuple of values and derivatives for each parameter field - * specified in the `DependsOn<...>` argument. - * - * @note The actual types of these arguments passed will be `double`, `tensor` or tuples thereof - * when doing direct evaluation. When differentiating with respect to one of the inputs, its stored - * values will change to `dual` numbers rather than `double`. (e.g. `tensor` becomes `tensor, - * 3>`) - * - * @param qdata the buffer of material internal variables at each quadrature point - * - * @pre MaterialType must have a public method `density` which can take FiniteElementState parameter inputs - * @pre MaterialType must have a public method 'pkStress' which returns the first Piola-Kirchhoff stress - * - */ - template - void setMaterial(DependsOn, std::string /*body_name*/, const MaterialType& material, - const mfem::IntegrationRule& displacement_ir) - { - mfem::future::tuple qfunction_inputs{ - // TODO: figure out how to pass in dt - mfem::future::Gradient{}, mfem::future::Gradient{}, - mfem::future::Gradient{}, mfem::future::Gradient{}, - mfem::future::Weight{}, mfem::future::Value{}...}; - mfem::future::tuple qfunction_outputs{mfem::future::Gradient{}}; - // TODO: find out the right attributes from the body name - mfem::Array solid_domain_attributes(DfemResidual2::mesh_->mfemParMesh().attributes.Max()); - auto stress_div_qf = StressDivQFunction{.material = material}; - DfemResidual2::addBodyIntegral(stress_div_qf, solutionFieldDescriptors(), parameter_fields_, qfunction_inputs, - qfunction_outputs, displacement_ir, solid_domain_attributes, - std::integer_sequence{}); - } - - protected: - /** - * @brief Construct the field descriptors for the differentiable operator - * - * @param parameter_fe_spaces The parameter finite element spaces - * @return A vector of field descriptors for the differentiable operator - */ - static std::vector parameterFieldDescriptors( - std::vector parameter_fe_spaces) - { - std::vector field_descriptors; - field_descriptors.reserve(parameter_fe_spaces.size()); - for (const auto& space : parameter_fe_spaces) { - field_descriptors.emplace_back(NUM_STATE_VARS + field_descriptors.size(), space); - } - return field_descriptors; - } - /** - * @brief Construct the field descriptors for the solution fields - * - */ - std::vector solutionFieldDescriptors() - { - return { - {DISP, &test_space_}, - {VELO, &test_space_}, - {ACCEL, &test_space_}, - {COORD, static_cast(DfemResidual2::mesh_->mfemParMesh().GetNodalFESpace())}, - }; - } - const mfem::ParFiniteElementSpace& test_space_; - std::vector parameter_fields_; -}; - -} // namespace serac From 11a7b922138ed418cd20b84139da552a4caf41b6 Mon Sep 17 00:00:00 2001 From: "Eric B. Chin" Date: Thu, 26 Jun 2025 09:51:36 -0700 Subject: [PATCH 09/71] get it building for now --- src/serac/physics/dfem_residual.hpp | 211 ++++++++------- src/serac/physics/solid_dfem_residual.hpp | 78 +++--- src/serac/physics/tests/CMakeLists.txt | 1 + .../physics/tests/test_dfem_vs_functional.cpp | 251 +++++++++--------- 4 files changed, 285 insertions(+), 256 deletions(-) diff --git a/src/serac/physics/dfem_residual.hpp b/src/serac/physics/dfem_residual.hpp index 21d08c7cb..ab8793452 100644 --- a/src/serac/physics/dfem_residual.hpp +++ b/src/serac/physics/dfem_residual.hpp @@ -26,6 +26,7 @@ namespace serac { * stiffness matrices based on body and boundary integrals. * */ +template class DfemResidual : public Residual { public: using SpacesT = std::vector; ///< typedef @@ -37,7 +38,23 @@ class DfemResidual : public Residual { * @param mesh The serac mesh * @param diff_op A differentiable operator that computes the residual and jacobian */ - DfemResidual(std::string physics_name, std::shared_ptr mesh) : Residual(physics_name), mesh_(mesh) {} + DfemResidual(std::string physics_name, std::shared_ptr mesh, + const mfem::ParFiniteElementSpace& output_mfem_space, + const std::vector& input_mfem_spaces) + : Residual(physics_name), + mesh_(mesh), + differentiable_operators_(input_mfem_spaces.size() + 1), + input_mfem_spaces_(input_mfem_spaces) + { + std::vector parameter_spaces; + parameter_spaces.reserve(input_mfem_spaces.size()); + for (size_t i = 1; i < input_mfem_spaces.size(); ++i) { + parameter_spaces.push_back(input_mfem_spaces[i]); + } + parameter_spaces.push_back(&output_mfem_space); + differentiable_operators_[0] = std::make_unique( + makeFieldDescriptors({input_mfem_spaces[0]}), makeFieldDescriptors(parameter_spaces, 1), mesh->mfemParMesh()); + } /** * @brief Add a body integral contribution to the residual @@ -59,56 +76,30 @@ class DfemResidual : public Residual { * 3>`) * */ - template - void addBodyIntegral(QFunctionType qfunction, const std::vector& solution_fields, - const std::vector& param_fields, InputType qfunction_inputs, - OutputType qfunction_outputs, const mfem::IntegrationRule& integration_rule, - mfem::Array domain_attributes, std::integer_sequence derivative_ids) + template + void addBodyIntegral(DependsOn, mfem::Array domain_attributes, BodyIntegralType body_integral, + const mfem::IntegrationRule& integration_rule, std::integer_sequence derivative_ids) { - // get field IDs - std::vector field_ids; - field_ids.reserve(solution_fields.size() + param_fields.size()); - for (const auto& field : solution_fields) { - field_ids.push_back(field.id); - } - for (const auto& field : param_fields) { - field_ids.push_back(field.id); - } - // just keep the first input field as a solution field, make the rest parameter fields to simplify the Mult() call - std::vector other_fields; - other_fields.reserve(param_fields.size() + solution_fields.size() - 1); - for (size_t i = 1; i < solution_fields.size(); ++i) { - other_fields.push_back(solution_fields[i]); - } - for (const auto& field : param_fields) { - other_fields.push_back(field); - } - residual_terms_with_field_ids_.emplace_back( - std::make_tuple(mfem::future::DifferentiableOperator({solution_fields[0]}, other_fields, mesh_->mfemParMesh()), - std::move(field_ids))); - std::get<0>(residual_terms_with_field_ids_.back()) - .AddDomainIntegrator(qfunction, qfunction_inputs, qfunction_outputs, integration_rule, domain_attributes, - derivative_ids); + InputType integral_inputs; + OutputType integral_outputs; + // ParameterReducer reduced_integral(body_integral); + differentiable_operators_[0]->AddDomainIntegrator(body_integral, integral_inputs, integral_outputs, + integration_rule, domain_attributes, derivative_ids); } /// @overload mfem::Vector residual(double, double, const std::vector& fields, int block_row = 0) const override { SLIC_ERROR_ROOT_IF(block_row != 0, "Invalid block row and column requested in fieldJacobian for DfemResidual"); - mfem::Vector resid(fields[0]->space().GetVSize()); + mfem::Vector resid(fields[0]->space().GetTrueVSize()); resid = 0.0; - for (auto& residual_term_with_field_ids : residual_terms_with_field_ids_) { - std::vector other_fields; - auto& field_ids = std::get<1>(residual_term_with_field_ids); - other_fields.reserve(field_ids.size() - 1); - for (size_t i = 1; i < field_ids.size(); ++i) { - other_fields.push_back(&fields[field_ids[i]]->gridFunction()); - } - std::get<0>(residual_term_with_field_ids).SetParameters(other_fields); - mfem::Vector term_resid(fields[0]->space().GetTrueVSize()); - std::get<0>(residual_term_with_field_ids).Mult(*fields[0], term_resid); - resid += term_resid; + std::vector other_fields; + other_fields.reserve(fields.size() - 1); + for (size_t i = 1; i < fields.size(); ++i) { + other_fields.push_back(&fields[i]->gridFunction()); } + differentiable_operators_[0]->SetParameters(other_fields); + differentiable_operators_[0]->Mult(*fields[0], resid); return resid; } @@ -128,52 +119,52 @@ class DfemResidual : public Residual { "Invalid number of field sensitivities relative to the number of fields"); SLIC_ERROR_ROOT_IF(jvp_reactions.size() != 1, "DfemResidual nonlinear systems only supports 1 output residual"); - *jvp_reactions[0] = 0.0; - - for (size_t input_col = 0; input_col < fields.size(); ++input_col) { - if (v_fields[input_col] == nullptr) { - continue; - } - for (auto& residual_term_with_field_ids : residual_terms_with_field_ids_) { - const auto& field_ids = std::get<1>(residual_term_with_field_ids); - auto field_id_it = std::find(field_ids.begin(), field_ids.end(), input_col); - // this residual term does not depend on this input field - if (field_id_it == field_ids.end()) { - continue; - } - // if (field_ids[0] != input_col) { - // continue; - // } - // auto& residual_term_ = std::get<0>(residual_term_with_field_ids); - // std::vector other_fields; - // other_fields.reserve(field_ids.size() - 1); - // for (size_t i = 1; i < field_ids.size(); ++i) { - // other_fields.push_back(&fields[field_ids[i]]->gridFunction()); - // } - // residual_term_.SetParameters(other_fields); - // mfem::Vector K(fields[0]->space().GetTrueVSize()); - // residual_term_.AddMult(*v_fields[input_col], K); - // (*jvp_reactions[0]) += K; - // } - // // find the residual term that corresponds to this input field - // auto it = std::find_if(residual_terms_with_field_ids_.begin(), residual_terms_with_field_ids_.end(), - // [&](const auto& term_and_ids) { - // const auto& field_ids = std::get<1>(term_and_ids); - // return field_ids[0] == input_col; - // }); - // if (it == residual_terms_with_field_ids_.end()) { - // continue; - // } - // const auto& residual_term_ = std::get<0>(*it); - // std::vector other_fields; - // const auto& field_ids = std::get<1>(*it); - // other_fields.reserve(field_ids.size() - 1); - // for (size_t i = 1; i < field_ids.size(); ++i) { - // other_fields.push_back(&fields[field_ids[i]]->gridFunction()); - // } - // residual_term_.SetParameters(other_fields); - } - } + // *jvp_reactions[0] = 0.0; + + // for (size_t input_col = 0; input_col < fields.size(); ++input_col) { + // if (v_fields[input_col] == nullptr) { + // continue; + // } + // for (auto& residual_term_with_field_ids : residual_terms_with_field_ids_) { + // const auto& field_ids = std::get<1>(residual_term_with_field_ids); + // auto field_id_it = std::find(field_ids.begin(), field_ids.end(), input_col); + // // this residual term does not depend on this input field + // if (field_id_it == field_ids.end()) { + // continue; + // } + // if (field_ids[0] != input_col) { + // continue; + // } + // auto& residual_term_ = std::get<0>(residual_term_with_field_ids); + // std::vector other_fields; + // other_fields.reserve(field_ids.size() - 1); + // for (size_t i = 1; i < field_ids.size(); ++i) { + // other_fields.push_back(&fields[field_ids[i]]->gridFunction()); + // } + // residual_term_.SetParameters(other_fields); + // mfem::Vector K(fields[0]->space().GetTrueVSize()); + // residual_term_.AddMult(*v_fields[input_col], K); + // (*jvp_reactions[0]) += K; + // } + // // find the residual term that corresponds to this input field + // auto it = std::find_if(residual_terms_with_field_ids_.begin(), residual_terms_with_field_ids_.end(), + // [&](const auto& term_and_ids) { + // const auto& field_ids = std::get<1>(term_and_ids); + // return field_ids[0] == input_col; + // }); + // if (it == residual_terms_with_field_ids_.end()) { + // continue; + // } + // const auto& residual_term_ = std::get<0>(*it); + // std::vector other_fields; + // const auto& field_ids = std::get<1>(*it); + // other_fields.reserve(field_ids.size() - 1); + // for (size_t i = 1; i < field_ids.size(); ++i) { + // other_fields.push_back(&fields[field_ids[i]]->gridFunction()); + // } + // residual_term_.SetParameters(other_fields); + // } + // } } /// @overload @@ -184,16 +175,52 @@ class DfemResidual : public Residual { "Invalid number of field sensitivities relative to the number of fields"); SLIC_ERROR_ROOT_IF(v_fields.size() != 1, "FunctionalResidual nonlinear systems only supports 1 output residual"); - auto grad_op = residual_term_.GetDerivative(0, {fields[0]}, {mesh_nodes_}); - grad_op->AddMultTranspose(*v_fields[0], *vjp_sensitivities[0]); + // auto grad_op = residual_term_.GetDerivative(0, {fields[0]}, {mesh_nodes_}); + // grad_op->AddMultTranspose(*v_fields[0], *vjp_sensitivities[0]); } protected: + static std::vector makeFieldDescriptors( + const std::vector& spaces, size_t offset = 0) + { + std::vector field_descriptors; + field_descriptors.reserve(spaces.size()); + for (size_t i = 0; i < spaces.size(); ++i) { + field_descriptors.emplace_back(i + offset, spaces[i]); + } + return field_descriptors; + } + /// @brief primary mesh std::shared_ptr mesh_; + std::vector> differentiable_operators_; + std::vector input_mfem_spaces_; + + private: + template + struct ParameterReducer { + ParameterReducer(Integrand integrand) : integrand_(integrand) {} + + template + SERAC_HOST_DEVICE inline auto operator()(Args&&... args) const + { + auto arg_tuple = mfem::future::make_tuple(std::forward(args)...); + return integrand_(std::get(arg_tuple)...); + } + + SERAC_HOST_DEVICE inline auto operator()(const mfem::future::tensor& du_dxi, + const mfem::future::tensor& dv_dxi, + const mfem::future::tensor& da_dxi, + const mfem::future::tensor& dX_dxi, + mfem::real_t weight, double K) const + { + auto arg_tuple = mfem::future::make_tuple(du_dxi, dv_dxi, da_dxi, dX_dxi, weight, K); + return integrand_(std::get(arg_tuple)...); + } - mutable std::vector>> - residual_terms_with_field_ids_; + private: + Integrand integrand_; + }; }; } // namespace serac diff --git a/src/serac/physics/solid_dfem_residual.hpp b/src/serac/physics/solid_dfem_residual.hpp index af9d7e467..3613b229c 100644 --- a/src/serac/physics/solid_dfem_residual.hpp +++ b/src/serac/physics/solid_dfem_residual.hpp @@ -31,7 +31,7 @@ struct StressDivQFunction { auto du_dX = mfem::future::dot(du_dxi, dxi_dX); auto dv_dX = mfem::future::dot(dv_dxi, dxi_dX); double dt = 1.0; // TODO: figure out how to pass this in to the qfunction - auto P = mfem::future::get<0>(material.pkStress(dt, du_dX, dv_dxi, params...)); + auto P = mfem::future::get<0>(material.pkStress(dt, du_dX, dv_dX, params...)); auto JxW = mfem::future::det(dX_dxi) * weight * mfem::future::transpose(dxi_dX); return mfem::future::tuple{P * JxW}; } @@ -48,8 +48,12 @@ struct StressDivQFunction { * @tparam order The order of the discretization of the displacement and velocity fields * @tparam dim The spatial dimension of the mesh */ -template -class SolidDfemResidual2 : public DfemResidual { +template +class SolidDfemResidual + : public DfemResidual< + mfem::future::tuple, mfem::future::Gradient<1>, mfem::future::Gradient<2>, + mfem::future::Gradient<3>, mfem::future::Weight, ParameterTypes...>, + mfem::future::tuple>> { public: /// @brief disp, velo, accel static constexpr int NUM_STATE_VARS = 4; @@ -62,21 +66,21 @@ class SolidDfemResidual2 : public DfemResidual { COORD ///< coordinates }; + using BaseClassT = DfemResidual< + mfem::future::tuple, mfem::future::Gradient, mfem::future::Gradient, + mfem::future::Gradient, mfem::future::Weight, ParameterTypes...>, + mfem::future::tuple>>; + /** * @brief Construct a new SolidResidual object * * @param physics_name A name for the physics module instance * @param mesh The serac Mesh */ - SolidDfemResidual2(std::string physics_name, std::shared_ptr mesh, - const mfem::ParFiniteElementSpace& test_space, - const std::vector& parameter_fe_spaces = {}) - : DfemResidual(physics_name, mesh), test_space_(test_space) + SolidDfemResidual(std::string physics_name, std::shared_ptr mesh, const mfem::ParFiniteElementSpace& test_space, + const std::vector& parameter_fe_spaces = {}) + : BaseClassT(physics_name, mesh, test_space, makeInputSpaces(test_space, mesh, parameter_fe_spaces)) { - parameter_fields_.reserve(parameter_fe_spaces.size()); - for (auto space : parameter_fe_spaces) { - parameter_fields_.emplace_back(NUM_STATE_VARS + parameter_fields_.size(), space); - } } /** @@ -103,26 +107,13 @@ class SolidDfemResidual2 : public DfemResidual { * @pre MaterialType must have a public method 'pkStress' which returns the first Piola-Kirchhoff stress * */ - template - void setMaterial(DependsOn, const mfem::Array& domain_attributes, - const MaterialType& material, const mfem::IntegrationRule& displacement_ir) + template + void setMaterial(DependsOn, const mfem::Array& domain_attributes, const MaterialType& material, + const mfem::IntegrationRule& displacement_ir) { - std::vector active_parameter_fields({parameter_fields_[active_parameters]...}); - // NOTE: this is used to size all of the arguments. for scalars, the qfunction argument can just be a double, but - // for tensors, it has to be a mfem::future::tensor - mfem::future::tuple qfunction_inputs{ - // TODO: figure out how to pass in dt - mfem::future::Gradient{}, mfem::future::Gradient{}, - mfem::future::Gradient{}, mfem::future::Gradient{}, - mfem::future::Weight{}, mfem::future::Value{}...}; - mfem::future::tuple qfunction_outputs{mfem::future::Gradient{}}; - // TODO: use serac::Domain instead of mfem attributes - // NOTE: MaterialParameters is used to set the qfunction parameters. TODO: potentially embed this into the material - // itself through e.g. a tuple? - auto stress_div_qfunction = StressDivQFunction{.material = material}; - DfemResidual::addBodyIntegral(stress_div_qfunction, solutionFieldDescriptors(), active_parameter_fields, - qfunction_inputs, qfunction_outputs, displacement_ir, domain_attributes, - std::integer_sequence{}); + auto stress_div_integral = StressDivQFunction{.material = material}; + BaseClassT::addBodyIntegral(DependsOn{}, domain_attributes, stress_div_integral, displacement_ir, + std::integer_sequence{}); } /// @overload @@ -134,21 +125,22 @@ class SolidDfemResidual2 : public DfemResidual { } protected: - /** - * @brief Construct the field descriptors for the solution fields - * - */ - std::vector solutionFieldDescriptors() + private: + std::vector makeInputSpaces( + const mfem::ParFiniteElementSpace& test_space, const std::shared_ptr& mesh, + const std::vector& parameter_fe_spaces) { - return { - {DISP, &test_space_}, - {VELO, &test_space_}, - {ACCEL, &test_space_}, - {COORD, static_cast(DfemResidual::mesh_->mfemParMesh().GetNodalFESpace())}, - }; + std::vector input_spaces; + input_spaces.reserve(4 + parameter_fe_spaces.size()); + for (int i = 0; i < 3; ++i) { + input_spaces.push_back(&test_space); + } + input_spaces.push_back(static_cast(mesh->mfemParMesh().GetNodalFESpace())); + for (auto space : parameter_fe_spaces) { + input_spaces.push_back(space); + } + return input_spaces; } - const mfem::ParFiniteElementSpace& test_space_; - std::vector parameter_fields_; }; } // namespace serac diff --git a/src/serac/physics/tests/CMakeLists.txt b/src/serac/physics/tests/CMakeLists.txt index afb44edce..512947114 100644 --- a/src/serac/physics/tests/CMakeLists.txt +++ b/src/serac/physics/tests/CMakeLists.txt @@ -29,6 +29,7 @@ set(physics_serial_test_sources finite_element_vector_set_over_domain.cpp solid_multi_material.cpp test_solid_residual.cpp + test_dfem_vs_functional.cpp test_heat_residual.cpp test_functional_residual.cpp thermomech_finite_diff.cpp diff --git a/src/serac/physics/tests/test_dfem_vs_functional.cpp b/src/serac/physics/tests/test_dfem_vs_functional.cpp index b9b74ae9c..10aed75b3 100644 --- a/src/serac/physics/tests/test_dfem_vs_functional.cpp +++ b/src/serac/physics/tests/test_dfem_vs_functional.cpp @@ -11,64 +11,63 @@ #include "serac/physics/materials/solid_material.hpp" #include "serac/physics/mesh.hpp" #include "serac/physics/state/state_manager.hpp" + +#include "serac/physics/tests/physics_test_utils.hpp" #include "serac/physics/solid_residual.hpp" +#include "serac/physics/solid_dfem_residual.hpp" auto element_shape = mfem::Element::QUADRILATERAL; -template -auto getPointers(std::vector& values) -{ - std::vector pointers; - for (auto& t : values) { - pointers.push_back(&t); - } - return pointers; -} +namespace serac { -template -auto getPointers(std::vector& states, std::vector& params) -{ - assert(params.size() > 0); - std::vector pointers; - for (auto& t : states) { - pointers.push_back(&t); - } - for (auto& t : params) { - pointers.push_back(&t); - } - return pointers; -} +template