Compare commits

..
Author SHA1 Message Date
camierjs 9bef1cc28f Cleanup pass 2026-03-01 15:30:39 -08:00
camierjs 6b45d9a69b Merge branch 'master' into checkpointing_examples_gpu 2026-03-01 14:38:23 -08:00
camierjs db32ae9807 Cleanup 2026-03-01 14:34:24 -08:00
camierjs 545a257ab5 LE style, compile fix 2026-02-10 12:43:39 -08:00
blaz 4b4aff736d adjoint struct 2026-02-09 00:23:04 -08:00
blaz 727339e0ae check for objective gradients 2026-02-08 22:07:05 -08:00
blaz d8d4d5b874 Merge branch 'dfem-transpose' into checkpointing_examples 2026-02-08 20:24:10 -08:00
blaz 5a7a017a4e objective integradn 2026-02-08 20:21:33 -08:00
blaz 224aff5ac9 forward part of the adjoint computations 2026-02-07 21:09:30 -08:00
blaz e42762365f adjoint initial version 2026-02-07 20:09:35 -08:00
blaz 765fb4a608 added spectral matrix preconditioner 2026-02-07 19:38:59 -08:00
John Camier 591156be2f Merge branch 'master' into dfem-transpose 2026-01-30 18:42:29 -08:00
blaz 8e2478f964 Checkpointing examples 2026-01-25 19:45:26 -08:00
blaz 06c485bb3e updates 2026-01-21 11:03:43 -08:00
bslazarov 1cd27605f0 memory transfer to Device for q-functions 2026-01-20 14:43:06 -08:00
blaz 547fd02911 anisotrpic elasticity structure 2026-01-20 09:57:29 -08:00
blaz aa92498033 explicit linear elasticity 2026-01-20 09:51:15 -08:00
blaz 024576620a initial REVOLVE 2026-01-16 11:29:25 -08:00
blaz 5b670bda3f file and segement storage 2026-01-15 22:59:32 -08:00
blaz 6b0e135bdb rename files 2026-01-15 22:18:33 -08:00
blaz 204b9df2d6 checkpinting into file - multiple files - one file per snapshot 2026-01-15 22:16:31 -08:00
blaz 8b49379f87 fixed slot in-memory or file storage for checkpointing 2026-01-15 22:15:09 -08:00
blaz d32a00fb2f CMake list 2026-01-15 19:57:10 -08:00
blaz 59c9ca70e4 initial version 2026-01-15 19:49:02 -08:00
camierjs 9646ce4d1a Fix dFEM linearized sections leaks 2026-01-12 17:10:40 -08:00
Julian Andrej f46a81bf3a clearer comment 2026-01-12 13:09:39 -08:00
Julian Andrej 460e5ba483 correct sparsematrix test 2026-01-12 13:09:04 -08:00
Julian Andrej 34593dccac properly initialize matrix pointer 2026-01-09 16:37:16 -08:00
Julian Andrej 248ab78e77 missing dual types 2026-01-09 16:04:01 -08:00
Julian Andrej 6b9bdb11df fix hypreparmatrix assemble with multiple integrators 2026-01-09 13:18:37 -08:00
Julian Andrej 505d2e5123 fix multiple integrators 2026-01-09 12:59:25 -08:00
Julian Andrej 47cff96e2e documentation 2026-01-08 10:57:34 -08:00
Julian Andrej 010f455382 documentation 2026-01-06 08:51:46 -08:00
Julian Andrej acebc322b8 cleanup 2026-01-06 08:43:02 -08:00
Julian Andrej 25ec78df9a make tests more strict 2026-01-05 16:21:58 -08:00
Julian Andrej 8c78387136 remove periodic tests for transposes 2025-12-30 13:40:59 -08:00
Julian Andrej 8684d9d5d1 working 2025-12-30 12:41:29 -08:00
Julian Andrej 2375953e95 debugging state 2025-12-29 13:35:24 -08:00
Julian Andrej ecbb196551 bug 2025-12-23 08:27:01 -08:00
Julian Andrej 90ca758a98 again 2025-12-23 07:40:55 -08:00
Julian Andrej 93d674a7a2 argument order 2025-12-23 07:31:27 -08:00
Julian Andrej 03ec3f739b correct size in parallel for ldof vectors 2025-12-23 07:27:36 -08:00
Julian Andrej a3eb9c8e6b doxygen 2025-12-22 16:30:39 -08:00
Julian Andrej 30f3e123cd more warnings 2025-12-22 16:27:42 -08:00
Julian Andrej 4d4c2325f9 warnings 2025-12-22 16:25:39 -08:00
Julian Andrej db84abedfc make dual numbers work 2025-12-22 16:21:48 -08:00
Julian Andrej 2609e3af6e transpose action using cache 2025-12-22 09:05:50 -08:00
129 changed files with 13603 additions and 2502 deletions
+1 -1
View File
@@ -25,7 +25,7 @@ runs:
steps:
- uses: ./.github/actions/sanitize/config
- uses: actions/cache@v5
- uses: actions/cache@v4
if: ${{env.DEBUG == 'true'}}
id: debug
with:
+1 -1
View File
@@ -36,7 +36,7 @@ runs:
steps:
- uses: ./.github/actions/sanitize/config
- uses: actions/cache@v5
- uses: actions/cache@v4
if: ${{env.DEBUG == 'true' && inputs.cache-skip != 'true'}}
id: debug
with:
+5 -5
View File
@@ -23,7 +23,7 @@ inputs:
runs:
using: 'composite'
steps:
- uses: actions/cache/restore@v5 # Cache for LLVM libcxx
- uses: actions/cache/restore@v4 # Cache for LLVM libcxx
with:
path: ${{env.LLVM_DIR}}
fail-on-cache-miss: true
@@ -32,14 +32,14 @@ runs:
- uses: ./.github/actions/sanitize/mpi
if: ${{inputs.par == 'true'}}
- uses: actions/cache/restore@v5 # Cache for Hypre
- uses: actions/cache/restore@v4 # Cache for Hypre
if: ${{inputs.par == 'true'}}
with:
path: ${{env.HYPRE_DIR}}
fail-on-cache-miss: true
key: ${{runner.os}}-ompi-build-${{env.HYPRE_DIR}}-int32-fp64-v2.5
- uses: actions/cache/restore@v5 # Cache for Metis
- uses: actions/cache/restore@v4 # Cache for Metis
if: ${{inputs.par == 'true'}}
with:
path: ${{env.METIS_DIR}}
@@ -51,13 +51,13 @@ runs:
run: ln -s -f ${{env.HYPRE_DIR}} hypre && ln -s -f ${{env.METIS_DIR}} metis-4.0
shell: bash
- uses: actions/cache/restore@v5 # Cache for LSAN suppression file
- uses: actions/cache/restore@v4 # Cache for LSAN suppression file
with:
path: ${{env.LSAN_DIR}}
fail-on-cache-miss: true
key: build-lsan-suppression-file
- uses: actions/checkout@v6 # Checkout the repository
- uses: actions/checkout@v4 # Checkout the repository
with:
path: mfem
# ref: ${{env.BRANCH}}
+1 -1
View File
@@ -43,7 +43,7 @@ jobs:
remove-docker-images: 'true'
- name: Checkout
uses: actions/checkout@v6
uses: actions/checkout@v4
# It's easier to reference named variables than indexes of the matrix
- name: Set Environment
+5 -6
View File
@@ -153,7 +153,7 @@ jobs:
# /home/runner/work/mfem/mfem/mfem
# Note: Done now to access "install-hypre" and "install-metis" actions.
- name: checkout mfem
uses: actions/checkout@v6
uses: actions/checkout@v4
with:
path: ${{ env.MFEM_TOP_DIR }}
# Fetch the complete history for codecov to access commits ID
@@ -225,7 +225,7 @@ jobs:
- name: cache hypre
id: hypre-cache
if: matrix.mpi == 'par'
uses: actions/cache@v5
uses: actions/cache@v4
with:
path: ${{ env.HYPRE_TOP_DIR }}
key: ${{ runner.os }}-ompi-build-${{ env.HYPRE_TOP_DIR }}-${{ matrix.hypre-target }}-${{ matrix.precision }}-v2.5
@@ -255,7 +255,7 @@ jobs:
- name: cache metis
id: metis-cache
if: matrix.mpi == 'par' && matrix.os != 'windows-latest'
uses: actions/cache@v5
uses: actions/cache@v4
with:
path: ${{ env.METIS_TOP_DIR }}
key: ${{ runner.os }}-build-${{ env.METIS_TOP_DIR }}-v2.5
@@ -270,7 +270,7 @@ jobs:
- name: cache vcpkg (Windows)
id: vcpkg-cache
if: matrix.os == 'windows-latest'
uses: actions/cache@v5
uses: actions/cache@v4
with:
path: vcpkg_cache
key: ${{ runner.os }}-${{ matrix.mpi }}-vcpkg-v1
@@ -295,8 +295,7 @@ jobs:
export HOMEBREW_NO_INSTALL_CLEANUP=1
brew update
brew install enzyme
ENZYME_LLVM=$(brew info enzyme | sed -n 's/^Required.*:.*\(llvm[^ ]*\).*/\1/p')
echo "ENZYME_LLVM=$ENZYME_LLVM"
ENZYME_LLVM=$(brew info enzyme | sed -n 's/^Required:.*\(llvm[^ ]*\).*/\1/p')
LLVM_PREFIX=$(brew --prefix $ENZYME_LLVM)
echo "LLVM_PREFIX=$LLVM_PREFIX" >> $GITHUB_ENV
echo "OMPI_CC=$LLVM_PREFIX/bin/clang" >> $GITHUB_ENV
+4 -4
View File
@@ -40,11 +40,11 @@ jobs:
steps:
- name: Checkout repository
uses: actions/checkout@v6
uses: actions/checkout@v4
# Initializes the CodeQL tools for scanning.
- name: Initialize CodeQL
uses: github/codeql-action/init@v4
uses: github/codeql-action/init@v2
with:
languages: ${{ matrix.language }}
# If you wish to specify custom queries, you can do so here or in a config file.
@@ -57,7 +57,7 @@ jobs:
# Autobuild attempts to build any compiled languages (C/C++, C#, or Java).
# If this step fails, then you should remove it and run the build manually (see below)
- name: Autobuild
uses: github/codeql-action/autobuild@v4
uses: github/codeql-action/autobuild@v2
# ️ Command-line programs to run using the OS shell.
# 📚 See https://docs.github.com/en/actions/using-workflows/workflow-syntax-for-github-actions#jobsjob_idstepsrun
@@ -70,4 +70,4 @@ jobs:
# ./location_of_script_within_repo/buildscript.sh
- name: Perform CodeQL Analysis
uses: github/codeql-action/analyze@v4
uses: github/codeql-action/analyze@v2
+3 -3
View File
@@ -39,7 +39,7 @@ jobs:
steps:
- name: checkout MFEM
uses: actions/checkout@v6
uses: actions/checkout@v4
with:
path: mfem
@@ -50,7 +50,7 @@ jobs:
- name: Cache Hypre Install
id: hypre-cache
uses: actions/cache@v5
uses: actions/cache@v4
with:
path: ${{ env.HYPRE_TOP_DIR }}
key: ${{ runner.os }}-ompi-build-${{ env.HYPRE_TOP_DIR }}-v2.5
@@ -65,7 +65,7 @@ jobs:
- name: Cache Metis Install
id: metis-cache
uses: actions/cache@v5
uses: actions/cache@v4
with:
path: ${{ env.METIS_TOP_DIR }}
key: ${{ runner.os }}-build-${{ env.METIS_TOP_DIR }}-v2.5
+4 -4
View File
@@ -38,7 +38,7 @@ jobs:
github.event.pull_request.head.repo.full_name != github.repository)
steps:
- name: checkout mfem
uses: actions/checkout@v6
uses: actions/checkout@v4
- name: copyright check
id: copyright
@@ -93,7 +93,7 @@ jobs:
github.event.pull_request.head.repo.full_name != github.repository)
steps:
- name: checkout mfem
uses: actions/checkout@v6
uses: actions/checkout@v4
- name: get astyle
run: |
@@ -110,7 +110,7 @@ jobs:
github.event.pull_request.head.repo.full_name != github.repository)
steps:
- name: checkout mfem
uses: actions/checkout@v6
uses: actions/checkout@v4
- name: get doxygen and graphviz
run: |
@@ -135,7 +135,7 @@ jobs:
runs-on: ubuntu-latest
steps:
- name: checkout mfem
uses: actions/checkout@v6
uses: actions/checkout@v4
with:
fetch-depth: 0
+2 -2
View File
@@ -17,11 +17,11 @@ jobs:
runs-on: ubuntu-latest
name: 2.19.0
steps:
- uses: actions/checkout@v6
- uses: actions/checkout@v4
- uses: ./.github/actions/sanitize/config
- name: Cache
id: cache
uses: actions/cache@v5
uses: actions/cache@v4
with:
path: ${{env.HYPRE_DIR}}
key: ${{runner.os}}-ompi-build-${{env.HYPRE_DIR}}-int32-fp64-v2.5
+2 -2
View File
@@ -27,13 +27,13 @@ jobs:
llvm_use_sanitizer: "Undefined"
name: ${{matrix.sanitizer}}
steps:
- uses: actions/checkout@v6
- uses: actions/checkout@v4
- uses: ./.github/actions/sanitize/config
with:
NO_FLAGS: true
- name: Cache
id: cache
uses: actions/cache@v5
uses: actions/cache@v4
with:
path: ${{env.LLVM_DIR}}
key: build-libcxx-${{env.LLVM_VER}}-${{matrix.sanitizer}}
+2 -2
View File
@@ -17,11 +17,11 @@ jobs:
runs-on: ubuntu-latest
name: lsan.supp
steps:
- uses: actions/checkout@v6
- uses: actions/checkout@v4
- uses: ./.github/actions/sanitize/config
- name: Cache
id: cache
uses: actions/cache@v5
uses: actions/cache@v4
with:
path: ${{env.LSAN_DIR}}
key: build-lsan-suppression-file
+2 -2
View File
@@ -17,11 +17,11 @@ jobs:
runs-on: ubuntu-latest
name: 4.0.3
steps:
- uses: actions/checkout@v6
- uses: actions/checkout@v4
- uses: ./.github/actions/sanitize/config
- name: Cache
id: cache
uses: actions/cache@v5
uses: actions/cache@v4
with:
path: ${{env.METIS_DIR}}
key: ${{runner.os}}-build-${{env.METIS_DIR}}-v2.5
+7 -7
View File
@@ -28,7 +28,7 @@ jobs:
build:
runs-on: ubuntu-latest
steps:
- uses: actions/checkout@v6
- uses: actions/checkout@v4
- uses: ./.github/actions/sanitize/mfem
with:
par: ${{inputs.par}}
@@ -40,7 +40,7 @@ jobs:
env:
ex: ${{inputs.par && 'ex1p' || 'ex1'}}
steps:
- uses: actions/checkout@v6
- uses: actions/checkout@v4
- uses: ./.github/actions/sanitize/restore
id: restore
with:
@@ -58,7 +58,7 @@ jobs:
env:
exclude: ${{inputs.par && '-E "_ser"' || ''}}
steps:
- uses: actions/checkout@v6
- uses: actions/checkout@v4
- uses: ./.github/actions/sanitize/restore
id: restore
with:
@@ -82,7 +82,7 @@ jobs:
env:
exclude: ${{inputs.par && '-E "_ser"' || ''}}
steps:
- uses: actions/checkout@v6
- uses: actions/checkout@v4
- uses: ./.github/actions/sanitize/restore
id: restore
with:
@@ -107,7 +107,7 @@ jobs:
run: ${{inputs.par && '-R "_cpu_np"' || ''}}
exclude: ${{inputs.par && '"unit_tests|debug"' || '"^unit_tests$|debug"'}}
steps:
- uses: actions/checkout@v6
- uses: actions/checkout@v4
- uses: ./.github/actions/sanitize/restore
id: restore
with:
@@ -131,7 +131,7 @@ jobs:
env:
unit_tests: ${{inputs.par && 'punit_tests' || 'unit_tests'}}
steps:
- uses: actions/checkout@v6
- uses: actions/checkout@v4
- uses: ./.github/actions/sanitize/restore
id: restore
with:
@@ -165,7 +165,7 @@ jobs:
unit_tests: ${{inputs.par && 'punit_tests' || 'unit_tests'}}
np: ${{inputs.par && '_np=2' || ''}}
steps:
- uses: actions/checkout@v6
- uses: actions/checkout@v4
- uses: ./.github/actions/sanitize/restore
id: restore
with:
-20
View File
@@ -8,21 +8,6 @@
https://mfem.org
Version 4.10 (development)
==========================
Discretization improvements
---------------------------
- Replaced legacy simplex quadrature rules with symmetric positive-weight
rules for triangles (orders 0-25) and tetrahedra (orders 0-20). These
rules guarantee all-positive weights and interior quadrature points,
improving numerical stability. Higher orders fall back to Grundmann-Moller.
Triangle rules: Witherden & Vincent, Comput. Math. Appl. 69(10):1232-1241,
2015.
Tet rules (d=1-13): Witherden & Vincent (ibid).
Tet rules (d=14-20): Chuluunbaatar et al., Comput. Math. Appl. 124:89-97,
2022.
Version 4.9.1 (development)
===========================
@@ -45,11 +30,6 @@ New and updated examples and miniapps
- Electromagnetics/lorentz miniapp has been updated to leverage the ParticleSet
capability.
Miscellaneous
-------------
- Removes the SidreDataCollection class from MFEM in favor of the
MFEMSidreDataCollection class in the Axom library (https://github.com/llnl/axom).
Version 4.9, released on Dec 11, 2025
=====================================
+10 -7
View File
@@ -75,10 +75,12 @@ set(XSDK_ENABLE_Fortran OFF)
# Check if we need to enable C or Fortran.
if (MFEM_USE_CONDUIT OR
MFEM_USE_SIDRE OR
MFEM_USE_PETSC)
# This seems to be needed by:
# * find_package(BLAS REQUIRED) and
# * find_package(HDF5 REQUIRED) and
# * find_package(HDF5 REQUIRED) needed, in turn, by:
# - find_package(AXOM REQUIRED)
# * find_package(PETSc REQUIRED)
set(XSDK_ENABLE_C ON)
endif()
@@ -476,6 +478,11 @@ if (MFEM_USE_FMS)
find_package(FMS REQUIRED fms)
endif()
# Axom/Sidre
if (MFEM_USE_SIDRE)
find_package(Axom REQUIRED Axom)
endif()
# PUMI
if (MFEM_USE_PUMI)
# If PUMI_DIR was specified, only link to that directory,
@@ -622,7 +629,7 @@ find_package(Threads REQUIRED)
# integers, the METIS header (with 32-bit indices, as used by mfem) needs to
# be before SuiteSparse.
set(MFEM_TPLS OPENMP HYPRE LAPACK BLAS SuperLUDist STRUMPACK METIS SuiteSparse
SUNDIALS PETSC SLEPC MUMPS FMS CONDUIT Ginkgo GNUTLS GSLIB HDF5
SUNDIALS PETSC SLEPC MUMPS AXOM FMS CONDUIT Ginkgo GNUTLS GSLIB HDF5
NETCDF MPFR PUMI HIOP POSIXCLOCKS MFEMBacktrace ZLIB OCCA CEED RAJA UMPIRE
ADIOS2 MKL_CPARDISO MKL_PARDISO AMGX MAGMA CUSPARSE CUBLAS CALIPER CODIPACK
BENCHMARK PARELAG TRIBOL MPI_CXX HIP HIPBLAS HIPSPARSE MOONOLITH BLITZ
@@ -645,8 +652,6 @@ foreach(TPL IN LISTS MFEM_TPLS)
endif()
endforeach(TPL)
# reverse to remove the first instance of entries in TPL_LIBRARIES
# so later duplicates are kept (for dependency ordering)
list(REVERSE TPL_LIBRARIES)
list(REMOVE_DUPLICATES TPL_LIBRARIES)
list(REVERSE TPL_LIBRARIES)
@@ -1010,7 +1015,5 @@ install(DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/data
# Create 'config.mk' from 'config.mk.in' for the build and install locations and
# define install rules for 'config.mk' and 'test.mk'
#-------------------------------------------------------------------------------
if (MFEM_USE_CUDA OR MFEM_USE_HIP)
option(MFEM_EXPORT_GPU_CONFIG "Export config.mk for GPU-enabled downstream packages" ON)
endif()
mfem_export_mk_files()
+17
View File
@@ -452,6 +452,13 @@ MFEM_USE_MPFR = YES/NO
quadrature rules. When enabled, this option uses the MPFR_* library options,
see below.
MFEM_USE_SIDRE = YES/NO
Sidre is a component of LLNL's axom project, https://github.com/LLNL/axom,
that provides an HDF5-based file format for visualization or restart
capability following the Conduit (https://github.com/LLNL/conduit) mesh
blueprint specification. When enabled, this option requires installation of
HDF5 (see also MFEM_USE_NETCDF), Conduit and LLNL's axom project.
MFEM_USE_SIMD = YES/NO
Enables the high performance templated classes to use architecture dependent
SIMD intrinsics instead of the generic implementation of class AutoSIMD in
@@ -771,6 +778,14 @@ The specific libraries and their options are:
Options: SLEPC_OPT, SLEPC_LIB.
Versions: SLEPc >= 3.8.0.
- Sidre (optional), part of LLNL's axom project, used when MFEM_USE_SIDRE = YES.
Starting with MFEM v4.1, Axom version 0.3.1 or later is required.
URL: https://github.com/LLNL/axom
https://github.com/LLNL/conduit (Conduit)
https://support.hdfgroup.org/HDF5 (HDF5)
Options: SIDRE_OPT, SIDRE_LIB.
Versions: Axom >= 0.3.1.
- Conduit (optional), used when MFEM_USE_CONDUIT = YES. Conduit Mesh Blueprint
support requires Conduit >= v0.3.1 and VisIt >= v2.13.1 to read the output.
URL: https://github.com/LLNL/conduit (Conduit)
@@ -1054,6 +1069,7 @@ MFEM_USE_OCCA
MFEM_USE_CEED
MFEM_USE_RAJA
MFEM_USE_UMPIRE
MFEM_USE_SIDRE
MFEM_USE_MOONOLITH
MFEM_USE_CALIPER
MFEM_USE_FMS
@@ -1117,6 +1133,7 @@ The CMake build system adds auto-detection for the following packages/libraries:
- OCCA
- RAJA
- UMPIRE
- AXOM - Used when MFEM_USE_SIDRE is enabled
- MOONOLITH
- CALIPER
- FMS
+4
View File
@@ -248,6 +248,10 @@ IF (DEFINED TPL_ENABLE_MPFR)
SET(MFEM_USE_MPFR ${TPL_ENABLE_MPFR} CACHE BOOL "Enable MPFR usage." FORCE)
ENDIF()
IF (DEFINED TPL_ENABLE_SIDRE)
SET(MFEM_USE_SIDRE ${TPL_ENABLE_SIDRE} CACHE BOOL "Enable Axom/Sidre usage" FORCE)
ENDIF()
IF (DEFINED TPL_ENABLE_FMS)
SET(MFEM_USE_FMS ${TPL_ENABLE_FMS} CACHE BOOL "Enable FMS usage" FORCE)
ENDIF()
+1
View File
@@ -46,6 +46,7 @@ set(MFEM_USE_NETCDF @MFEM_USE_NETCDF@)
set(MFEM_USE_PETSC @MFEM_USE_PETSC@)
set(MFEM_USE_SLEPC @MFEM_USE_SLEPC@)
set(MFEM_USE_MPFR @MFEM_USE_MPFR@)
set(MFEM_USE_SIDRE @MFEM_USE_SIDRE@)
set(MFEM_USE_FMS @MFEM_USE_FMS@)
set(MFEM_USE_CONDUIT @MFEM_USE_CONDUIT@)
set(MFEM_USE_PUMI @MFEM_USE_PUMI@)
+3
View File
@@ -120,6 +120,9 @@
// Enable secure socket streams based on the GNUTLS library.
#cmakedefine MFEM_USE_GNUTLS
// Enable Sidre support.
#cmakedefine MFEM_USE_SIDRE
// Enable the use of SIMD in the high performance templated classes.
#cmakedefine MFEM_USE_SIMD
-24
View File
@@ -1,24 +0,0 @@
# Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
# LICENSE and NOTICE for details. LLNL-CODE-806117.
#
# This file is part of the MFEM library. For more information and source code
# availability visit https://mfem.org.
#
# MFEM is free software; you can redistribute it and/or modify it under the
# terms of the BSD-3 license. We welcome feedback and contributions, see file
# CONTRIBUTING.md for details.
# Defines the following variables:
# - ADIAK_FOUND
# - ADIAK_LIBRARIES
# - ADIAK_INCLUDE_DIRS
include(MfemCmakeUtilities)
mfem_find_package(Adiak ADIAK ADIAK_DIR
"include" "adiak.h"
"lib" "adiak"
"Paths to headers required by Adiak."
"Libraries required by Adiak.")
-3
View File
@@ -13,9 +13,6 @@
# - AXOM_FOUND
# - AXOM_LIBRARIES
# - AXOM_INCLUDE_DIRS
#
# MFEM itself does not depend on Axom, however Tribol does. This module exists
# to support MFEM's Tribol integration (e.g. the contact miniapp).
include(MfemCmakeUtilities)
# Note: components are enabled based on the find_package() parameters.
-36
View File
@@ -1,36 +0,0 @@
# Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
# LICENSE and NOTICE for details. LLNL-CODE-806117.
#
# This file is part of the MFEM library. For more information and source code
# availability visit https://mfem.org.
#
# MFEM is free software; you can redistribute it and/or modify it under the
# terms of the BSD-3 license. We welcome feedback and contributions, see file
# CONTRIBUTING.md for details.
# Defines the following variables:
# - CAMP_FOUND
# - CAMP_LIBRARIES
# - CAMP_INCLUDE_DIRS
include(MfemCmakeUtilities)
mfem_find_package(CAMP CAMP CAMP_DIR
"include" "camp/camp.hpp"
"lib" "camp"
"Paths to headers required by CAMP."
"Libraries required by CAMP.")
# RAJA commonly lists "camp" in INTERFACE_LINK_LIBRARIES. If there is no CMake
# target named "camp", CMake treats it as a bare library name (-lcamp).
if (CAMP_FOUND AND NOT TARGET camp)
list(GET CAMP_LIBRARIES 0 _camp_lib0)
add_library(camp UNKNOWN IMPORTED)
set_target_properties(camp PROPERTIES
IMPORTED_LOCATION "${_camp_lib0}"
INTERFACE_INCLUDE_DIRECTORIES "${CAMP_INCLUDE_DIRS}")
set(CAMP_LIBRARIES "camp" CACHE STRING "CAMP imported target." FORCE)
unset(_camp_lib0)
endif()
-15
View File
@@ -21,21 +21,6 @@ mfem_find_package(Caliper CALIPER CALIPER_DIR
"Paths to headers required by Caliper."
"Libraries required by Caliper.")
# Some downstream CMake packages (notably RAJA) may list "caliper" in their
# INTERFACE_LINK_LIBRARIES. If there is no CMake target named "caliper", CMake
# treats it as a bare library name and will pass -lcaliper to the linker.
# Create a minimal imported target when we only located the library by path.
if (CALIPER_FOUND AND NOT TARGET caliper)
list(GET CALIPER_LIBRARIES 0 _caliper_lib0)
add_library(caliper UNKNOWN IMPORTED)
set_target_properties(caliper PROPERTIES
IMPORTED_LOCATION "${_caliper_lib0}"
INTERFACE_INCLUDE_DIRECTORIES "${CALIPER_INCLUDE_DIRS}")
# Prefer linking via the target.
set(CALIPER_LIBRARIES "caliper" CACHE STRING "Caliper imported target." FORCE)
unset(_caliper_lib0)
endif()
# Append adiak path/lib if the user provided ADIAK_DIR
if(ADIAK_DIR AND EXISTS ${ADIAK_DIR})
find_package(adiak NO_DEFAULT_PATH REQUIRED PATHS ${ADIAK_DIR}/lib/cmake/adiak ${ADIAK_DIR})
+18 -90
View File
@@ -701,6 +701,7 @@ endfunction(mfem_find_library)
# Extract compile and link options needed by the given target.
#
function(mfem_get_target_options Target CompileOptsVar LinkOptsVar)
if (NOT TARGET ${Target})
return()
endif()
@@ -798,12 +799,7 @@ function(mfem_get_target_options Target CompileOptsVar LinkOptsVar)
# message(STATUS "Lib = ${Lib}")
# Filter-out generator expressions
if (NOT ("${Lib}" MATCHES "^\\$"))
if(NOT ("${Lib}" STREQUAL "dl"))
list(APPEND LinkOpts "${Lib}")
else()
# for some reason libdl doesn't include the "-l"
list(APPEND LinkOpts "-ldl")
endif()
list(APPEND LinkOpts "${Lib}")
endif()
else()
mfem_get_target_options(${Lib} COpts LOpts)
@@ -878,7 +874,7 @@ function(mfem_export_mk_files)
MFEM_USE_SUITESPARSE MFEM_USE_SUPERLU MFEM_USE_SUPERLU5 MFEM_USE_MUMPS
MFEM_USE_STRUMPACK MFEM_USE_GINKGO MFEM_USE_AMGX MFEM_USE_MAGMA
MFEM_USE_GNUTLS MFEM_USE_NETCDF MFEM_USE_PETSC MFEM_USE_SLEPC
MFEM_USE_MPFR MFEM_USE_FMS MFEM_USE_CONDUIT MFEM_USE_PUMI
MFEM_USE_MPFR MFEM_USE_SIDRE MFEM_USE_FMS MFEM_USE_CONDUIT MFEM_USE_PUMI
MFEM_USE_HIOP MFEM_USE_GSLIB MFEM_USE_CUDA MFEM_USE_HIP MFEM_USE_RAJA
MFEM_USE_OCCA MFEM_USE_CEED MFEM_USE_CALIPER MFEM_USE_UMPIRE MFEM_USE_SIMD
MFEM_USE_ADIOS2 MFEM_USE_MKL_CPARDISO MFEM_USE_MKL_PARDISO
@@ -892,18 +888,9 @@ function(mfem_export_mk_files)
set(${var} NO)
endif()
endforeach()
if (MFEM_USE_CUDA AND MFEM_EXPORT_GPU_CONFIG)
set(MFEM_CXX ${CMAKE_CUDA_COMPILER})
if(MFEM_CUDA_COMPILER_IS_NVCC)
set(MFEM_HOST_CXX ${CMAKE_CUDA_HOST_COMPILER})
else()
set(MFEM_HOST_CXX ${CMAKE_CXX_COMPILER})
endif()
else()
# mfem doesn't use enable_language(HIP)
set(MFEM_CXX ${CMAKE_CXX_COMPILER})
set(MFEM_HOST_CXX ${CMAKE_CXX_COMPILER})
endif()
# TODO: Add support for MFEM_USE_CUDA=YES
set(MFEM_CXX ${CMAKE_CXX_COMPILER})
set(MFEM_HOST_CXX ${MFEM_CXX})
set(MFEM_CPPFLAGS "")
get_target_property(cxx_std mfem CXX_STANDARD)
# For now, we ignore the setting of the CXX_EXTENSIONS property. If this
@@ -913,50 +900,6 @@ function(mfem_export_mk_files)
string(STRIP
"${cxx_std_flag} ${CMAKE_CXX_FLAGS_${BUILD_TYPE}} ${CMAKE_CXX_FLAGS}"
MFEM_CXXFLAGS)
if(MFEM_EXPORT_GPU_CONFIG)
if (MFEM_USE_CUDA)
set(MFEM_CXXFLAGS "${MFEM_CXXFLAGS} ${CMAKE_CUDA_FLAGS}")
if (MFEM_CUDA_COMPILER_IS_NVCC)
set(MFEM_CXXFLAGS "-x=cu ${MFEM_CXXFLAGS} -ccbin ${CMAKE_CXX_COMPILER} --forward-unknown-to-host-compiler")
# The following intentionally hides CUDA deprecation warnings
foreach(ENTRY IN LISTS CUDAToolkit_INCLUDE_DIRS)
set(MFEM_CXXFLAGS "${MFEM_CXXFLAGS} -isystem ${ENTRY}")
endforeach()
if (CMAKE_VERSION VERSION_GREATER_EQUAL 3.18.0)
# architecture flags not part of CMAKE_CUDA_FLAGS
if ("all" STREQUAL "${CMAKE_CUDA_ARCHITECTURES}"
OR "native" STREQUAL "${CMAKE_CUDA_ARCHITECTURES}"
OR "all-major" STREQUAL "${CMAKE_CUDA_ARCHITECTURES}")
set(MFEM_CXXFLAGS "${MFEM_CXXFLAGS} -arch=${CMAKE_CUDA_ARCHITECTURES}")
else()
foreach (ENTRY IN LISTS CMAKE_CUDA_ARCHITECTURES)
set(MFEM_CXXFLAGS
"${MFEM_CXXFLAGS} -gencode arch=compute_${ENTRY},code=sm_${ENTRY}")
endforeach()
endif()
endif()
else()
set(MFEM_CXXFLAGS "${MFEM_CXXFLAGS} -xcuda --cuda-path=${CUDAToolkit_LIBRARY_ROOT}")
if (CMAKE_VERSION VERSION_GREATER_EQUAL 3.18.0)
# architecture flags not part of CMAKE_CUDA_FLAGS
if ("all" STREQUAL "${CMAKE_CUDA_ARCHITECTURES}"
OR "native" STREQUAL "${CMAKE_CUDA_ARCHITECTURES}"
OR "all-major" STREQUAL "${CMAKE_CUDA_ARCHITECTURES}")
# TODO: not supported
else()
foreach(ENTRY IN LISTS CMAKE_CUDA_ARCHITECTURES)
set(MFEM_CXXFLAGS "-cuda-gpu-arch=sm_${ENTRY} ${MFEM_CXXFLAGS}")
endforeach()
endif()
endif()
endif()
elseif (MFEM_USE_HIP)
set(MFEM_CXXFLAGS "${MFEM_CXXFLAGS} -xhip")
foreach(ENTRY IN LISTS CMAKE_HIP_ARCHITECTURES)
set(MFEM_CXXFLAGS "--offload-arch=${ENTRY} ${MFEM_CXXFLAGS}")
endforeach()
endif()
endif()
set(MFEM_TPLFLAGS "")
foreach(dir ${TPL_INCLUDE_DIRS})
set(MFEM_TPLFLAGS "${MFEM_TPLFLAGS} -I${dir}")
@@ -987,9 +930,6 @@ function(mfem_export_mk_files)
set(MFEM_SHARED NO)
set(MFEM_STATIC YES)
endif()
if (MFEM_USE_CUDA)
set(MFEM_EXT_LIBS "${MFEM_EXT_LIBS} -lcudart")
endif()
set(MFEM_BUILD_TAG "${CMAKE_SYSTEM}")
set(MFEM_PREFIX "${CMAKE_INSTALL_PREFIX}")
# For the next 4 variables, these are the values for the build-tree version of
@@ -998,15 +938,8 @@ function(mfem_export_mk_files)
set(MFEM_LIB_DIR "${PROJECT_BINARY_DIR}")
set(MFEM_TEST_MK "${PROJECT_SOURCE_DIR}/config/test.mk")
set(MFEM_CONFIG_EXTRA "MFEM_BUILD_DIR ?= ${PROJECT_BINARY_DIR}")
if (MFEM_USE_CUDA AND MFEM_EXPORT_GPU_CONFIG)
if (MFEM_CUDA_COMPILER_IS_NVCC)
set(MFEM_XLINKER "-Xlinker=")
else()
set(MFEM_XLINKER "${CMAKE_CUDA_LINKER_WRAPPER_FLAG}")
endif()
else()
set(MFEM_XLINKER "${CMAKE_CXX_LINKER_WRAPPER_FLAG}")
endif()
# TODO: CUDA/HIP support:
set(MFEM_XLINKER "${CMAKE_CXX_LINKER_WRAPPER_FLAG}")
set(MFEM_MPIEXEC ${MPIEXEC})
if (NOT MFEM_MPIEXEC)
set(MFEM_MPIEXEC "mpirun")
@@ -1054,21 +987,16 @@ function(mfem_export_mk_files)
# handle interfaces (e.g., SCOREC::apf)
if ("${lib}" MATCHES "SCOREC::.*" OR "${lib}" MATCHES "Ginkgo::.*" OR "${lib}" MATCHES "ParMoonolith::.*")
elseif (TARGET "${lib}")
mfem_get_target_options(${lib} CompileOpts2 LinkOpts2)
# remove generator expressions
string(GENEX_STRIP "${CompileOpts2}" CompileOpts)
string(GENEX_STRIP "${LinkOpts2}" LinkOpts)
mfem_get_target_options(${lib} CompileOpts LinkOpts)
# Removing duplicates may lead to issues:
# list(REMOVE_DUPLICATES CompileOpts)
# list(REMOVE_DUPLICATES LinkOpts)
# message(WARNING "${lib}[LinkOpts]: ${LinkOpts}")
# message(WARNING "${lib}[CompileOpts]: ${CompileOpts}")
foreach(LOpt IN LISTS LinkOpts)
set(MFEM_EXT_LIBS "${MFEM_EXT_LIBS} ${LOpt}")
endforeach()
foreach(COpt IN LISTS CompileOpts)
set(MFEM_TPLFLAGS "${MFEM_TPLFLAGS} ${COpt}")
endforeach()
string(REPLACE ";" " " COpts "${CompileOpts}")
string(REPLACE ";" " " LOpts "${LinkOpts}")
# message(STATUS "${lib}[COpts]: '${COpts}'")
# message(STATUS "${lib}[LOpts]: '${LOpts}'")
set(MFEM_TPLFLAGS "${MFEM_TPLFLAGS} ${COpts}")
set(MFEM_EXT_LIBS "${MFEM_EXT_LIBS} ${LOpts}")
# message(FATAL_ERROR "***** interface lib found ... exiting *****")
# handle static and shared libs
elseif ("${suffix}" STREQUAL "${CMAKE_SHARED_LIBRARY_SUFFIX}")
@@ -1076,7 +1004,7 @@ function(mfem_export_mk_files)
get_filename_component(fullLibName ${lib} NAME_WE)
string(REGEX REPLACE "^lib" "" libname ${fullLibName})
set(MFEM_EXT_LIBS
"${MFEM_EXT_LIBS} ${shared_link_flag}${dir} -L${dir} -l${libname}")
"${MFEM_EXT_LIBS} ${shared_link_flag}${dir} -L${dir} -l${libname}")
else()
set(MFEM_EXT_LIBS "${MFEM_EXT_LIBS} ${lib}")
endif()
@@ -1085,7 +1013,7 @@ function(mfem_export_mk_files)
# Create the build-tree version of 'config.mk'
configure_file(
"${PROJECT_SOURCE_DIR}/config/config.mk.in"
"${PROJECT_BINARY_DIR}/config/config.mk" @ONLY)
"${PROJECT_BINARY_DIR}/config/config.mk")
# Copy 'test.mk' from the source-tree to the build-tree
configure_file(
"${PROJECT_SOURCE_DIR}/config/test.mk"
@@ -1103,7 +1031,7 @@ function(mfem_export_mk_files)
# Create the install-tree version of 'config.mk'
configure_file(
"${PROJECT_SOURCE_DIR}/config/config.mk.in"
"${PROJECT_BINARY_DIR}/config/config-install.mk" @ONLY)
"${PROJECT_BINARY_DIR}/config/config-install.mk")
# Install rules for 'config.mk' and 'test.mk'
install(FILES ${PROJECT_SOURCE_DIR}/config/test.mk
+3
View File
@@ -120,6 +120,9 @@
// Enable secure socket streams based on the GNUTLS library.
// #define MFEM_USE_GNUTLS
// Enable Sidre support.
// #define MFEM_USE_SIDRE
// Enable the use of SIMD in the high performance templated classes.
// #define MFEM_USE_SIMD
+1
View File
@@ -45,6 +45,7 @@ MFEM_USE_NETCDF = @MFEM_USE_NETCDF@
MFEM_USE_PETSC = @MFEM_USE_PETSC@
MFEM_USE_SLEPC = @MFEM_USE_SLEPC@
MFEM_USE_MPFR = @MFEM_USE_MPFR@
MFEM_USE_SIDRE = @MFEM_USE_SIDRE@
MFEM_USE_FMS = @MFEM_USE_FMS@
MFEM_USE_CONDUIT = @MFEM_USE_CONDUIT@
MFEM_USE_PUMI = @MFEM_USE_PUMI@
+14 -9
View File
@@ -48,6 +48,7 @@ option(MFEM_USE_NETCDF "Enable NETCDF usage" OFF)
option(MFEM_USE_PETSC "Enable PETSc support." OFF)
option(MFEM_USE_SLEPC "Enable SLEPc support." OFF)
option(MFEM_USE_MPFR "Enable MPFR usage." OFF)
option(MFEM_USE_SIDRE "Enable Axom/Sidre usage" OFF)
option(MFEM_USE_FMS "Enable FMS usage" OFF)
option(MFEM_USE_CONDUIT "Enable Conduit usage" OFF)
option(MFEM_USE_PUMI "Enable PUMI" OFF)
@@ -223,8 +224,17 @@ set(FMS_DIR "${MFEM_DIR}/../fms" CACHE PATH
set(CONDUIT_DIR "${MFEM_DIR}/../conduit" CACHE PATH
"Path to the Conduit library.")
set(AXOM_DIR "${MFEM_DIR}/../axom" CACHE PATH
"Path to the Axom library (required by Tribol for the contact mini-app).")
set(AXOM_DIR "${MFEM_DIR}/../axom" CACHE PATH "Path to the Axom library.")
# May need to add "Boost" as requirement.
if (MFEM_USE_SIDRE)
if (MFEM_USE_MPI)
set(Axom_REQUIRED_PACKAGES "Conduit/blueprint/blueprint_mpi/relay/relay_mpi" CACHE STRING
"Additional packages required by Axom.")
elseif()
set(Axom_REQUIRED_PACKAGES "Conduit/blueprint/relay" CACHE STRING
"Additional packages required by Axom.")
endif()
endif()
set(PUMI_DIR "${MFEM_DIR}/../pumi-2.1.0" CACHE STRING
"Directory where PUMI is installed")
@@ -242,7 +252,6 @@ set(MKL_PARDISO_DIR "" CACHE STRING "MKL installation path.")
set(OCCA_DIR "${MFEM_DIR}/../occa" CACHE PATH "Path to OCCA")
set(RAJA_DIR "${MFEM_DIR}/../raja" CACHE PATH "Path to RAJA")
set(CAMP_DIR "${MFEM_DIR}/../camp" CACHE PATH "Path to CAMP (required by RAJA/Umpire)")
set(CEED_DIR "${MFEM_DIR}/../libCEED" CACHE PATH "Path to libCEED")
set(UMPIRE_DIR "${MFEM_DIR}/../umpire" CACHE PATH "Path to Umpire")
set(CALIPER_DIR "${MFEM_DIR}/../caliper" CACHE PATH "Path to Caliper")
@@ -263,12 +272,8 @@ set(PARELAG_LIBRARIES "${PARELAG_DIR}/build/src/libParELAG.a" CACHE STRING
"The ParELAG library.")
set(TRIBOL_DIR "${MFEM_DIR}/../tribol" CACHE PATH "Path to Tribol")
# Tribol requires Axom. Many Tribol builds also enable optional TPLs like
# RAJA/UMPIRE/Caliper, and may pull additional Axom components (e.g. quest,
# lumberjack) via its exported targets.
set(Tribol_REQUIRED_PACKAGES
"REQUIRED:;Axom/core/primal/mint/slam/slic/quest/lumberjack;OPTIONAL:;Adiak;CAMP;RAJA;UMPIRE;Caliper"
CACHE STRING "Additional packages required by Tribol")
set(Tribol_REQUIRED_PACKAGES "Axom/core/mint/slam/slic" CACHE STRING
"Additional packages required by Tribol")
set(ENZYME_DIR "${MFEM_DIR}/../enzyme" CACHE PATH "Path to Enzyme")
+15 -78
View File
@@ -162,6 +162,7 @@ MFEM_USE_NETCDF = NO
MFEM_USE_PETSC = NO
MFEM_USE_SLEPC = NO
MFEM_USE_MPFR = NO
MFEM_USE_SIDRE = NO
MFEM_USE_FMS = NO
MFEM_USE_CONDUIT = NO
MFEM_USE_PUMI = NO
@@ -248,15 +249,6 @@ endif
# METIS library configuration
ifeq ($(MFEM_USE_SUPERLU)$(MFEM_USE_STRUMPACK)$(MFEM_USE_MUMPS),NONONO)
# MFEM_USE_METIS_5: when the user supplies METIS_DIR, try to auto-detect
# METIS 5 installs that follow the common <prefix>/{include,lib,lib64} layout.
ifeq ($(MFEM_USE_METIS_5),NO)
ifneq ($(wildcard $(METIS_DIR)/include/metis.h),)
ifneq ($(wildcard $(METIS_DIR)/lib/libmetis.* $(METIS_DIR)/lib64/libmetis.*),)
MFEM_USE_METIS_5 = YES
endif
endif
endif
ifeq ($(MFEM_USE_METIS_5),NO)
METIS_DIR = @MFEM_DIR@/../metis-4.0
METIS_OPT =
@@ -495,6 +487,17 @@ ifneq (,$(wildcard $(CONDUIT_HDF5_HEADER)))
-lhdf5 $(ZLIB_LIB)
endif
# Sidre and required libraries configuration
# Be sure to check the HDF5_DIR (set above) is correct
SIDRE_DIR = @MFEM_DIR@/../axom
SIDRE_OPT = -I$(SIDRE_DIR)/include -I$(CONDUIT_DIR)/include/conduit\
-I$(HDF5_DIR)/include
SIDRE_LIB = \
$(XLINKER)-rpath,$(SIDRE_DIR)/lib -L$(SIDRE_DIR)/lib \
$(XLINKER)-rpath,$(CONDUIT_DIR)/lib -L$(CONDUIT_DIR)/lib \
$(XLINKER)-rpath,$(HDF5_DIR)/lib -L$(HDF5_DIR)/lib \
-laxom -lconduit -lconduit_relay -lconduit_blueprint -lhdf5 $(ZLIB_LIB) -ldl
# PUMI
# Note that PUMI_DIR is needed -- it is used to check for gmi_sim.h
PUMI_DIR = @MFEM_DIR@/../pumi-2.1.0
@@ -576,13 +579,7 @@ ifdef CUB_DIR
RAJA_OPT += -I$(CUB_DIR)
endif
# CAMP library configuration (required by RAJA/Umpire for most installs)
CAMP_LIB = -lcamp
# If the common sibling layout exists, use it as a default (handles versioned
# directories like camp-<hash>).
ifneq ($(wildcard $(RAJA_DIR)/../camp*/include/camp/camp.hpp),)
CAMP_DIR ?= $(patsubst %/include/camp/camp.hpp,%,$(firstword $(wildcard $(RAJA_DIR)/../camp*/include/camp/camp.hpp)))
endif
ifdef CAMP_DIR
RAJA_OPT += -I$(CAMP_DIR)/include
CAMP_LIB = $(XLINKER)-rpath,$(CAMP_DIR)/lib -L$(CAMP_DIR)/lib -lcamp
@@ -592,12 +589,7 @@ RAJA_LIB = $(XLINKER)-rpath,$(RAJA_DIR)/lib -L$(RAJA_DIR)/lib -lRAJA $(CAMP_LIB)
# UMPIRE library configuration
UMPIRE_DIR = @MFEM_DIR@/../umpire
UMPIRE_OPT = -I$(UMPIRE_DIR)/include $(if $(CAMP_DIR), -I$(CAMP_DIR)/include)
UMPIRE_LIB = -L$(UMPIRE_DIR)/lib -L$(UMPIRE_DIR)/lib64 -lumpire $(CAMP_LIB) -lpthread
# If the common sibling layout exists, use it as a default (handles versioned
# directories like fmt-<hash>).
ifneq ($(wildcard $(UMPIRE_DIR)/../fmt*/include/fmt/format.h),)
FMT_DIR ?= $(patsubst %/include/fmt/format.h,%,$(firstword $(wildcard $(UMPIRE_DIR)/../fmt*/include/fmt/format.h)))
endif
UMPIRE_LIB = -L$(UMPIRE_DIR)/lib -L$(UMPIRE_DIR)/lib64 -lumpire $(CAMP_LIB)
ifdef FMT_DIR
UMPIRE_OPT += -I$(FMT_DIR)/include
UMPIRE_LIB += -L$(FMT_DIR)/lib -L$(FMT_DIR)/lib64 -lfmt
@@ -629,63 +621,8 @@ PARELAG_LIB = -L$(PARELAG_DIR)/build/src -lParELAG
AXOM_DIR = @MFEM_DIR@/../axom
TRIBOL_DIR = @MFEM_DIR@/../tribol
TRIBOL_OPT = -I$(TRIBOL_DIR)/include -I$(AXOM_DIR)/include
# Tribol may be built with optional dependencies (e.g. RAJA/UMPIRE/CALIPER).
# Add those options only when the corresponding headers/libraries exist.
ifneq ($(wildcard $(RAJA_DIR)/include/RAJA/RAJA.hpp),)
TRIBOL_OPT += $(RAJA_OPT)
endif
ifneq ($(wildcard $(UMPIRE_DIR)/include/umpire/Umpire.hpp),)
TRIBOL_OPT += $(UMPIRE_OPT)
endif
ifneq ($(wildcard $(CALIPER_DIR)/include/caliper/cali.h),)
TRIBOL_OPT += $(CALIPER_OPT)
endif
TRIBOL_LIB = -L$(TRIBOL_DIR)/lib -L$(TRIBOL_DIR)/lib64
ifneq ($(wildcard $(TRIBOL_DIR)/lib/libtribol.* $(TRIBOL_DIR)/lib64/libtribol.*),)
TRIBOL_LIB += -ltribol
endif
ifneq ($(wildcard $(TRIBOL_DIR)/lib/libtribol_shared.* $(TRIBOL_DIR)/lib64/libtribol_shared.*),)
TRIBOL_LIB += -ltribol_shared
endif
ifneq ($(wildcard $(TRIBOL_DIR)/lib/libredecomp.* $(TRIBOL_DIR)/lib64/libredecomp.*),)
TRIBOL_LIB += -lredecomp
endif
TRIBOL_LIB += -L$(AXOM_DIR)/lib -L$(AXOM_DIR)/lib64
ifneq ($(wildcard $(AXOM_DIR)/lib/libaxom_quest.* $(AXOM_DIR)/lib64/libaxom_quest.*),)
TRIBOL_LIB += -laxom_quest
endif
ifneq ($(wildcard $(AXOM_DIR)/lib/libaxom_mint.* $(AXOM_DIR)/lib64/libaxom_mint.*),)
TRIBOL_LIB += -laxom_mint
endif
ifneq ($(wildcard $(AXOM_DIR)/lib/libaxom_slam.* $(AXOM_DIR)/lib64/libaxom_slam.*),)
TRIBOL_LIB += -laxom_slam
endif
ifneq ($(wildcard $(AXOM_DIR)/lib/libaxom_slic.* $(AXOM_DIR)/lib64/libaxom_slic.*),)
TRIBOL_LIB += -laxom_slic
endif
ifneq ($(wildcard $(AXOM_DIR)/lib/libaxom_lumberjack.* $(AXOM_DIR)/lib64/libaxom_lumberjack.*),)
TRIBOL_LIB += -laxom_lumberjack
endif
ifneq ($(wildcard $(AXOM_DIR)/lib/libaxom_core.* $(AXOM_DIR)/lib64/libaxom_core.*),)
TRIBOL_LIB += -laxom_core
endif
# Add common optional Tribol TPLs when their libraries are present.
ifneq ($(wildcard $(ADIAK_DIR)/lib/libadiak.* $(ADIAK_DIR)/lib64/libadiak.*),)
TRIBOL_LIB += $(XLINKER)-rpath,$(ADIAK_DIR)/lib64 $(XLINKER)-rpath,$(ADIAK_DIR)/lib \
-L$(ADIAK_DIR)/lib64 -L$(ADIAK_DIR)/lib -ladiak -ldl
endif
ifneq ($(wildcard $(UMPIRE_DIR)/lib/libumpire.* $(UMPIRE_DIR)/lib64/libumpire.*),)
TRIBOL_LIB += $(UMPIRE_LIB)
endif
ifneq ($(wildcard $(RAJA_DIR)/lib/libRAJA.* $(RAJA_DIR)/lib64/libRAJA.*),)
TRIBOL_LIB += $(RAJA_LIB)
endif
ifneq ($(wildcard $(CALIPER_DIR)/lib/libcaliper.* $(CALIPER_DIR)/lib64/libcaliper.*),)
TRIBOL_LIB += $(CALIPER_LIB)
endif
TRIBOL_LIB = -L$(TRIBOL_DIR)/lib -ltribol -lredecomp -L$(AXOM_DIR)/lib -laxom_mint\
-laxom_slam -laxom_slic -laxom_core
# Enzyme configuration
ENZYME_DIR = @MFEM_DIR@/../enzyme
+2 -2
View File
@@ -5,9 +5,9 @@
// Sample runs:
// mpirun -np 4 ex12p -m ../data/beam-tri.mesh
// mpirun -np 4 ex12p -m ../data/beam-quad.mesh
// mpirun -np 4 ex12p -m ../data/beam-tet.mesh -s 464 -n 10 -o 2 -elast
// mpirun -np 4 ex12p -m ../data/beam-tet.mesh -s 462 -n 10 -o 2 -elast
// mpirun -np 4 ex12p -m ../data/beam-hex.mesh -s 3878
// mpirun -np 4 ex12p -m ../data/beam-wedge.mesh -s 82
// mpirun -np 4 ex12p -m ../data/beam-wedge.mesh -s 81
// mpirun -np 4 ex12p -m ../data/beam-tri.mesh -s 3877 -o 2 -sys
// mpirun -np 4 ex12p -m ../data/beam-quad.mesh -s 4544 -n 6 -o 3 -elast
// mpirun -np 4 ex12p -m ../data/beam-quad-nurbs.mesh
+9 -27
View File
@@ -302,21 +302,15 @@ int main(int argc, char *argv[])
char vishost[] = "localhost";
int visport = 19916;
socketstream sol_sock_r(vishost, visport);
socketstream sol_sock_i(vishost, visport);
sol_sock_r << "parallel " << num_procs << " " << myid << "\n";
sol_sock_i << "parallel " << num_procs << " " << myid << "\n";
sol_sock_r.precision(8);
sol_sock_i.precision(8);
sol_sock_r << "solution\n" << *pmesh << u_exact->real()
<< "window_title 'Exact: Real Part'" << flush;
// Make sure all ranks have sent their real solution before initiating
// another set of GLVis connections (one from each rank):
MPI_Barrier(pmesh->GetComm());
socketstream sol_sock_i(vishost, visport);
sol_sock_i << "parallel " << num_procs << " " << myid << "\n";
sol_sock_i.precision(8);
sol_sock_i << "solution\n" << *pmesh << u_exact->imag()
<< "window_title 'Exact: Imaginary Part'" << flush;
// Make sure all ranks have sent their imaginary solution before initiating
// another set of GLVis connections (one from each rank):
MPI_Barrier(pmesh->GetComm());
}
// 11. Set up the parallel sesquilinear form a(.,.) on the finite element
@@ -540,21 +534,15 @@ int main(int argc, char *argv[])
char vishost[] = "localhost";
int visport = 19916;
socketstream sol_sock_r(vishost, visport);
socketstream sol_sock_i(vishost, visport);
sol_sock_r << "parallel " << num_procs << " " << myid << "\n";
sol_sock_i << "parallel " << num_procs << " " << myid << "\n";
sol_sock_r.precision(8);
sol_sock_i.precision(8);
sol_sock_r << "solution\n" << *pmesh << u.real()
<< "window_title 'Solution: Real Part'" << flush;
// Make sure all ranks have sent their real solution before initiating
// another set of GLVis connections (one from each rank):
MPI_Barrier(pmesh->GetComm());
socketstream sol_sock_i(vishost, visport);
sol_sock_i << "parallel " << num_procs << " " << myid << "\n";
sol_sock_i.precision(8);
sol_sock_i << "solution\n" << *pmesh << u.imag()
<< "window_title 'Solution: Imaginary Part'" << flush;
// Make sure all ranks have sent their imaginary solution before initiating
// another set of GLVis connections (one from each rank):
MPI_Barrier(pmesh->GetComm());
}
if (visualization && exact_sol)
{
@@ -563,21 +551,15 @@ int main(int argc, char *argv[])
char vishost[] = "localhost";
int visport = 19916;
socketstream sol_sock_r(vishost, visport);
socketstream sol_sock_i(vishost, visport);
sol_sock_r << "parallel " << num_procs << " " << myid << "\n";
sol_sock_i << "parallel " << num_procs << " " << myid << "\n";
sol_sock_r.precision(8);
sol_sock_i.precision(8);
sol_sock_r << "solution\n" << *pmesh << u_exact->real()
<< "window_title 'Error: Real Part'" << flush;
// Make sure all ranks have sent their real solution before initiating
// another set of GLVis connections (one from each rank):
MPI_Barrier(pmesh->GetComm());
socketstream sol_sock_i(vishost, visport);
sol_sock_i << "parallel " << num_procs << " " << myid << "\n";
sol_sock_i.precision(8);
sol_sock_i << "solution\n" << *pmesh << u_exact->imag()
<< "window_title 'Error: Imaginary Part'" << flush;
// Make sure all ranks have sent their imaginary solution before initiating
// another set of GLVis connections (one from each rank):
MPI_Barrier(pmesh->GetComm());
}
if (visualization)
{
+4 -6
View File
@@ -433,18 +433,16 @@ int main(int argc, char *argv[])
u.ProjectCoefficient(*u0);
// Create data collection for solution output: either VisItDataCollection for
// ascii data files, or ConduitDataCollection for binary data files.
// ascii data files, or SidreDataCollection for binary data files.
DataCollection *dc = NULL;
if (visit)
{
if (binary)
{
#ifdef MFEM_USE_CONDUIT
auto conduit_dc = new ConduitDataCollection("Example41", &mesh);
conduit_dc->SetProtocol("hdf5");
dc = conduit_dc;
#ifdef MFEM_USE_SIDRE
dc = new SidreDataCollection("Example41", &mesh);
#else
MFEM_ABORT("Must build with MFEM_USE_CONDUIT=YES for binary output.");
MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
#endif
}
else
+3 -5
View File
@@ -518,12 +518,10 @@ int main(int argc, char *argv[])
{
if (binary)
{
#ifdef MFEM_USE_CONDUIT
auto conduit_dc = new ConduitDataCollection("Example41-Parallel", pmesh);
conduit_dc->SetProtocol("hdf5");
dc = conduit_dc;
#ifdef MFEM_USE_SIDRE
dc = new SidreDataCollection("Example41-Parallel", pmesh);
#else
MFEM_ABORT("Must build with MFEM_USE_CONDUIT=YES for binary output.");
MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
#endif
}
else
+4 -6
View File
@@ -305,18 +305,16 @@ int main(int argc, char *argv[])
}
// Create data collection for solution output: either VisItDataCollection for
// ascii data files, or ConduitDataCollection for binary data files.
// ascii data files, or SidreDataCollection for binary data files.
DataCollection *dc = NULL;
if (visit)
{
if (binary)
{
#ifdef MFEM_USE_CONDUIT
auto conduit_dc = new ConduitDataCollection("Example9", &mesh);
conduit_dc->SetProtocol("hdf5");
dc = conduit_dc;
#ifdef MFEM_USE_SIDRE
dc = new SidreDataCollection("Example9", &mesh);
#else
MFEM_ABORT("Must build with MFEM_USE_CONDUIT=YES for binary output.");
MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
#endif
}
else
+4 -6
View File
@@ -441,18 +441,16 @@ int main(int argc, char *argv[])
}
// Create data collection for solution output: either VisItDataCollection for
// ascii data files, or ConduitDataCollection for binary data files.
// ascii data files, or SidreDataCollection for binary data files.
DataCollection *dc = NULL;
if (visit)
{
if (binary)
{
#ifdef MFEM_USE_CONDUIT
auto conduit_dc = new ConduitDataCollection("Example9-Parallel", pmesh);
conduit_dc->SetProtocol("hdf5");
dc = conduit_dc;
#ifdef MFEM_USE_SIDRE
dc = new SidreDataCollection("Example9-Parallel", pmesh);
#else
MFEM_ABORT("Must build with MFEM_USE_CONDUIT=YES for binary output.");
MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
#endif
}
else
+4 -6
View File
@@ -354,18 +354,16 @@ int main(int argc, char *argv[])
}
// Create data collection for solution output: either VisItDataCollection for
// ascii data files, or ConduitDataCollection for binary data files.
// ascii data files, or SidreDataCollection for binary data files.
DataCollection *dc = NULL;
if (visit)
{
if (binary)
{
#ifdef MFEM_USE_CONDUIT
auto conduit_dc = new ConduitDataCollection("Example9", mesh);
conduit_dc->SetProtocol("hdf5");
dc = conduit_dc;
#ifdef MFEM_USE_SIDRE
dc = new SidreDataCollection("Example9", mesh);
#else
MFEM_ABORT("Must build with MFEM_USE_CONDUIT=YES for binary output.");
MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
#endif
}
else
+4 -6
View File
@@ -414,18 +414,16 @@ int main(int argc, char *argv[])
}
// Create data collection for solution output: either VisItDataCollection for
// ascii data files, or ConduitDataCollection for binary data files.
// ascii data files, or SidreDataCollection for binary data files.
DataCollection *dc = NULL;
if (visit)
{
if (binary)
{
#ifdef MFEM_USE_CONDUIT
auto conduit_dc = new ConduitDataCollection("Example9-Parallel", pmesh);
conduit_dc->SetProtocol("hdf5");
dc = conduit_dc;
#ifdef MFEM_USE_SIDRE
dc = new SidreDataCollection("Example9-Parallel", pmesh);
#else
MFEM_ABORT("Must build with MFEM_USE_CONDUIT=YES for binary output.");
MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
#endif
}
else
+4 -6
View File
@@ -368,18 +368,16 @@ int main(int argc, char *argv[])
}
// Create data collection for solution output: either VisItDataCollection for
// ascii data files, or ConduitDataCollection for binary data files.
// ascii data files, or SidreDataCollection for binary data files.
DataCollection *dc = NULL;
if (visit)
{
if (binary)
{
#ifdef MFEM_USE_CONDUIT
auto conduit_dc = new ConduitDataCollection("Example9-Parallel", pmesh);
conduit_dc->SetProtocol("hdf5");
dc = conduit_dc;
#ifdef MFEM_USE_SIDRE
dc = new SidreDataCollection("Example9-Parallel", pmesh);
#else
MFEM_ABORT("Must build with MFEM_USE_CONDUIT=YES for binary output.");
MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
#endif
}
else
+4 -6
View File
@@ -316,18 +316,16 @@ int main(int argc, char *argv[])
}
// Create data collection for solution output: either VisItDataCollection for
// ascii data files, or ConduitDataCollection for binary data files.
// ascii data files, or SidreDataCollection for binary data files.
DataCollection *dc = NULL;
if (visit)
{
if (binary)
{
#ifdef MFEM_USE_CONDUIT
auto conduit_dc = new ConduitDataCollection("Example9", &mesh);
conduit_dc->SetProtocol("hdf5");
dc = conduit_dc;
#ifdef MFEM_USE_SIDRE
dc = new SidreDataCollection("Example9", &mesh);
#else
MFEM_ABORT("Must build with MFEM_USE_CONDUIT=YES for binary output.");
MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
#endif
}
else
+4 -6
View File
@@ -453,18 +453,16 @@ int main(int argc, char *argv[])
}
// Create data collection for solution output: either VisItDataCollection for
// ascii data files, or ConduitDataCollection for binary data files.
// ascii data files, or SidreDataCollection for binary data files.
DataCollection *dc = NULL;
if (visit)
{
if (binary)
{
#ifdef MFEM_USE_CONDUIT
auto conduit_dc = new ConduitDataCollection("Example9-Parallel", pmesh);
conduit_dc->SetProtocol("hdf5");
dc = conduit_dc;
#ifdef MFEM_USE_SIDRE
dc = new SidreDataCollection("Example9-Parallel", pmesh);
#else
MFEM_ABORT("Must build with MFEM_USE_CONDUIT=YES for binary output.");
MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
#endif
}
else
+5
View File
@@ -312,6 +312,11 @@ set(HDRS
particleset.hpp
)
if (MFEM_USE_SIDRE)
list(APPEND SRCS sidredatacollection.cpp)
list(APPEND HDRS sidredatacollection.hpp)
endif()
if (MFEM_USE_CONDUIT)
list(APPEND SRCS conduitdatacollection.cpp)
list(APPEND HDRS conduitdatacollection.hpp)
+3 -7
View File
@@ -729,8 +729,7 @@ void BilinearForm::Assemble(int skip_zeros)
tr = mesh -> GetBdrFaceTransformations (i);
if (tr != NULL)
{
mfem::DofTransformation doftrans;
fes -> GetElementVDofs (tr -> Elem1No, vdofs, doftrans);
fes -> GetElementVDofs (tr -> Elem1No, vdofs);
fe1 = fes -> GetFE (tr -> Elem1No);
// The fe2 object is really a dummy and not used on the boundaries,
// but we can't dereference a NULL pointer, and we don't want to
@@ -744,7 +743,6 @@ void BilinearForm::Assemble(int skip_zeros)
boundary_face_integs[k] -> AssembleFaceMatrix (*fe1, *fe2, *tr,
elemmat);
doftrans.TransformDual(elemmat);
mat -> AddSubMatrix (vdofs, vdofs, elemmat, skip_zeros);
}
}
@@ -1725,7 +1723,6 @@ void MixedBilinearForm::Assemble(int skip_zeros)
}
}
DofTransformation dom_dof_trans, ran_dof_trans;
for (int i = 0; i < trial_fes -> GetNBE(); i++)
{
const int bdr_attr = mesh->GetBdrAttribute(i);
@@ -1734,8 +1731,8 @@ void MixedBilinearForm::Assemble(int skip_zeros)
ftr = mesh -> GetBdrFaceTransformations (i);
if (ftr != NULL)
{
trial_fes->GetElementVDofs(ftr->Elem1No, trial_vdofs, dom_dof_trans);
test_fes->GetElementVDofs(ftr->Elem1No, test_vdofs, ran_dof_trans);
trial_fes->GetElementVDofs(ftr->Elem1No, trial_vdofs);
test_fes->GetElementVDofs(ftr->Elem1No, test_vdofs);
trial_fe1 = trial_fes->GetFE(ftr->Elem1No);
test_fe1 = test_fes->GetFE(ftr->Elem1No);
// The test_fe2 object is really a dummy and not used on the
@@ -1751,7 +1748,6 @@ void MixedBilinearForm::Assemble(int skip_zeros)
boundary_face_integs[k]->AssembleFaceMatrix(*trial_fe1, *test_fe1, *trial_fe2,
*test_fe2,
*ftr, elemmat);
TransformDual(ran_dof_trans, dom_dof_trans, elemmat);
mat->AddSubMatrix(test_vdofs, trial_vdofs, elemmat, skip_zeros);
}
}
+1 -1
View File
@@ -2710,7 +2710,7 @@ public:
/** Integrator for $(-Q u, \nabla v)$ for Nedelec ($u$) and $H^1$ ($v$) elements.
This is equivalent to a weak divergence of the $H(curl)$ basis functions. */
This is equivalent to a weak divergence of the $H(curl$ basis functions. */
class VectorFEWeakDivergenceIntegrator: public BilinearFormIntegrator
{
protected:
-19
View File
@@ -82,25 +82,6 @@ public:
/// underlying #fes
int VectorDim() const;
/// Copy assignment. Only the data of the base class Vector is copied.
/** It is assumed that this object and @a rhs use FiniteElementSpace%s that
have the same size.
@note Defining this method overwrites the implicitly defined copy
assignment operator. */
ComplexGridFunction &operator=(const ComplexGridFunction &rhs)
{ return operator=((const Vector &)rhs); }
/// Copy the data from @a v.
/** The size of @a v must be equal to double of the size of the associated
FiniteElementSpace #fes. */
ComplexGridFunction &operator=(const Vector &v)
{
MFEM_ASSERT(fes && v.Size() == 2*fes->GetVSize(), "");
Vector::operator=(v);
return *this;
}
/// Assign constant values to the ComplexGridFunction data.
ComplexGridFunction &operator=(const std::complex<real_t> & value)
{ *gfr = value.real(); *gfi = value.imag(); return *this; }
+2
View File
@@ -1453,6 +1453,8 @@ ConduitDataCollection::LoadMeshAndFields(int domain_id,
std::string
ConduitDataCollection::ElementTypeToShapeName(Element::Type element_type)
{
// Adapted from SidreDataCollection
// Note -- the mapping from Element::Type to string is based on
// enum Element::Type { POINT, SEGMENT, TRIANGLE, QUADRILATERAL,
// TETRAHEDRON, HEXAHEDRON };
+4 -4
View File
@@ -34,10 +34,10 @@ namespace mfem
- HDF5 library, https://support.hdfgroup.org/HDF5
@note The ConduitDataCollection only wraps the mfem objects to save them and
creates them on load, Conduit does not own any of the data.
The MFEMSidreDataCollection in the Axom package (https://github.com/LLNL/axom)
derives from mfem::DataCollection and provides more features, for example
it allocates and will own the data backing the mfem objects in the data collection.
creates them on load, Conduit does not own any of the data. The
SidreDataCollection provides more features, for example the
SidreDataCollection allocates and will own the data backing the mfem objects
in the data collection.
This class also provides public static methods that convert between MFEM
Meshes and GridFunctions and Conduit Mesh Blueprint descriptions.
+215 -41
View File
@@ -84,29 +84,33 @@ public:
const std::vector<derivative_action_t> &derivative_actions,
const FieldDescriptor &direction,
const int &daction_l_size,
const std::vector<derivative_action_t> &derivative_actions_transpose,
const int &derivative_action_tr_l_size,
const std::vector<derivative_action_t> &derivative_tr_actions,
const FieldDescriptor &transpose_direction,
const int &daction_transpose_l_size,
const std::vector<Vector *> &solutions_l,
const std::vector<Vector *> &parameters_l,
const restriction_callback_t &restriction_callback,
const std::function<void(Vector &, Vector &)> &prolongation_transpose,
const std::function<void(Vector &, Vector &)> &tr_prolongation_transpose,
const std::vector<assemble_derivative_sparsematrix_callback_t>
&assemble_derivative_sparsematrix_callbacks,
const std::vector<assemble_derivative_hypreparmatrix_callback_t>
&assemble_derivative_hypreparmatrix_callbacks) :
const assemble_derivative_hypreparmatrix_callback_t
&assemble_derivative_hypreparmatrix_callback) :
Operator(height, width),
derivative_actions(derivative_actions),
direction(direction),
daction_l(daction_l_size),
daction_l_size(daction_l_size),
derivative_actions_transpose(derivative_actions_transpose),
derivative_action_tr_l_size(derivative_action_tr_l_size),
derivative_tr_actions(derivative_tr_actions),
transpose_direction(transpose_direction),
prolongation_transpose(prolongation_transpose),
tr_prolongation_transpose(tr_prolongation_transpose),
assemble_derivative_sparsematrix_callbacks(
assemble_derivative_sparsematrix_callbacks),
assemble_derivative_hypreparmatrix_callbacks(
assemble_derivative_hypreparmatrix_callbacks)
assemble_derivative_hypreparmatrix_callback(
assemble_derivative_hypreparmatrix_callback)
{
std::vector<Vector> s_l(solutions_l.size());
for (size_t i = 0; i < s_l.size(); i++)
@@ -156,18 +160,18 @@ public:
/// direction_t on T-dofs.
void MultTranspose(const Vector &direction_t, Vector &result_t) const override
{
MFEM_ASSERT(!derivative_actions_transpose.empty(),
MFEM_ASSERT(!derivative_tr_actions.empty(),
"derivative can't be used to be multiplied in transpose mode");
daction_l.SetSize(width);
daction_l.SetSize(derivative_action_tr_l_size);
daction_l = 0.0;
prolongation(transpose_direction, direction_t, direction_l);
for (const auto &f : derivative_actions_transpose)
for (const auto &f : derivative_tr_actions)
{
f(fields_e, direction_l, daction_l);
}
prolongation_transpose(daction_l, result_t);
tr_prolongation_transpose(daction_l, result_t);
};
/// @brief Assemble the derivative operator into a SparseMatrix.
@@ -183,6 +187,10 @@ public:
{
f(fields_e, A);
}
// SparseMatrix A is finalized after all callbacks have contributed to
// it.
A->Finalize();
}
/// @brief Assemble the derivative operator into a HypreParMatrix.
@@ -191,13 +199,7 @@ public:
/// be an uninitialized object.
void Assemble(HypreParMatrix *&A)
{
MFEM_ASSERT(!assemble_derivative_hypreparmatrix_callbacks.empty(),
"derivative can't be assembled into a HypreParMatrix");
for (const auto &f : assemble_derivative_hypreparmatrix_callbacks)
{
f(fields_e, A);
}
assemble_derivative_hypreparmatrix_callback(fields_e, A);
}
private:
@@ -212,10 +214,12 @@ private:
const int daction_l_size;
const int derivative_action_tr_l_size;
/// Transpose Derivative action callbacks. Depending on the requested
/// derivatives in DifferentiableOperator the callbacks represent certain
/// combinations of actions of derivatives of the forward operator.
std::vector<derivative_action_t> derivative_actions_transpose;
std::vector<derivative_action_t> derivative_tr_actions;
FieldDescriptor transpose_direction;
@@ -225,13 +229,15 @@ private:
std::function<void(Vector &, Vector &)> prolongation_transpose;
std::function<void(Vector &, Vector &)> tr_prolongation_transpose;
/// Callbacks that assemble derivatives into a SparseMatrix.
std::vector<assemble_derivative_sparsematrix_callback_t>
assemble_derivative_sparsematrix_callbacks;
/// Callbacks that assemble derivatives into a HypreParMatrix.
std::vector<assemble_derivative_hypreparmatrix_callback_t>
assemble_derivative_hypreparmatrix_callbacks;
assemble_derivative_hypreparmatrix_callback_t
assemble_derivative_hypreparmatrix_callback;
};
/// Class representing a differentiable operator which acts on solution and
@@ -457,7 +463,10 @@ public:
dir_l = s_l[derivative_idx];
}
derivative_setup_callbacks[derivative_id][0](fields_e, dir_l);
for (size_t i = 0; i < derivative_setup_callbacks[derivative_id].size(); i++)
{
derivative_setup_callbacks[derivative_id][i](fields_e, dir_l);
}
return std::make_shared<DerivativeOperator>(
height,
@@ -465,15 +474,17 @@ public:
derivative_action_callbacks[derivative_id],
fields[derivative_idx],
residual_l.Size(),
daction_transpose_callbacks[derivative_id],
derivative_action_tr_l_size[derivative_id],
derivative_action_tr_callbacks[derivative_id],
fields[test_space_field_idx],
GetVSize(fields[test_space_field_idx]),
sol_l,
par_l,
restriction_callback,
prolongation_transpose,
derivative_tr_prolongation_transpose[derivative_id],
assemble_derivative_sparsematrix_callbacks[derivative_id],
assemble_derivative_hypreparmatrix_callbacks[derivative_id]);
assemble_derivative_hypreparmatrix_callback[derivative_id]);
}
private:
@@ -486,13 +497,15 @@ private:
std::map<size_t,
std::vector<derivative_action_t>> derivative_action_callbacks;
std::map<size_t,
std::vector<derivative_action_t>> daction_transpose_callbacks;
std::vector<derivative_action_t>> derivative_action_tr_callbacks;
std::map<size_t,
std::function<void(Vector &, Vector &)>> derivative_tr_prolongation_transpose;
std::map<size_t, int> derivative_action_tr_l_size;
std::map<size_t,
std::vector<assemble_derivative_sparsematrix_callback_t>>
assemble_derivative_sparsematrix_callbacks;
std::map<size_t,
std::vector<assemble_derivative_hypreparmatrix_callback_t>>
assemble_derivative_hypreparmatrix_callbacks;
std::map<size_t, assemble_derivative_hypreparmatrix_callback_t>
assemble_derivative_hypreparmatrix_callback;
std::vector<FieldDescriptor> solutions;
std::vector<FieldDescriptor> parameters;
@@ -510,7 +523,7 @@ private:
std::function<void(Vector &, Vector &)> output_restriction_transpose;
restriction_callback_t restriction_callback;
std::map<size_t, Vector> derivative_qp_caches;
std::map<size_t, std::vector<Vector>> derivative_qp_caches;
std::map<size_t, size_t> assembled_vector_sizes;
@@ -769,9 +782,10 @@ void DifferentiableOperator::AddIntegrator(
auto input_size_on_qp =
get_input_size_on_qp(inputs, std::make_index_sequence<num_inputs> {});
// printf("calculate shmem action info\n");
auto action_shmem_info =
get_shmem_info<entity_t, num_fields, num_inputs, num_outputs>
(input_dtq_maps, output_dtq_maps, fields, num_entities, inputs, num_qp,
(input_dtq_maps, output_dtq_maps, fields, num_entities, num_qp,
input_size_on_qp, residual_size_on_qp, element_dof_ordering);
Vector shmem_cache(action_shmem_info.total_size);
@@ -894,9 +908,10 @@ void DifferentiableOperator::AddIntegrator(
const int da_size_on_qp =
GetSizeOnQP<entity_t>(output_fop, fields[test_space_field_idx]);
// printf("calculate shmem derivative action info\n");
auto shmem_info =
get_shmem_info<entity_t, num_fields, num_inputs, num_outputs>(
input_dtq_maps, output_dtq_maps, fields, num_entities, inputs,
input_dtq_maps, output_dtq_maps, fields, num_entities,
num_qp, input_size_on_qp, residual_size_on_qp,
element_dof_ordering, d_field_idx);
@@ -968,9 +983,12 @@ void DifferentiableOperator::AddIntegrator(
// Quadrature point local derivative cache for each element, with data
// layout:
// [test_vdim, test_op_dim, trial_vdim, trial_op_dim, qp, num_entities].
derivative_qp_caches[derivative_id] = Vector(test_vdim * test_op_dim *
trial_vdim *
total_trial_op_dim * num_qp * num_entities);
derivative_qp_caches[derivative_id].push_back(
Vector(test_vdim * test_op_dim * trial_vdim * total_trial_op_dim * num_qp *
num_entities));
const int cache_index = this->derivative_qp_caches[derivative_id].size() - 1;
// Create local references for MSVC lambda capture compatibility
auto& fields_ref = this->fields;
auto& derivative_qp_caches_ref = this->derivative_qp_caches[derivative_id];
@@ -1012,6 +1030,8 @@ void DifferentiableOperator::AddIntegrator(
trial_vdim,
inputs_trial_op_dim,
qpdc_idx = cache_index,
// capture by ref:
&qpdc_mem = derivative_qp_caches_ref
](std::vector<Vector> &f_e, const Vector &dir_l) mutable
@@ -1024,7 +1044,7 @@ void DifferentiableOperator::AddIntegrator(
shmem_info.direction_size,
num_entities);
auto qpdc = Reshape(qpdc_mem.ReadWrite(), test_vdim, test_op_dim,
auto qpdc = Reshape(qpdc_mem[qpdc_idx].ReadWrite(), test_vdim, test_op_dim,
trial_vdim, total_trial_op_dim, num_qp, num_entities);
auto itod = Reshape(inputs_trial_op_dim.Read(), num_inputs);
@@ -1095,6 +1115,7 @@ void DifferentiableOperator::AddIntegrator(
inputs_trial_op_dim,
total_trial_op_dim,
trial_vdim,
qpdc_idx = cache_index,
// capture by ref:
&qpdc_mem = derivative_qp_caches_ref,
&or_transpose
@@ -1112,7 +1133,7 @@ void DifferentiableOperator::AddIntegrator(
shmem_info.direction_size,
num_entities);
auto qpdc = Reshape(qpdc_mem.Read(), test_vdim, test_op_dim,
auto qpdc = Reshape(qpdc_mem[qpdc_idx].Read(), test_vdim, test_op_dim,
trial_vdim, total_trial_op_dim, num_qp, num_entities);
auto itod = Reshape(inputs_trial_op_dim.Read(), num_inputs);
@@ -1157,6 +1178,154 @@ void DifferentiableOperator::AddIntegrator(
or_transpose(derivative_action_e, der_action_l);
});
// This prevents Sum/Identity. These are invalid
// as input FieldOperators anyways.
constexpr auto dummy_fop = Value<0> {};
auto [input_rt,
input_e_sz] = get_restriction_transpose<entity_t>
(fields[d_field_idx],
element_dof_ordering, dummy_fop);
const auto input_restriction_transpose = input_rt;
derivative_tr_prolongation_transpose[derivative_id] =
get_prolongation_transpose(
fields[d_field_idx], dummy_fop, mesh.GetComm());
const auto d_tr_field_idx = test_space_field_idx;
const auto direction_tr = fields[d_tr_field_idx];
auto output_size_on_qp =
get_input_size_on_qp(outputs, std::make_index_sequence<num_outputs> {});
const int residual_tr_size_on_qp = trial_vdim * total_trial_op_dim;
auto shmem_tr_info =
get_shmem_info<entity_t, num_fields, num_outputs, num_inputs>(
output_dtq_maps, input_dtq_maps, fields, num_entities,
num_qp, output_size_on_qp, residual_tr_size_on_qp,
element_dof_ordering, test_space_field_idx);
// print_shared_memory_info(shmem_tr_info);
// TODO: this is a hack to extend the shared memory with a known
// offset for a temp variable
Vector shmem_tr_cache(shmem_tr_info.total_size + residual_tr_size_on_qp *
num_qp);
Vector direction_tr_e(get_restriction<entity_t>(
fields[test_space_field_idx],
element_dof_ordering)->Height());
derivative_action_tr_l_size[derivative_id] =
get_restriction<entity_t>(fields[d_field_idx],
element_dof_ordering)->Width();
Vector derivative_action_tr_e(input_e_sz);
derivative_action_tr_e = 0.0;
derivative_action_tr_callbacks[derivative_id].push_back(
[
// capture by copy:
dimension, // int
num_entities, // int
num_trial_dof, // int
num_qp, // int
q1d, // int
test_vdim, // int (= output_fop.vdim)
test_op_dim, // int (derived from output_fop)
inputs, // mfem::future::tuple
outputs, // mfem::future::tuple
attributes, // Array<int>
ir_weights, // DeviceTensor
use_sum_factorization, // bool
input_dtq_maps, // std::array<DofToQuadMap, num_fields>
output_dtq_maps, // std::array<DofToQuadMap, num_fields>
// output_fop, // class derived from FieldOperator
thread_blocks, // ThreadBlocks
shmem_tr_cache, // Vector (local)
shmem_tr_info, // SharedMemoryInfo
// TODO: make this Array<int> a member of the DifferentiableOperator
// and capture it by ref.
elem_attributes, // Array<int>
input_is_dependent,
direction_tr, // FieldDescriptor
direction_tr_e, // Vector
derivative_action_tr_e, // Vector
element_dof_ordering, // ElementDofOrdering
inputs_trial_op_dim,
total_trial_op_dim,
trial_vdim,
input_restriction_transpose,
qpdc_idx = cache_index,
// capture by ref:
&qpdc_mem = derivative_qp_caches_ref
](
std::vector<Vector> &f_e, const Vector &dir_tr_l,
Vector &derivative_action_tr_l) mutable
{
restriction<entity_t>(direction_tr, dir_tr_l, direction_tr_e,
element_dof_ordering);
auto ye = Reshape(derivative_action_tr_e.ReadWrite(), num_trial_dof,
trial_vdim, num_entities);
auto wrapped_fields_e = wrap_fields(f_e, shmem_tr_info.field_sizes,
num_entities);
auto wrapped_direction_e = Reshape(direction_tr_e.ReadWrite(),
shmem_tr_info.direction_size,
num_entities);
auto qpdc = Reshape(qpdc_mem[qpdc_idx].Read(), test_vdim, test_op_dim,
trial_vdim, total_trial_op_dim, num_qp, num_entities);
auto itod = Reshape(inputs_trial_op_dim.Read(), num_inputs);
const bool has_attr = attributes.Size() > 0;
const auto d_attr = attributes.Read();
const auto d_elem_attr = elem_attributes->Read();
derivative_action_tr_e = 0.0;
forall([=] MFEM_HOST_DEVICE (int e, real_t *shmem)
{
if (has_attr && !d_attr[d_elem_attr[e] - 1]) { return; }
auto [output_dtq_shmem, input_dtq_shmem, fields_shmem,
direction_shmem, input_shmem,
shadow_shmem_, residual_shmem,
scratch_shmem] =
unpack_shmem(shmem, shmem_tr_info, output_dtq_maps, input_dtq_maps,
wrapped_fields_e, wrapped_direction_e, num_qp, e);
auto &shadow_shmem = shadow_shmem_;
std::array<bool, num_outputs> all_true{true};
map_direction_to_quadrature_data_conditional(
shadow_shmem, direction_shmem, output_dtq_shmem, outputs,
ir_weights, scratch_shmem, all_true, dimension,
use_sum_factorization);
auto fhat = Reshape(&residual_shmem(0, 0), trial_vdim,
total_trial_op_dim, num_qp);
auto qpdce = Reshape(&qpdc(0, 0, 0, 0, 0, e), test_vdim, test_op_dim,
trial_vdim, total_trial_op_dim, num_qp);
constexpr bool transpose = true;
apply_qpdc(fhat, shadow_shmem, qpdce, itod, q1d, dimension,
use_sum_factorization, transpose);
auto y = Reshape(&ye(0, 0, e), num_trial_dof, trial_vdim);
auto fi_shmem = Reshape(shmem + shmem_tr_info.total_size, trial_vdim,
total_trial_op_dim, num_qp);
map_quadrature_data_to_fields_conditional(
y, fhat, inputs, itod, input_dtq_shmem, scratch_shmem, fi_shmem,
input_is_dependent, dimension, use_sum_factorization);
}, num_entities, thread_blocks, shmem_tr_info.total_size,
shmem_tr_cache.ReadWrite());
input_restriction_transpose(derivative_action_tr_e, derivative_action_tr_l);
});
assemble_derivative_sparsematrix_callbacks[derivative_id].push_back(
[
// capture by copy:
@@ -1190,7 +1359,7 @@ void DifferentiableOperator::AddIntegrator(
inputs_trial_op_dim,
Ae_mem,
output_to_field,
qpdc_idx = cache_index,
// capture by ref:
&qpdc_mem = derivative_qp_caches_ref,
&fields = fields_ref
@@ -1202,7 +1371,7 @@ void DifferentiableOperator::AddIntegrator(
shmem_info.direction_size,
num_entities);
auto qpdc = Reshape(qpdc_mem.Read(), test_vdim, test_op_dim,
auto qpdc = Reshape(qpdc_mem[qpdc_idx].Read(), test_vdim, test_op_dim,
trial_vdim, total_trial_op_dim, num_qp, num_entities);
auto itod = Reshape(inputs_trial_op_dim.Read(), num_inputs);
@@ -1250,7 +1419,10 @@ void DifferentiableOperator::AddIntegrator(
auto test_fes = *std::get_if<const ParFiniteElementSpace *>
(&fields[output_to_field[0]].data);
A = new SparseMatrix(test_fes->GetVSize(), trial_fes->GetVSize());
if (A == nullptr)
{
A = new SparseMatrix(test_fes->GetVSize(), trial_fes->GetVSize());
}
auto tmp = Reshape(Ae_mem.HostReadWrite(), num_test_dof * test_vdim,
num_trial_dof * trial_vdim, num_entities);
@@ -1321,14 +1493,15 @@ void DifferentiableOperator::AddIntegrator(
A->AddSubMatrix(test_vdofs, trial_vdofs, Aee, 1);
}
}
A->Finalize();
// Don't finalize here since multiple callbacks might contribute to the same matrix
// A->Finalize() will be called after all callbacks have contributed
});
// Create local references for MSVC lambda capture compatibility
auto& assemble_derivative_sparsematrix_callbacks_ref =
this->assemble_derivative_sparsematrix_callbacks[derivative_id];
assemble_derivative_hypreparmatrix_callbacks[derivative_id].push_back(
assemble_derivative_hypreparmatrix_callback[derivative_id] =
[
input_is_dependent,
input_to_field,
@@ -1342,6 +1515,7 @@ void DifferentiableOperator::AddIntegrator(
{
f(f_e, spmat);
}
spmat->Finalize();
if (spmat == nullptr)
{
@@ -1395,7 +1569,7 @@ void DifferentiableOperator::AddIntegrator(
trial_fes->Dof_TrueDof_Matrix());
}
delete spmat;
});
};
}, derivative_ids);
}
}
+69 -11
View File
@@ -90,8 +90,8 @@ void map_quadrature_data_to_fields_impl(
}
else
{
MFEM_ABORT_KERNEL("quadrature data mapping to field is not implemented"
" for this field descriptor");
MFEM_ABORT("quadrature data mapping to field is not implemented for"
" this field descriptor");
}
}
@@ -169,9 +169,8 @@ void map_quadrature_data_to_fields_tensor_impl_1d(
}
else
{
MFEM_ABORT_KERNEL("quadrature data mapping to field is not implemented"
"for this field descriptor with sum factorization on"
" tensor product elements");
MFEM_ABORT("quadrature data mapping to field is not implemented for"
" this field descriptor with sum factorization on tensor product elements");
}
}
@@ -307,9 +306,8 @@ void map_quadrature_data_to_fields_tensor_impl_2d(
}
else
{
MFEM_ABORT_KERNEL("quadrature data mapping to field is not implemented"
" for this field descriptor with sum factorization on"
" tensor product elements");
MFEM_ABORT("quadrature data mapping to field is not implemented for"
" this field descriptor with sum factorization on tensor product elements");
}
}
@@ -494,9 +492,8 @@ void map_quadrature_data_to_fields_tensor_impl_3d(
}
else
{
MFEM_ABORT_KERNEL("quadrature data mapping to field is not implemented"
" for this field descriptor with sum factorization on"
" tensor product elements");
MFEM_ABORT("quadrature data mapping to field is not implemented for"
" this field descriptor with sum factorization on tensor product elements");
}
}
@@ -533,4 +530,65 @@ void map_quadrature_data_to_fields(
}
}
template <size_t N, typename field_operator_ts>
MFEM_HOST_DEVICE
void map_quadrature_data_to_fields_conditional(
DeviceTensor<2, real_t> &y,
const DeviceTensor<3, real_t> &f,
const field_operator_ts &fops,
const DeviceTensor<1, const real_t> &op_dims,
const std::array<DofToQuadMap, N> &dtqmaps,
std::array<DeviceTensor<1>, 6> &scratch_mem,
const DeviceTensor<3> &fi_shmem,
const std::array<bool, N> &conditions,
const int &dimension,
const bool &use_sum_factorization)
{
int offset = 0;
for_constexpr<N>([&](auto i)
{
if (conditions[i])
{
[[maybe_unused]] const auto [K, unused, M] = f.GetShape();
const int L = static_cast<int>(op_dims(static_cast<size_t>(i)));
auto fi = Reshape(&fi_shmem(0, 0, 0), K, L, M);
for (int k = 0; k < K; k++)
{
for (int l = 0; l < L; l++)
{
for (int m = 0; m < M; m++)
{
fi(k, l, m) = f(k, l + offset, m);
}
}
}
if (use_sum_factorization)
{
if (dimension == 1)
{
map_quadrature_data_to_fields_tensor_impl_1d(
y, fi, get<i>(fops), dtqmaps[i], scratch_mem);
}
else if (dimension == 2)
{
map_quadrature_data_to_fields_tensor_impl_2d(
y, fi, get<i>(fops), dtqmaps[i], scratch_mem);
}
else if (dimension == 3)
{
map_quadrature_data_to_fields_tensor_impl_3d(
y, fi, get<i>(fops), dtqmaps[i], scratch_mem);
}
else { MFEM_ABORT_KERNEL("dimension not supported"); }
}
else
{
map_quadrature_data_to_fields_impl(y, fi, get<i>(fops), dtqmaps[i]);
}
offset += L;
}
});
}
} // namespace mfem::future
+12 -17
View File
@@ -505,13 +505,13 @@ void map_field_to_quadrature_data(
}
}
template <typename field_operator_ts, size_t num_inputs, size_t num_fields>
template <typename field_operator_ts, size_t N, size_t M>
MFEM_HOST_DEVICE inline
void map_fields_to_quadrature_data(
std::array<DeviceTensor<2>, num_inputs> &fields_qp,
const std::array<DeviceTensor<1>, num_fields> &fields_e,
const std::array<DofToQuadMap, num_inputs> &dtqmaps,
const std::array<size_t, num_inputs> &input_to_field,
std::array<DeviceTensor<2>, N> &fields_qp,
const std::array<DeviceTensor<1>, M> &fields_e,
const std::array<DofToQuadMap, N> &dtqmaps,
const std::array<size_t, N> &input_to_field,
const field_operator_ts &fops,
const DeviceTensor<1, const real_t> &integration_weights,
const std::array<DeviceTensor<1>, 6> &scratch_mem,
@@ -523,7 +523,7 @@ void map_fields_to_quadrature_data(
// attached to them and we create a dummy field which is not accessed
// inside the functions it is passed to.
const auto dummy_field_weight = DeviceTensor<1>(nullptr, 0);
for_constexpr<num_inputs>([&](auto i)
for_constexpr<N>([&](auto i)
{
const DeviceTensor<1> &field_e =
(input_to_field[i] == SIZE_MAX) ? dummy_field_weight :
@@ -549,12 +549,7 @@ void map_fields_to_quadrature_data(
fields_qp[i], dtqmaps[i], field_e, get<i>(fops),
integration_weights, scratch_mem);
}
else
{
#if !(defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP))
MFEM_ABORT("unsupported dimension");
#endif
}
else { MFEM_ABORT_KERNEL("unsupported dimension"); }
}
else
{
@@ -627,20 +622,20 @@ void map_fields_to_quadrature_data_conditional(
});
}
template <size_t num_inputs, typename field_operator_ts>
template <size_t N, typename field_operator_ts>
MFEM_HOST_DEVICE
void map_direction_to_quadrature_data_conditional(
std::array<DeviceTensor<2>, num_inputs> &directions_qp,
std::array<DeviceTensor<2>, N> &directions_qp,
const DeviceTensor<1> &direction_e,
const std::array<DofToQuadMap, num_inputs> &dtqmaps,
const std::array<DofToQuadMap, N> &dtqmaps,
field_operator_ts fops,
const DeviceTensor<1, const real_t> &integration_weights,
const std::array<DeviceTensor<1>, 6> &scratch_mem,
const std::array<bool, num_inputs> &conditions,
const std::array<bool, N> &conditions,
const int &dimension,
const bool &use_sum_factorization)
{
for_constexpr<num_inputs>([&](auto i)
for_constexpr<N>([&](auto i)
{
if (conditions[i])
{
+1 -1
View File
@@ -82,7 +82,7 @@ protected:
};
/// @brief Uniform parameter space
class UniformParameterSpace : public ParameterSpace
class UniformParameterSpace final : public ParameterSpace
{
public:
/// @brief Constructor for a uniform parameter space
+88 -38
View File
@@ -379,54 +379,102 @@ namespace detail
/// @param shadow_shmem the shadow shared memory.
/// @param qpdc the quadrature point data cache holding the resulting
/// Jacobians on each quadrature point.
/// @param itod inputs trial operator dimension.
/// If input is dependent the value corresponds to the spatial dimension, otherwise
/// a zero indicates non-dependence on the variable.
/// @param op_dims operator dimensions.
/// If an operator is dependent, the value corresponds to the spatial dimension.
/// Otherwise a zero indicates indepence on the variable.
/// @param q the current quadrature point index.
template <size_t num_fields>
/// @param transpose switch to use transpose action.
template <size_t N>
MFEM_HOST_DEVICE inline
void apply_qpdc(
DeviceTensor<3> &fhat,
const std::array<DeviceTensor<2>, num_fields> &shadow_shmem,
const std::array<DeviceTensor<2>, N> &shadow_shmem,
const DeviceTensor<5, const real_t> &qpdc,
const DeviceTensor<1, const real_t> &itod,
const int &q)
const DeviceTensor<1, const real_t> &op_dims,
const int &q,
bool transpose)
{
const size_t num_ops = op_dims.GetShape()[0];
const int test_vdim = qpdc.GetShape()[0];
const int test_op_dim = qpdc.GetShape()[1];
const int trial_vdim = qpdc.GetShape()[2];
const int num_qp = qpdc.GetShape()[4];
const size_t num_inputs = itod.GetShape()[0];
const int total_trial_op_dim = qpdc.GetShape()[3];
for (int i = 0; i < test_vdim; i++)
const int num_qp = qpdc.GetShape()[4];
if (transpose)
{
for (int k = 0; k < test_op_dim; k++)
for (int j = 0; j < trial_vdim; j++)
{
real_t sum = 0.0;
int m_offset = 0;
for (size_t s = 0; s < num_inputs; s++)
for (int m = 0; m < total_trial_op_dim; m++)
{
const int trial_op_dim = static_cast<int>(itod(s));
if (trial_op_dim == 0)
{
continue;
}
const auto d_qp =
Reshape(&(shadow_shmem[s])[0], trial_vdim, trial_op_dim, num_qp);
for (int j = 0; j < trial_vdim; j++)
{
for (int m = 0; m < trial_op_dim; m++)
{
sum += qpdc(i, k, j, m + m_offset, q) * d_qp(j, m, q);
}
}
m_offset += trial_op_dim;
fhat(j, m, q) = 0.0;
}
}
// Since we don't support more than output space right now
// shadow_shmem will always be of size 1.
constexpr int shadow_idx_tr = 0;
auto d_qp = Reshape(&(shadow_shmem[shadow_idx_tr])[0], test_vdim, test_op_dim,
num_qp);
int m_offset = 0;
for (size_t s = 0; s < num_ops; s++)
{
const int trial_op_dim = static_cast<int>(op_dims(s));
if (trial_op_dim == 0) { continue; }
for (int j = 0; j < trial_vdim; j++)
{
for (int m = 0; m < trial_op_dim; m++)
{
real_t sum = 0.0;
for (int i = 0; i < test_vdim; i++)
{
for (int k = 0; k < test_op_dim; k++)
{
const real_t contrib = qpdc(i, k, j, m + m_offset, q) * d_qp(i, k, q);
sum += contrib;
}
}
fhat(j, m + m_offset, q) += sum;
}
}
m_offset += trial_op_dim;
}
}
else
{
for (int i = 0; i < test_vdim; i++)
{
for (int k = 0; k < test_op_dim; k++)
{
real_t sum = 0.0;
int m_offset = 0;
for (size_t s = 0; s < num_ops; s++)
{
const int trial_op_dim = static_cast<int>(op_dims(s));
if (trial_op_dim == 0) { continue; }
const auto d_qp =
Reshape(&(shadow_shmem[s])[0], trial_vdim, trial_op_dim, num_qp);
for (int j = 0; j < trial_vdim; j++)
{
for (int m = 0; m < trial_op_dim; m++)
{
sum += qpdc(i, k, j, m + m_offset, q) * d_qp(j, m, q);
}
}
m_offset += trial_op_dim;
}
fhat(i, k, q) = sum;
}
fhat(i, k, q) = sum;
}
}
}
}
} // namespace detail
/// @brief Apply the quadrature point data cache (qpdc) to a vector
/// (usually a direction).
@@ -445,16 +493,18 @@ void apply_qpdc(
/// @param q1d number of quadrature points in 1D.
/// @param dimension spatial dimension.
/// @param use_sum_factorization whether to use sum factorization.
template <size_t num_fields>
/// @param T switch to use transpose application.
template <size_t N>
MFEM_HOST_DEVICE inline
void apply_qpdc(
DeviceTensor<3> &fhat,
const std::array<DeviceTensor<2>, num_fields> &shadow_shmem,
const std::array<DeviceTensor<2>, N> &shadow_shmem,
const DeviceTensor<5, const real_t> &qpdc,
const DeviceTensor<1, const real_t> &itod,
const int &q1d,
const int &dimension,
const bool &use_sum_factorization)
const bool &use_sum_factorization,
const bool T = false)
{
if (use_sum_factorization)
{
@@ -462,7 +512,7 @@ void apply_qpdc(
{
MFEM_FOREACH_THREAD_DIRECT(q, x, q1d)
{
detail::apply_qpdc(fhat, shadow_shmem, qpdc, itod, q);
detail::apply_qpdc(fhat, shadow_shmem, qpdc, itod, q, T);
}
}
else if (dimension == 2)
@@ -472,7 +522,7 @@ void apply_qpdc(
MFEM_FOREACH_THREAD_DIRECT(qy, y, q1d)
{
const int q = qx + q1d * qy;
detail::apply_qpdc(fhat, shadow_shmem, qpdc, itod, q);
detail::apply_qpdc(fhat, shadow_shmem, qpdc, itod, q, T);
}
}
}
@@ -485,7 +535,7 @@ void apply_qpdc(
MFEM_FOREACH_THREAD_DIRECT(qz, z, q1d)
{
const int q = qx + q1d * (qy + q1d * qz);
detail::apply_qpdc(fhat, shadow_shmem, qpdc, itod, q);
detail::apply_qpdc(fhat, shadow_shmem, qpdc, itod, q, T);
}
}
}
@@ -500,7 +550,7 @@ void apply_qpdc(
const int num_qp = qpdc.GetShape()[4];
MFEM_FOREACH_THREAD_DIRECT(q, x, num_qp)
{
detail::apply_qpdc(fhat, shadow_shmem, qpdc, itod, q);
detail::apply_qpdc(fhat, shadow_shmem, qpdc, itod, q, T);
}
}
}
+29 -16
View File
@@ -243,6 +243,35 @@ void process_qf_arg(
}
}
template <typename T, int n>
MFEM_HOST_DEVICE inline
void process_qf_arg(
const DeviceTensor<1, T> &u,
const DeviceTensor<1, T> &v,
tensor<T, n> &arg)
{
for (int i = 0; i < n; i++)
{
arg(i) = u(i);
}
}
template <typename T, int n, int m>
MFEM_HOST_DEVICE inline
void process_qf_arg(
const DeviceTensor<1, T> &u,
const DeviceTensor<1, T> &v,
tensor<T, n, m> &arg)
{
for (int i = 0; i < m; i++)
{
for (int j = 0; j < n; j++)
{
arg(j, i) = u((i * n) + j);
}
}
}
template <typename arg_type>
MFEM_HOST_DEVICE inline
void process_qf_arg(const DeviceTensor<2> &u, arg_type &arg, int qp)
@@ -327,20 +356,4 @@ void process_qf_result(
}
}
template <typename T, int n, int m>
MFEM_HOST_DEVICE inline
void process_qf_arg(
const DeviceTensor<1, T> &u,
const DeviceTensor<1, T> &v,
tensor<T, n, m> &arg)
{
for (int i = 0; i < m; i++)
{
for (int j = 0; j < n; j++)
{
arg(j, i) = u((i * n) + j);
}
}
}
} // namespace mfem::future
+9 -3
View File
@@ -1202,7 +1202,14 @@ std::function<void(const Vector&, Vector&)> get_prolongation_transpose(
const Operator *P = get_prolongation(f);
auto PT = [=](const Vector &r_local, Vector &y)
{
P->MultTranspose(r_local, y);
if (P)
{
P->MultTranspose(r_local, y);
}
else
{
y = r_local;
}
};
return PT;
}
@@ -1580,14 +1587,13 @@ struct SharedMemoryInfo
std::array<int, 6> temp_sizes;
};
template <typename entity_t, std::size_t num_fields, std::size_t num_inputs, std::size_t num_outputs, typename input_t>
template <typename entity_t, std::size_t num_fields, std::size_t num_inputs, std::size_t num_outputs>
SharedMemoryInfo<num_fields, num_inputs, num_outputs>
get_shmem_info(
const std::array<DofToQuadMap, num_inputs> &input_dtq_maps,
const std::array<DofToQuadMap, num_outputs> &output_dtq_maps,
const std::vector<FieldDescriptor> &fields,
const int &num_entities,
const input_t &inputs,
const int &num_qp,
const std::vector<int> &input_size_on_qp,
const int &residual_size_on_qp,
+5 -82
View File
@@ -1044,50 +1044,9 @@ void VectorFiniteElement::SetDerivMembers()
switch (map_type)
{
case H_DIV:
switch (dim)
{
case 3: // div: 3D H_DIV -> 3D INTEGRAL
deriv_type = DIV;
deriv_range_type = SCALAR;
deriv_map_type = INTEGRAL;
break;
case 2: // div: 2D H_DIV -> 2D INTEGRAL
deriv_type = DIV;
deriv_range_type = SCALAR;
deriv_map_type = INTEGRAL;
break;
default:
MFEM_ABORT("Invalid dimension, Dim = " << dim);
}
break;
case H_DIV_R2D:
switch (dim)
{
case 2: // div: 2D H_DIV_R2D -> 2D INTEGRAL
deriv_type = DIV;
deriv_range_type = SCALAR;
deriv_map_type = INTEGRAL;
break;
case 1: // div: 1D H_DIV_R2D -> 1D INTEGRAL
deriv_type = DIV;
deriv_range_type = SCALAR;
deriv_map_type = INTEGRAL;
break;
default:
MFEM_ABORT("Invalid dimension, Dim = " << dim);
}
break;
case H_DIV_R1D:
switch (dim)
{
case 1: // div: 1D H_DIV_R1D -> 1D INTEGRAL
deriv_type = DIV;
deriv_range_type = SCALAR;
deriv_map_type = INTEGRAL;
break;
default:
MFEM_ABORT("Invalid dimension, Dim = " << dim);
}
deriv_type = DIV;
deriv_range_type = SCALAR;
deriv_map_type = INTEGRAL;
break;
case H_CURL:
switch (dim)
@@ -1105,49 +1064,13 @@ void VectorFiniteElement::SetDerivMembers()
break;
case 1:
deriv_type = NONE;
deriv_range_type = UNKNOWN_RANGE_TYPE;
deriv_map_type = UNKNOWN_MAP_TYPE;
deriv_range_type = SCALAR;
deriv_map_type = INTEGRAL;
break;
default:
MFEM_ABORT("Invalid dimension, Dim = " << dim);
}
break;
case H_CURL_R2D:
switch (dim)
{
case 2:
// curl: 2D H_CURL_R2D -> H_DIV_R2D
deriv_type = CURL;
deriv_range_type = VECTOR;
deriv_map_type = H_DIV_R2D;
break;
case 1:
// curl: 1D H_CURL_R2D -> H_DIV_R2D
deriv_type = CURL;
deriv_range_type = VECTOR;
deriv_map_type = H_DIV_R2D;
break;
default:
MFEM_ABORT("Invalid dimension, Dim = " << dim);
}
break;
case H_CURL_R1D:
switch (dim)
{
case 1:
// curl: 1D H_CURL_R1D -> H_DIV_R1D
deriv_type = CURL;
deriv_range_type = VECTOR;
deriv_map_type = H_DIV_R1D;
break;
case 0:
deriv_type = NONE;
deriv_range_type = UNKNOWN_RANGE_TYPE;
deriv_map_type = UNKNOWN_MAP_TYPE;
default:
MFEM_ABORT("Invalid dimension, Dim = " << dim);
}
break;
default:
MFEM_ABORT("Invalid MapType = " << map_type);
}
+3 -31
View File
@@ -295,20 +295,10 @@ public:
$ u(x) = (1/w) \hat u(\hat x) $ */
H_DIV, /**< For vector fields; preserves surface integrals of the
normal component $ u(x) = (J/w) \hat u(\hat x) $ */
H_CURL, /**< For vector fields; preserves line integrals of the
H_CURL /**< For vector fields; preserves line integrals of the
tangential component
$ u(x) = J^{-t} \hat u(\hat x) $ (square J),
$ u(x) = J(J^t J)^{-1} \hat u(\hat x) $ (general J) */
H_DIV_R2D, /**< For 3-component vector fields in 2D; equivalent to a
direct sum of an H_DIV basis and an INTEGRAL basis */
H_CURL_R2D,/**< For 3-component vector fields in 2D; equivalent to a
direct sum of an H_CURL basis and a VALUE basis */
H_DIV_R1D, /**< For 3-component vector fields in 1D; equivalent to a
direct sum of a VALUE basis and a pair of INTEGRAL
bases */
H_CURL_R1D /**< For 3-component vector fields in 1D; equivalent to a
direct sum of an INTEGRAL basis and a pair of VALUE
bases */
};
/** @brief Enumeration for DerivType: defines which derivative method
@@ -340,28 +330,12 @@ public:
int GetDim() const { return dim; }
/** @brief Returns the vector dimension for vector-valued finite elements,
which is also the dimension of the interpolation operation and the
width of the DenseMatrix argument in
CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape). */
which is also the dimension of the interpolation operation. */
int GetRangeDim() const { return vdim; }
/** @brief Returns the vector dimension, in physical space, for
vector-valued finite elements, which is also the width of the
DenseMatrix argument in
CalcPhysVShape(ElementTransformation &Trans, DenseMatrix &shape). */
virtual int GetPhysRangeDim(int /* space_dim */) const { return vdim; }
/** Returns the dimension of the curl for vector-valued finite elements,
which is also the width of the DenseMatrix argument in
CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape). */
/// Returns the dimension of the curl for vector-valued finite elements.
int GetCurlDim() const { return cdim; }
/** Returns the dimension, in physical space, of the curl for vector-valued
finite elements, which is also the width of the DenseMatrix argument in
CalcPhysCurlShape(ElementTransformation &Trans, DenseMatrix &curl_shape).
*/
virtual int GetPhysCurlDim(int /* space_dim */) const { return cdim; }
/// Returns the Geometry::Type of the reference element.
Geometry::Type GetGeomType() const { return geom_type; }
@@ -1016,8 +990,6 @@ protected:
public:
VectorFiniteElement(int D, Geometry::Type G, int Do, int O, int M,
int F = FunctionSpace::Pk);
int GetPhysRangeDim(int space_dim) const override { return space_dim; }
};
/// @brief Class for computing 1D special polynomials and their associated basis
+4 -4
View File
@@ -2531,7 +2531,7 @@ void ND_FuentesPyramidElement::calcCurlBasis(const int p,
ND_R1D_PointElement::ND_R1D_PointElement(int p)
: VectorFiniteElement(1, Geometry::POINT, 2, p,
H_CURL_R1D, FunctionSpace::Pk)
H_CURL, FunctionSpace::Pk)
{
// VectorFiniteElement::SetDerivMembers doesn't support 0D H_CURL elements
// so we mimic a 1D element and then correct the dimension here.
@@ -2562,7 +2562,7 @@ ND_R1D_SegmentElement::ND_R1D_SegmentElement(const int p,
const int cb_type,
const int ob_type)
: VectorFiniteElement(1, Geometry::SEGMENT, 3 * p + 2, p,
H_CURL_R1D, FunctionSpace::Pk),
H_CURL, FunctionSpace::Pk),
dof2tk(dof),
cbasis1d(poly1d.GetBasis(p, VerifyClosed(cb_type))),
obasis1d(poly1d.GetBasis(p - 1, VerifyOpen(ob_type)))
@@ -2839,7 +2839,7 @@ ND_R2D_SegmentElement::ND_R2D_SegmentElement(const int p,
const int cb_type,
const int ob_type)
: VectorFiniteElement(1, Geometry::SEGMENT, 2 * p + 1, p,
H_CURL_R2D, FunctionSpace::Pk),
H_CURL, FunctionSpace::Pk),
dof2tk(dof),
cbasis1d(poly1d.GetBasis(p, VerifyClosed(cb_type))),
obasis1d(poly1d.GetBasis(p - 1, VerifyOpen(ob_type)))
@@ -3023,7 +3023,7 @@ void ND_R2D_SegmentElement::Project(VectorCoefficient &vc,
ND_R2D_FiniteElement::ND_R2D_FiniteElement(int p, Geometry::Type G, int Do,
const real_t *tk_fe)
: VectorFiniteElement(2, G, Do, p,
H_CURL_R2D, FunctionSpace::Pk),
H_CURL, FunctionSpace::Pk),
tk(tk_fe),
dof_map(dof),
dof2tk(dof)
-6
View File
@@ -663,9 +663,6 @@ public:
const int cb_type = BasisType::GaussLobatto,
const int ob_type = BasisType::GaussLegendre);
int GetPhysRangeDim(int space_dim) const override { return 2; }
int GetPhysCurlDim(int space_dim) const override { return 1; }
void CalcVShape(const IntegrationPoint &ip,
DenseMatrix &shape) const override;
@@ -708,9 +705,6 @@ private:
DenseMatrix &I) const;
public:
int GetPhysRangeDim(int space_dim) const override { return 3; }
int GetPhysCurlDim(int space_dim) const override { return 3; }
using FiniteElement::CalcVShape;
using FiniteElement::CalcPhysCurlShape;
+3 -3
View File
@@ -2006,7 +2006,7 @@ RT_R1D_SegmentElement::RT_R1D_SegmentElement(const int p,
const int cb_type,
const int ob_type)
: VectorFiniteElement(1, Geometry::SEGMENT, 3 * p + 4, p + 1,
H_DIV_R1D, FunctionSpace::Pk),
H_DIV, FunctionSpace::Pk),
dof2nk(dof),
cbasis1d(poly1d.GetBasis(p + 1, VerifyClosed(cb_type))),
obasis1d(poly1d.GetBasis(p, VerifyOpen(ob_type)))
@@ -2281,7 +2281,7 @@ const real_t RT_R2D_SegmentElement::nk[2] = { 0.,1.};
RT_R2D_SegmentElement::RT_R2D_SegmentElement(const int p,
const int ob_type)
: VectorFiniteElement(1, Geometry::SEGMENT, p + 1, p + 1,
H_DIV_R2D, FunctionSpace::Pk),
H_DIV, FunctionSpace::Pk),
dof2nk(dof),
obasis1d(poly1d.GetBasis(p, VerifyOpen(ob_type)))
{
@@ -2392,7 +2392,7 @@ void RT_R2D_SegmentElement::LocalInterpolation(const VectorFiniteElement &cfe,
RT_R2D_FiniteElement::RT_R2D_FiniteElement(int p, Geometry::Type G, int Do,
const real_t *nk_fe)
: VectorFiniteElement(2, G, Do, p + 1,
H_DIV_R2D, FunctionSpace::Pk),
H_DIV, FunctionSpace::Pk),
nk(nk_fe),
dof_map(dof),
dof2nk(dof)
-6
View File
@@ -510,9 +510,6 @@ public:
RT_R2D_SegmentElement(const int p,
const int ob_type = BasisType::GaussLegendre);
int GetPhysRangeDim(int space_dim) const override { return 2; }
int GetPhysCurlDim(int space_dim) const override { return 0; }
void CalcVShape(const IntegrationPoint &ip,
DenseMatrix &shape) const override;
@@ -550,9 +547,6 @@ private:
DenseMatrix &I) const;
public:
int GetPhysRangeDim(int space_dim) const override { return 3; }
int GetPhysCurlDim(int space_dim) const override { return 0; }
using FiniteElement::CalcVShape;
void CalcVShape(ElementTransformation &Trans,
+4
View File
@@ -62,6 +62,10 @@
#include "pnonlinearform.hpp"
#endif
#ifdef MFEM_USE_SIDRE
#include "sidredatacollection.hpp"
#endif
#ifdef MFEM_USE_CONDUIT
#include "conduitdatacollection.hpp"
#endif
-15
View File
@@ -3934,16 +3934,6 @@ const FiniteElement *FiniteElementSpace::GetBE(int i) const
return BE;
}
const FiniteElement *FiniteElementSpace::GetTypicalBE() const
{
if (mesh->GetNBE() > 0) { return GetBE(0); }
Geometry::Type geom = mesh->GetTypicalFaceGeometry();
const FiniteElement *be = fec->FiniteElementForGeometry(geom);
MFEM_VERIFY(be != nullptr, "Could not determine a typical BE!");
return be;
}
const FiniteElement *FiniteElementSpace::GetFaceElement(int i) const
{
MFEM_VERIFY(!IsVariableOrder(), "not implemented");
@@ -3974,11 +3964,6 @@ const FiniteElement *FiniteElementSpace::GetFaceElement(int i) const
return fe;
}
const FiniteElement *FiniteElementSpace::GetTypicalFaceElement() const
{
return fec->FiniteElementForGeometry(mesh->GetTypicalFaceGeometry());
}
const FiniteElement *FiniteElementSpace::GetEdgeElement(int i,
int variant) const
{
+1 -13
View File
@@ -839,7 +839,7 @@ public:
Note: For vector-valued elements, the results pads up the range dimension
to the spatial dimension. E.g., consider a stack of 5 vector-valued
elements each representing 2D vectors, living in a 3 dimensional space.
Then this function would give 15, not 10.
Then this fucntion would give 15, not 10.
*/
int GetVectorDim() const;
@@ -1323,24 +1323,12 @@ public:
associated with i'th boundary face in the mesh object. */
const FiniteElement *GetBE(int i) const;
/// @brief Return a "typical" boundary element.
///
/// This can be used in situations where the local mesh partition may be
/// empty.
const FiniteElement *GetTypicalBE() const;
/** @brief Returns pointer to the FiniteElement in the FiniteElementCollection
associated with i'th face in the mesh object. Faces in this case refer
to the MESHDIM-1 primitive so in 2D they are segments and in 1D they are
points.*/
const FiniteElement *GetFaceElement(int i) const;
/// @brief Return a "typical" face element.
///
/// This can be used in situations where the local mesh partition may be
/// empty.
const FiniteElement *GetTypicalFaceElement() const;
/** @brief Returns pointer to the FiniteElement in the FiniteElementCollection
associated with i'th edge in the mesh object. */
const FiniteElement *GetEdgeElement(int i, int variant = 0) const;
+66 -76
View File
@@ -345,6 +345,27 @@ void GridFunction::ComputeFlux(BilinearFormIntegrator &blfi,
}
}
int GridFunction::VectorDim() const
{
const FiniteElement *fe = fes->GetTypicalFE();
if (!fe || fe->GetRangeType() == FiniteElement::SCALAR)
{
return fes->GetVDim();
}
return fes->GetVDim()*std::max(fes->GetMesh()->SpaceDimension(),
fe->GetRangeDim());
}
int GridFunction::CurlDim() const
{
const FiniteElement *fe = fes->GetTypicalFE();
if (!fe || fe->GetRangeType() == FiniteElement::SCALAR)
{
return 2 * fes->GetMesh()->SpaceDimension() - 3;
}
return fes->GetVDim()*fe->GetCurlDim();
}
void GridFunction::GetTrueDofs(Vector &tv) const
{
const SparseMatrix *R = fes->GetRestrictionMatrix();
@@ -2029,18 +2050,6 @@ void GridFunction::AccumulateAndCountBdrValues(
Coefficient *coeff[], VectorCoefficient *vcoeff, const Array<int> &attr,
Array<int> &values_counter)
{
if (vcoeff)
{
MFEM_VERIFY(fes->GetVDim() == vcoeff->GetVDim(),
"vcoeff vdim != fes VDim");
MFEM_VERIFY(fes->GetTypicalBE()->GetMapType() == FiniteElement::VALUE &&
fes->GetTypicalBE()->GetRangeType() ==
FiniteElement::SCALAR,
"Can only call ProjectBdrCoefficient on scalar value-type "
"boundary elements. "
"Did you intended to call ProjectBdrCoefficientNormal or "
"ProjectBdrCoefficientTangent for vector finite elements?");
}
Array<int> vdofs;
Vector vc;
@@ -2193,9 +2202,6 @@ void GridFunction::AccumulateAndCountBdrTangentValues(
VectorCoefficient &vcoeff, const Array<int> &bdr_attr,
Array<int> &values_counter)
{
MFEM_VERIFY(fes->GetTypicalBE()->GetPhysRangeDim(
fes->GetMesh()->SpaceDimension()) == vcoeff.GetVDim(),
"vcoeff vdim != PhysRangeDim");
const FiniteElement *fe;
ElementTransformation *T;
Array<int> dofs;
@@ -2349,9 +2355,6 @@ void GridFunction::ProjectDeltaCoefficient(DeltaCoefficient &delta_coeff,
void GridFunction::ProjectCoefficient(Coefficient &coeff, ProjectType type)
{
MFEM_VERIFY(
VectorDim() == 1,
"Cannot project scalar Coefficient onto vector GridFunction");
DeltaCoefficient *delta_c = dynamic_cast<DeltaCoefficient *>(&coeff);
DofTransformation doftrans;
Array<int> vdofs;
@@ -2627,7 +2630,6 @@ void GridFunction::ProjectCoefficient(
void GridFunction::ProjectCoefficient(VectorCoefficient &vcoeff,
ProjectType type)
{
MFEM_VERIFY(VectorDim() == vcoeff.GetVDim(), "vcoeff vdim != VectorDim()");
Array<int> vdofs;
Vector vals;
DofTransformation doftrans;
@@ -2943,7 +2945,6 @@ void GridFunction::ProjectCoefficientElementL2(VectorCoefficient &vcoeff)
void GridFunction::ProjectCoefficient(
VectorCoefficient &vcoeff, Array<int> &dofs)
{
MFEM_VERIFY(VectorDim() == vcoeff.GetVDim(), "vcoeff vdim != VectorDim()");
int el = -1;
ElementTransformation *T = NULL;
const FiniteElement *fe = NULL;
@@ -2973,7 +2974,6 @@ void GridFunction::ProjectCoefficient(
void GridFunction::ProjectCoefficient(VectorCoefficient &vcoeff, int attribute)
{
MFEM_VERIFY(VectorDim() == vcoeff.GetVDim(), "vcoeff vdim != VectorDim()");
int i;
Array<int> vdofs;
Vector vals;
@@ -3033,7 +3033,6 @@ void GridFunction::ProjectCoefficient(Coefficient *coeff[])
void GridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff,
Array<int> &dof_attr)
{
MFEM_VERIFY(VectorDim() == coeff.GetVDim(), "coeff vdim != VectorDim()");
Array<int> vdofs;
Vector vals;
@@ -3065,7 +3064,6 @@ void GridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff,
void GridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff)
{
MFEM_VERIFY(VectorDim() == coeff.GetVDim(), "coeff vdim != VectorDim()");
Array<int> dof_attr;
ProjectDiscCoefficient(coeff, dof_attr);
}
@@ -3075,10 +3073,6 @@ void GridFunction::ProjectDiscCoefficient(Coefficient &coeff, AvgType type)
// Harmonic (x1 ... xn) = [ (1/x1 + ... + 1/xn) / n ]^-1.
// Arithmetic(x1 ... xn) = (x1 + ... + xn) / n.
MFEM_VERIFY(
VectorDim() == 1,
"Cannot project a scalar coefficient onto a vector GridFunction");
Array<int> zones_per_vdof;
AccumulateAndCountZones(coeff, type, zones_per_vdof);
@@ -3088,7 +3082,6 @@ void GridFunction::ProjectDiscCoefficient(Coefficient &coeff, AvgType type)
void GridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff,
AvgType type)
{
MFEM_VERIFY(VectorDim() == coeff.GetVDim(), "coeff vdim != VectorDim()");
Array<int> zones_per_vdof;
AccumulateAndCountZones(coeff, type, zones_per_vdof);
@@ -3144,33 +3137,52 @@ void GridFunction::ProjectBdrCoefficient(Coefficient *coeff[],
}
void GridFunction::ProjectBdrCoefficientNormal(
Coefficient *coeff, VectorCoefficient *vcoeff, const Array<int> &bdr_attr)
VectorCoefficient &vcoeff, const Array<int> &bdr_attr)
{
MFEM_VERIFY(fes->GetVDim() == 1, "fespace VDim != 1");
MFEM_VERIFY(fes->GetTypicalBE()->GetRangeType() == FiniteElement::SCALAR &&
fes->GetTypicalBE()->GetMapType() == FiniteElement::INTEGRAL,
"Not an RT FE space!");
if (vcoeff)
{
MFEM_VERIFY(vcoeff->GetVDim() == fes->GetMesh()->SpaceDimension(),
"vcoeff vdim (" << vcoeff->GetVDim()
<< ") != SpaceDimension ("
<< fes->GetMesh()->SpaceDimension() << ")");
}
#if 0
// implementation for the case when the face dofs are integrals of the
// normal component.
const FiniteElement *fe;
ElementTransformation *T;
Array<int> dofs;
int dim = vcoeff.GetVDim();
Vector vc(dim), nor(dim), lvec, shape;
for (int i = 0; i < fes->GetNBE(); i++)
{
if (bdr_attr[fes->GetBdrAttribute(i)-1] == 0)
{
continue;
}
fe = fes->GetBE(i);
T = fes->GetBdrElementTransformation(i);
int intorder = 2*fe->GetOrder(); // !!!
const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(), intorder);
int nd = fe->GetDof();
lvec.SetSize(nd);
shape.SetSize(nd);
lvec = 0.0;
for (int j = 0; j < ir.GetNPoints(); j++)
{
const IntegrationPoint &ip = ir.IntPoint(j);
T->SetIntPoint(&ip);
vcoeff.Eval(vc, *T, ip);
CalcOrtho(T->Jacobian(), nor);
fe->CalcShape(ip, shape);
lvec.Add(ip.weight * (vc * nor), shape);
}
fes->GetBdrElementDofs(i, dofs);
SetSubVector(dofs, lvec);
}
#else
// implementation for the case when the face dofs are scaled point
// values of the normal component.
const FiniteElement *fe;
ElementTransformation *T;
Array<int> dofs;
Vector vc, nor, lvec;
int dim = vcoeff.GetVDim();
Vector vc(dim), nor(dim), lvec;
DofTransformation doftrans;
if (vcoeff)
{
const int dim = vcoeff->GetVDim();
vc.SetSize(dim);
nor.SetSize(dim);
}
for (int i = 0; i < fes->GetNBE(); i++)
{
@@ -3186,22 +3198,15 @@ void GridFunction::ProjectBdrCoefficientNormal(
{
const IntegrationPoint &ip = ir.IntPoint(j);
T->SetIntPoint(&ip);
if (coeff)
{
const real_t c = coeff->Eval(*T, ip);
lvec(j) = c * T->Weight();
}
else if (vcoeff)
{
vcoeff->Eval(vc, *T, ip);
CalcOrtho(T->Jacobian(), nor);
lvec(j) = (vc * nor);
}
vcoeff.Eval(vc, *T, ip);
CalcOrtho(T->Jacobian(), nor);
lvec(j) = (vc * nor);
}
fes->GetBdrElementDofs(i, dofs, doftrans);
doftrans.TransformPrimal(lvec);
SetSubVector(dofs, lvec);
}
#endif
}
void GridFunction::ProjectBdrCoefficientTangent(
@@ -5002,14 +5007,6 @@ real_t ExtrudeCoefficient::Eval(ElementTransformation &T,
return sol_in.Eval(*T_in, ip);
}
void VectorExtrudeCoefficient::Eval(Vector &v, ElementTransformation &T,
const IntegrationPoint &ip)
{
ElementTransformation *T_in =
mesh_in->GetElementTransformation(T.ElementNo / n);
T_in->SetIntPoint(&ip);
sol_in.Eval(v, *T_in, ip);
}
GridFunction *Extrude1DGridFunction(Mesh *mesh, Mesh *mesh2d,
GridFunction *sol, const int ny)
@@ -5060,17 +5057,10 @@ GridFunction *Extrude1DGridFunction(Mesh *mesh, Mesh *mesh2d,
return NULL;
}
FiniteElementSpace *solfes2d;
const int vdim = sol->FESpace()->GetVDim();
solfes2d = new FiniteElementSpace(mesh2d, solfec2d, vdim);
// assuming sol is scalar
solfes2d = new FiniteElementSpace(mesh2d, solfec2d);
sol2d = new GridFunction(solfes2d);
sol2d->MakeOwner(solfec2d);
if (vdim > 1)
{
VectorGridFunctionCoefficient vcsol(sol);
VectorExtrudeCoefficient vc2d(mesh, vcsol, ny);
sol2d->ProjectCoefficient(vc2d);
}
else
{
GridFunctionCoefficient csol(sol);
ExtrudeCoefficient c2d(mesh, csol, ny);
@@ -5768,4 +5758,4 @@ std::pair<real_t, real_t> GridFunction::EstimateFunctionMaximum(
return std::make_pair(global_max_lower, global_max_upper);
}
}
}
+13 -69
View File
@@ -150,13 +150,11 @@ public:
FiniteElementCollection *OwnFEC() { return fec_owned; }
/** @brief Shortcut for calling FiniteElementSpace::GetVectorDim() on the
underlying #fes */
int VectorDim() const { return fes->GetVectorDim(); }
/// Shortcut for calling FiniteElementSpace::GetVectorDim() on the underlying #fes
int VectorDim() const;
/** @brief Shortcut for calling FiniteElementSpace::GetCurlDim() on the
underlying #fes */
int CurlDim() const { return fes->GetCurlDim(); }
/// Shortcut for calling FiniteElementSpace::GetCurlDim() on the underlying #fes
int CurlDim() const;
/// Read only access to the (optional) internal true-dof Vector.
const Vector &GetTrueVector() const
@@ -534,9 +532,6 @@ public:
std::unique_ptr<GridFunction> ProlongateToMaxOrder() const;
protected:
void ProjectBdrCoefficientNormal(Coefficient *coeff, VectorCoefficient *vcoeff,
const Array<int> &attr);
/** @brief Accumulates (depending on @a type) the values of @a coeff at all
shared vdofs and counts in how many zones each vdof appears. */
void AccumulateAndCountZones(Coefficient &coeff, AvgType type,
@@ -661,26 +656,15 @@ public:
virtual void ProjectBdrCoefficient(Coefficient *coeff[],
const Array<int> &attr);
/** @brief Project the normal component of the given VectorCoefficient on
the boundary. */
/** Only boundary attributes that are marked in @a bdr_attr are
projected. Assumes RT-type vector finite element GridFunction. */
/** Project the normal component of the given VectorCoefficient on
the boundary. Only boundary attributes that are marked in
'bdr_attr' are projected. Assumes RT-type VectorFE GridFunction. */
void ProjectBdrCoefficientNormal(VectorCoefficient &vcoeff,
const Array<int> &bdr_attr)
{ ProjectBdrCoefficientNormal(NULL, &vcoeff, bdr_attr); }
/** @brief Project the given Coefficient in the normal direction on the
boundary. */
/** Only boundary attributes that are marked in @a bdr_attr are projected.
Assumes RT-type vector finite element GridFunction. */
void ProjectBdrCoefficientNormal(Coefficient &coeff,
const Array<int> &bdr_attr)
{ ProjectBdrCoefficientNormal(&coeff, NULL, bdr_attr); }
const Array<int> &bdr_attr);
/** @brief Project the tangential components of the given VectorCoefficient
on the boundary. */
/** Only boundary attributes that are marked in @a bdr_attr
are projected. Assumes ND-type vector finite element GridFunction. */
on the boundary. Only boundary attributes that are marked in @a bdr_attr
are projected. Assumes ND-type VectorFE GridFunction. */
virtual void ProjectBdrCoefficientTangent(VectorCoefficient &vcoeff,
const Array<int> &bdr_attr);
@@ -1930,7 +1914,7 @@ real_t ComputeElementLpDistance(real_t p, int i,
GridFunction& gf1, GridFunction& gf2);
/// Class used for extruding a scalar coefficient
/// Class used for extruding scalar GridFunctions
class ExtrudeCoefficient : public Coefficient
{
private:
@@ -1938,53 +1922,13 @@ private:
Mesh *mesh_in;
Coefficient &sol_in;
public:
/// Constructs an instance of VectorExtrudeCoefficient
/**
* @param m 1D mesh
* @param s 1D vector coefficient
* @param n_ number of transverse elements of the extruded mesh
*/
ExtrudeCoefficient(Mesh *m, Coefficient &s, int n_)
: n(n_), mesh_in(m), sol_in(s)
{ MFEM_VERIFY(n > 0, "Number of transverse elements must be positive!"); }
: n(n_), mesh_in(m), sol_in(s) { }
real_t Eval(ElementTransformation &T, const IntegrationPoint &ip) override;
virtual ~ExtrudeCoefficient() { }
};
/// Class used for extruding a vector coefficient
class VectorExtrudeCoefficient : public VectorCoefficient
{
private:
int n;
Mesh *mesh_in;
VectorCoefficient &sol_in;
public:
/// Constructs an instance of VectorExtrudeCoefficient
/**
* @param m 1D mesh
* @param s 1D vector coefficient
* @param n_ number of transverse elements of the extruded mesh
*/
VectorExtrudeCoefficient(Mesh *m, VectorCoefficient &s, int n_)
: VectorCoefficient(s.GetVDim()), n(n_), mesh_in(m), sol_in(s)
{ MFEM_VERIFY(n > 0, "Number of transverse elements must be positive!"); }
void Eval(Vector &v, ElementTransformation &T,
const IntegrationPoint &ip) override;
using VectorCoefficient::Eval;
virtual ~VectorExtrudeCoefficient() { }
};
/// Extrude a 1D GridFunction, after extruding the mesh with Extrude1D()
/**
* @param mesh 1D mesh
* @param mesh2d extruded mesh
* @param sol grid function
* @param ny number of transverse elements of the extruded mesh
*/
/// Extrude a scalar 1D GridFunction, after extruding the mesh with Extrude1D.
GridFunction *Extrude1DGridFunction(Mesh *mesh, Mesh *mesh2d,
GridFunction *sol, const int ny);
+8 -18
View File
@@ -197,21 +197,15 @@ static void EAHdivAssemble3D(const int NE,
// Assemble (one row per thread)
MFEM_FOREACH_THREAD(idx_i, x, NDOF)
{
// NOTE: due to an llvm backend bug, usage of the modulus operator
// has been removed from this foreach section.
const int ic = idx_i / NDOF_C;
const int idx_ii = idx_i - ic * NDOF_C; // idx_i % NDOF_C
const int idx_ii = idx_i % NDOF_C;
const int nx_i = (ic == 0) ? D1D : D1D-1;
const int ny_i = (ic == 1) ? D1D : D1D-1;
const int qx_i = idx_ii / nx_i;
const int ix = idx_ii - qx_i * nx_i; // idx_ii % nx_i
const int qy_i = qx_i / ny_i;
const int iy = qx_i - qy_i * ny_i; // (idx_ii / nx_i) % ny_i
const int iz = qy_i; // (idx_ii / nx_i) / ny_i
const int ix = idx_ii % nx_i;
const int iy = (idx_ii / nx_i) % ny_i;
const int iz = (idx_ii / nx_i) / ny_i;
const real_t (&Bi1)[MQ1][MD1] = (ic == 0) ? r_Bc : r_Bo;
const real_t (&Bi2)[MQ1][MD1] = (ic == 1) ? r_Bc : r_Bo;
@@ -220,18 +214,14 @@ static void EAHdivAssemble3D(const int NE,
for (int idx_j = 0; idx_j < NDOF; ++idx_j)
{
const int jc = idx_j / NDOF_C;
const int idx_jj = idx_j - jc * NDOF_C; // idx_j % NDOF_C
const int idx_jj = idx_j % NDOF_C;
const int nx_j = (jc == 0) ? D1D : D1D-1;
const int ny_j = (jc == 1) ? D1D : D1D-1;
const int qx_j = idx_jj / nx_j;
const int jx = idx_jj - qx_j * nx_j; // idx_jj % nx_j
const int qy_j = qx_j / ny_j;
const int jy = qx_j - qy_j * ny_j; // (idx_jj / nx_j) % ny_j
const int jz = qy_j; // (idx_jj / nx_j) / ny_j
const int jx = idx_jj % nx_j;
const int jy = (idx_jj / nx_j) % ny_j;
const int jz = (idx_jj / nx_j) / ny_j;
const real_t (&Bj1)[MQ1][MD1] = (jc == 0) ? r_Bc : r_Bo;
const real_t (&Bj2)[MQ1][MD1] = (jc == 1) ? r_Bc : r_Bo;
+327 -811
View File
File diff suppressed because it is too large Load Diff
+27 -30
View File
@@ -125,6 +125,18 @@ private:
void AddTriPoints3b(const int off, const real_t b, const real_t weight)
{ AddTriPoints3(off, (1. - b)/2., b, weight); }
void AddTriPoints3R(const int off, const real_t a, const real_t b,
const real_t c, const real_t weight)
{
IntPoint(off + 0).Set2w(a, b, weight);
IntPoint(off + 1).Set2w(c, a, weight);
IntPoint(off + 2).Set2w(b, c, weight);
}
void AddTriPoints3R(const int off, const real_t a, const real_t b,
const real_t weight)
{ AddTriPoints3R(off, a, b, 1. - a - b, weight); }
void AddTriPoints6(const int off, const real_t a, const real_t b,
const real_t c, const real_t weight)
{
@@ -171,6 +183,14 @@ private:
AddTetPoints3(off + 1, a, 1. - 3.*a, weight);
}
// given b, add the permutations of (a,a,a,b), where 3*a + b = 1
void AddTetPoints4b(const int off, const real_t b, const real_t weight)
{
const real_t a = (1. - b)/3.;
IntPoint(off).Set(a, a, a, weight);
AddTetPoints3(off + 1, a, b, weight);
}
// add the permutations of (a,a,b,b), 2*(a + b) = 1
void AddTetPoints6(const int off, const real_t a, const real_t weight)
{
@@ -189,37 +209,14 @@ private:
AddTetPoints6(off + 6, a, bc, cb, weight);
}
// add all 24 permutations of (a,b,c,d) where a+b+c+d = 1, all distinct
void AddTetPoints24(const int off, const real_t a, const real_t b,
const real_t c, const real_t weight)
// given (b,c), add the permutations of (a,a,b,c), 2*a + b + c = 1
void AddTetPoints12bc(const int off, const real_t b, const real_t c,
const real_t weight)
{
const real_t d = 1. - a - b - c;
// all 24 permutations of 4 distinct barycentric coordinates
// permuting which coordinate goes to x, y, z (4th is 1-x-y-z)
IntPoint(off + 0).Set(a, b, c, weight);
IntPoint(off + 1).Set(a, b, d, weight);
IntPoint(off + 2).Set(a, c, b, weight);
IntPoint(off + 3).Set(a, c, d, weight);
IntPoint(off + 4).Set(a, d, b, weight);
IntPoint(off + 5).Set(a, d, c, weight);
IntPoint(off + 6).Set(b, a, c, weight);
IntPoint(off + 7).Set(b, a, d, weight);
IntPoint(off + 8).Set(b, c, a, weight);
IntPoint(off + 9).Set(b, c, d, weight);
IntPoint(off + 10).Set(b, d, a, weight);
IntPoint(off + 11).Set(b, d, c, weight);
IntPoint(off + 12).Set(c, a, b, weight);
IntPoint(off + 13).Set(c, a, d, weight);
IntPoint(off + 14).Set(c, b, a, weight);
IntPoint(off + 15).Set(c, b, d, weight);
IntPoint(off + 16).Set(c, d, a, weight);
IntPoint(off + 17).Set(c, d, b, weight);
IntPoint(off + 18).Set(d, a, b, weight);
IntPoint(off + 19).Set(d, a, c, weight);
IntPoint(off + 20).Set(d, b, a, weight);
IntPoint(off + 21).Set(d, b, c, weight);
IntPoint(off + 22).Set(d, c, a, weight);
IntPoint(off + 23).Set(d, c, b, weight);
const real_t a = (1. - b - c)/2.;
AddTetPoints3(off, a, b, weight);
AddTetPoints3(off + 3, a, c, weight);
AddTetPoints6(off + 6, a, b, c, weight);
}
public:
+1 -3
View File
@@ -297,8 +297,7 @@ void LinearForm::Assemble()
tr = mesh->GetBdrFaceTransformations(i);
if (tr != NULL)
{
mfem::DofTransformation doftrans;
fes -> GetElementVDofs (tr -> Elem1No, vdofs, doftrans);
fes -> GetElementVDofs (tr -> Elem1No, vdofs);
for (int k = 0; k < boundary_face_integs.Size(); k++)
{
if (boundary_face_integs_marker[k] &&
@@ -308,7 +307,6 @@ void LinearForm::Assemble()
boundary_face_integs[k]->
AssembleRHSElementVect(*fes->GetFE(tr->Elem1No),
*tr, elemvect);
doftrans.TransformDual(elemvect);
AddElementVector (vdofs, elemvect);
}
}
-14
View File
@@ -545,8 +545,6 @@ void ParGridFunction::GetElementDofValues(int el, Vector &dof_vals) const
void ParGridFunction::ProjectCoefficient(Coefficient &coeff, ProjectType type)
{
MFEM_VERIFY(VectorDim() == 1,
"Cannot project scalar coefficient onto vector ParGridFunction");
DeltaCoefficient *delta_c = dynamic_cast<DeltaCoefficient *>(&coeff);
if (delta_c == NULL)
@@ -719,7 +717,6 @@ void ParGridFunction::ProjectCoefficientElementL2(VectorCoefficient &vcoeff)
void ParGridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff)
{
MFEM_VERIFY(VectorDim() == coeff.GetVDim(), "coeff vdim != VectorDim()");
// local maximal element attribute for each dof
Array<int> ldof_attr;
@@ -764,9 +761,6 @@ void ParGridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff)
void ParGridFunction::ProjectDiscCoefficient(Coefficient &coeff, AvgType type)
{
MFEM_VERIFY(
VectorDim() == 1,
"Cannot project scalar coefficient onto a vector ParGridFunction");
// Harmonic (x1 ... xn) = [ (1/x1 + ... + 1/xn) / n ]^-1.
// Arithmetic(x1 ... xn) = (x1 + ... + xn) / n.
@@ -792,8 +786,6 @@ void ParGridFunction::ProjectDiscCoefficient(VectorCoefficient &vcoeff,
// Harmonic (x1 ... xn) = [ (1/x1 + ... + 1/xn) / n ]^-1.
// Arithmetic(x1 ... xn) = (x1 + ... + xn) / n.
MFEM_VERIFY(VectorDim() == vcoeff.GetVDim(), "vcoeff vdim != VectorDim()");
// Number of zones that contain a given dof.
Array<int> zones_per_vdof;
AccumulateAndCountZones(vcoeff, type, zones_per_vdof);
@@ -866,12 +858,6 @@ void ParGridFunction::ProjectBdrCoefficient(
#endif
}
void ParGridFunction::ProjectBdrCoefficient(VectorCoefficient &vcoeff,
const Array<int> &attr)
{
ProjectBdrCoefficient(NULL, &vcoeff, attr);
}
void ParGridFunction::ProjectBdrCoefficientTangent(VectorCoefficient &vcoeff,
const Array<int> &bdr_attr)
{
+2 -1
View File
@@ -280,7 +280,8 @@ public:
using GridFunction::ProjectBdrCoefficient;
void ProjectBdrCoefficient(VectorCoefficient &vcoeff,
const Array<int> &attr) override;
const Array<int> &attr) override
{ ProjectBdrCoefficient(NULL, &vcoeff, attr); }
void ProjectBdrCoefficient(Coefficient *coeff[],
const Array<int> &attr) override
+5 -11
View File
@@ -321,17 +321,12 @@ void ParL2FaceRestriction::DoubleValuedConformingMult(
const int vd = vdim;
const bool t = byvdim;
const int threshold = ndofs;
const int nsdofs = pfes.GetFaceNbrVSize() / vd;
const int nsdofs = pfes.GetFaceNbrVSize();
auto d_indices1 = scatter_indices1.Read();
auto d_indices2 = scatter_indices2.Read();
auto d_x = Reshape(x.Read(), t?vd:ndofs, t?ndofs:vd);
const int ne_shared = nsdofs / elem_dofs;
const int nedof = elem_dofs;
// Note: the shape of face_nbr_data, as determined by
// ParFiniteElementSpace::ExchangeFaceNbrData, is (elem_dofs, vdim,
// ne_shared), independent of the ordering (byNODES or byVDIM) of the finite
// element space.
auto d_x_shared = Reshape(face_nbr_data.Read(), elem_dofs, vd, ne_shared);
auto d_x_shared = Reshape(face_nbr_data.Read(),
t?vd:nsdofs, t?nsdofs:vd);
auto d_y = Reshape(y.Write(), nface_dofs, vd, 2, nf);
mfem::forall(nfdofs, [=] MFEM_HOST_DEVICE (int i)
{
@@ -351,9 +346,8 @@ void ParL2FaceRestriction::DoubleValuedConformingMult(
}
else if (idx2>=threshold) // shared boundary
{
const int e_shared = (idx2 - threshold) / nedof;
const int i_shared = (idx2 - threshold) % nedof;
d_y(dof, c, 1, face) = d_x_shared(i_shared,c,e_shared);
d_y(dof, c, 1, face) = d_x_shared(t?c:(idx2-threshold),
t?(idx2-threshold):c);
}
else // true boundary
{
+6 -3
View File
@@ -1398,17 +1398,20 @@ void L2FaceRestriction::PermuteAndSetSharedFaceDofsScatterIndices2(
const int dim = fes.GetMesh()->Dimension();
const int dof1d = fes.GetTypicalFE()->GetOrder()+1;
fes.GetTypicalFE()->GetFaceMap(face_id2, face_map);
Array<int> face_nbr_dofs;
const ParFiniteElementSpace &pfes =
static_cast<const ParFiniteElementSpace&>(this->fes);
pfes.GetFaceNbrElementVDofs(elem_index, face_nbr_dofs);
for (int face_dof_elem1 = 0; face_dof_elem1 < face_dofs; ++face_dof_elem1)
{
const int face_dof_elem2 = PermuteFaceL2(dim, face_id1, face_id2,
orientation, dof1d, face_dof_elem1);
const int volume_dof_elem2 = face_map[face_dof_elem2];
// Encode the volume DOF index and element index
const int global_dof_elem2 = elem_index*elem_dofs + volume_dof_elem2;
const int global_dof_elem2 = face_nbr_dofs[volume_dof_elem2];
const int restriction_dof_elem2 = face_dofs*face_index + face_dof_elem1;
// Trick to differentiate dof location inter/shared
scatter_indices2[restriction_dof_elem2] = ndofs + global_dof_elem2;
scatter_indices2[restriction_dof_elem2] = ndofs+global_dof_elem2;
}
#endif
}
File diff suppressed because it is too large Load Diff
+539
View File
@@ -0,0 +1,539 @@
// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
// LICENSE and NOTICE for details. LLNL-CODE-806117.
//
// This file is part of the MFEM library. For more information and source code
// availability visit https://mfem.org.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the BSD-3 license. We welcome feedback and contributions, see file
// CONTRIBUTING.md for details.
#ifndef MFEM_SIDREDATACOLLECTION
#define MFEM_SIDREDATACOLLECTION
#include "../config/config.hpp"
#ifdef MFEM_USE_SIDRE
#include "datacollection.hpp"
// Ignore warnings from the axom/sidre header (GCC + Clang versions)
#ifdef MFEM_HAVE_GCC_PRAGMA_DIAGNOSTIC
# pragma GCC diagnostic push
# if defined(__clang__)
# pragma GCC diagnostic ignored "-Wextra-semi"
# else // real GCC?
# pragma GCC diagnostic ignored "-Wpedantic"
# endif
#endif
#include <axom/sidre.hpp>
#ifdef MFEM_HAVE_GCC_PRAGMA_DIAGNOSTIC
# pragma GCC diagnostic pop
#endif
namespace mfem
{
/** @brief Data collection with Sidre routines following the Conduit mesh
blueprint specification. */
/** SidreDataCollection provides an HDF5-based file format for visualization or
restart capability. This functionality is aimed primarily at customers of
LLNL's axom project that run problems at extreme scales.
For more information, see:
- Sidre component of LLNL's axom project (to be open-sourced), http://goo.gl/cZyJdn
- LLNL conduit/blueprint library, https://github.com/LLNL/conduit
- HDF5 library, https://support.hdfgroup.org/HDF5
The layout created in the Sidre DataStore is: (`"──"` denote groups,
`"─•"` denote views, `"─>"` denote links, i.e. shallow-copy view)
<root>
<collection-name>_global (global group)
blueprint_index
<collection-name> (bp_index group)
state
cycle
time
number_of_domains = <mesh-mpi-comm-size>
coordsets
coords
path = "<bp-path>/coordsets/coords"
type > <bp-grp>/coordsets/coords/type = "explicit"
coord_system = "x"|"xy"|"xyz"
topologies
mesh
path = "<bp-path>/topologies/mesh"
type > <bp-grp>/topologies/mesh/type = "unstructured"
coordset > <bp-grp>/topologies/mesh/coordset = "coords"
grid_function > <bp-grp>/topologies/mesh/grid_function = "<nodes-field-name>"
boundary_topology > <bp-grp>/topologies/mesh/boundary_topology = "boundary"
boundary
path = "<bp-path>/topologies/mesh"
type > <bp-grp>/topologies/boundary/type = "unstructured"
coordset > <bp-grp>/topologies/boundary/coordset = "coords"
fields
mesh_material_attribute
path = "<bp-path>/fields/mesh_material_attribute"
association > <bp-grp>/fields/mesh_material_attribute/association = "element"
topology > <bp-grp>/fields/mesh_material_attribute/topology = "mesh"
number_of_components = 1
boundary_material_attribute
path = "<bp-path>/fields/boundary_material_attribute"
association > <bp-grp>/fields/boundary_material_attribute/association = "element"
topology > <bp-grp>/fields/boundary_material_attribute/topology = "boundary"
number_of_components = 1
grid-function-1
path = "<bp-path>/fields/grid-function-1"
basis > <bp-grp>/fields/grid-function-1/basis = "<fe-coll-name>"
topology > <bp-grp>/fields/grid-function-1/topology = "mesh"
number_of_components = gf1->VectorDim()
grid-function-2
path = "<bp-path>/fields/grid-function-2"
basis > <bp-grp>/fields/grid-function-2/basis = "<fe-coll-name>"
topology > <bp-grp>/fields/grid-function-2/topology = "mesh"
number_of_components = gf2->VectorDim()
...
...
<collection-name> (domain group)
blueprint (blueprint group)
state
cycle
time
domain = <mesh-mpi-rank>
time_step
coordsets
coords
type = "explicit"
values
x = view in <vertex-coords-buffer>/<ext-double-data>
y = view in <vertex-coords-buffer>/<ext-double-data>
z = view in <vertex-coords-buffer>/<ext-double-data>
topologies
mesh
type = "unstructured"
elements
shape = "points"|"lines"|...
connectivity = <vert-idx-array>
coordset = "coords"
grid_function = "<nodes-field-name>"
boundary_topology = "boundary"
boundary
type = "unstructured"
elements
shape = "points"|"lines"|...
connectivity = <vert-idx-array>
coordset = "coords"
fields
mesh_material_attribute
association = "element"
topology = "mesh"
values = <attr-array>
boundary_material_attribute
association = "element"
topology = "boundary"
values = <attr-array>
grid-function-1 (name can include path)
basis = "<fe-coll-name>"
topology = "mesh"
values = <ext-double-array>/<named-buffer> (vdim == 1)
grid-function-2 (name can include path)
basis = "<fe-coll-name>"
topology = "mesh"
values (vdim > 1)
x0 = view into <ext-double-array>/<named-buffer>
x1 = view into <ext-double-array>/<named-buffer>
x2 = view into <ext-double-array>/<named-buffer>
...
...
named_buffers (named_buffers group)
vertex_coords = <double-array>
grid-function-1 = <double-array>
grid-function-2 = <double-array>
...
@note blueprint_index is used both in serial and in parallel. In parallel,
only rank 0 will add entries to the blueprint index.
@note QuadratureFunction%s (q-fields) are not supported.
@note SidreDataCollection does not manage the FiniteElementSpace%s and
FiniteElementCollection%s associated with registered GridFunction%s.
Therefore, field registration is left to the user of SidreDataCollection and
there are no methods that automatically register GridFunction%s using just
the content of the Sidre DataStore. Such capabilities can be implemented in
a derived class, adding any desired object management routines.
@warning This class is still _experimental_, meaning that in future
releases, it may not be backward compatible, and the output files generated
by the current version may become unreadable.
*/
class SidreDataCollection : public DataCollection
{
public:
typedef NamedFieldsMap< Array<int> > AttributeFieldMap;
AttributeFieldMap attr_map;
public:
/// Constructor that allocates and initializes a Sidre DataStore.
/**
@param[in] collection_name Name of the collection used as a file name
when saving
@param[in] the_mesh Mesh shared by all grid functions in the
collection (can be NULL)
@param[in] owns_mesh_data Does the SidreDC own the mesh vertices?
With this constructor, the SidreDataCollection owns the allocated Sidre
DataStore.
*/
SidreDataCollection(const std::string& collection_name,
Mesh *the_mesh = NULL,
bool owns_mesh_data = false);
/// Constructor that links to an external Sidre DataStore.
/** Specifically, the global and domain groups can be at arbitrary paths.
@param[in] collection_name Name of the collection used as a file name
when saving
@param[in] bp_index_grp Pointer to the blueprint index group in the
datastore, see the above schematic
@param[in] domain_grp Pointer to the domain group in the datastore,
see the above schematic
@param[in] owns_mesh_data Does the SidreDC own the mesh vertices?
With this constructor, the SidreDataCollection does not own the Sidre
DataStore.
@note No mesh or fields are read from the given Groups. The mesh has
to be set with SetMesh() and fields registered with RegisterField().
*/
SidreDataCollection(const std::string& collection_name,
axom::sidre::Group * bp_index_grp,
axom::sidre::Group * domain_grp,
bool owns_mesh_data = false);
#ifdef MFEM_USE_MPI
/// Associate an MPI communicator with the collection.
/** If no mesh was associated with the collection, this method should be
called before using any of the Load() methods to read parallel data. */
void SetComm(MPI_Comm comm);
#endif
/// Register a GridFunction in the Sidre DataStore.
/** This method is a shortcut for the call
`RegisterField(field_name, gf, field_name, 0)`.
*/
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
{
RegisterField(field_name, gf, field_name, 0);
}
/// Register a GridFunction in the Sidre DataStore.
/** The registration procedure is as follows:
- if (@a gf's data is NULL), allocate named buffer with the name
@a buffer_name with size _offset + gf->FESpace()->GetVSize()_ and use
its data (plus the given @a offset) to set @a gf's data;
- else, if (DataStore has a named buffer @a buffer_name), replace @a gf's
data array with that named buffer plus the given @a offset;
- else, use @a gf's data as external data associated with @a field_name
in the DataStore;
- register @a field_name in #field_map.
Both the @a field_name and @a buffer_name can contain a path prefix.
@note If @a field_name or @a buffer_name is empty, the method does
nothing.
@note If the GridFunction pointer @a gf or it's FiniteElementSpace
pointer are NULL, the method does nothing.
*/
void RegisterField(const std::string &field_name, GridFunction *gf,
const std::string &buffer_name,
axom::sidre::IndexType offset);
/// Registers an attribute field in the Sidre DataStore
/** The registration process is similar to that of RegisterField()
The attribute field is associated with the elements of the mesh
when @a is_bdry is false, and with the boundary elements, when
@a is_bdry is true.
@sa RegisterField() */
void RegisterAttributeField(const std::string& name, bool is_bdry);
void DeregisterAttributeField(const std::string& name);
/** Returns a pointer to the attribute field associated with
@a field_name, or NULL when there is no associated field */
Array<int>* GetAttributeField(const std::string& field_name) const
{ return attr_map.Get(field_name); }
/** Checks if there is an attribute field associated with @a field_name */
bool HasAttributeField(const std::string& field_name) const
{ return attr_map.Has(field_name); }
/** Checks if any rank in the mesh has boundary elements */
bool HasBoundaryMesh() const;
/// Set the name of the mesh nodes field.
/** This name will be used by SetMesh() to register the mesh nodes, if not
already registered. Also, this method should be called if the mesh nodes
GridFunction was or will be registered directly by the user. The default
value for the name is "mesh_nodes". */
void SetMeshNodesName(const std::string &nodes_name)
{
if (!nodes_name.empty()) { m_meshNodesGFName = nodes_name; }
}
/// De-register @a field_name from the SidreDataCollection.
/** The field is removed from the #field_map and the DataStore, including
deleting it from the named_buffers group, if allocated. */
virtual void DeregisterField(const std::string& field_name);
/// Delete all owned data.
virtual ~SidreDataCollection();
/// Set/change the mesh associated with the collection
/** Uses the field name "mesh_nodes" or the value set by SetMeshNodesName()
to register the mesh nodes GridFunction, if the mesh uses nodes. */
virtual void SetMesh(Mesh *new_mesh);
#ifdef MFEM_USE_MPI
/// Set/change the mesh associated with the collection
/** Uses the field name "mesh_nodes" or the value set by SetMeshNodesName()
to register the mesh nodes GridFunction, if the mesh uses nodes. */
virtual void SetMesh(MPI_Comm comm, Mesh *new_mesh);
#endif
/// Reset the domain and global datastore group pointers.
/** These are set in the constructor, but if a host code changes the
datastore contents ( such as wiping out the datastore and loading in new
contents from a file, i.e. a restart ) these pointers will need to be
reset to valid groups in the datastore.
@sa Load(const std::string &path, const std::string &protocol).
*/
void SetGroupPointers(axom::sidre::Group * global_grp,
axom::sidre::Group * domain_grp);
axom::sidre::Group * GetBPGroup() { return m_bp_grp; }
axom::sidre::Group * GetBPIndexGroup() { return m_bp_index_grp; }
/// Prepare the DataStore for writing
virtual void PrepareToSave();
/// Save the collection to file.
/** This method calls `Save(collection_name, "sidre_hdf5")`. */
virtual void Save();
/// Save the collection to @a filename.
/** The collection path prefix is prepended to the @a filename and the
current cycle is appended, if cycle >= 0. */
void Save(const std::string& filename, const std::string& protocol);
/// Load the Sidre DataStore from file.
/** No mesh or fields are read from the loaded DataStore.
If the data collection created the datastore, it knows the layout of
where the domain and global groups are, and can restore them after the
Load().
If, however, the data collection does not own the datastore (e.g. it did
not create the datastore), the host code must reset these pointers after
the load operation, using SetGroupPointers(), and also reset the state
variables, using UpdateStateFromDS().
*/
void Load(const std::string& path, const std::string& protocol);
/// Load SidreDataCollection from file.
/** The used file path is based on the current prefix path, collection name,
and the given @a cycle_. The protocol is "sidre_hdf5".
@sa Load(const std::string &path, const std::string &protocol).
*/
virtual void Load(int cycle_ = 0)
{
SetCycle(cycle_);
Load(get_file_path(name), "sidre_hdf5");
}
/// Load external data after registering externally owned fields.
void LoadExternalData(const std::string& path);
/** @brief Updates the DataCollection's cycle, time, and time-step variables
with the values from the data store. */
void UpdateStateFromDS();
/** @brief Updates the data store's cycle, time, and time-step variables with
the values from the SidreDataCollection. */
void UpdateStateToDS();
/** @name Methods for named buffer access and manipulation. */
///@{
/** @brief Get a pointer to the sidre::View holding the named buffer for
@a buffer_name. */
/** If such named buffer is not allocated, the method returns NULL.
@note To access the underlying pointer, use View::getData().
@note To query the size of the buffer, use View::getNumElements().
*/
axom::sidre::View *
GetNamedBuffer(const std::string& buffer_name) const
{
return named_buffers_grp()->hasView(buffer_name)
? named_buffers_grp()->getView(buffer_name)
: NULL;
}
/// Return newly allocated or existing named buffer for @a buffer_name.
/** The buffer is stored in the named_buffers group. If the currently
allocated buffer size is smaller than @a sz, then the buffer is
reallocated with size @a sz, destroying its contents.
@note To access the underlying pointer, use View::getData().
*/
axom::sidre::View *
AllocNamedBuffer(const std::string& buffer_name,
axom::sidre::IndexType sz,
axom::sidre::TypeID type =
axom::sidre::DOUBLE_ID);
/// Deallocate the named buffer @a buffer_name.
void FreeNamedBuffer(const std::string& buffer_name)
{ named_buffers_grp()->destroyViewAndData(buffer_name); }
///@}
private:
// Used if the Sidre data collection is providing the datastore itself.
const bool m_owns_datastore;
// TODO - Need to evaluate if this bool member can be combined with own_data
// in parent data collection class. m_owns_mesh_data indicates whether the
// Sidre dc owns the mesh element data and node positions gf. The DC base
// class own_data indicates if the dc owns the mesh object pointer itself and
// GF objects. Can we use one flag and just have DC own all objects vs none?
const bool m_owns_mesh_data;
// Name to be used for registering the mesh nodes in the SidreDataCollection.
// This name is used by SetMesh() and can be overwritten by the method
// SetMeshNodesName().
// Default value: "mesh_nodes".
std::string m_meshNodesGFName;
// If the data collection owns the datastore, it will store a pointer to it.
// Otherwise, this pointer is NULL.
axom::sidre::DataStore * m_datastore_ptr;
protected:
axom::sidre::Group *named_buffers_grp() const;
axom::sidre::View *
alloc_view(axom::sidre::Group *grp,
const std::string &view_name);
axom::sidre::View *
alloc_view(axom::sidre::Group *grp,
const std::string &view_name,
const axom::sidre::DataType &dtype);
axom::sidre::Group *
alloc_group(axom::sidre::Group *grp,
const std::string &group_name);
// return the filename based on prefix_path, collection name and cycle.
std::string get_file_path(const std::string &filename) const;
private:
// If the data collection does not own the datastore, it will need pointers
// to the blueprint and blueprint index group to use.
axom::sidre::Group * m_bp_grp;
axom::sidre::Group * m_bp_index_grp;
// This is stored for convenience.
axom::sidre::Group * m_named_bufs_grp;
// Private helper functions
void RegisterFieldInBPIndex(const std::string& field_name,
GridFunction *gf);
void DeregisterFieldInBPIndex(const std::string & field_name);
void RegisterAttributeFieldInBPIndex(const std::string& attr_name);
void DeregisterAttributeFieldInBPIndex(const std::string& attr_name);
/** @brief Return a string with the conduit blueprint name for the given
Element::Type. */
std::string getElementName( Element::Type elementEnum );
/**
* \brief A private helper function to set up the views associated with the
data of a scalar valued grid function in the blueprint style.
* \pre gf is not null
* \note This function is expected to be called by RegisterField()
* \note Handles cases where hierarchy is already set up,
* where the data was allocated by this data collection
* and where the grid function data is external to Sidre
*/
void addScalarBasedGridFunction(const std::string& field_name,
GridFunction* gf,
const std::string &buffer_name,
axom::sidre::IndexType offset);
/**
* \brief A private helper function to set up the views associated with the
data of a vector valued grid function in the blueprint style.
* \pre gf is not null
* \note This function is expected to be called by RegisterField()
* \note Handles cases where hierarchy is already set up,
* where the data was allocated by this data collection
* and where the grid function data is external to Sidre
*/
void addVectorBasedGridFunction(const std::string& field_name,
GridFunction* gf,
const std::string &buffer_name,
axom::sidre::IndexType offset);
/** @brief A private helper function to set up the Views associated with
attribute field named @a field_name */
void addIntegerAttributeField(const std::string& field_name, bool is_bdry);
/// Sets up the four main mesh blueprint groups.
/**
* \param hasBP Indicates whether the blueprint has already been set up.
*/
void createMeshBlueprintStubs(bool hasBP);
/// Sets up the mesh blueprint 'state' group.
/**
* \param hasBP Indicates whether the blueprint has already been set up.
*/
void createMeshBlueprintState(bool hasBP);
/// Sets up the mesh blueprint 'coordsets' group.
/**
* \param hasBP Indicates whether the blueprint has already been set up.
*/
void createMeshBlueprintCoordset(bool hasBP);
/// Sets up the mesh blueprint 'topologies' group.
/**
* This method is called from SetMesh().
* \param hasBP Indicates whether the blueprint has already been set up.
* \param mesh_name The name of the topology.
* \note Valid values for @a mesh_name are "mesh" and "boundary" and the
former has to be created with this method before the latter.
*/
void createMeshBlueprintTopologies(bool hasBP, const std::string& mesh_name);
#ifdef MFEM_USE_MPI
/// Sets up the mesh blueprint 'adjacencies' group.
/**
* \param hasBP Indicates whether the blueprint has already been set up.
* \note Only valid when using parallel meshes
*/
void createMeshBlueprintAdjacencies(bool hasBP);
#endif
/// Verifies that the contents of the mesh blueprint data is valid.
void verifyMeshBlueprint();
};
} // end namespace mfem
#endif
#endif
+3 -12
View File
@@ -278,18 +278,9 @@ void ArraysByName<T>::Load(std::istream &in)
q1 = ArrayLine.find(' ');
ArrayName = ArrayLine.substr(0,q1-1);
}
if (q1+2 < ArrayLine.size())
{
// Read the remainder of the line which contains the array data
std::istringstream ArrayDataStream(ArrayLine.substr(q1+2,
ArrayLine.size()));
data[ArrayName].Load(ArrayDataStream, 0);
}
else
{
// Read the array data starting on the next line
data[ArrayName].Load(in, 0);
}
// Ignore the remainder of the line which may contain explanatory comments
data[ArrayName].Load(in, 0);
}
}
+4 -4
View File
@@ -726,16 +726,16 @@ std::string Device::GetUUID(const int device_id)
MFEM_GPU_CHECK(cudaGetDeviceProperties(&prop, device_id));
for (int i = 0; i < 16; ++i)
{
const unsigned b = static_cast<unsigned char>(prop.uuid.bytes[i]);
res << std::setfill('0') << std::setw(2) << std::hex << b;
res << std::setfill('0') << std::setw(2) << std::hex
<< static_cast<unsigned>(prop.uuid.bytes[i]);
}
#elif defined(MFEM_USE_HIP)
hipUUID uuid;
MFEM_GPU_CHECK(hipDeviceGetUuid(&uuid, device_id));
for (int i = 0; i < 16; ++i)
{
const unsigned b = static_cast<unsigned char>(uuid.bytes[i]);
res << std::setfill('0') << std::setw(2) << std::hex << b;
res << std::setfill('0') << std::setw(2) << std::hex
<< static_cast<unsigned>(uuid.bytes[i]);
}
#endif
return res.str();
+14 -1
View File
@@ -22,6 +22,14 @@
//#define _WIN32
//#define _aligned_malloc(s,a) malloc(s)
#ifdef NVTX_DEBUG_HPP
#undef NVTX_COLOR
#define NVTX_COLOR ::nvtx::kGold
#include NVTX_DEBUG_HPP
#else
#define dbg(...)
#endif
#ifndef _WIN32
#include <unistd.h>
#include <signal.h>
@@ -760,7 +768,12 @@ private:
switch (mt)
{
case MT::HOST_DEBUG:
if (GetEnv("MFEM_MMU_STD")) { return new StdHostMemorySpace(); }
if (GetEnv("MFEM_MMU_STD"))
{
dbg("Using STD memory space for debug device!");
return new StdHostMemorySpace();
}
dbg("Using MMU memory space for debug device!");
return new MmuHostMemorySpace();
#ifdef MFEM_USE_UMPIRE
case MT::HOST_UMPIRE:
+8
View File
@@ -200,6 +200,10 @@ void OptionsParser::Parse()
isValid = isValidAsInt(argv[i]);
*(int *)(options[j].var_ptr) = atoi(argv[i++]);
break;
case LONG:
isValid = isValidAsInt(argv[i]);
*(long long *)(options[j].var_ptr) = atoi(argv[i++]);
break;
case DOUBLE:
isValid = isValidAsDouble(argv[i]);
*(real_t *)(options[j].var_ptr) = atof(argv[i++]);
@@ -278,6 +282,10 @@ void OptionsParser::WriteValue(const Option &opt, std::ostream &os)
case INT:
os << *(int *)(opt.var_ptr);
break;
case LONG:
os << *(long long *)(opt.var_ptr);
break;
case DOUBLE:
os << *(real_t *)(opt.var_ptr);
+9 -1
View File
@@ -31,7 +31,7 @@ class Vector;
class OptionsParser
{
public:
enum OptionType { INT, DOUBLE, STRING, STD_STRING, ENABLE, DISABLE, ARRAY, VECTOR };
enum OptionType { INT, DOUBLE, STRING, STD_STRING, ENABLE, DISABLE, ARRAY, VECTOR , LONG };
private:
struct Option
@@ -98,6 +98,14 @@ public:
required));
}
/// Add a long integer option and set 'var' to receive the value.
void AddOption(long long *var, const char *short_name, const char *long_name,
const char *description, bool required = false)
{
options.Append(Option(LONG, var, short_name, long_name, description,
required));
}
/// Add a double option and set 'var' to receive the value.
void AddOption(real_t *var, const char *short_name, const char *long_name,
const char *description, bool required = false)
+3
View File
@@ -160,6 +160,9 @@ const char *GetConfigStr()
#ifdef MFEM_USE_RAJA
"MFEM_USE_RAJA\n"
#endif
#ifdef MFEM_USE_SIDRE
"MFEM_USE_SIDRE\n"
#endif
#ifdef MFEM_USE_SIMD
"MFEM_USE_SIMD\n"
#endif
-3
View File
@@ -317,9 +317,6 @@ void HypreParVector::WrapHypreParVector(hypre_ParVector *y, bool owner)
Vector * HypreParVector::GlobalVector() const
{
MFEM_VERIFY(size > 0,
"GlobalVector method can only be called on vectors wherein each "
"process owns one or more entries");
hypre_Vector *hv = hypre_ParVectorToVectorAll(*this);
Vector *v = new Vector(hv->data, internal::to_int(hv->size));
v->MakeDataOwner();
+81 -58
View File
@@ -38,13 +38,6 @@
#if PETSC_VERSION_LT(3,19,0)
#define PETSC_SUCCESS 0
#endif
#if PETSC_VERSION_LT(3,23,0)
#define PetscContainerSetCtxDestroy(A,B) PetscContainerSetUserDestroy(A,B)
typedef PetscErrorCode (PetscCtxDestroyFn)(void**);
#endif
#if PETSC_VERSION_LT(3,24,0)
typedef PetscErrorCode KSPMonitorFn(KSP,PetscInt,PetscReal,void*);
#endif
#include <fstream>
#include <iomanip>
@@ -84,17 +77,13 @@ static PetscErrorCode __mfem_mat_shell_apply_transpose(Mat,Vec,Vec);
static PetscErrorCode __mfem_mat_shell_destroy(Mat);
static PetscErrorCode __mfem_mat_shell_copy(Mat,Mat,MatStructure);
#if PETSC_VERSION_LT(3,23,0)
typedef void *PetscCtxRt;
#elif PETSC_VERSION_LT(3,25,0)
typedef void **PetscCtxRt;
#endif
static PetscErrorCode __mfem_array_container_destroy(PetscCtxRt);
static PetscErrorCode __mfem_matarray_container_destroy(PetscCtxRt);
#if PETSC_VERSION_LT(3,23,0)
static PetscErrorCode __mfem_monitor_ctx_destroy(void**);
static PetscErrorCode __mfem_array_container_destroy(void*);
static PetscErrorCode __mfem_matarray_container_destroy(void *);
#else
static PetscErrorCode __mfem_monitor_ctx_destroy(PetscCtxRt);
static PetscErrorCode __mfem_array_container_destroy(void**);
static PetscErrorCode __mfem_matarray_container_destroy(void**);
#endif
static PetscErrorCode __mfem_monitor_ctx_destroy(void**);
// auxiliary functions
static PetscErrorCode Convert_Array_IS(MPI_Comm,bool,const mfem::Array<int>*,
@@ -1328,7 +1317,11 @@ BlockDiagonalConstructor(MPI_Comm comm,
ierr = PetscContainerCreate(comm,&c); CCHKERRQ(comm,ierr);
ierr = PetscContainerSetPointer(c,ptrs[i]); CCHKERRQ(comm,ierr);
#if PETSC_VERSION_LT(3,23,0)
ierr = PetscContainerSetUserDestroy(c,__mfem_array_container_destroy);
#else
ierr = PetscContainerSetCtxDestroy(c,__mfem_array_container_destroy);
#endif
CCHKERRQ(comm,ierr);
ierr = PetscObjectCompose((PetscObject)A,names[i],(PetscObject)c);
CCHKERRQ(comm,ierr);
@@ -1655,7 +1648,11 @@ void PetscParMatrix::ConvertOperator(MPI_Comm comm, const Operator &op, Mat* A,
PetscContainer c;
ierr = PetscContainerCreate(comm,&c); CCHKERRQ(comm,ierr);
ierr = PetscContainerSetPointer(c,vmatsl2l); PCHKERRQ(c,ierr);
#if PETSC_VERSION_LT(3,23,0)
ierr = PetscContainerSetUserDestroy(c,__mfem_matarray_container_destroy);
#else
ierr = PetscContainerSetCtxDestroy(c,__mfem_matarray_container_destroy);
#endif
PCHKERRQ(c,ierr);
ierr = PetscObjectCompose((PetscObject)(*A),"_MatIS_PtAP_l2l",(PetscObject)c);
PCHKERRQ((*A),ierr);
@@ -1751,7 +1748,11 @@ void PetscParMatrix::ConvertOperator(MPI_Comm comm, const Operator &op, Mat* A,
ierr = PetscContainerCreate(PETSC_COMM_SELF,&c); PCHKERRQ(B,ierr);
ierr = PetscContainerSetPointer(c,ptrs[i]); PCHKERRQ(B,ierr);
#if PETSC_VERSION_LT(3,23,0)
ierr = PetscContainerSetUserDestroy(c,__mfem_array_container_destroy);
#else
ierr = PetscContainerSetCtxDestroy(c,__mfem_array_container_destroy);
#endif
PCHKERRQ(B,ierr);
ierr = PetscObjectCompose((PetscObject)(B),names[i],(PetscObject)c);
PCHKERRQ(B,ierr);
@@ -2197,7 +2198,11 @@ PetscParMatrix * RAP(PetscParMatrix *Rt, PetscParMatrix *A, PetscParMatrix *P)
ierr = PetscContainerCreate(PetscObjectComm((PetscObject)B),&c);
PCHKERRQ(B,ierr);
ierr = PetscContainerSetPointer(c,vmatsl2l); PCHKERRQ(c,ierr);
#if PETSC_VERSION_LT(3,23,0)
ierr = PetscContainerSetUserDestroy(c,__mfem_matarray_container_destroy);
#else
ierr = PetscContainerSetCtxDestroy(c,__mfem_matarray_container_destroy);
#endif
PCHKERRQ(c,ierr);
ierr = PetscObjectCompose((PetscObject)B,"_MatIS_PtAP_l2l",(PetscObject)c);
PCHKERRQ(B,ierr);
@@ -2480,6 +2485,7 @@ void PetscSolver::SetMaxIter(int max_iter)
void PetscSolver::SetPrintLevel(int plev)
{
typedef PetscErrorCode (*myPetscFunc)(void**);
PetscViewerAndFormat *vf = NULL;
PetscViewer viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm(obj));
@@ -2492,6 +2498,7 @@ void PetscSolver::SetPrintLevel(int plev)
{
// there are many other options, see the function KSPSetFromOptions() in
// src/ksp/ksp/interface/itcl.c
typedef PetscErrorCode (*myMonitor)(KSP,PetscInt,PetscReal,void*);
KSP ksp = (KSP)obj;
if (plev >= 0)
{
@@ -2500,29 +2507,29 @@ void PetscSolver::SetPrintLevel(int plev)
if (plev == 1)
{
#if PETSC_VERSION_LT(3,15,0)
ierr = KSPMonitorSet(ksp,(KSPMonitorFn *)KSPMonitorDefault,vf,
ierr = KSPMonitorSet(ksp,(myMonitor)KSPMonitorDefault,vf,
#else
ierr = KSPMonitorSet(ksp,(KSPMonitorFn *)KSPMonitorResidual,vf,
ierr = KSPMonitorSet(ksp,(myMonitor)KSPMonitorResidual,vf,
#endif
(PetscCtxDestroyFn *)PetscViewerAndFormatDestroy);
(myPetscFunc)PetscViewerAndFormatDestroy);
PCHKERRQ(ksp,ierr);
}
else if (plev > 1)
{
ierr = KSPSetComputeSingularValues(ksp,PETSC_TRUE); PCHKERRQ(ksp,ierr);
ierr = KSPMonitorSet(ksp,(KSPMonitorFn *)KSPMonitorSingularValue,vf,
(PetscCtxDestroyFn *)PetscViewerAndFormatDestroy);
ierr = KSPMonitorSet(ksp,(myMonitor)KSPMonitorSingularValue,vf,
(myPetscFunc)PetscViewerAndFormatDestroy);
PCHKERRQ(ksp,ierr);
if (plev > 2)
{
ierr = PetscViewerAndFormatCreate(viewer,PETSC_VIEWER_DEFAULT,&vf);
PCHKERRQ(viewer,ierr);
#if PETSC_VERSION_LT(3,15,0)
ierr = KSPMonitorSet(ksp,(KSPMonitorFn *)KSPMonitorTrueResidualNorm,vf,
ierr = KSPMonitorSet(ksp,(myMonitor)KSPMonitorTrueResidualNorm,vf,
#else
ierr = KSPMonitorSet(ksp,(KSPMonitorFn *)KSPMonitorTrueResidual,vf,
ierr = KSPMonitorSet(ksp,(myMonitor)KSPMonitorTrueResidual,vf,
#endif
(PetscCtxDestroyFn *)PetscViewerAndFormatDestroy);
(myPetscFunc)PetscViewerAndFormatDestroy);
PCHKERRQ(ksp,ierr);
}
}
@@ -2538,7 +2545,7 @@ void PetscSolver::SetPrintLevel(int plev)
if (plev > 0)
{
ierr = SNESMonitorSet(snes,(myMonitor)SNESMonitorDefault,vf,
(PetscCtxDestroyFn *)PetscViewerAndFormatDestroy);
(myPetscFunc)PetscViewerAndFormatDestroy);
PCHKERRQ(snes,ierr);
}
}
@@ -4156,31 +4163,20 @@ void PetscNonlinearSolver::SetUpdate(void (*update)(Operator *,int,
void PetscNonlinearSolver::Mult(const Vector &b, Vector &x) const
{
SNES snes = (SNES)obj;
MPI_Comm comm = PetscObjectComm(obj);
// Reduction needed: some processes may have null local size while others don't,
// and VecPlaceArray (used by PlaceMemory) is a logically collective operation.
PetscBool b_nonempty = b.Size() ? PETSC_TRUE : PETSC_FALSE;
#if PETSC_VERSION_LT(3,24,0)
mpiierr = MPI_Allreduce(MPI_IN_PLACE,&b_nonempty,1,MPIU_BOOL,MPI_LOR,comm);
#else
mpiierr = MPI_Allreduce(MPI_IN_PLACE,&b_nonempty,1,MPI_C_BOOL,MPI_LOR,comm);
#endif
CCHKERRQ(comm,mpiierr);
// Always create B with allocate=false so that PlaceMemory can be called on
// it regardless of whether b was empty on a previous call.
if (!B) { B = new PetscParVector(comm, *this, true, false); }
if (!X) { X = new PetscParVector(comm, *this, false, false); }
bool b_nonempty = b.Size();
if (!B) { B = new PetscParVector(PetscObjectComm(obj), *this, true); }
if (!X) { X = new PetscParVector(PetscObjectComm(obj), *this, false, false); }
X->PlaceMemory(x.GetMemory(),iterative_mode);
if (b_nonempty) { B->PlaceMemory(b.GetMemory()); }
else { *B = 0.0; }
Customize();
if (!iterative_mode) { *X = 0.; }
// Solve the system. Pass nullptr for b when empty (PETSc treats it as zero RHS).
ierr = SNESSolve(snes, b_nonempty ? B->x : nullptr, X->x); PCHKERRQ(snes, ierr);
// Solve the system.
ierr = SNESSolve(snes, B->x, X->x); PCHKERRQ(snes, ierr);
X->ResetMemory();
if (b_nonempty) { B->ResetMemory(); }
}
@@ -5333,27 +5329,21 @@ static PetscErrorCode __mfem_pc_shell_destroy(PC pc)
PetscFunctionReturn(PETSC_SUCCESS);
}
static PetscErrorCode __mfem_array_container_destroy(PetscCtxRt ptr)
#if PETSC_VERSION_LT(3,23,0)
static PetscErrorCode __mfem_array_container_destroy(void *ptr)
{
PetscErrorCode ierr;
PetscFunctionBeginUser;
#if PETSC_VERSION_LT(3,23,0)
ierr = PetscFree(ptr); CHKERRQ(ierr);
#else
ierr = PetscFree(*(void**)ptr); CHKERRQ(ierr);
#endif
PetscFunctionReturn(PETSC_SUCCESS);
}
static PetscErrorCode __mfem_matarray_container_destroy(PetscCtxRt ptr)
static PetscErrorCode __mfem_matarray_container_destroy(void *ptr)
{
#if PETSC_VERSION_LT(3,23,0)
mfem::Array<Mat> *a = (mfem::Array<Mat>*)ptr;
#else
mfem::Array<Mat> *a = *(mfem::Array<Mat>**)ptr;
#endif
PetscErrorCode ierr;
PetscErrorCode ierr;
PetscFunctionBeginUser;
for (int i=0; i<a->Size(); i++)
@@ -5366,16 +5356,41 @@ static PetscErrorCode __mfem_matarray_container_destroy(PetscCtxRt ptr)
PetscFunctionReturn(PETSC_SUCCESS);
}
#if PETSC_VERSION_LT(3,23,0)
static PetscErrorCode __mfem_monitor_ctx_destroy(void **ctx)
#else
static PetscErrorCode __mfem_monitor_ctx_destroy(PetscCtxRt ctx)
#endif
static PetscErrorCode __mfem_array_container_destroy(void **ptr)
{
PetscErrorCode ierr;
PetscFunctionBeginUser;
ierr = PetscFree(*(void**)ctx); CHKERRQ(ierr);
ierr = PetscFree(*ptr); CHKERRQ(ierr);
PetscFunctionReturn(PETSC_SUCCESS);
}
static PetscErrorCode __mfem_matarray_container_destroy(void **ptr)
{
mfem::Array<Mat> *a = (mfem::Array<Mat>*)*ptr;
PetscErrorCode ierr;
PetscFunctionBeginUser;
for (int i=0; i<a->Size(); i++)
{
Mat M = (*a)[i];
MPI_Comm comm = PetscObjectComm((PetscObject)M);
ierr = MatDestroy(&M); CCHKERRQ(comm,ierr);
}
delete a;
PetscFunctionReturn(PETSC_SUCCESS);
}
#endif
static PetscErrorCode __mfem_monitor_ctx_destroy(void **ctx)
{
PetscErrorCode ierr;
PetscFunctionBeginUser;
ierr = PetscFree(*ctx); CHKERRQ(ierr);
PetscFunctionReturn(PETSC_SUCCESS);
}
@@ -5620,7 +5635,11 @@ static PetscErrorCode MatConvert_hypreParCSR_AIJ(hypre_ParCSRMatrix* hA,Mat* pA)
ierr = PetscContainerCreate(comm,&c); CHKERRQ(ierr);
ierr = PetscContainerSetPointer(c,ptrs[i]); CHKERRQ(ierr);
#if PETSC_VERSION_LT(3,23,0)
ierr = PetscContainerSetUserDestroy(c,__mfem_array_container_destroy);
#else
ierr = PetscContainerSetCtxDestroy(c,__mfem_array_container_destroy);
#endif
CHKERRQ(ierr);
ierr = PetscObjectCompose((PetscObject)(*pA),names[i],(PetscObject)c);
CHKERRQ(ierr);
@@ -5714,7 +5733,11 @@ static PetscErrorCode MatConvert_hypreParCSR_IS(hypre_ParCSRMatrix* hA,Mat* pA)
ierr = PetscContainerCreate(PETSC_COMM_SELF,&c); CHKERRQ(ierr);
ierr = PetscContainerSetPointer(c,ptrs[i]); CHKERRQ(ierr);
#if PETSC_VERSION_LT(3,23,0)
ierr = PetscContainerSetUserDestroy(c,__mfem_array_container_destroy);
#else
ierr = PetscContainerSetCtxDestroy(c,__mfem_array_container_destroy);
#endif
CHKERRQ(ierr);
ierr = PetscObjectCompose((PetscObject)lA,names[i],(PetscObject)c);
CHKERRQ(ierr);
+5 -4
View File
@@ -126,11 +126,11 @@ EXAMPLE_TEST_DIRS := examples
MINIAPP_SUBDIRS = common electromagnetics meshing performance tools \
toys nurbs gslib adjoint solvers shifted mtop parelag tribol autodiff dfem \
hooke multidomain dpg hdiv-linear-solver spde diag-smoothers contact \
fluids/navier fluids/schrodinger-flow plasma
fluids/navier fluids/schrodinger-flow
MINIAPP_DIRS := $(addprefix miniapps/,$(MINIAPP_SUBDIRS))
MINIAPP_TEST_DIRS := $(filter-out %/common,$(MINIAPP_DIRS))
MINIAPP_USE_COMMON := $(addprefix miniapps/,electromagnetics meshing tools \
toys shifted dpg diag-smoothers fluids/navier plasma)
toys shifted dpg diag-smoothers fluids/navier)
EM_DIRS = $(EXAMPLE_DIRS) $(MINIAPP_DIRS)
@@ -299,7 +299,7 @@ ifeq ($(MFEM_USE_LEGACY_OPENMP),YES)
endif
# List of MFEM dependencies, that require the *_LIB variable to be non-empty
MFEM_REQ_LIB_DEPS = SUPERLU MUMPS METIS FMS CONDUIT LAPACK SUNDIALS\
MFEM_REQ_LIB_DEPS = SUPERLU MUMPS METIS FMS CONDUIT SIDRE LAPACK SUNDIALS\
SUITESPARSE STRUMPACK GINKGO GNUTLS HDF5 NETCDF SLEPC PETSC MPFR PUMI HIOP\
GSLIB OCCA CEED RAJA UMPIRE MKL_CPARDISO MKL_PARDISO AMGX MAGMA CALIPER PARELAG\
TRIBOL BENCHMARK MOONOLITH ALGOIM
@@ -365,7 +365,7 @@ MFEM_DEFINES = MFEM_VERSION MFEM_VERSION_STRING MFEM_GIT_STRING MFEM_USE_MPI\
MFEM_USE_LEGACY_OPENMP MFEM_USE_MEMALLOC MFEM_TIMER_TYPE MFEM_USE_SUNDIALS\
MFEM_USE_SUITESPARSE MFEM_USE_GINKGO MFEM_USE_SUPERLU MFEM_USE_SUPERLU5\
MFEM_USE_STRUMPACK MFEM_USE_GNUTLS MFEM_USE_HDF5 MFEM_USE_NETCDF MFEM_USE_PETSC\
MFEM_USE_SLEPC MFEM_USE_MPFR MFEM_USE_FMS MFEM_USE_CONDUIT\
MFEM_USE_SLEPC MFEM_USE_MPFR MFEM_USE_SIDRE MFEM_USE_FMS MFEM_USE_CONDUIT\
MFEM_USE_PUMI MFEM_USE_HIOP MFEM_USE_GSLIB MFEM_USE_CUDA MFEM_USE_HIP\
MFEM_USE_OCCA MFEM_USE_MOONOLITH MFEM_USE_CEED MFEM_USE_RAJA MFEM_USE_UMPIRE\
MFEM_USE_SIMD MFEM_USE_ADIOS2 MFEM_USE_MKL_CPARDISO MFEM_USE_MKL_PARDISO MFEM_USE_AMGX\
@@ -746,6 +746,7 @@ status info:
$(info MFEM_USE_PETSC = $(MFEM_USE_PETSC))
$(info MFEM_USE_SLEPC = $(MFEM_USE_SLEPC))
$(info MFEM_USE_MPFR = $(MFEM_USE_MPFR))
$(info MFEM_USE_SIDRE = $(MFEM_USE_SIDRE))
$(info MFEM_USE_FMS = $(MFEM_USE_FMS))
$(info MFEM_USE_CONDUIT = $(MFEM_USE_CONDUIT))
$(info MFEM_USE_PUMI = $(MFEM_USE_PUMI))
-12
View File
@@ -3206,22 +3206,10 @@ public:
/// Extrude a 1D mesh
/**
* @param mesh 1D mesh
* @param ny number of transverse elements of the extruded mesh
* @param sy physical size in the direction of extrusion
* @param closed if false, only the original boundaries are extruded,
* otherwise boundaries are generated all around the domain
*/
Mesh *Extrude1D(Mesh *mesh, const int ny, const real_t sy,
const bool closed = false);
/// Extrude a 2D mesh
/**
* @param mesh 2D mesh
* @param nz number of transverse elements of the extruded mesh
* @param sz physical size in the direction of extrusion
*/
Mesh *Extrude2D(Mesh *mesh, const int nz, const real_t sz);
/** @brief Constructs the smallest possible [0,1]^dim serial mesh that can be
+3 -6
View File
@@ -1516,15 +1516,12 @@ void Mesh::ReadInlineMesh(std::istream &input, bool generate_edges)
void Mesh::ReadGmshMesh(std::istream &input, int &curved, int &read_gf)
{
string buff;
string version;
real_t version;
int binary, dsize;
input >> version >> binary >> dsize;
if (version != "2.2")
if (version < 2.2)
{
MFEM_ABORT("Gmsh file version must be 2.2, found version "
<< version << ".\n"
"To convert your mesh to the required format, use:\n"
" gmsh -format msh22 -save -o output.msh input.msh");
MFEM_ABORT("Gmsh file version < 2.2");
}
if (dsize != sizeof(double))
{
-6
View File
@@ -5639,12 +5639,6 @@ Mesh ParMesh::GetSerialMesh(int save_rank) const
}
}
if (MyRank == save_rank)
{
attribute_sets.Copy(serialmesh.attribute_sets);
bdr_attribute_sets.Copy(serialmesh.bdr_attribute_sets);
}
MPI_Barrier(MyComm);
return serialmesh;
}
+2 -1
View File
@@ -31,11 +31,12 @@ add_subdirectory(hdiv-linear-solver)
add_subdirectory(hooke)
add_subdirectory(meshing)
add_subdirectory(mtop)
add_subdirectory(mtop/chpt)
add_subdirectory(mtop/examples)
add_subdirectory(multidomain)
add_subdirectory(nurbs)
add_subdirectory(parelag)
add_subdirectory(performance)
add_subdirectory(plasma)
add_subdirectory(shifted)
add_subdirectory(solvers)
add_subdirectory(spde)
+1 -38
View File
@@ -82,43 +82,6 @@ Although Tribol can be built automatically via **uberenv** and **Spack**,
for this miniapp it is simpler to build **Axom** and **MFEM** manually and
point Tribol to them. The steps are as follows:
### Using pre-built Tribol/Axom installs
If you already have compatible installs of Tribol and Axom, point MFEM to the install prefixes.
- Hypre install prefix: `<path/to/hypre>`
- METIS install prefix: `<path/to/metis>`
- Axom install prefix: `<path/to/axom>`
- Tribol install prefix: `<path/to/tribol>`
**MFEM make build (configure):**
```bash
make config MFEM_USE_MPI=YES MFEM_USE_METIS=YES MFEM_USE_TRIBOL=YES \
HYPRE_DIR=<path/to/hypre> METIS_DIR=<path/to/metis> \
AXOM_DIR=<path/to/axom> TRIBOL_DIR=<path/to/tribol> ADIAK_DIR=<path/to/adiak> CAMP_DIR=<path/to/camp> RAJA_DIR=<path/to/raja> \
UMPIRE_DIR=<path/to/umpire> FMT_DIR=<path/to/fmt> CALIPER_DIR=<path/to/caliper>
```
**MFEM CMake build (configure):**
```bash
cmake -S . -B <mfem-build-dir> -DMFEM_USE_MPI=YES -DMFEM_USE_METIS=YES -DMFEM_USE_TRIBOL=YES \
HYPRE_DIR=<path/to/hypre> METIS_DIR=<path/to/metis> \
AXOM_DIR=<path/to/axom> TRIBOL_DIR=<path/to/tribol> ADIAK_DIR=<path/to/adiak> CAMP_DIR=<path/to/camp> RAJA_DIR=<path/to/raja> \
UMPIRE_DIR=<path/to/umpire> FMT_DIR=<path/to/fmt> CALIPER_DIR=<path/to/caliper>
```
Note: RAJA/UMPIRE/CALIPER are optional for MFEM itself, but many Tribol builds
enable them. If your Tribol install does not depend on them, you can omit the
corresponding `*_DIR` entries above.
Note: `FMT_DIR` only needs to be added for the make-based build (and only when
the Umpire install uses `fmt`). If `FMT_DIR` is not set and a sibling `fmt-*`
directory exists next to your `UMPIRE_DIR`, MFEM's make configuration will try
to pick it up automatically.
Note: when using pre-built Tribol/Axom, you typically need to use a compatible
compiler/MPI wrapper (same C++ standard library ABI).
### Manual Build Steps
1. Pull axom and tribol (starting from the mfem folder):
@@ -136,7 +99,7 @@ compiler/MPI wrapper (same C++ standard library ABI).
TRIBOL_DIR = @MFEM_DIR@/../tribol-repo/tribol
TRIBOL_OPT = -I$(TRIBOL_DIR)/include -I$(AXOM_DIR)/include
TRIBOL_LIB = -L$(TRIBOL_DIR)/lib -ltribol -lredecomp -L$(AXOM_DIR)/lib \
-laxom_quest -laxom_mint -laxom_slam -laxom_slic -laxom_lumberjack -laxom_core
-laxom_mint -laxom_slam -laxom_slic -laxom_core
```
3. [**Axom:**](https://github.com/LLNL/axom.git) Starting from the MFEM root
directory (we assume this directory is named mfem):
@@ -9,17 +9,24 @@
# terms of the BSD-3 license. We welcome feedback and contributions, see file
# CONTRIBUTING.md for details.
if (MFEM_USE_MPI)
list(APPEND PLASMA_COMMON_SOURCES)
if(MFEM_USE_MPI)
list(APPEND PLASMA_COMMON_HEADERS
plasma.hpp)
list(APPEND AN_SOURCES linear_anisotropic_elasticity.hpp)
list(APPEND AN_HEADERS linear_anisotropic_elasticity.cpp)
convert_filenames_to_full_paths(PLASMA_COMMON_SOURCES)
convert_filenames_to_full_paths(PLASMA_COMMON_HEADERS)
convert_filenames_to_full_paths(AN_SOURCES)
convert_filenames_to_full_paths(AN_HEADERS)
set(EX_COMMON_FILES
EXTRA_SOURCES ${AN_SOURCES}
EXTRA_HEADERS ${AN_HEADERS})
add_mfem_miniapp(mtop_test_anisotropic_le
MAIN test_anisotropic_le.cpp
${EX_COMMON_FILES}
LIBRARIES mfem)
endif (MFEM_USE_MPI)
set(PLASMA_COMMON_FILES
EXTRA_SOURCES ${PLASMA_COMMON_SOURCES}
EXTRA_HEADERS ${PLASMA_COMMON_HEADERS})
endif()
File diff suppressed because it is too large Load Diff
+65
View File
@@ -0,0 +1,65 @@
# Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
# LICENSE and NOTICE for details. LLNL-CODE-806117.
#
# This file is part of the MFEM library. For more information and source code
# availability visit https://mfem.org.
#
# MFEM is free software; you can redistribute it and/or modify it under the
# terms of the BSD-3 license. We welcome feedback and contributions, see file
# CONTRIBUTING.md for details.
list(APPEND CHPT_SOURCES )
list(APPEND CHPT_HEADERS )
convert_filenames_to_full_paths(CHPT_SOURCES)
convert_filenames_to_full_paths(CHPT_HEADERS)
set(CHPT_COMMON_FILES
EXTRA_SOURCES ${CHPT_SOURCES}
EXTRA_HEADERS ${CHPT_HEADERS})
add_mfem_miniapp(mtop_test_scalar_chpt
MAIN tst_scalar_field.cpp
${CHPT_COMMON_FILES}
LIBRARIES mfem)
add_mfem_miniapp(mtop_test_vector_chpt
MAIN tst_vector_field.cpp
${CHPT_COMMON_FILES}
LIBRARIES mfem)
add_mfem_miniapp(mtop_test_fixedslot
MAIN tst_fixedslot_scalar_vector.cpp
${CHPT_COMMON_FILES}
LIBRARIES mfem)
add_mfem_miniapp(mtop_test_segmentio
MAIN tst_segmentio_scalar_vector_mfem.cpp
${CHPT_COMMON_FILES}
LIBRARIES mfem)
add_mfem_miniapp(mtop_test_file
MAIN tst_fileio_scalar_vector_mfem.cpp
${CHPT_COMMON_FILES}
LIBRARIES mfem)
add_mfem_miniapp(mtop_test_revolve
MAIN revolve.cpp
${CHPT_COMMON_FILES}
LIBRARIES mfem)
add_mfem_miniapp(mtop_test_dyn
MAIN tst_dynamic_chpt.cpp
${CHPT_COMMON_FILES}
LIBRARIES mfem)
add_mfem_miniapp(mtop_test_rev_fixedslot
MAIN tst_revolve_chpt_fixedslot.cpp
${CHPT_COMMON_FILES}
LIBRARIES mfem)
@@ -0,0 +1,309 @@
#ifndef MFEM_DYNAMIC_CHECKPOINTING_HPP
#define MFEM_DYNAMIC_CHECKPOINTING_HPP
#include "mfem.hpp"
#include <map>
#include <memory>
#include <limits>
#include <vector>
#include <utility>
namespace mfem
{
template <typename Snapshot>
class InMemoryCheckpointStorage
{
public:
using Handle = int;
Handle InvalidHandle() const { return -1; }
bool IsValid(const Handle &h) const { return h >= 0; }
Handle Store(Snapshot &&snap)
{
Handle id = InvalidHandle();
if (!free_.empty())
{
id = free_.back();
free_.pop_back();
if (slots_[id])
{
// reuse existing allocation
*slots_[id] = std::move(snap);
}
else
{
slots_[id].reset(new Snapshot(std::move(snap)));
}
}
else
{
id = (Handle) slots_.size();
slots_.push_back(std::unique_ptr<Snapshot>(new Snapshot(std::move(snap))));
}
return id;
}
template <typename Func>
void Read(const Handle &h, Func &&f) const
{
MFEM_VERIFY(IsValid(h), "InMemoryCheckpointStorage: invalid handle.");
MFEM_VERIFY(h < (Handle) slots_.size(), "InMemoryCheckpointStorage: handle out of range.");
MFEM_VERIFY(slots_[h].get() != nullptr, "InMemoryCheckpointStorage: empty slot.");
f(*slots_[h]);
}
void Erase(Handle &h)
{
if (!IsValid(h)) { h = InvalidHandle(); return; }
MFEM_VERIFY(h < (Handle) slots_.size(), "InMemoryCheckpointStorage: handle out of range.");
slots_[h].reset();
free_.push_back(h);
h = InvalidHandle();
}
private:
std::vector<std::unique_ptr<Snapshot>> slots_;
std::vector<Handle> free_;
};
/**
* Dynamic checkpointing manager (WangMoinIaccarino 2009), with pluggable storage.
*
* Snapshot: checkpointed object (often the primal State itself).
* Storage : must provide:
* using Handle = ...
* Handle InvalidHandle() const;
* bool IsValid(const Handle&) const;
* Handle Store(Snapshot&&);
* template<class F> void Read(const Handle&, F&&) const; // calls f(const Snapshot&)
* void Erase(Handle&);
*/
template <typename Snapshot,
typename Storage = InMemoryCheckpointStorage<Snapshot>>
class DynamicCheckpointing
{
public:
using Step = long long;
using Handle = typename Storage::Handle;
struct Checkpoint
{
int level = 0;
Handle h; // InvalidHandle => placeholder
Checkpoint() = default;
Checkpoint(int lvl, const Handle &inv) : level(lvl), h(inv) {}
};
explicit DynamicCheckpointing(int s)
: s_(s),
owned_storage_(new Storage()),
storage_(owned_storage_.get())
{
MFEM_VERIFY(s_ > 0, "DynamicCheckpointing: s must be > 0.");
Reset();
}
DynamicCheckpointing(int s, Storage &external_storage)
: s_(s),
owned_storage_(nullptr),
storage_(&external_storage)
{
MFEM_VERIFY(s_ > 0, "DynamicCheckpointing: s must be > 0.");
Reset();
}
~DynamicCheckpointing() { ReleaseAll(); }
void Reset()
{
ReleaseAll();
cps_.emplace(Step(0), Checkpoint(InfLevel(), storage_->InvalidHandle()));
}
struct CheckpointInfo
{
Step step;
int level;
bool stored;
};
std::vector<CheckpointInfo> GetCheckpointInfo() const
{
std::vector<CheckpointInfo> out;
out.reserve(cps_.size());
for (const auto &kv : cps_)
{
out.push_back({kv.first, kv.second.level, storage_->IsValid(kv.second.h)});
}
return out;
}
Step GetMaxStep() const
{
MFEM_VERIFY(!cps_.empty(), "DynamicCheckpointing: checkpoint map is empty.");
return cps_.rbegin()->first;
}
template <typename State, typename PrimalStepFn, typename MakeSnapshotFn>
void ForwardStep(const Step i,
State &u_i_inout,
PrimalStepFn &&primal_step,
MakeSnapshotFn &&make_snapshot)
{
AllocateCheckpointForNextStep(i);
auto it = cps_.find(i);
if (it != cps_.end())
{
storage_->Erase(it->second.h);
Snapshot snap = make_snapshot(u_i_inout);
it->second.h = storage_->Store(std::move(snap));
}
primal_step(u_i_inout, i);
}
template <typename State,
typename AdjState,
typename PrimalStepFn,
typename AdjointStepFn,
typename MakeSnapshotFn,
typename RestoreSnapshotFn>
void BackwardStep(const Step i,
AdjState &q_ip1_inout,
State &u_work_inout,
PrimalStepFn &&primal_step,
AdjointStepFn &&adjoint_step,
MakeSnapshotFn &&make_snapshot,
RestoreSnapshotFn &&restore_snapshot)
{
// remove placeholder at i+1
const Step ph = i + 1;
auto it_ph = cps_.find(ph);
MFEM_VERIFY(it_ph != cps_.end(),
"DynamicCheckpointing: expected checkpoint at i+1 before BackwardStep.");
storage_->Erase(it_ph->second.h);
cps_.erase(it_ph);
MFEM_ASSERT(GetMaxStep() <= i,
"DynamicCheckpointing: found a checkpoint beyond current adjoint step.");
auto restore_from_handle = [&](Handle &h)
{
storage_->Read(h, [&](const Snapshot &snap)
{
restore_snapshot(snap, u_work_inout);
});
storage_->Erase(h); // retrieved => placeholder (Algorithm 4 semantics)
};
if (GetMaxStep() == i)
{
Handle h = TakeHandleMakePlaceholder(i);
restore_from_handle(h);
}
else
{
const Step k = GetMaxStep();
Handle hk = TakeHandleMakePlaceholder(k);
restore_from_handle(hk);
for (Step t = k; t < i; ++t)
{
ForwardStep(t, u_work_inout, primal_step, make_snapshot);
}
}
adjoint_step(q_ip1_inout, u_work_inout, i);
}
private:
int s_ = 0;
std::unique_ptr<Storage> owned_storage_;
Storage *storage_ = nullptr;
std::map<Step, Checkpoint> cps_;
static int InfLevel() { return std::numeric_limits<int>::max(); }
void ReleaseAll()
{
if (!storage_) { cps_.clear(); return; }
for (auto &kv : cps_) { storage_->Erase(kv.second.h); }
cps_.clear();
}
bool FindDispensableLargestStep(Step &out_step) const
{
int max_level_seen = std::numeric_limits<int>::min();
for (auto it = cps_.rbegin(); it != cps_.rend(); ++it)
{
const Step step = it->first;
const int lvl = it->second.level;
if (max_level_seen > lvl) { out_step = step; return true; }
max_level_seen = (lvl > max_level_seen) ? lvl : max_level_seen;
}
return false;
}
void AllocateCheckpointForNextStep(const Step i)
{
const Step new_step = i + 1;
MFEM_VERIFY(cps_.find(new_step) == cps_.end(),
"DynamicCheckpointing: checkpoint at i+1 already exists.");
const Handle inv = storage_->InvalidHandle();
// allow growth to s+1 entries (incl. placeholder)
if ((int)cps_.size() <= s_)
{
cps_.emplace(new_step, Checkpoint(0, inv));
return;
}
Step disp = -1;
if (FindDispensableLargestStep(disp))
{
auto it = cps_.find(disp);
MFEM_ASSERT(it != cps_.end(), "Internal error: dispensable checkpoint not found.");
storage_->Erase(it->second.h);
cps_.erase(it);
cps_.emplace(new_step, Checkpoint(0, inv));
return;
}
auto it_i = cps_.find(i);
MFEM_VERIFY(it_i != cps_.end(),
"DynamicCheckpointing: promotion expected checkpoint at step i but none found.");
MFEM_VERIFY(i != 0, "DynamicCheckpointing: attempted to remove step 0 checkpoint.");
const int l = it_i->second.level;
storage_->Erase(it_i->second.h);
cps_.erase(it_i);
cps_.emplace(new_step, Checkpoint(l + 1, inv));
}
Handle TakeHandleMakePlaceholder(const Step i)
{
auto it = cps_.find(i);
MFEM_VERIFY(it != cps_.end(),
"DynamicCheckpointing: TakeHandle requested a non-existent checkpoint.");
MFEM_VERIFY(storage_->IsValid(it->second.h),
"DynamicCheckpointing: TakeHandle requested a checkpoint with no snapshot.");
Handle h = std::move(it->second.h);
it->second.h = storage_->InvalidHandle();
return h;
}
};
} // namespace mfem
#endif // MFEM_DYNAMIC_CHECKPOINTING_HPP
@@ -0,0 +1,291 @@
#ifndef MFEM_FILE_CHECKPOINT_STORAGE_HPP
#define MFEM_FILE_CHECKPOINT_STORAGE_HPP
#include "mfem.hpp"
#include <cstdint>
#include <cstdio> // std::remove, std::rename
#include <cstring>
#include <string>
#include <vector>
#include <sstream>
#include <iomanip>
#include <fstream>
#include <type_traits>
#if __has_include(<filesystem>)
#include <filesystem>
namespace mfem_fs = std::filesystem;
#define MFEM_HAVE_FILESYSTEM 1
#else
#define MFEM_HAVE_FILESYSTEM 0
#endif
namespace mfem
{
/**
* @brief Default binary serializer for checkpoint snapshots.
*
* Contract:
* - static void Write(std::ostream&, const Snapshot&)
* - static Snapshot Read (std::istream&)
*
* Default implementation supports trivially-copyable POD types.
* Specialization for mfem::Vector is provided below.
*
* Portability note:
* This binary format is NOT portable across endianness or differing sizeof(real_t).
* If you need portability, provide a custom SnapshotIO (e.g., text, XDR, HDF5).
*/
template <typename Snapshot, typename Enable = void>
struct DefaultCheckpointBinaryIO
{
static void Write(std::ostream &, const Snapshot &)
{
static_assert(sizeof(Snapshot) == 0,
"DefaultCheckpointBinaryIO: no implementation for this Snapshot type. "
"Provide a custom SnapshotIO with Write/Read.");
}
static Snapshot Read(std::istream &)
{
static_assert(sizeof(Snapshot) == 0,
"DefaultCheckpointBinaryIO: no implementation for this Snapshot type. "
"Provide a custom SnapshotIO with Write/Read.");
return Snapshot();
}
};
// POD / trivially-copyable types (double, int, structs of POD, etc.)
template <typename Snapshot>
struct DefaultCheckpointBinaryIO<
Snapshot,
typename std::enable_if<std::is_trivially_copyable<Snapshot>::value>::type>
{
static void Write(std::ostream &os, const Snapshot &x)
{
os.write(reinterpret_cast<const char*>(&x), sizeof(Snapshot));
MFEM_VERIFY(os.good(), "DefaultCheckpointBinaryIO: failed to write POD snapshot.");
}
static Snapshot Read(std::istream &is)
{
Snapshot x;
is.read(reinterpret_cast<char*>(&x), sizeof(Snapshot));
MFEM_VERIFY(is.good(), "DefaultCheckpointBinaryIO: failed to read POD snapshot.");
return x;
}
};
// Specialization for mfem::Vector
template <>
struct DefaultCheckpointBinaryIO<mfem::Vector, void>
{
static void Write(std::ostream &os, const mfem::Vector &v)
{
const std::int64_t n = (std::int64_t) v.Size();
os.write(reinterpret_cast<const char*>(&n), sizeof(n));
MFEM_VERIFY(os.good(), "VectorBinaryIO: failed to write vector size.");
if (n > 0)
{
const mfem::real_t *data = v.GetData();
os.write(reinterpret_cast<const char*>(data),
(std::streamsize)(n * (std::int64_t)sizeof(mfem::real_t)));
MFEM_VERIFY(os.good(), "VectorBinaryIO: failed to write vector data.");
}
}
static mfem::Vector Read(std::istream &is)
{
std::int64_t n = 0;
is.read(reinterpret_cast<char*>(&n), sizeof(n));
MFEM_VERIFY(is.good(), "VectorBinaryIO: failed to read vector size.");
MFEM_VERIFY(n >= 0, "VectorBinaryIO: invalid negative vector size.");
mfem::Vector v((int)n);
if (n > 0)
{
mfem::real_t *data = v.GetData();
is.read(reinterpret_cast<char*>(data),
(std::streamsize)(n * (std::int64_t)sizeof(mfem::real_t)));
MFEM_VERIFY(is.good(), "VectorBinaryIO: failed to read vector data.");
}
return v;
}
};
/**
* @brief File-based checkpoint storage backend.
*
* Stores each snapshot in a separate file:
* <directory>/<prefix><id><extension>
*
* Handle is a monotonically-increasing integer id (reused if keep_files==false).
*
* Template parameters:
* - Snapshot : stored snapshot type
* - SnapshotIO : serializer with static Write/Read methods (see DefaultCheckpointBinaryIO)
*
* Threading:
* - Not thread-safe. Typical adjoint/checkpointing usage is single-threaded control flow.
*/
template <typename Snapshot,
typename SnapshotIO = DefaultCheckpointBinaryIO<Snapshot>>
class FileCheckpointStorage
{
public:
using Handle = std::int64_t;
/**
* @param directory directory where checkpoint files live (created if create_dir==true)
* @param prefix filename prefix (e.g. "ckpt_")
* @param extension filename extension (e.g. ".bin")
* @param create_dir create directory if missing (requires <filesystem>)
* @param keep_files if true, Erase() will NOT delete files (useful for debugging),
* and ids are NOT reused.
*/
FileCheckpointStorage(const std::string &directory,
const std::string &prefix = "ckpt_",
const std::string &extension = ".bin",
bool create_dir = true,
bool keep_files = false)
: dir_(directory),
prefix_(prefix),
ext_(extension),
keep_files_(keep_files)
{
MFEM_VERIFY(!dir_.empty(), "FileCheckpointStorage: directory must be non-empty.");
MFEM_VERIFY(!prefix_.empty(), "FileCheckpointStorage: prefix must be non-empty.");
MFEM_VERIFY(!ext_.empty(), "FileCheckpointStorage: extension must be non-empty.");
if (create_dir)
{
#if MFEM_HAVE_FILESYSTEM
std::error_code ec;
mfem_fs::create_directories(mfem_fs::path(dir_), ec);
MFEM_VERIFY(!ec, "FileCheckpointStorage: failed to create directory: " << dir_);
#else
MFEM_ABORT("FileCheckpointStorage: create_dir=true requires <filesystem> support.");
#endif
}
}
Handle InvalidHandle() const { return (Handle)-1; }
bool IsValid(const Handle &h) const { return h >= 0; }
/**
* @brief Store snapshot to file, return handle.
*
* Uses atomic-ish pattern:
* - write to "<path>.tmp"
* - rename to "<path>"
*/
Handle Store(Snapshot &&snap)
{
const Handle id = AllocateId_();
const std::string path = Path_(id);
const std::string tmp = path + ".tmp";
{
std::ofstream os(tmp, std::ios::binary | std::ios::trunc);
MFEM_VERIFY(os.is_open(), "FileCheckpointStorage: failed to open for write: " << tmp);
// Write payload
SnapshotIO::Write(os, snap);
MFEM_VERIFY(os.good(), "FileCheckpointStorage: write failed for: " << tmp);
}
// Rename tmp -> final
const int rc = std::rename(tmp.c_str(), path.c_str());
MFEM_VERIFY(rc == 0, "FileCheckpointStorage: rename failed: " << tmp << " -> " << path);
return id;
}
/**
* @brief Read snapshot from file and pass it to callback f(const Snapshot&).
*
* The reference passed to f is valid only during the call.
*/
template <typename Func>
void Read(const Handle &h, Func &&f) const
{
MFEM_VERIFY(IsValid(h), "FileCheckpointStorage: Read called with invalid handle.");
const std::string path = Path_(h);
std::ifstream is(path, std::ios::binary);
MFEM_VERIFY(is.is_open(), "FileCheckpointStorage: failed to open for read: " << path);
Snapshot snap = SnapshotIO::Read(is);
MFEM_VERIFY(is.good(), "FileCheckpointStorage: read failed for: " << path);
f(snap);
}
/**
* @brief Erase snapshot (delete file unless keep_files==true), set handle invalid.
*
* Id reuse policy:
* - if keep_files_ == false: deleted ids are reused to avoid creating huge numbers of files
* - if keep_files_ == true : ids are not reused (avoid overwriting old debug files)
*/
void Erase(Handle &h)
{
if (!IsValid(h)) { h = InvalidHandle(); return; }
if (!keep_files_)
{
const std::string path = Path_(h);
const int rc = std::remove(path.c_str());
MFEM_VERIFY(rc == 0, "FileCheckpointStorage: failed to remove file: " << path);
free_.push_back(h);
}
h = InvalidHandle();
}
private:
std::string dir_;
std::string prefix_;
std::string ext_;
bool keep_files_ = false;
mutable Handle next_id_ = 0;
mutable std::vector<Handle> free_;
Handle AllocateId_()
{
if (!keep_files_ && !free_.empty())
{
const Handle id = free_.back();
free_.pop_back();
return id;
}
return next_id_++;
}
std::string Path_(Handle id) const
{
MFEM_VERIFY(IsValid(id), "FileCheckpointStorage: Path_ called with invalid id.");
std::ostringstream oss;
oss << prefix_ << std::setw(12) << std::setfill('0') << id << ext_;
#if MFEM_HAVE_FILESYSTEM
mfem_fs::path p = mfem_fs::path(dir_) / mfem_fs::path(oss.str());
return p.string();
#else
// Fallback: simple concatenation; assumes dir_ ends without trailing slash if needed.
return dir_ + "/" + oss.str();
#endif
}
};
} // namespace mfem
#endif // MFEM_FILE_CHECKPOINT_STORAGE_HPP
@@ -0,0 +1,385 @@
#ifndef MFEM_FIXED_SLOT_CHECKPOINT_STORAGE_HPP
#define MFEM_FIXED_SLOT_CHECKPOINT_STORAGE_HPP
#include "mfem.hpp"
#include <cstddef>
#include <cstdint>
#include <cstring>
#include <string>
#include <vector>
#include <fstream>
#include <limits>
#include <type_traits>
namespace mfem
{
// ------------------------------------------------------------
// Packer 1: trivially-copyable snapshots (memcpy)
// ------------------------------------------------------------
template <typename Snapshot>
class TrivialFixedPacker
{
public:
TrivialFixedPacker()
{
static_assert(std::is_trivially_copyable<Snapshot>::value,
"TrivialFixedPacker requires Snapshot to be trivially copyable.");
}
std::size_t SlotBytes() const { return sizeof(Snapshot); }
void Pack(const Snapshot &snap, void *dst) const
{
std::memcpy(dst, &snap, sizeof(Snapshot));
}
void Unpack(const void *src, Snapshot &out) const
{
std::memcpy(&out, src, sizeof(Snapshot));
}
};
// ------------------------------------------------------------
// Packer 2: mfem::Vector of fixed length n
// (Serialized size is fixed: n*sizeof(real_t))
// ------------------------------------------------------------
class FixedVectorPacker
{
public:
explicit FixedVectorPacker(int n) : n_(n)
{
MFEM_VERIFY(n_ > 0, "FixedVectorPacker: n must be > 0.");
}
std::size_t SlotBytes() const
{
return (std::size_t)n_ * sizeof(mfem::real_t);
}
void Pack(const mfem::Vector &v, void *dst) const
{
MFEM_VERIFY(v.Size() == n_, "FixedVectorPacker: vector size mismatch.");
std::memcpy(dst, v.GetData(), SlotBytes());
}
void Unpack(const void *src, mfem::Vector &out) const
{
out.SetSize(n_);
std::memcpy(out.GetData(), src, SlotBytes());
}
int Size() const { return n_; }
private:
int n_ = 0;
};
// ------------------------------------------------------------
// Fixed-slot MEMORY storage: one big RAM block
// ------------------------------------------------------------
template <typename Snapshot, typename Packer = TrivialFixedPacker<Snapshot>>
class FixedSlotMemoryCheckpointStorage
{
public:
using Handle = int;
FixedSlotMemoryCheckpointStorage(int max_slots, const Packer &packer = Packer())
: max_slots_(max_slots), packer_(packer)
{
MFEM_VERIFY(max_slots_ > 0, "FixedSlotMemoryCheckpointStorage: max_slots must be > 0.");
slot_bytes_ = packer_.SlotBytes();
MFEM_VERIFY(slot_bytes_ > 0, "FixedSlotMemoryCheckpointStorage: SlotBytes must be > 0.");
// Single contiguous block
bytes_.resize((std::size_t)max_slots_ * slot_bytes_);
in_use_.assign(max_slots_, 0);
free_.reserve(max_slots_);
for (int i = 0; i < max_slots_; ++i) { free_.push_back(i); }
}
Handle InvalidHandle() const { return -1; }
bool IsValid(const Handle &h) const { return h >= 0; }
int MaxSlots() const { return max_slots_; }
std::size_t SlotBytes() const { return slot_bytes_; }
Handle Store(Snapshot &&snap)
{
MFEM_VERIFY(!free_.empty(),
"FixedSlotMemoryCheckpointStorage: out of slots (increase max_slots).");
const int slot = free_.back();
free_.pop_back();
in_use_[slot] = 1;
void *dst = SlotPtr_(slot);
packer_.Pack(snap, dst);
return slot;
}
template <typename Func>
void Read(const Handle &h, Func &&f) const
{
MFEM_VERIFY(IsValid(h), "FixedSlotMemoryCheckpointStorage: Read invalid handle.");
MFEM_VERIFY(h < max_slots_, "FixedSlotMemoryCheckpointStorage: Read handle out of range.");
MFEM_VERIFY(in_use_[h] == 1, "FixedSlotMemoryCheckpointStorage: Read from free slot.");
Snapshot tmp;
const void *src = SlotPtrConst_(h);
packer_.Unpack(src, tmp);
f(tmp);
}
void Erase(Handle &h)
{
if (!IsValid(h)) { h = InvalidHandle(); return; }
MFEM_VERIFY(h < max_slots_, "FixedSlotMemoryCheckpointStorage: Erase handle out of range.");
MFEM_VERIFY(in_use_[h] == 1, "FixedSlotMemoryCheckpointStorage: double-free / invalid erase.");
in_use_[h] = 0;
free_.push_back(h);
h = InvalidHandle();
}
/// Optional: return all slots to the free list (does not zero memory).
void Reset()
{
free_.clear();
for (int i = 0; i < max_slots_; ++i) { in_use_[i] = 0; free_.push_back(i); }
}
private:
int max_slots_ = 0;
std::size_t slot_bytes_ = 0;
Packer packer_;
std::vector<unsigned char> bytes_; // single block
std::vector<unsigned char> in_use_;
std::vector<int> free_;
void *SlotPtr_(int slot)
{
return (void*)(&bytes_[(std::size_t)slot * slot_bytes_]);
}
const void *SlotPtrConst_(int slot) const
{
return (const void*)(&bytes_[(std::size_t)slot * slot_bytes_]);
}
};
// ------------------------------------------------------------
// Fixed-slot FILE storage: one single pre-sized file
// ------------------------------------------------------------
template <typename Snapshot, typename Packer = TrivialFixedPacker<Snapshot>>
class FixedSlotFileCheckpointStorage
{
public:
using Handle = int;
struct Header
{
char magic[8]; // "MFCKPTFS"
std::uint64_t version; // 1
std::uint64_t slot_bytes;
std::uint64_t max_slots;
std::uint64_t reserved[4]; // future use / padding
};
FixedSlotFileCheckpointStorage(const std::string &path,
int max_slots,
const Packer &packer = Packer(),
bool truncate = true,
bool flush_on_store = false)
: path_(path),
max_slots_(max_slots),
packer_(packer),
flush_on_store_(flush_on_store)
{
MFEM_VERIFY(!path_.empty(), "FixedSlotFileCheckpointStorage: empty file path.");
MFEM_VERIFY(max_slots_ > 0, "FixedSlotFileCheckpointStorage: max_slots must be > 0.");
slot_bytes_ = packer_.SlotBytes();
MFEM_VERIFY(slot_bytes_ > 0, "FixedSlotFileCheckpointStorage: SlotBytes must be > 0.");
Open_(truncate);
in_use_.assign(max_slots_, 0);
free_.reserve(max_slots_);
for (int i = 0; i < max_slots_; ++i) { free_.push_back(i); }
scratch_.resize(slot_bytes_);
}
~FixedSlotFileCheckpointStorage()
{
if (file_.is_open()) { file_.close(); }
}
Handle InvalidHandle() const { return -1; }
bool IsValid(const Handle &h) const { return h >= 0; }
int MaxSlots() const { return max_slots_; }
std::size_t SlotBytes() const { return slot_bytes_; }
const std::string &Path() const { return path_; }
Handle Store(Snapshot &&snap)
{
MFEM_VERIFY(!free_.empty(),
"FixedSlotFileCheckpointStorage: out of slots (increase max_slots).");
const int slot = free_.back();
free_.pop_back();
in_use_[slot] = 1;
// Pack into scratch buffer then write into fixed slot offset
packer_.Pack(snap, scratch_.data());
const std::uint64_t off = SlotOffset_(slot);
file_.seekp((std::streamoff)off, std::ios::beg);
MFEM_VERIFY(file_.good(), "FixedSlotFileCheckpointStorage: seekp failed.");
file_.write(reinterpret_cast<const char*>(scratch_.data()),
(std::streamsize)slot_bytes_);
MFEM_VERIFY(file_.good(), "FixedSlotFileCheckpointStorage: write failed.");
if (flush_on_store_) { file_.flush(); }
return slot;
}
template <typename Func>
void Read(const Handle &h, Func &&f) const
{
MFEM_VERIFY(IsValid(h), "FixedSlotFileCheckpointStorage: Read invalid handle.");
MFEM_VERIFY(h < max_slots_, "FixedSlotFileCheckpointStorage: Read handle out of range.");
MFEM_VERIFY(in_use_[h] == 1, "FixedSlotFileCheckpointStorage: Read from free slot.");
const std::uint64_t off = SlotOffset_(h);
file_.seekg((std::streamoff)off, std::ios::beg);
MFEM_VERIFY(file_.good(), "FixedSlotFileCheckpointStorage: seekg failed.");
file_.read(reinterpret_cast<char*>(scratch_.data()),
(std::streamsize)slot_bytes_);
MFEM_VERIFY(file_.good(), "FixedSlotFileCheckpointStorage: read failed.");
Snapshot tmp;
packer_.Unpack(scratch_.data(), tmp);
f(tmp);
}
void Erase(Handle &h)
{
if (!IsValid(h)) { h = InvalidHandle(); return; }
MFEM_VERIFY(h < max_slots_, "FixedSlotFileCheckpointStorage: Erase handle out of range.");
MFEM_VERIFY(in_use_[h] == 1, "FixedSlotFileCheckpointStorage: double-free / invalid erase.");
// No file deletion; just return slot to free list.
in_use_[h] = 0;
free_.push_back(h);
h = InvalidHandle();
}
/// Optional: return all slots to free list (file contents remain).
void Reset()
{
free_.clear();
for (int i = 0; i < max_slots_; ++i) { in_use_[i] = 0; free_.push_back(i); }
}
private:
std::string path_;
int max_slots_ = 0;
std::size_t slot_bytes_ = 0;
Packer packer_;
bool flush_on_store_ = false;
// mutable because Read() is const but needs to seek/read
mutable std::fstream file_;
mutable std::vector<unsigned char> scratch_;
std::vector<unsigned char> in_use_;
std::vector<int> free_;
static Header MakeHeader_(std::uint64_t slot_bytes, std::uint64_t max_slots)
{
Header h;
std::memset(&h, 0, sizeof(h));
h.magic[0] = 'M'; h.magic[1] = 'F'; h.magic[2] = 'C'; h.magic[3] = 'K';
h.magic[4] = 'P'; h.magic[5] = 'T'; h.magic[6] = 'F'; h.magic[7] = 'S';
h.version = 1;
h.slot_bytes = slot_bytes;
h.max_slots = max_slots;
return h;
}
void Open_(bool truncate)
{
const std::ios::openmode mode =
std::ios::binary | std::ios::in | std::ios::out | (truncate ? std::ios::trunc : (std::ios::openmode)0);
file_.open(path_.c_str(), mode);
MFEM_VERIFY(file_.is_open(), "FixedSlotFileCheckpointStorage: failed to open file.");
const Header expected = MakeHeader_((std::uint64_t)slot_bytes_, (std::uint64_t)max_slots_);
if (truncate)
{
// Write header
file_.seekp(0, std::ios::beg);
file_.write(reinterpret_cast<const char*>(&expected), sizeof(expected));
MFEM_VERIFY(file_.good(), "FixedSlotFileCheckpointStorage: header write failed.");
// Pre-size file to: header + max_slots*slot_bytes
const std::uint64_t total = (std::uint64_t)sizeof(Header)
+ (std::uint64_t)max_slots_ * (std::uint64_t)slot_bytes_;
MFEM_VERIFY(total > 0, "FixedSlotFileCheckpointStorage: invalid total file size.");
file_.seekp((std::streamoff)(total - 1), std::ios::beg);
MFEM_VERIFY(file_.good(), "FixedSlotFileCheckpointStorage: seekp for resize failed.");
const char zero = 0;
file_.write(&zero, 1);
MFEM_VERIFY(file_.good(), "FixedSlotFileCheckpointStorage: resize write failed.");
file_.flush();
}
else
{
// Validate existing header
Header got;
file_.seekg(0, std::ios::beg);
file_.read(reinterpret_cast<char*>(&got), sizeof(got));
MFEM_VERIFY(file_.good(), "FixedSlotFileCheckpointStorage: header read failed.");
MFEM_VERIFY(std::memcmp(got.magic, expected.magic, 8) == 0,
"FixedSlotFileCheckpointStorage: magic mismatch.");
MFEM_VERIFY(got.version == expected.version,
"FixedSlotFileCheckpointStorage: version mismatch.");
MFEM_VERIFY(got.slot_bytes == expected.slot_bytes,
"FixedSlotFileCheckpointStorage: slot_bytes mismatch.");
MFEM_VERIFY(got.max_slots == expected.max_slots,
"FixedSlotFileCheckpointStorage: max_slots mismatch.");
}
}
std::uint64_t SlotOffset_(int slot) const
{
return (std::uint64_t)sizeof(Header) + (std::uint64_t)slot * (std::uint64_t)slot_bytes_;
}
};
} // namespace mfem
#endif // MFEM_FIXED_SLOT_CHECKPOINT_STORAGE_HPP
+694
View File
@@ -0,0 +1,694 @@
#include "mfem.hpp"
#include <cstdint>
#include <cstring>
#include <fstream>
#include <string>
#include <vector>
#include <algorithm>
using mfem::out;
// -------------------------------
// Fixed-slot storage: Memory
// -------------------------------
class FixedSlotMemoryStorage
{
private:
int max_slots_ = 0;
size_t slot_bytes_ = 0;
std::vector<uint8_t> data_;
public:
FixedSlotMemoryStorage(int max_slots, size_t slot_bytes)
: max_slots_(max_slots), slot_bytes_(slot_bytes),
data_(size_t(max_slots)*slot_bytes, 0)
{
MFEM_VERIFY(max_slots_ > 0, "max_slots must be > 0");
MFEM_VERIFY(slot_bytes_ > 0, "slot_bytes must be > 0");
}
int MaxSlots() const { return max_slots_; }
size_t SlotBytes() const { return slot_bytes_; }
void Save(int slot, const uint8_t *src, size_t bytes)
{
MFEM_VERIFY(0 <= slot && slot < max_slots_, "Save: slot out of range");
MFEM_VERIFY(bytes == slot_bytes_, "Save: bytes mismatch");
MFEM_VERIFY(src != nullptr, "Save: src is null");
std::memcpy(data_.data() + size_t(slot)*slot_bytes_, src, slot_bytes_);
}
void Load(int slot, uint8_t *dst, size_t bytes) const
{
MFEM_VERIFY(0 <= slot && slot < max_slots_, "Load: slot out of range");
MFEM_VERIFY(bytes == slot_bytes_, "Load: bytes mismatch");
MFEM_VERIFY(dst != nullptr, "Load: dst is null");
std::memcpy(dst, data_.data() + size_t(slot)*slot_bytes_, slot_bytes_);
}
};
// -------------------------------
// Fixed-slot storage: Single file
// -------------------------------
class FixedSlotFileStorage
{
private:
int max_slots_ = 0;
size_t slot_bytes_ = 0;
std::string filename_;
mutable std::fstream file_;
std::streamoff Offset(int slot) const
{
return std::streamoff(size_t(slot) * slot_bytes_);
}
public:
FixedSlotFileStorage(const std::string &filename,
int max_slots,
size_t slot_bytes)
: max_slots_(max_slots), slot_bytes_(slot_bytes), filename_(filename)
{
MFEM_VERIFY(max_slots_ > 0, "max_slots must be > 0");
MFEM_VERIFY(slot_bytes_ > 0, "slot_bytes must be > 0");
MFEM_VERIFY(!filename_.empty(), "filename must not be empty");
// Create/truncate file and size it.
{
std::ofstream ofs(filename_, std::ios::binary | std::ios::trunc);
MFEM_VERIFY(ofs.good(), "Failed to create checkpoint file");
const size_t total = size_t(max_slots_) * slot_bytes_;
if (total > 0)
{
ofs.seekp(std::streamoff(total - 1));
char zero = 0;
ofs.write(&zero, 1);
}
}
file_.open(filename_, std::ios::binary | std::ios::in | std::ios::out);
MFEM_VERIFY(file_.good(), "Failed to open checkpoint file");
}
~FixedSlotFileStorage()
{
if (file_.is_open()) { file_.close(); }
}
int MaxSlots() const { return max_slots_; }
size_t SlotBytes() const { return slot_bytes_; }
void Save(int slot, const uint8_t *src, size_t bytes)
{
MFEM_VERIFY(0 <= slot && slot < max_slots_, "Save: slot out of range");
MFEM_VERIFY(bytes == slot_bytes_, "Save: bytes mismatch");
MFEM_VERIFY(src != nullptr, "Save: src is null");
file_.seekp(Offset(slot));
MFEM_VERIFY(file_.good(), "Save: seekp failed");
file_.write(reinterpret_cast<const char*>(src), std::streamsize(slot_bytes_));
MFEM_VERIFY(file_.good(), "Save: write failed");
file_.flush();
}
void Load(int slot, uint8_t *dst, size_t bytes) const
{
MFEM_VERIFY(0 <= slot && slot < max_slots_, "Load: slot out of range");
MFEM_VERIFY(bytes == slot_bytes_, "Load: bytes mismatch");
MFEM_VERIFY(dst != nullptr, "Load: dst is null");
file_.seekg(Offset(slot));
MFEM_VERIFY(file_.good(), "Load: seekg failed");
file_.read(reinterpret_cast<char*>(dst), std::streamsize(slot_bytes_));
MFEM_VERIFY(file_.good(), "Load: read failed");
}
};
// -------------------------------
// REVOLVE controller (actions)
// (Transcribed from revolve.c used in ADOL-C; Algorithm 799.)
// -------------------------------
enum class RevolveAction
{
advance,
takeshot,
restore,
firsturn,
youturn,
terminate
};
class RevolveController
{
private:
int snaps_ = 0;
int check_ = -1;
int capo_ = 0;
int fine_ = 0;
int turn_ = 0;
int oldfine_ = 0;
int oldsnaps_ = 0;
// Stack of checkpoint times, indexed by 'check_'.
std::vector<int> ch_;
public:
RevolveController() = default;
RevolveController(int snaps, int capo0, int fine0)
: snaps_(snaps), check_(-1), capo_(capo0), fine_(fine0),
turn_(0), oldfine_(fine0), oldsnaps_(snaps), ch_(snaps, 0)
{
MFEM_VERIFY(snaps_ > 0, "REVOLVE snaps must be > 0");
MFEM_VERIFY(capo_ <= fine_, "REVOLVE: capo must be <= fine");
// Match revolve.c initialization behavior.
if (check_ == -1 && capo_ < fine_)
{
turn_ = 0;
ch_[0] = capo_ - 1;
}
}
int Snaps() const { return snaps_; }
int Check() const { return check_; }
int Capo() const { return capo_; }
int Fine() const { return fine_; }
const std::vector<int>& CheckpointTimes() const { return ch_; }
RevolveAction Next()
{
MFEM_VERIFY(!(check_ < -1), "REVOLVE: check < -1");
MFEM_VERIFY(!(capo_ > fine_), "REVOLVE: capo > fine");
if ((check_ == -1) && (capo_ < fine_))
{
turn_ = 0;
ch_[0] = capo_ - 1;
}
const int diff = fine_ - capo_;
switch (diff)
{
case 0:
{
// Terminate or restore to next checkpoint on the stack.
if (check_ == -1 || capo_ == ch_[0])
{
check_ -= 1; // mirror revolve.c behavior
return RevolveAction::terminate;
}
else
{
capo_ = ch_[check_];
oldfine_ = fine_;
return RevolveAction::restore;
}
}
case 1:
{
// One adjoint step available.
fine_ -= 1;
if (check_ >= 0 && ch_[check_] == capo_) { check_ -= 1; }
if (turn_ == 0)
{
turn_ = 1;
oldfine_ = fine_;
return RevolveAction::firsturn;
}
else
{
oldfine_ = fine_;
return RevolveAction::youturn;
}
}
default:
{
// diff > 1
if (check_ == -1 || ch_[check_] != capo_)
{
// Take a new checkpoint at current capo.
check_ += 1;
MFEM_VERIFY(check_ + 1 <= snaps_, "REVOLVE: exceeded snaps");
ch_[check_] = capo_;
oldfine_ = fine_;
return RevolveAction::takeshot;
}
else
{
// Advance capo forward within (capo, fine).
// This follows the binomial logic in revolve.c.
MFEM_VERIFY(!((oldfine_ < fine_) && (snaps_ == check_ + 1)),
"REVOLVE: fine increased unexpectedly with full stack");
const int oldcapo = capo_;
const int ds = snaps_ - check_;
MFEM_VERIFY(ds >= 1, "REVOLVE: ds < 1");
int reps = 0;
long long range = 1;
while (range < (fine_ - capo_))
{
reps += 1;
// range = range*(reps+ds)/reps (integer arithmetic)
range = range * (reps + ds) / reps;
}
MFEM_VERIFY(reps >= 1, "REVOLVE: reps < 1");
// Binomial helper values (integer)
const long long bino1 = range * reps / (ds + reps);
const long long bino2 = (ds > 1) ? (bino1 * ds / (ds + reps - 1)) : 1;
const long long bino3 =
(ds == 1) ? 0 :
(ds > 2) ? (bino2 * (ds - 1) / (ds + reps - 2)) : 1;
const long long bino4 = bino2 * (reps - 1) / ds;
const long long bino5 =
(ds < 3) ? 0 :
(ds > 3) ? (bino3 * (ds - 2) / reps) : 1;
// Kowarz "new version": keep l^ as small as possible
const long long bino6 = bino1 * ds / reps;
const long long gap = fine_ - capo_;
if (gap <= bino1 + bino3)
{
capo_ += int(bino4);
}
else if (gap < bino1 + bino2)
{
capo_ = fine_ - int(bino2 + bino3);
}
else if (gap <= bino1 + bino2 + bino5)
{
capo_ += int(bino1 - bino3);
}
else
{
capo_ = fine_ - int(bino6);
}
if (capo_ == oldcapo) { capo_ = oldcapo + 1; }
oldfine_ = fine_;
return RevolveAction::advance;
}
}
}
}
};
// -------------------------------
// Fixed-step REVOLVE checkpointing
// -------------------------------
template <typename StorageT>
class FixedStepRevolveCheckpointing
{
public:
struct Shot
{
int time = 0; // state index
int slot = 0; // checkpoint slot index
};
private:
int num_steps_ = 0; // total number of primal steps (0..num_steps)
int num_checkpoints_ = 0; // REVOLVE "snaps"
size_t snapshot_bytes_ = 0;
StorageT *storage_ = nullptr;
std::vector<Shot> forward_shots_;
int forward_shot_cursor_ = 0;
// Controller state at the beginning of the reverse sweep (pre-firsturn).
RevolveController ctrl_init_;
RevolveController ctrl_;
// Two scratch buffers for (de)serialization.
std::vector<uint8_t> io_buf_;
std::vector<uint8_t> prefinal_buf_;
bool prefinal_valid_ = false;
// Reverse sweep bookkeeping.
bool reverse_started_ = false;
int u_work_time_ = -1;
public:
FixedStepRevolveCheckpointing(int num_steps,
int num_checkpoints,
size_t snapshot_bytes,
StorageT &storage)
: num_steps_(num_steps),
num_checkpoints_(num_checkpoints),
snapshot_bytes_(snapshot_bytes),
storage_(&storage),
io_buf_(snapshot_bytes, 0),
prefinal_buf_(snapshot_bytes, 0)
{
MFEM_VERIFY(num_steps_ >= 0, "num_steps must be >= 0");
MFEM_VERIFY(num_checkpoints_ > 0, "num_checkpoints must be > 0");
MFEM_VERIFY(snapshot_bytes_ > 0, "snapshot_bytes must be > 0");
MFEM_VERIFY(storage_->MaxSlots() == num_checkpoints_,
"Storage MaxSlots() must match num_checkpoints");
MFEM_VERIFY(storage_->SlotBytes() == snapshot_bytes_,
"Storage SlotBytes() must match snapshot_bytes");
BuildForwardPlanAndInitialControllerState();
Reset();
}
void Reset()
{
forward_shot_cursor_ = 0;
prefinal_valid_ = false;
reverse_started_ = false;
u_work_time_ = -1;
ctrl_ = ctrl_init_;
}
int NumSteps() const { return num_steps_; }
int NumCheckpoints() const { return num_checkpoints_; }
size_t SnapshotBytes() const { return snapshot_bytes_; }
// ForwardStep: called for i=0..num_steps-1
template <typename State, typename PrimalStep, typename MakeSnapshot>
void ForwardStep(int i,
State &u,
PrimalStep &&primal_step,
MakeSnapshot &&make_snapshot)
{
MFEM_VERIFY(0 <= i && i < num_steps_, "ForwardStep: i out of range");
// Take any planned shots at time i (before advancing).
while (forward_shot_cursor_ < (int)forward_shots_.size() &&
forward_shots_[forward_shot_cursor_].time == i)
{
const int slot = forward_shots_[forward_shot_cursor_].slot;
make_snapshot(u, io_buf_.data(), snapshot_bytes_);
storage_->Save(slot, io_buf_.data(), snapshot_bytes_);
forward_shot_cursor_++;
}
// Cache u_{num_steps-1} so reverse can start there (REVOLVE expects capo=num_steps-1).
if (i == num_steps_ - 1)
{
make_snapshot(u, prefinal_buf_.data(), snapshot_bytes_);
prefinal_valid_ = true;
}
// Advance one step.
primal_step(i, u);
}
// BackwardStep: called for i=num_steps-1..0
template <typename State, typename AdjointState,
typename PrimalStep, typename AdjointStep,
typename MakeSnapshot, typename RestoreSnapshot>
void BackwardStep(int i,
AdjointState &lambda,
State &u_work,
PrimalStep &&primal_step,
AdjointStep &&adjoint_step,
MakeSnapshot &&make_snapshot,
RestoreSnapshot &&restore_snapshot)
{
MFEM_VERIFY(0 <= i && i < num_steps_, "BackwardStep: i out of range");
if (!reverse_started_)
{
MFEM_VERIFY(prefinal_valid_ || num_steps_ == 0,
"Reverse started but prefinal state was not captured. "
"Did you run the forward loop through i=num_steps-1?");
ctrl_ = ctrl_init_;
reverse_started_ = true;
if (num_steps_ > 0)
{
// Restore u_{num_steps-1} into u_work.
restore_snapshot(u_work, prefinal_buf_.data(), snapshot_bytes_);
u_work_time_ = num_steps_ - 1;
}
}
// Execute controller actions until we perform exactly one adjoint step.
for (;;)
{
RevolveAction act = ctrl_.Next();
switch (act)
{
case RevolveAction::takeshot:
{
const int slot = ctrl_.Check();
MFEM_VERIFY(u_work_time_ == ctrl_.Capo(),
"takeshot: u_work_time must equal capo");
make_snapshot(u_work, io_buf_.data(), snapshot_bytes_);
storage_->Save(slot, io_buf_.data(), snapshot_bytes_);
break;
}
case RevolveAction::restore:
{
const int slot = ctrl_.Check();
storage_->Load(slot, io_buf_.data(), snapshot_bytes_);
restore_snapshot(u_work, io_buf_.data(), snapshot_bytes_);
u_work_time_ = ctrl_.Capo();
break;
}
case RevolveAction::advance:
{
const int target = ctrl_.Capo();
MFEM_VERIFY(u_work_time_ >= 0, "advance: u_work_time not initialized");
MFEM_VERIFY(target >= u_work_time_, "advance: target < current time");
for (int t = u_work_time_; t < target; ++t)
{
primal_step(t, u_work);
}
u_work_time_ = target;
break;
}
case RevolveAction::firsturn:
case RevolveAction::youturn:
{
// After firsturn/youturn, ctrl_.Fine() has been decremented and equals ctrl_.Capo().
const int step = ctrl_.Fine();
MFEM_VERIFY(step == ctrl_.Capo(), "youturn: fine != capo");
MFEM_VERIFY(step == u_work_time_, "youturn: u_work_time != step");
MFEM_VERIFY(step == i, "BackwardStep called with i that doesn't match REVOLVE schedule");
// One adjoint step.
adjoint_step(step, u_work, lambda);
return;
}
case RevolveAction::terminate:
{
MFEM_ABORT("REVOLVE terminated early: BackwardStep called after completion?");
break;
}
}
}
}
private:
void BuildForwardPlanAndInitialControllerState()
{
forward_shots_.clear();
RevolveController sim(num_checkpoints_, /*capo=*/0, /*fine=*/num_steps_);
// Forward-plan phase: only TAKESHOT/ADVANCE should occur while fine-capo>1.
while (sim.Fine() - sim.Capo() > 1)
{
RevolveAction a = sim.Next();
if (a == RevolveAction::takeshot)
{
forward_shots_.push_back({sim.Capo(), sim.Check()});
}
else if (a == RevolveAction::advance)
{
// nothing to record; capo moved forward inside sim
}
else
{
MFEM_ABORT("Unexpected REVOLVE action during forward planning phase");
}
}
// At this point, sim is in the pre-firsturn state (capo = num_steps-1, fine = num_steps)
// for num_steps>=1. For num_steps<=1, fine-capo<=1 from the start.
ctrl_init_ = sim;
}
};
// ============================================================
// Mini example 1: scalar state (double)
// u_{n+1} = factor * u_n
// J = 0.5 (u_N - target)^2 => lambda_N = (u_N - target)
// lambda_n = factor * lambda_{n+1}
// ============================================================
static void RunScalarExample(int Nsteps, int Ncheck, double factor)
{
out << "\n--- Scalar REVOLVE example ---\n";
out << "Nsteps=" << Nsteps << ", Ncheckpoints=" << Ncheck << ", factor=" << factor << "\n";
const size_t snap_bytes = sizeof(double);
FixedSlotMemoryStorage storage(Ncheck, snap_bytes);
FixedStepRevolveCheckpointing<FixedSlotMemoryStorage> ckpt(Nsteps, Ncheck, snap_bytes, storage);
auto make_snapshot = [](const double &u, uint8_t *outb, size_t bytes)
{
MFEM_VERIFY(bytes == sizeof(double), "scalar snapshot size mismatch");
std::memcpy(outb, &u, sizeof(double));
};
auto restore_snapshot = [](double &u, const uint8_t *inb, size_t bytes)
{
MFEM_VERIFY(bytes == sizeof(double), "scalar snapshot size mismatch");
std::memcpy(&u, inb, sizeof(double));
};
auto primal_step = [factor](int /*step*/, double &u)
{
u *= factor;
};
auto adjoint_step = [factor](int /*step*/, const double & /*u_step*/, double &lambda)
{
lambda *= factor;
};
const double target = 2.0;
// Forward
double u = 1.0;
for (int i = 0; i < Nsteps; ++i)
{
ckpt.ForwardStep(i, u, primal_step, make_snapshot);
}
const double uN = u;
// Init adjoint at final state
double lambda = (uN - target);
// Reverse
double u_work = 0.0; // will be overwritten by restore_snapshot
for (int i = Nsteps - 1; i >= 0; --i)
{
ckpt.BackwardStep(i, lambda, u_work,
primal_step, adjoint_step,
make_snapshot, restore_snapshot);
}
// Analytic lambda_0 = factor^Nsteps * (uN - target)
double factorN = 1.0;
for (int k = 0; k < Nsteps; ++k) { factorN *= factor; }
const double lambda0_exact = factorN * (uN - target);
out << "uN = " << uN << "\n";
out << "lambda0 (computed) = " << lambda << "\n";
out << "lambda0 (exact) = " << lambda0_exact << "\n";
}
// ============================================================
// Mini example 2: mfem::Vector state
// u_{n+1} = factor * u_n (elementwise scalar multiply)
// J = 0.5 ||u_N - target||^2 => lambda_N = (u_N - target)
// lambda_n = factor * lambda_{n+1}
// ============================================================
static void RunVectorExample(int Nsteps, int Ncheck, int dim, double factor)
{
out << "\n--- mfem::Vector REVOLVE example ---\n";
out << "Nsteps=" << Nsteps << ", Ncheckpoints=" << Ncheck
<< ", dim=" << dim << ", factor=" << factor << "\n";
MFEM_VERIFY(dim > 0, "dim must be > 0");
const size_t snap_bytes = sizeof(double) * size_t(dim);
FixedSlotMemoryStorage storage(Ncheck, snap_bytes);
FixedStepRevolveCheckpointing<FixedSlotMemoryStorage> ckpt(Nsteps, Ncheck, snap_bytes, storage);
auto make_snapshot = [](const mfem::Vector &u, uint8_t *outb, size_t bytes)
{
MFEM_VERIFY(bytes == sizeof(double) * size_t(u.Size()), "Vector snapshot size mismatch");
std::memcpy(outb, u.GetData(), bytes);
};
auto restore_snapshot = [](mfem::Vector &u, const uint8_t *inb, size_t bytes)
{
MFEM_VERIFY(bytes == sizeof(double) * size_t(u.Size()), "Vector snapshot size mismatch");
std::memcpy(u.GetData(), inb, bytes);
};
auto primal_step = [factor](int /*step*/, mfem::Vector &u)
{
u *= factor;
};
auto adjoint_step = [factor](int /*step*/, const mfem::Vector & /*u_step*/, mfem::Vector &lambda)
{
lambda *= factor;
};
mfem::Vector u(dim);
for (int k = 0; k < dim; ++k) { u[k] = 1.0 + 0.1*k; }
mfem::Vector target(dim);
target = 2.0;
// Forward
for (int i = 0; i < Nsteps; ++i)
{
ckpt.ForwardStep(i, u, primal_step, make_snapshot);
}
mfem::Vector uN(u); // copy for reporting
// Init adjoint: lambda_N = uN - target
mfem::Vector lambda(dim);
lambda = u;
lambda -= target;
// Reverse
mfem::Vector u_work(dim);
u_work = 0.0;
for (int i = Nsteps - 1; i >= 0; --i)
{
ckpt.BackwardStep(i, lambda, u_work,
primal_step, adjoint_step,
make_snapshot, restore_snapshot);
}
out << "||uN||_2 = " << uN.Norml2() << "\n";
out << "||lambda0||_2 = " << lambda.Norml2() << "\n";
}
int main(int argc, char *argv[])
{
mfem::Device device("cpu");
device.Print();
int Nsteps = 20;
int Ncheck = 3;
int dim = 5;
double factor = 1.05;
mfem::OptionsParser args(argc, argv);
args.AddOption(&Nsteps, "-n", "--num-steps", "Number of primal steps.");
args.AddOption(&Ncheck, "-s", "--num-checkpoints", "Number of checkpoints (snaps).");
args.AddOption(&dim, "-d", "--dim", "Vector dimension for the mfem::Vector example.");
args.AddOption(&factor, "-f", "--factor", "Scalar factor in the toy update u_{n+1}=f*u_n.");
args.Parse();
if (!args.Good())
{
args.PrintUsage(out);
return 1;
}
args.PrintOptions(out);
RunScalarExample(Nsteps, Ncheck, factor);
RunVectorExample(Nsteps, Ncheck, dim, factor);
return 0;
}
@@ -0,0 +1,535 @@
#ifndef MFEM_REVOLVE_CHECKPOINTING_HPP
#define MFEM_REVOLVE_CHECKPOINTING_HPP
#include "mfem.hpp"
#include <cstdint>
#include <cstring>
#include <fstream>
#include <string>
#include <vector>
#include <algorithm>
using mfem::out;
// -------------------------------
// Fixed-slot storage: Memory
// -------------------------------
class FixedSlotMemoryStorage
{
private:
int max_slots_ = 0;
size_t slot_bytes_ = 0;
std::vector<uint8_t> data_;
public:
FixedSlotMemoryStorage(int max_slots, size_t slot_bytes)
: max_slots_(max_slots), slot_bytes_(slot_bytes),
data_(size_t(max_slots)*slot_bytes, 0)
{
MFEM_VERIFY(max_slots_ > 0, "max_slots must be > 0");
MFEM_VERIFY(slot_bytes_ > 0, "slot_bytes must be > 0");
}
int MaxSlots() const { return max_slots_; }
size_t SlotBytes() const { return slot_bytes_; }
void Save(int slot, const uint8_t *src, size_t bytes)
{
MFEM_VERIFY(0 <= slot && slot < max_slots_, "Save: slot out of range");
MFEM_VERIFY(bytes == slot_bytes_, "Save: bytes mismatch");
MFEM_VERIFY(src != nullptr, "Save: src is null");
std::memcpy(data_.data() + size_t(slot)*slot_bytes_, src, slot_bytes_);
}
void Load(int slot, uint8_t *dst, size_t bytes) const
{
MFEM_VERIFY(0 <= slot && slot < max_slots_, "Load: slot out of range");
MFEM_VERIFY(bytes == slot_bytes_, "Load: bytes mismatch");
MFEM_VERIFY(dst != nullptr, "Load: dst is null");
std::memcpy(dst, data_.data() + size_t(slot)*slot_bytes_, slot_bytes_);
}
};
// -------------------------------
// Fixed-slot storage: Single file
// -------------------------------
class FixedSlotFileStorage
{
private:
int max_slots_ = 0;
size_t slot_bytes_ = 0;
std::string filename_;
mutable std::fstream file_;
std::streamoff Offset(int slot) const
{
return std::streamoff(size_t(slot) * slot_bytes_);
}
public:
FixedSlotFileStorage(const std::string &filename,
int max_slots,
size_t slot_bytes)
: max_slots_(max_slots), slot_bytes_(slot_bytes), filename_(filename)
{
MFEM_VERIFY(max_slots_ > 0, "max_slots must be > 0");
MFEM_VERIFY(slot_bytes_ > 0, "slot_bytes must be > 0");
MFEM_VERIFY(!filename_.empty(), "filename must not be empty");
// Create/truncate file and size it.
{
std::ofstream ofs(filename_, std::ios::binary | std::ios::trunc);
MFEM_VERIFY(ofs.good(), "Failed to create checkpoint file");
const size_t total = size_t(max_slots_) * slot_bytes_;
if (total > 0)
{
ofs.seekp(std::streamoff(total - 1));
char zero = 0;
ofs.write(&zero, 1);
}
}
file_.open(filename_, std::ios::binary | std::ios::in | std::ios::out);
MFEM_VERIFY(file_.good(), "Failed to open checkpoint file");
}
~FixedSlotFileStorage()
{
if (file_.is_open()) { file_.close(); }
}
int MaxSlots() const { return max_slots_; }
size_t SlotBytes() const { return slot_bytes_; }
void Save(int slot, const uint8_t *src, size_t bytes)
{
MFEM_VERIFY(0 <= slot && slot < max_slots_, "Save: slot out of range");
MFEM_VERIFY(bytes == slot_bytes_, "Save: bytes mismatch");
MFEM_VERIFY(src != nullptr, "Save: src is null");
file_.seekp(Offset(slot));
MFEM_VERIFY(file_.good(), "Save: seekp failed");
file_.write(reinterpret_cast<const char*>(src), std::streamsize(slot_bytes_));
MFEM_VERIFY(file_.good(), "Save: write failed");
file_.flush();
}
void Load(int slot, uint8_t *dst, size_t bytes) const
{
MFEM_VERIFY(0 <= slot && slot < max_slots_, "Load: slot out of range");
MFEM_VERIFY(bytes == slot_bytes_, "Load: bytes mismatch");
MFEM_VERIFY(dst != nullptr, "Load: dst is null");
file_.seekg(Offset(slot));
MFEM_VERIFY(file_.good(), "Load: seekg failed");
file_.read(reinterpret_cast<char*>(dst), std::streamsize(slot_bytes_));
MFEM_VERIFY(file_.good(), "Load: read failed");
}
};
// -------------------------------
// REVOLVE controller (actions)
// (Transcribed from revolve.c used in ADOL-C; Algorithm 799.)
// -------------------------------
enum class RevolveAction
{
advance,
takeshot,
restore,
firsturn,
youturn,
terminate
};
class RevolveController
{
private:
int snaps_ = 0;
int check_ = -1;
int capo_ = 0;
int fine_ = 0;
int turn_ = 0;
int oldfine_ = 0;
int oldsnaps_ = 0;
// Stack of checkpoint times, indexed by 'check_'.
std::vector<int> ch_;
public:
RevolveController() = default;
RevolveController(int snaps, int capo0, int fine0)
: snaps_(snaps), check_(-1), capo_(capo0), fine_(fine0),
turn_(0), oldfine_(fine0), oldsnaps_(snaps), ch_(snaps, 0)
{
MFEM_VERIFY(snaps_ > 0, "REVOLVE snaps must be > 0");
MFEM_VERIFY(capo_ <= fine_, "REVOLVE: capo must be <= fine");
// Match revolve.c initialization behavior.
if (check_ == -1 && capo_ < fine_)
{
turn_ = 0;
ch_[0] = capo_ - 1;
}
}
int Snaps() const { return snaps_; }
int Check() const { return check_; }
int Capo() const { return capo_; }
int Fine() const { return fine_; }
const std::vector<int>& CheckpointTimes() const { return ch_; }
RevolveAction Next()
{
MFEM_VERIFY(!(check_ < -1), "REVOLVE: check < -1");
MFEM_VERIFY(!(capo_ > fine_), "REVOLVE: capo > fine");
if ((check_ == -1) && (capo_ < fine_))
{
turn_ = 0;
ch_[0] = capo_ - 1;
}
const int diff = fine_ - capo_;
switch (diff)
{
case 0:
{
// Terminate or restore to next checkpoint on the stack.
if (check_ == -1 || capo_ == ch_[0])
{
check_ -= 1; // mirror revolve.c behavior
return RevolveAction::terminate;
}
else
{
capo_ = ch_[check_];
oldfine_ = fine_;
return RevolveAction::restore;
}
}
case 1:
{
// One adjoint step available.
fine_ -= 1;
if (check_ >= 0 && ch_[check_] == capo_) { check_ -= 1; }
if (turn_ == 0)
{
turn_ = 1;
oldfine_ = fine_;
return RevolveAction::firsturn;
}
else
{
oldfine_ = fine_;
return RevolveAction::youturn;
}
}
default:
{
// diff > 1
if (check_ == -1 || ch_[check_] != capo_)
{
// Take a new checkpoint at current capo.
check_ += 1;
MFEM_VERIFY(check_ + 1 <= snaps_, "REVOLVE: exceeded snaps");
ch_[check_] = capo_;
oldfine_ = fine_;
return RevolveAction::takeshot;
}
else
{
// Advance capo forward within (capo, fine).
// This follows the binomial logic in revolve.c.
MFEM_VERIFY(!((oldfine_ < fine_) && (snaps_ == check_ + 1)),
"REVOLVE: fine increased unexpectedly with full stack");
const int oldcapo = capo_;
const int ds = snaps_ - check_;
MFEM_VERIFY(ds >= 1, "REVOLVE: ds < 1");
int reps = 0;
long long range = 1;
while (range < (fine_ - capo_))
{
reps += 1;
// range = range*(reps+ds)/reps (integer arithmetic)
range = range * (reps + ds) / reps;
}
MFEM_VERIFY(reps >= 1, "REVOLVE: reps < 1");
// Binomial helper values (integer)
const long long bino1 = range * reps / (ds + reps);
const long long bino2 = (ds > 1) ? (bino1 * ds / (ds + reps - 1)) : 1;
const long long bino3 =
(ds == 1) ? 0 :
(ds > 2) ? (bino2 * (ds - 1) / (ds + reps - 2)) : 1;
const long long bino4 = bino2 * (reps - 1) / ds;
const long long bino5 =
(ds < 3) ? 0 :
(ds > 3) ? (bino3 * (ds - 2) / reps) : 1;
// Kowarz "new version": keep l^ as small as possible
const long long bino6 = bino1 * ds / reps;
const long long gap = fine_ - capo_;
if (gap <= bino1 + bino3)
{
capo_ += int(bino4);
}
else if (gap < bino1 + bino2)
{
capo_ = fine_ - int(bino2 + bino3);
}
else if (gap <= bino1 + bino2 + bino5)
{
capo_ += int(bino1 - bino3);
}
else
{
capo_ = fine_ - int(bino6);
}
if (capo_ == oldcapo) { capo_ = oldcapo + 1; }
oldfine_ = fine_;
return RevolveAction::advance;
}
}
}
}
};
// -------------------------------
// Fixed-step REVOLVE checkpointing
// -------------------------------
template <typename StorageT>
class FixedStepRevolveCheckpointing
{
public:
struct Shot
{
int time = 0; // state index
int slot = 0; // checkpoint slot index
};
private:
int num_steps_ = 0; // total number of primal steps (0..num_steps)
int num_checkpoints_ = 0; // REVOLVE "snaps"
size_t snapshot_bytes_ = 0;
StorageT *storage_ = nullptr;
std::vector<Shot> forward_shots_;
int forward_shot_cursor_ = 0;
// Controller state at the beginning of the reverse sweep (pre-firsturn).
RevolveController ctrl_init_;
RevolveController ctrl_;
// Two scratch buffers for (de)serialization.
std::vector<uint8_t> io_buf_;
std::vector<uint8_t> prefinal_buf_;
bool prefinal_valid_ = false;
// Reverse sweep bookkeeping.
bool reverse_started_ = false;
int u_work_time_ = -1;
public:
FixedStepRevolveCheckpointing(int num_steps,
int num_checkpoints,
size_t snapshot_bytes,
StorageT &storage)
: num_steps_(num_steps),
num_checkpoints_(num_checkpoints),
snapshot_bytes_(snapshot_bytes),
storage_(&storage),
io_buf_(snapshot_bytes, 0),
prefinal_buf_(snapshot_bytes, 0)
{
MFEM_VERIFY(num_steps_ >= 0, "num_steps must be >= 0");
MFEM_VERIFY(num_checkpoints_ > 0, "num_checkpoints must be > 0");
MFEM_VERIFY(snapshot_bytes_ > 0, "snapshot_bytes must be > 0");
MFEM_VERIFY(storage_->MaxSlots() == num_checkpoints_,
"Storage MaxSlots() must match num_checkpoints");
MFEM_VERIFY(storage_->SlotBytes() == snapshot_bytes_,
"Storage SlotBytes() must match snapshot_bytes");
BuildForwardPlanAndInitialControllerState();
Reset();
}
void Reset()
{
forward_shot_cursor_ = 0;
prefinal_valid_ = false;
reverse_started_ = false;
u_work_time_ = -1;
ctrl_ = ctrl_init_;
}
int NumSteps() const { return num_steps_; }
int NumCheckpoints() const { return num_checkpoints_; }
size_t SnapshotBytes() const { return snapshot_bytes_; }
// ForwardStep: called for i=0..num_steps-1
template <typename State, typename PrimalStep, typename MakeSnapshot>
void ForwardStep(int i,
State &u,
PrimalStep &&primal_step,
MakeSnapshot &&make_snapshot)
{
MFEM_VERIFY(0 <= i && i < num_steps_, "ForwardStep: i out of range");
// Take any planned shots at time i (before advancing).
while (forward_shot_cursor_ < (int)forward_shots_.size() &&
forward_shots_[forward_shot_cursor_].time == i)
{
const int slot = forward_shots_[forward_shot_cursor_].slot;
make_snapshot(u, io_buf_.data(), snapshot_bytes_);
storage_->Save(slot, io_buf_.data(), snapshot_bytes_);
forward_shot_cursor_++;
}
// Cache u_{num_steps-1} so reverse can start there (REVOLVE expects capo=num_steps-1).
if (i == num_steps_ - 1)
{
make_snapshot(u, prefinal_buf_.data(), snapshot_bytes_);
prefinal_valid_ = true;
}
// Advance one step.
primal_step(i, u);
}
// BackwardStep: called for i=num_steps-1..0
template <typename State, typename AdjointState,
typename PrimalStep, typename AdjointStep,
typename MakeSnapshot, typename RestoreSnapshot>
void BackwardStep(int i,
AdjointState &lambda,
State &u_work,
PrimalStep &&primal_step,
AdjointStep &&adjoint_step,
MakeSnapshot &&make_snapshot,
RestoreSnapshot &&restore_snapshot)
{
MFEM_VERIFY(0 <= i && i < num_steps_, "BackwardStep: i out of range");
if (!reverse_started_)
{
MFEM_VERIFY(prefinal_valid_ || num_steps_ == 0,
"Reverse started but prefinal state was not captured. "
"Did you run the forward loop through i=num_steps-1?");
ctrl_ = ctrl_init_;
reverse_started_ = true;
if (num_steps_ > 0)
{
// Restore u_{num_steps-1} into u_work.
restore_snapshot(u_work, prefinal_buf_.data(), snapshot_bytes_);
u_work_time_ = num_steps_ - 1;
}
}
// Execute controller actions until we perform exactly one adjoint step.
for (;;)
{
RevolveAction act = ctrl_.Next();
switch (act)
{
case RevolveAction::takeshot:
{
const int slot = ctrl_.Check();
MFEM_VERIFY(u_work_time_ == ctrl_.Capo(),
"takeshot: u_work_time must equal capo");
make_snapshot(u_work, io_buf_.data(), snapshot_bytes_);
storage_->Save(slot, io_buf_.data(), snapshot_bytes_);
break;
}
case RevolveAction::restore:
{
const int slot = ctrl_.Check();
storage_->Load(slot, io_buf_.data(), snapshot_bytes_);
restore_snapshot(u_work, io_buf_.data(), snapshot_bytes_);
u_work_time_ = ctrl_.Capo();
break;
}
case RevolveAction::advance:
{
const int target = ctrl_.Capo();
MFEM_VERIFY(u_work_time_ >= 0, "advance: u_work_time not initialized");
MFEM_VERIFY(target >= u_work_time_, "advance: target < current time");
for (int t = u_work_time_; t < target; ++t)
{
primal_step(t, u_work);
}
u_work_time_ = target;
break;
}
case RevolveAction::firsturn:
case RevolveAction::youturn:
{
// After firsturn/youturn, ctrl_.Fine() has been decremented and equals ctrl_.Capo().
const int step = ctrl_.Fine();
MFEM_VERIFY(step == ctrl_.Capo(), "youturn: fine != capo");
MFEM_VERIFY(step == u_work_time_, "youturn: u_work_time != step");
MFEM_VERIFY(step == i, "BackwardStep called with i that doesn't match REVOLVE schedule");
// One adjoint step.
adjoint_step(step, u_work, lambda);
return;
}
case RevolveAction::terminate:
{
MFEM_ABORT("REVOLVE terminated early: BackwardStep called after completion?");
break;
}
}
}
}
private:
void BuildForwardPlanAndInitialControllerState()
{
forward_shots_.clear();
RevolveController sim(num_checkpoints_, /*capo=*/0, /*fine=*/num_steps_);
// Forward-plan phase: only TAKESHOT/ADVANCE should occur while fine-capo>1.
while (sim.Fine() - sim.Capo() > 1)
{
RevolveAction a = sim.Next();
if (a == RevolveAction::takeshot)
{
forward_shots_.push_back({sim.Capo(), sim.Check()});
}
else if (a == RevolveAction::advance)
{
// nothing to record; capo moved forward inside sim
}
else
{
MFEM_ABORT("Unexpected REVOLVE action during forward planning phase");
}
}
// At this point, sim is in the pre-firsturn state (capo = num_steps-1, fine = num_steps)
// for num_steps>=1. For num_steps<=1, fine-capo<=1 from the start.
ctrl_init_ = sim;
}
};
#endif //MFEM_REVOLVE_CHECKPOINTING_HPP
+357
View File
@@ -0,0 +1,357 @@
#ifndef MFEM_REVOLVE_SEPARATE_FILE_STORAGE_HPP
#define MFEM_REVOLVE_SEPARATE_FILE_STORAGE_HPP
#include "mfem.hpp"
#include <cstdint>
#include <cstdio> // std::remove
#include <fstream>
#include <iomanip>
#include <limits>
#include <sstream>
#include <string>
#include <vector>
#if __cplusplus >= 201703L
#include <filesystem>
#include <system_error>
#else
#if defined(_WIN32)
#include <direct.h> // _mkdir
#else
#include <sys/stat.h> // mkdir
#include <sys/types.h>
#include <errno.h>
#endif
#endif
namespace mfem
{
/**
* @brief Separate-file checkpoint storage for REVOLVE: one file per checkpoint slot.
*
* This storage is "slot indexed":
* slot = 0..max_slots-1
*
* File naming:
* <dir>/<prefix><slot as zero-padded integer><ext>
*
* Example:
* dir="ckpt", prefix="rev_", ext=".bin", max_slots=8
* -> ckpt/rev_00.bin ... ckpt/rev_07.bin
*
* Intended REVOLVE interface:
* int MaxSlots() const;
* size_t SlotBytes() const;
* void Save(int slot, const uint8_t *src, size_t bytes);
* void Load(int slot, uint8_t *dst, size_t bytes) const;
*
* Notes:
* - Slot files are pre-created and resized on construction if truncate_files=true.
* - Save() overwrites the full slot file content from offset 0.
* - Load() reads exactly SlotBytes() from offset 0.
* - If keep_files=false, destructor deletes the slot files.
* - If keep_open=true, keeps N files open (faster, but uses file descriptors).
*/
class RevolveSeparateFileCheckpointStorage
{
public:
RevolveSeparateFileCheckpointStorage(const std::string &directory,
const std::string &prefix,
const std::string &ext,
int max_slots,
std::size_t slot_bytes,
bool create_dir = true,
bool truncate_files = true,
bool keep_files = true,
bool keep_open = false,
bool flush_on_save = false)
: dir_(directory),
prefix_(prefix),
ext_(ext),
max_slots_(max_slots),
slot_bytes_(slot_bytes),
keep_files_(keep_files),
keep_open_(keep_open),
flush_on_save_(flush_on_save)
{
MFEM_VERIFY(max_slots_ > 0, "RevolveSeparateFileCheckpointStorage: max_slots must be > 0.");
MFEM_VERIFY(slot_bytes_ > 0, "RevolveSeparateFileCheckpointStorage: slot_bytes must be > 0.");
MFEM_VERIFY(!prefix_.empty(), "RevolveSeparateFileCheckpointStorage: prefix must not be empty.");
if (create_dir) { EnsureDirectory_(dir_); }
// Precompute slot paths.
slot_paths_.resize((std::size_t)max_slots_);
for (int s = 0; s < max_slots_; ++s)
{
slot_paths_[(std::size_t)s] = MakeSlotPath_(s);
}
if (truncate_files)
{
PrecreateAll_();
}
else
{
// Optional light sanity check: try opening one file. (User may want to reuse existing.)
// We'll rely on Save/Load verification otherwise.
}
if (keep_open_)
{
OpenAll_();
}
}
~RevolveSeparateFileCheckpointStorage()
{
CloseAll_();
if (!keep_files_)
{
RemoveAllFiles_();
}
}
int MaxSlots() const { return max_slots_; }
std::size_t SlotBytes() const { return slot_bytes_; }
/// Return the full path for a given slot (useful for debugging).
const std::string &SlotPath(int slot) const
{
MFEM_VERIFY(0 <= slot && slot < max_slots_, "SlotPath: slot out of range.");
return slot_paths_[(std::size_t)slot];
}
/**
* @brief Save a checkpoint image into slot file.
*
* Requirements:
* - slot in [0, MaxSlots())
* - bytes == SlotBytes()
* - src != nullptr
*/
void Save(int slot, const std::uint8_t *src, std::size_t bytes)
{
MFEM_VERIFY(0 <= slot && slot < max_slots_, "Save: slot out of range.");
MFEM_VERIFY(bytes == slot_bytes_, "Save: bytes mismatch.");
MFEM_VERIFY(src != nullptr, "Save: src is null.");
if (keep_open_)
{
std::fstream &f = files_[(std::size_t)slot];
MFEM_VERIFY(f.is_open(), "Save: file not open (keep_open).");
f.clear();
f.seekp(0, std::ios::beg);
MFEM_VERIFY(f.good(), "Save: seekp failed (keep_open).");
f.write(reinterpret_cast<const char*>(src), (std::streamsize)slot_bytes_);
MFEM_VERIFY(f.good(), "Save: write failed (keep_open).");
if (flush_on_save_) { f.flush(); }
return;
}
// Open on demand
std::fstream f(SlotPath(slot).c_str(),
std::ios::binary | std::ios::in | std::ios::out);
if (!f.is_open())
{
// If missing, create it sized correctly, then reopen.
PrecreateOne_(slot);
f.open(SlotPath(slot).c_str(), std::ios::binary | std::ios::in | std::ios::out);
}
MFEM_VERIFY(f.is_open(), "Save: failed to open slot file.");
f.seekp(0, std::ios::beg);
MFEM_VERIFY(f.good(), "Save: seekp failed.");
f.write(reinterpret_cast<const char*>(src), (std::streamsize)slot_bytes_);
MFEM_VERIFY(f.good(), "Save: write failed.");
if (flush_on_save_) { f.flush(); }
}
/**
* @brief Load a checkpoint image from slot file.
*
* Requirements:
* - slot in [0, MaxSlots())
* - bytes == SlotBytes()
* - dst != nullptr
*/
void Load(int slot, std::uint8_t *dst, std::size_t bytes) const
{
MFEM_VERIFY(0 <= slot && slot < max_slots_, "Load: slot out of range.");
MFEM_VERIFY(bytes == slot_bytes_, "Load: bytes mismatch.");
MFEM_VERIFY(dst != nullptr, "Load: dst is null.");
if (keep_open_)
{
std::fstream &f = files_[(std::size_t)slot];
MFEM_VERIFY(f.is_open(), "Load: file not open (keep_open).");
f.clear();
f.seekg(0, std::ios::beg);
MFEM_VERIFY(f.good(), "Load: seekg failed (keep_open).");
f.read(reinterpret_cast<char*>(dst), (std::streamsize)slot_bytes_);
MFEM_VERIFY(f.good(), "Load: read failed (keep_open).");
return;
}
std::ifstream f(SlotPath(slot).c_str(), std::ios::binary);
MFEM_VERIFY(f.is_open(), "Load: failed to open slot file.");
f.read(reinterpret_cast<char*>(dst), (std::streamsize)slot_bytes_);
MFEM_VERIFY(f.good(), "Load: read failed.");
}
private:
std::string dir_;
std::string prefix_;
std::string ext_;
int max_slots_ = 0;
std::size_t slot_bytes_ = 0;
bool keep_files_ = true;
bool keep_open_ = false;
bool flush_on_save_ = false;
std::vector<std::string> slot_paths_;
mutable std::vector<std::fstream> files_; // only used if keep_open_==true
static std::string JoinPath_(const std::string &dir, const std::string &file)
{
if (dir.empty()) { return file; }
const char last = dir.back();
if (last == '/' || last == '\\') { return dir + file; }
return dir + "/" + file;
}
int SlotDigits_() const
{
int x = max_slots_ - 1;
int d = 1;
while (x >= 10) { x /= 10; ++d; }
return d;
}
std::string MakeSlotFilename_(int slot) const
{
std::ostringstream os;
os << prefix_
<< std::setw(SlotDigits_()) << std::setfill('0') << slot
<< ext_;
return os.str();
}
std::string MakeSlotPath_(int slot) const
{
return JoinPath_(dir_, MakeSlotFilename_(slot));
}
void PrecreateOne_(int slot) const
{
MFEM_VERIFY(0 <= slot && slot < max_slots_, "PrecreateOne: slot out of range.");
// Create/truncate and set file size to slot_bytes_ by writing last byte.
std::ofstream ofs(SlotPath(slot).c_str(), std::ios::binary | std::ios::trunc);
MFEM_VERIFY(ofs.is_open(), "PrecreateOne: failed to create slot file.");
const std::uint64_t sb = (std::uint64_t)slot_bytes_;
MFEM_VERIFY(sb <= (std::uint64_t)std::numeric_limits<std::streamoff>::max(),
"PrecreateOne: slot_bytes too large for streamoff.");
if (sb > 0)
{
ofs.seekp((std::streamoff)(sb - 1), std::ios::beg);
MFEM_VERIFY(ofs.good(), "PrecreateOne: seekp failed.");
const char zero = 0;
ofs.write(&zero, 1);
MFEM_VERIFY(ofs.good(), "PrecreateOne: size write failed.");
}
}
void PrecreateAll_() const
{
for (int slot = 0; slot < max_slots_; ++slot)
{
PrecreateOne_(slot);
}
}
void OpenAll_()
{
files_.resize((std::size_t)max_slots_);
for (int slot = 0; slot < max_slots_; ++slot)
{
std::fstream &f = files_[(std::size_t)slot];
f.open(SlotPath(slot).c_str(), std::ios::binary | std::ios::in | std::ios::out);
if (!f.is_open())
{
PrecreateOne_(slot);
f.open(SlotPath(slot).c_str(), std::ios::binary | std::ios::in | std::ios::out);
}
MFEM_VERIFY(f.is_open(), "OpenAll: failed to open slot file.");
}
}
void CloseAll_()
{
if (!files_.empty())
{
for (auto &f : files_) { if (f.is_open()) { f.close(); } }
files_.clear();
}
}
void RemoveAllFiles_() const
{
for (int slot = 0; slot < max_slots_; ++slot)
{
// Ignore remove errors (e.g., already removed), but you can tighten if desired.
std::remove(SlotPath(slot).c_str());
}
}
static void EnsureDirectory_(const std::string &dir)
{
if (dir.empty()) { return; }
#if __cplusplus >= 201703L
namespace fs = std::filesystem;
std::error_code ec;
if (!fs::exists(dir, ec))
{
fs::create_directories(dir, ec);
}
MFEM_VERIFY(!ec, "EnsureDirectory: failed to create directory: " + dir);
#else
#if defined(_WIN32)
const int rc = _mkdir(dir.c_str());
if (rc != 0)
{
// If directory already exists, _mkdir fails. We accept that.
// There's no reliable portable "exists" check in pre-C++17 without more code.
}
#else
const int rc = mkdir(dir.c_str(), 0755);
if (rc != 0 && errno != EEXIST)
{
MFEM_ABORT("EnsureDirectory: failed to create directory: " + dir);
}
#endif
#endif
}
};
} // namespace mfem
#endif // MFEM_REVOLVE_SEPARATE_FILE_STORAGE_HPP
@@ -0,0 +1,273 @@
#ifndef MFEM_SEGMENT_CHECKPOINT_STORAGE_HPP
#define MFEM_SEGMENT_CHECKPOINT_STORAGE_HPP
#include "mfem.hpp"
#include "file_checkpoint_storage.hpp" // reuses DefaultCheckpointBinaryIO
#include <cstdint>
#include <string>
#include <vector>
#include <fstream>
#include <sstream>
#include <iomanip>
#include <set>
#if __has_include(<filesystem>)
#include <filesystem>
namespace mfem_fs = std::filesystem;
#define MFEM_HAVE_FILESYSTEM 1
#else
#define MFEM_HAVE_FILESYSTEM 0
#endif
namespace mfem
{
/**
* @brief Segment/range file storage: one file per handle-id range + in-file offsets.
*
* Segment file = <dir>/<prefix><segment_id><ext>
* Where segment_id = handle / records_per_file.
*
* Pros:
* - Much fewer files than file-per-snapshot
* - No per-snapshot file create/delete
*
* Cons:
* - Append-only: Erase() does not reclaim file space (records remain)
* - Index is in-memory (not restartable across processes unless you persist it)
*/
template <typename Snapshot,
typename SnapshotIO = DefaultCheckpointBinaryIO<Snapshot> >
class SegmentedFileCheckpointStorage
{
public:
using Handle = std::int64_t;
SegmentedFileCheckpointStorage(const std::string &directory,
std::int64_t records_per_file = 4096,
const std::string &prefix = "seg_",
const std::string &extension = ".bin",
bool create_dir = true,
bool keep_files = false)
: dir_(directory),
prefix_(prefix),
ext_(extension),
keep_files_(keep_files),
records_per_file_(records_per_file)
{
MFEM_VERIFY(records_per_file_ > 0, "SegmentedFileCheckpointStorage: records_per_file must be > 0.");
MFEM_VERIFY(!dir_.empty(), "SegmentedFileCheckpointStorage: empty directory.");
if (create_dir)
{
#if MFEM_HAVE_FILESYSTEM
std::error_code ec;
mfem_fs::create_directories(mfem_fs::path(dir_), ec);
MFEM_VERIFY(!ec, "SegmentedFileCheckpointStorage: failed to create directory.");
#else
MFEM_ABORT("SegmentedFileCheckpointStorage: create_dir=true requires <filesystem> support.");
#endif
}
}
Handle InvalidHandle() const { return (Handle)-1; }
bool IsValid(const Handle &h) const { return h >= 0; }
Handle Store(Snapshot &&snap)
{
const Handle id = AllocateId_();
EnsureMetaSize_(id);
const std::int64_t seg = SegmentId_(id);
const std::string path = SegmentPath_(seg);
const std::uint64_t offset = AppendRecord_(path, snap);
const std::uint64_t bytes = last_payload_bytes_;
meta_[id].valid = true;
meta_[id].seg_id = seg;
meta_[id].offset = offset;
meta_[id].bytes = bytes;
touched_segments_.insert(seg);
return id;
}
template <typename Func>
void Read(const Handle &h, Func &&f) const
{
MFEM_VERIFY(IsValid(h), "SegmentedFileCheckpointStorage: Read invalid handle.");
MFEM_VERIFY(h < (Handle)meta_.size(), "SegmentedFileCheckpointStorage: Read out-of-range handle.");
MFEM_VERIFY(meta_[h].valid, "SegmentedFileCheckpointStorage: Read on erased handle.");
const std::string path = SegmentPath_(meta_[h].seg_id);
std::ifstream is(path, std::ios::binary);
MFEM_VERIFY(is.is_open(), "SegmentedFileCheckpointStorage: open-for-read failed.");
// Seek to record start
is.seekg((std::streamoff)meta_[h].offset, std::ios::beg);
MFEM_VERIFY(is.good(), "SegmentedFileCheckpointStorage: seekg failed.");
std::uint64_t payload_bytes = 0;
is.read(reinterpret_cast<char*>(&payload_bytes), sizeof(payload_bytes));
MFEM_VERIFY(is.good(), "SegmentedFileCheckpointStorage: read header failed.");
MFEM_VERIFY(payload_bytes == meta_[h].bytes, "SegmentedFileCheckpointStorage: payload size mismatch.");
Snapshot snap = SnapshotIO::Read(is);
MFEM_VERIFY(is.good(), "SegmentedFileCheckpointStorage: payload read failed.");
f(snap);
}
void Erase(Handle &h)
{
if (!IsValid(h)) { h = InvalidHandle(); return; }
MFEM_VERIFY(h < (Handle)meta_.size(), "SegmentedFileCheckpointStorage: Erase out-of-range handle.");
meta_[h].valid = false;
if (!keep_files_)
{
free_.push_back(h); // allow handle reuse
}
h = InvalidHandle();
}
/**
* @brief Optional cleanup helper (NOT part of the required storage interface).
*
* Deletes all segment files touched by this storage object.
*/
void PurgeAllFiles()
{
#if MFEM_HAVE_FILESYSTEM
for (auto seg : touched_segments_)
{
const std::string path = SegmentPath_(seg);
std::error_code ec;
mfem_fs::remove(mfem_fs::path(path), ec);
// Best-effort cleanup:
MFEM_VERIFY(!ec, "SegmentedFileCheckpointStorage: failed to remove segment file.");
}
touched_segments_.clear();
#else
MFEM_ABORT("SegmentedFileCheckpointStorage: PurgeAllFiles requires <filesystem> support.");
#endif
}
private:
struct Meta
{
bool valid = false;
std::int64_t seg_id = 0;
std::uint64_t offset = 0;
std::uint64_t bytes = 0;
};
std::string dir_, prefix_, ext_;
bool keep_files_ = false;
std::int64_t records_per_file_ = 4096;
mutable Handle next_id_ = 0;
mutable std::vector<Handle> free_;
std::vector<Meta> meta_;
mutable std::uint64_t last_payload_bytes_ = 0;
std::set<std::int64_t> touched_segments_;
void EnsureMetaSize_(Handle id)
{
if ((std::size_t)id >= meta_.size())
{
meta_.resize((std::size_t)id + 1);
}
}
Handle AllocateId_()
{
if (!keep_files_ && !free_.empty())
{
const Handle id = free_.back();
free_.pop_back();
return id;
}
return next_id_++;
}
std::int64_t SegmentId_(Handle id) const
{
return (std::int64_t)(id / records_per_file_);
}
std::string SegmentPath_(std::int64_t seg) const
{
std::ostringstream oss;
oss << prefix_ << std::setw(8) << std::setfill('0') << seg << ext_;
#if MFEM_HAVE_FILESYSTEM
mfem_fs::path p = mfem_fs::path(dir_) / mfem_fs::path(oss.str());
return p.string();
#else
return dir_ + "/" + oss.str();
#endif
}
static void EnsureFileExists_(const std::string &path)
{
// Try open for read/write; if missing, create.
std::fstream fs(path, std::ios::binary | std::ios::in | std::ios::out);
if (!fs.is_open())
{
std::ofstream create(path, std::ios::binary | std::ios::out);
MFEM_VERIFY(create.is_open(), "SegmentedFileCheckpointStorage: file create failed.");
}
}
// Append framed record to file. Returns offset of record start. Updates last_payload_bytes_.
std::uint64_t AppendRecord_(const std::string &path, const Snapshot &snap)
{
EnsureFileExists_(path);
std::fstream fs(path, std::ios::binary | std::ios::in | std::ios::out);
MFEM_VERIFY(fs.is_open(), "SegmentedFileCheckpointStorage: open-for-append failed.");
fs.seekp(0, std::ios::end);
MFEM_VERIFY(fs.good(), "SegmentedFileCheckpointStorage: seekp(end) failed.");
const std::streamoff begin = (std::streamoff)fs.tellp();
MFEM_VERIFY(begin >= 0, "SegmentedFileCheckpointStorage: tellp failed.");
// placeholder for payload bytes
std::uint64_t payload_bytes = 0;
fs.write(reinterpret_cast<const char*>(&payload_bytes), sizeof(payload_bytes));
MFEM_VERIFY(fs.good(), "SegmentedFileCheckpointStorage: write header failed.");
// payload
SnapshotIO::Write(fs, snap);
MFEM_VERIFY(fs.good(), "SegmentedFileCheckpointStorage: write payload failed.");
const std::streamoff end = (std::streamoff)fs.tellp();
MFEM_VERIFY(end >= begin, "SegmentedFileCheckpointStorage: tellp end failed.");
payload_bytes = (std::uint64_t)(end - begin - (std::streamoff)sizeof(std::uint64_t));
last_payload_bytes_ = payload_bytes;
// backpatch payload size
fs.seekp(begin, std::ios::beg);
MFEM_VERIFY(fs.good(), "SegmentedFileCheckpointStorage: seekp(begin) failed.");
fs.write(reinterpret_cast<const char*>(&payload_bytes), sizeof(payload_bytes));
MFEM_VERIFY(fs.good(), "SegmentedFileCheckpointStorage: backpatch failed.");
return (std::uint64_t)begin;
}
};
} // namespace mfem
#endif // MFEM_SEGMENT_CHECKPOINT_STORAGE_HPP
+701
View File
@@ -0,0 +1,701 @@
#include "mfem.hpp"
#include "dynamic_checkpointing.hpp"
#include "fixed_slot_checkpoint_storage.hpp"
#include "file_checkpoint_storage.hpp"
#include <cmath>
#include <iomanip>
#include <string>
using namespace mfem;
/************************************************************************
* The forward state consists of vector and several additional values.
* The goal is to demontrate how different storages can be used together
* with DynamicCheckpointing in order to avoid unnecessary memory
* allocations, data copies, and deallocations.
* *********************************************************************/
struct State
{
mfem::real_t time = 0.0;
mfem::real_t obj = 0.0;
mfem::Vector v;
};
// ---------------------------
// Snapshot type used by storage:
// a lightweight view (non-owning).
// ---------------------------
struct StateSnapshotView
{
mfem::real_t time = 0.0;
mfem::real_t obj = 0.0;
// Points to n*sizeof(real_t) bytes.
// - during Store(): points to State::v host data
// - during Read(): points into storage's internal scratch buffer
const unsigned char *v_bytes = nullptr;
};
// ---------------------------
// Packer for fixed-slot storage
// Layout in slot bytes:
// [ time | obj | v[0..n-1] ]
// ---------------------------
class StateSnapshotViewPacker
{
public:
explicit StateSnapshotViewPacker(int n) : n_(n)
{
MFEM_VERIFY(n_ > 0, "StateSnapshotViewPacker: n must be > 0.");
}
std::size_t SlotBytes() const
{
return (std::size_t)(2 + n_) * sizeof(mfem::real_t);
}
void Pack(const StateSnapshotView &s, void *dst) const
{
MFEM_VERIFY(dst != nullptr, "Pack: dst is null.");
MFEM_VERIFY(s.v_bytes != nullptr, "Pack: v_bytes is null.");
unsigned char *b = static_cast<unsigned char*>(dst);
std::memcpy(b + 0*sizeof(mfem::real_t), &s.time, sizeof(mfem::real_t));
std::memcpy(b + 1*sizeof(mfem::real_t), &s.obj, sizeof(mfem::real_t));
std::memcpy(b + 2*sizeof(mfem::real_t),
s.v_bytes,
(std::size_t)n_ * sizeof(mfem::real_t));
}
// IMPORTANT:
// Unpack produces a view into the provided src buffer.
// The resulting pointer is only valid as long as src remains unchanged.
void Unpack(const void *src, StateSnapshotView &out) const
{
MFEM_VERIFY(src != nullptr, "Unpack: src is null.");
const unsigned char *b = static_cast<const unsigned char*>(src);
std::memcpy(&out.time, b + 0*sizeof(mfem::real_t), sizeof(mfem::real_t));
std::memcpy(&out.obj, b + 1*sizeof(mfem::real_t), sizeof(mfem::real_t));
out.v_bytes = b + 2*sizeof(mfem::real_t);
}
int N() const { return n_; }
private:
int n_ = 0;
};
/*****************************************************************************
* Provides a recipe for runing dynamic checkpointing with memory storage
* alocated as one single big block.
* **************************************************************************/
void run_fixed_slot_memory_storage( mfem::real_t Tfinal /*Final time*/,
mfem::real_t dtime /*time step*/,
int s /* checkpoint budget (real checkpoints) */,
int n /*State vector size*/,
bool print=true)
{
StateSnapshotViewPacker packer(n);
using Storage = mfem::FixedSlotMemoryCheckpointStorage<StateSnapshotView,
StateSnapshotViewPacker>;
Storage storage(s, packer);
// Snapshot type is StateSnapshotView
using Checkpointing = mfem::DynamicCheckpointing<StateSnapshotView, Storage>;
Checkpointing ckpt(s, storage);
// Returns view of the State and avoids data transfer
auto make_snapshot = [&](const State &u) -> StateSnapshotView
{
MFEM_VERIFY(u.v.Size() == n, "make_snapshot: State.v size changed!");
const mfem::real_t *vh = u.v.HostRead(); // host pointer
StateSnapshotView snap;
snap.time = u.time;
snap.obj = u.obj;
snap.v_bytes = reinterpret_cast<const unsigned char*>(vh);
return snap;
};
//Transfers data from the snaphot view to the State u_out.
auto restore_snapshot = [&](const StateSnapshotView &snap, State &u_out)
{
u_out.time = snap.time;
u_out.obj = snap.obj;
if (u_out.v.Size() != n) { u_out.v.SetSize(n); }
mfem::real_t *vh = u_out.v.HostWrite();
std::memcpy(vh,
snap.v_bytes,
(std::size_t)n * sizeof(mfem::real_t));
};
using Step = Checkpointing::Step;
State u;
u.v.SetSize(n); u.v=0.0;
auto primal_step = [&](State &u_inout, Step i)
{
// advance u_inout.v, u_inout.time, u_inout.obj
if(print)
std::cout<<"Forward step: "<<i<<" time="<<u_inout.time<<" obj="<<u_inout.obj<<std::endl;
u_inout.time=i*0.1;
u_inout.obj=i*0.2;
};
struct AdjointState { /* ... */ };
AdjointState q;
auto adjoint_step = [&](AdjointState &q_ip1_inout, const State &u_i, Step i)
{
// use u_i.time/u_i.obj/u_i.v to update adjoint
if(print)
std::cout<<"Adjoint step: "<<i<<" time="<<u_i.time<<" obj="<<u_i.obj<<std::endl;
};
// Forward sweep (unknown number of steps)
mfem::real_t t=0.0;
Step i=0;
while(t<Tfinal)
{
// Store snapshot of u_i when manager decides; then primal_step u_i -> u_{i+1}
ckpt.ForwardStep(i, u, primal_step, make_snapshot);
t+=dtime;
++i;
}
mfem::out << "Total number of steps m="<<i-1<<"\n";
mfem::out << "\nBackward sweep.\n";
const Step m=i;
// Backward sweep
State u_work;
u_work.v.SetSize(n); // allocate once
for (Step i = m - 1; i >= 0; --i)
{
ckpt.BackwardStep(i, q, u_work,
primal_step, adjoint_step,
make_snapshot, restore_snapshot);
if (i == 0) { break; }
}
}
/****************************************************************************
* Provides a recipe for runing dynamic checkpointing with file storage
* alocated as one single big file.
****************************************************************************/
/****************************************************************************
* The dynamic algorithm conceptually maintains s + 1 checkpoints, but the
* last one is a placeholder checkpoint that stores no solution and takes
* little memory; only s are real stored solutions.
* So FixedSlotFileCheckpointStorage needs only s slots.
* File behavior for FixedSlotFileCheckpointStorage:
- Writes a small header + s * slot_bytes payload region,
- Erase(handle) just returns the slot to the free list
(file contents remain, will be overwritten later),
no filesystem metadata churn.
* Device memory: Runs on GPU -packing to a host file requires a host pointer.
Using:
u.v.HostRead() for packing,
u_out.v.HostWrite() for restore,
keeps memmory access correct.
* Snapshot pointer lifetime: In the read path, StateSnapshotViewPacker::Unpack
sets snap.v_bytes to point into the storages internal read buffer. That
pointer is only valid until the storage performs another Read/Store
(and in general, one should treat it as valid only during restore_snapshot).
*****************************************************************************/
void run_fixed_slot_file_storage(mfem::real_t Tfinal /*Final time*/,
mfem::real_t dtime /*time step*/,
int s /* checkpoint budget (real checkpoints) */,
int n /*State vector size*/,
const std::string &file_path,
bool print=true)
{
// s = number of REAL checkpoints (placeholder doesn't store a snapshot)
MFEM_VERIFY(s > 0, "Need s > 0.");
MFEM_VERIFY(n > 0, "Need n > 0.");
StateSnapshotViewPacker packer(n);
// Fixed-slot FILE storage:
// - single file (file_path)
// - max_slots = s (real checkpoints)
// - slot_bytes = packer.SlotBytes()
// - truncate=true to start fresh
// - flush_on_store=false for performance (set true if you want extra safety)
using Storage = mfem::FixedSlotFileCheckpointStorage<StateSnapshotView, StateSnapshotViewPacker>;
Storage storage(file_path,
/*max_slots=*/s,
/*packer=*/packer,
/*truncate=*/true,
/*flush_on_store=*/false);
// Dynamic checkpointing manager using this storage.
using Checkpointing = mfem::DynamicCheckpointing<StateSnapshotView, Storage>;
Checkpointing ckpt(s, storage);
// ---------------------------
// Callbacks
// ---------------------------
// Convert State -> SnapshotView (no allocation).
auto make_snapshot = [&](const State &u) -> StateSnapshotView
{
MFEM_VERIFY(u.v.Size() == n, "make_snapshot: State.v size changed!");
// Ensure host pointer valid even when MFEM is using device memory.
const mfem::real_t *vh = u.v.HostRead();
StateSnapshotView snap;
snap.time = u.time;
snap.obj = u.obj;
snap.v_bytes = reinterpret_cast<const unsigned char*>(vh);
return snap;
};
// Convert SnapshotView -> State (must COPY out of v_bytes).
auto restore_snapshot = [&](const StateSnapshotView &snap, State &u_out)
{
u_out.time = snap.time;
u_out.obj = snap.obj;
if (u_out.v.Size() != n) { u_out.v.SetSize(n); }
mfem::real_t *vh = u_out.v.HostWrite();
std::memcpy(vh,
snap.v_bytes,
(std::size_t)n * sizeof(mfem::real_t));
};
// Example primal/adjoint step signatures (replace with your own):
using Step = Checkpointing::Step;
auto primal_step = [&](State &u_inout, Step i)
{
if(print)
std::cout<<"Forward step: "<<i<<" time="<<u_inout.time<<" obj="<<u_inout.obj<<std::endl;
u_inout.time=i*0.1;
u_inout.obj=i*0.2;
// update u_inout.time, u_inout.obj, u_inout.v
};
struct AdjointState
{
// your adjoint variables, e.g., mfem::Vector lambda;
};
auto adjoint_step = [&](AdjointState &q_ip1_inout, const State &u_i, Step i)
{
// update adjoint using u_i
if(print)
std::cout<<"Adjoint step: "<<i<<" time="<<u_i.time<<" obj="<<u_i.obj<<std::endl;
};
// ---------------------------
// Forward/backward skeleton
// ---------------------------
State u;
u.v.SetSize(n);
u.time = 0.0;
u.obj = 0.0;
u.v = 0.0;
// Forward sweep (unknown number of steps)
mfem::real_t t=0.0;
Step i=0;
while(t<Tfinal)
{
ckpt.ForwardStep(i, u, primal_step, make_snapshot);
t+=dtime;
++i;
}
mfem::out << "Total number of steps m="<<i-1<<"\n";
mfem::out << "\nBackward sweep.\n";
// Backward sweep
AdjointState q;
const Step m=i;
State u_work;
u_work.v.SetSize(n);
for (Step i = m - 1; i >= 0; --i)
{
ckpt.BackwardStep(i, q, u_work,
primal_step, adjoint_step,
make_snapshot, restore_snapshot);
if (i == 0) { break; }
}
}
/******************************************************************************
The following classes implement the IO and the data necessery for dynamic
checkpointing with mfem::FileCheckpointStorage (one file per stored snapshot).
With FileCheckpointStorage, the Snapshot must own the data which will be
written, i.e. it cannot be a view containing pointers into State::v
as those pointers would be meaningless when read back.
*******************************************************************************/
// The StateSnapshot is just a copy of the State
using StateSnapshot=State;
struct StateSnapshotBinaryIO
{
static void Write(std::ostream &os, const StateSnapshot &s)
{
os.write(reinterpret_cast<const char*>(&s.time), sizeof(mfem::real_t));
os.write(reinterpret_cast<const char*>(&s.obj), sizeof(mfem::real_t));
MFEM_VERIFY(os.good(), "StateSnapshotBinaryIO: write time/obj failed.");
const std::int64_t n = (std::int64_t)s.v.Size();
os.write(reinterpret_cast<const char*>(&n), sizeof(n));
MFEM_VERIFY(os.good(), "StateSnapshotBinaryIO: write vector size failed.");
if (n > 0)
{
const mfem::real_t *vh = s.v.HostRead(); // ensure host pointer
os.write(reinterpret_cast<const char*>(vh),
(std::streamsize)(n * (std::int64_t)sizeof(mfem::real_t)));
MFEM_VERIFY(os.good(), "StateSnapshotBinaryIO: write vector payload failed.");
}
}
static StateSnapshot Read(std::istream &is)
{
StateSnapshot s;
is.read(reinterpret_cast<char*>(&s.time), sizeof(mfem::real_t));
is.read(reinterpret_cast<char*>(&s.obj), sizeof(mfem::real_t));
MFEM_VERIFY(is.good(), "StateSnapshotBinaryIO: read time/obj failed.");
std::int64_t n = 0;
is.read(reinterpret_cast<char*>(&n), sizeof(n));
MFEM_VERIFY(is.good(), "StateSnapshotBinaryIO: read vector size failed.");
MFEM_VERIFY(n >= 0, "StateSnapshotBinaryIO: invalid negative vector size.");
s.v.SetSize((int)n);
if (n > 0)
{
mfem::real_t *vh = s.v.HostWrite(); // ensure host pointer
is.read(reinterpret_cast<char*>(vh),
(std::streamsize)(n * (std::int64_t)sizeof(mfem::real_t)));
MFEM_VERIFY(is.good(), "StateSnapshotBinaryIO: read vector payload failed.");
}
return s;
}
};
/***************************************************************************
Important differences vs FixedSlotFileCheckpointStorage
Metadata overhead
FileCheckpointStorage typically:
*creates a new file for each Store() (plus a temp file rename),
*deletes a file on each Erase() (unless keep_files=true).
Dynamic checkpointing can perform many store/erase operations during
forward and during recomputation in reverse, so this can hammer filesystem
metadata on parallel filesystems.
If thats a concern, prefer:
*fixed-slot single-file storage, or
*segmented range file storage (few files + offsets).
Correctness note about pointers:
This is why we used an owning snapshot (StateSnapshot): storing a view
with v_bytes into file storage would store pointer values, which are
meaningless when read back.
*****************************************************************************/
void run_file_storage(mfem::real_t Tfinal /*Final time*/,
mfem::real_t dtime /*time step*/,
int s /* checkpoint budget (real checkpoints) */,
int n /*State vector size*/,
const std::string &directory,
bool print=true)
{
MFEM_VERIFY(s > 0, "Need s > 0.");
MFEM_VERIFY(n > 0, "Need n > 0.");
// One file per snapshot (create_dir=true). keep_files=false means Erase() removes files.
using Storage = mfem::FileCheckpointStorage<StateSnapshot, StateSnapshotBinaryIO>;
Storage storage(directory, "ckpt_", ".bin", /*create_dir=*/true, /*keep_files=*/false);
// Snapshot type = StateSnapshot (owning)
using Checkpointing = mfem::DynamicCheckpointing<StateSnapshot, Storage>;
Checkpointing ckpt(s, storage);
using Step = Checkpointing::Step;
// ---- Callbacks ----
// Make an owning snapshot from the current State (deep copy of vector payload)
auto make_snapshot = [&](const State &u) -> StateSnapshot
{
MFEM_VERIFY(u.v.Size() == n, "make_snapshot: State.v size changed!");
StateSnapshot snap;
snap.time = u.time;
snap.obj = u.obj;
snap.v.SetSize(n);
const mfem::real_t *src = u.v.HostRead();
mfem::real_t *dst = snap.v.HostWrite();
std::memcpy(dst, src, (std::size_t)n * sizeof(mfem::real_t));
return snap; // move into storage
};
// Restore from snapshot into an actual State (copy payload)
auto restore_snapshot = [&](const StateSnapshot &snap, State &u_out)
{
u_out.time = snap.time;
u_out.obj = snap.obj;
MFEM_VERIFY(snap.v.Size() == n, "restore_snapshot: snapshot vector size mismatch.");
if (u_out.v.Size() != n) { u_out.v.SetSize(n); }
const mfem::real_t *src = snap.v.HostRead();
mfem::real_t *dst = u_out.v.HostWrite();
std::memcpy(dst, src, (std::size_t)n * sizeof(mfem::real_t));
};
// Example primal step signature (replace with your integrator)
auto primal_step = [&](State &u_inout, Step i)
{
// advance u_inout.time, u_inout.obj, u_inout.v
if(print)
std::cout<<"Forward step: "<<i<<" time="<<u_inout.time<<" obj="<<u_inout.obj<<std::endl;
u_inout.time=i*0.1;
u_inout.obj=i*0.2;
};
// Example adjoint state and adjoint step signature (replace with yours)
struct AdjointState
{
// e.g., mfem::Vector lambda;
};
auto adjoint_step = [&](AdjointState &q_ip1_inout, const State &u_i, Step i)
{
// update adjoint using u_i
if(print)
std::cout<<"Adjoint step: "<<i<<" time="<<u_i.time<<" obj="<<u_i.obj<<std::endl;
};
// ---- Forward/backward skeleton ----
State u;
u.v.SetSize(n);
u.time = 0.0;
u.obj = 0.0;
u.v = 0.0;
// Forward sweep (unknown number of steps)
mfem::real_t t=0.0;
Step i=0;
while(t<Tfinal)
{
ckpt.ForwardStep(i, u, primal_step, make_snapshot);
t+=dtime;
++i;
}
mfem::out << "Total number of steps m="<<i-1<<"\n";
mfem::out << "\nBackward sweep.\n";
// backward
const Step m=i;
AdjointState q;
State u_work;
u_work.v.SetSize(n);
for (Step i = m - 1; i >= 0; --i)
{
ckpt.BackwardStep(i, q, u_work,
primal_step, adjoint_step,
make_snapshot, restore_snapshot);
if (i == 0) { break; }
}
}
/****************************************************************************
* Using standard in memory storage. Works with variable snaphot sizes.
* The size of the vector (n) is provided in order to compare the run to
* the other implementations.
***************************************************************************/
void run_in_memory_storage(mfem::real_t Tfinal /*Final time*/,
mfem::real_t dtime /*time step*/,
int s /* checkpoint budget (real checkpoints) */,
int n /*State vector size*/,
bool print=true)
{
InMemoryCheckpointStorage<StateSnapshot> storage;
using Checkpointing = mfem::DynamicCheckpointing<StateSnapshot,
InMemoryCheckpointStorage<StateSnapshot>>;
Checkpointing ckpt(s, storage);
auto make_snapshot = [&](const State &u) -> StateSnapshot
{
StateSnapshot snap;
snap.time = u.time;
snap.obj = u.obj;
snap.v = u.v; // deep copy (MFEM handles memory)
return snap; // moved into storage by Store()
};
auto restore_snapshot = [&](const StateSnapshot &snap, State &u_out)
{
u_out.time = snap.time;
u_out.obj = snap.obj;
u_out.v = snap.v; // deep copy back into working state
};
using Step = Checkpointing::Step;
// Primal step: u_i -> u_{i+1}
auto primal_step = [&](State &u, Step i)
{
// ... update u.time, u.obj, u.v ...
if(print)
std::cout<<"Forward step: "<<i<<" time="<<u.time<<" obj="<<u.obj<<std::endl;
u.time=i*0.1;
u.obj=i*0.2;
};
// Adjoint step: q_{i+1} -> q_i using u_i
struct AdjointState
{
// e.g. mfem::Vector lambda;
};
auto adjoint_step = [&](AdjointState &q, const State &u_i, Step i)
{
// ... update q using u_i ...
if(print)
std::cout<<"Adjoint step: "<<i<<" time="<<u_i.time<<" obj="<<u_i.obj<<std::endl;
};
State u;
u.v.SetSize(n);
u.v=0.0;
u.time = 0.0;
u.obj = 0.0;
// Forward sweep (unknown number of steps)
mfem::real_t t=0.0;
Step i=0;
while(t<Tfinal)
{
ckpt.ForwardStep(i, u, primal_step, make_snapshot);
t+=dtime;
++i;
}
mfem::out << "Total number of steps m="<<i-1<<"\n";
mfem::out << "\nBackward sweep.\n";
// Backward
const Step m=i;
AdjointState q;
State u_work; // used to restore/recompute primal state at step i
u_work.v.SetSize(u.v.Size());
for (Step i = m - 1; i >= 0; --i)
{
ckpt.BackwardStep(i, q, u_work,
primal_step, adjoint_step,
make_snapshot, restore_snapshot);
if (i == 0) { break; }
}
}
int main(int argc, char *argv[])
{
int s = 5/* checkpoint budget (real checkpoints) */;
int n = 30 /* fixed State.v size */;
std::string filepath="dynamic_ckpts.bin"; /*file name for fixed slot single file*/
std::string directory="dyn_ckpts";
mfem::real_t Tfinal=1;
mfem::real_t dt=0.1;
int print=true;
OptionsParser args(argc, argv);
args.AddOption(&s, "-s",
"--checkpoints", "Checkpoint budget s (real checkpoints).");
args.AddOption(&n, "-n",
"--size", "Size of the state vector.");
args.AddOption(&Tfinal, "-T",
"--tfinal", "Terminate when accumulated time reaches Tfinal.");
args.AddOption(&dt, "-dt", "--dt", "Time step.");
args.AddOption(&filepath, "-file", "--file",
"File name for storing the checkpoints.");
args.AddOption(&directory, "-dir", "--dir",
"Directory for storing the checkpoints.");
args.AddOption(&print, "-print", "--print",
"Turn on/off printing." );
args.Parse();
if (!args.Good())
{
args.PrintUsage(mfem::out);
return 1;
}
args.PrintOptions(mfem::out);
double run_time;
mfem::out<<"\nFixedSlotMemoryStorage Start\n";
mfem::tic();
run_fixed_slot_memory_storage(Tfinal, dt, s, n, print);
run_time=mfem::toc();
mfem::out<<"FixedSlotMemoryStorage time= "<<run_time<<std::endl;
mfem::out<<"\nFixedSlotFileStorage Start\n";
mfem::tic();
run_fixed_slot_file_storage(Tfinal,dt, s,n, filepath,print);
run_time=mfem::toc();
mfem::out<<"FixedSlotFileStorage time= "<<run_time<<std::endl;
mfem::out<<"\nFileStorage Start\n";
mfem::tic();
run_file_storage(Tfinal,dt, s,n, directory, print);
run_time=mfem::toc();
mfem::out<<"FileStorage time= "<<run_time<<std::endl;
mfem::out<<"\nInMemoryStorage\n";
mfem::tic();
run_in_memory_storage(Tfinal,dt, s,n, print);
run_time=mfem::toc();
mfem::out<<"InMemoryStorage time= "<<run_time<<std::endl;
mfem::out << "\nDone.\n";
return 0;
}

Some files were not shown because too many files have changed in this diff Show More