Merge branch 'master' into checkpointing

Add option to skip branch-history check

.github/workflows/builds-and-tests.yml

@@ -142,6 +142,10 @@ jobs:
 
     continue-on-error: ${{ matrix.enzyme && true || false }}
 
+    # Enable ccache for all jobs except Windows (would need sccache).
+    env:
+      USE_CCACHE: ${{ matrix.os != 'windows-latest' }}
+
     steps:
       # Fix 'No space left on device' errors for Ubuntu builds.
       - name: Run Actions Cleaner
@@ -290,6 +294,52 @@ jobs:
           echo "OMPI_CC=$LLVM_PREFIX/bin/clang" >> $GITHUB_ENV
           echo "OMPI_CXX=$LLVM_PREFIX/bin/clang++" >> $GITHUB_ENV
 
+      # Restore the compiler cache (ccache). The key embeds the run id, so new
+      # runs save a fresh snapshot; the restore-keys prefix warm-starts from the
+      # most recent prior run (incl. the base branch for PRs).
+      - name: cache ccache
+        if: ${{ env.USE_CCACHE == 'true' }}
+        uses: actions/cache@v5
+        with:
+          path: .ccache
+          key: ccache-${{ matrix.os }}-${{ matrix.build-system }}-${{ matrix.target }}-${{ matrix.mpi }}-${{ matrix.hypre-target }}-${{ matrix.precision }}${{ matrix.enzyme && '-enzyme' || '' }}-${{ github.run_id }}
+          restore-keys: |
+            ccache-${{ matrix.os }}-${{ matrix.build-system }}-${{ matrix.target }}-${{ matrix.mpi }}-${{ matrix.hypre-target }}-${{ matrix.precision }}${{ matrix.enzyme && '-enzyme' || '' }}-
+
+      # Configure ccache and select how it is injected into the MFEM build:
+      # - make:  set CXX="ccache g++"; for MPI, OMPI_CXX="ccache g++" so mpicxx
+      #          runs ccache around g++ (not ccache around the mpicxx wrapper).
+      # - cmake: set CMAKE_<LANG>_COMPILER_LAUNCHER=ccache.
+      # - enzyme: wrap the brew clang++ via OMPI_CXX.
+      # The chosen options are passed through build-mfem's 'config-options'
+      # input (see the build step below).
+      - name: configure ccache
+        if: ${{ env.USE_CCACHE == 'true' }}
+        run: |
+          command -v ccache >/dev/null 2>&1 || {
+            if [[ "${{ runner.os }}" == "Linux" ]]; then
+              sudo apt-get update && sudo apt-get install -y ccache
+            else
+              brew install ccache
+            fi
+          }
+          echo "CCACHE_DIR=${{ github.workspace }}/.ccache" >> $GITHUB_ENV
+          echo "CCACHE_MAXSIZE=1G" >> $GITHUB_ENV
+          echo "CCACHE_COMPILERCHECK=content" >> $GITHUB_ENV
+          # Ignore header timestamps (restamped by each checkout) so direct mode hits.
+          echo "CCACHE_SLOPPINESS=include_file_mtime,include_file_ctime,time_macros" >> $GITHUB_ENV
+          # Hash absolute paths relative to the workspace.
+          echo "CCACHE_BASEDIR=${{ github.workspace }}" >> $GITHUB_ENV
+          if [[ "${{ matrix.enzyme }}" == "true" ]]; then
+            echo "OMPI_CXX=ccache $LLVM_PREFIX/bin/clang++" >> $GITHUB_ENV
+          elif [[ "${{ matrix.build-system }}" == "cmake" ]]; then
+            echo 'CCACHE_CONFIG_OPTS=-DCMAKE_CXX_COMPILER_LAUNCHER=ccache -DCMAKE_C_COMPILER_LAUNCHER=ccache' >> $GITHUB_ENV
+          else
+            echo "OMPI_CXX=ccache g++" >> $GITHUB_ENV
+            echo 'CCACHE_CONFIG_OPTS=CXX="ccache g++" MPICXX="mpicxx"' >> $GITHUB_ENV
+          fi
+        shell: bash
+
       # MFEM build and test
       - name: build
         uses: mfem/github-actions/build-mfem@v2.7
@@ -305,9 +355,14 @@ jobs:
           metis-dir: ${{ env.METIS_TOP_DIR }}
           mfem-dir: ${{ env.MFEM_TOP_DIR }}
           precision: ${{ matrix.precision }}
-          config-options: ${{ matrix.config-opts }}
+          config-options: ${{ matrix.config-opts }} ${{ env.CCACHE_CONFIG_OPTS }}
           library-only: ${{ matrix.target == 'dbg' && matrix.os != 'ubuntu-latest' }}
 
+      - name: ccache stats
+        if: ${{ env.USE_CCACHE == 'true' }}
+        run: ccache -s
+        shell: bash
+
       # Run checks (and only checks) on debug targets
       - name: checks
         if: matrix.build-system == 'make' && matrix.target == 'dbg'

.github/workflows/cleanup-pr-caches.yml

@@ -0,0 +1,42 @@
+# 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.
+---
+# A closed PR's caches can never be restored again, so delete them to free
+# space against the 10 GB per-repo cache limit.
+name: Cleanup PR caches
+
+on:
+  pull_request:
+    types: [closed]
+
+permissions:
+  actions: write
+
+jobs:
+  cleanup:
+    runs-on: ubuntu-latest
+    steps:
+      - name: Delete caches for the closed PR
+        env:
+          GH_TOKEN: ${{ secrets.GITHUB_TOKEN }}
+          GH_REPO: ${{ github.repository }}
+          PR_REF: refs/pull/${{ github.event.pull_request.number }}/merge
+        run: |
+          echo "Deleting caches for $PR_REF"
+          while :; do
+            ids=$(gh cache list --ref "$PR_REF" --limit 100 --json id --jq '.[].id')
+            [ -n "$ids" ] || break
+            echo "$ids" | while read -r id; do
+              [ -n "$id" ] || continue
+              echo "Deleting cache $id"
+              gh cache delete "$id" || echo "  (already gone)"
+            done
+          done

.github/workflows/repo-check.yml

@@ -13,6 +13,7 @@ name: "Checks"
 
 permissions:
   actions: write
+  pull-requests: read
 
 on:
   push:
@@ -29,6 +30,11 @@ concurrency:
 # by checking if the workflow trigger is 'push' ("github.event_name == 'push'")
 # and if we are in a fork ("github.event.pull_request.head.repo.full_name !=
 # github.repository").
+#
+# The logic for the branch-history check is slightly different, since that check
+# also inspects the PR's labels to allow for overriding failures. In this case,
+# we run on all 'pull_request' triggers, but only run for 'push' triggers that
+# do not correspond to any open PRs.
 
 jobs:
   file-headers-check:
@@ -128,10 +134,7 @@ jobs:
 
   branch-history:
     if: |
-      github.ref != 'refs/heads/next' &&
-      github.ref != 'refs/heads/master' &&
-      (github.event_name == 'push' ||
-       github.event.pull_request.head.repo.full_name != github.repository)
+      github.ref != 'refs/heads/next' && github.ref != 'refs/heads/master'
     runs-on: ubuntu-latest
     steps:
       - name: checkout mfem
@@ -139,7 +142,27 @@ jobs:
         with:
           fetch-depth: 0
 
+      - name: check for pull request
+        id: check_pr
+        if: github.event_name == 'push'
+        env:
+          GH_TOKEN: ${{ github.token }}
+        run: |
+          pr_exists=$(gh pr list --repo "$GITHUB_REPOSITORY" \
+            --head "$GITHUB_REF_NAME" \
+            --state open \
+            --json number \
+            --jq 'length > 0')
+          echo "pr_exists=$pr_exists" >> "$GITHUB_OUTPUT"
+
       - name: branch-history
+        id: branch_history
+        if: |
+          (github.event_name == 'pull_request' ||
+           github.event_name == 'workflow_dispatch' ||
+           steps.check_pr.outputs.pr_exists == 'false')
+        continue-on-error: ${{ contains(github.event.pull_request.labels.*.name,
+                                        'branch-history-override') }}
         run: |
           # We override origin to make sure we point to the main repo.
           # This is to have consistent test results on PRs from forks.
@@ -147,3 +170,9 @@ jobs:
           git remote add origin https://github.com/mfem/mfem.git
           git checkout -b gh-actions-branch-history
           ./config/githooks/pre-push --history
+
+      - name: report branch-history override
+        if: steps.branch_history.outcome == 'failure'
+        run: |
+          echo "::warning::branch-history check failed, but the" \
+            "'branch-history-override' label is set."

config/githooks/README.md

@@ -39,3 +39,8 @@ when a picture was added for documentation.
 If that is the case, make sure the failure is indeed justified, and rerun the
 push command with the `--no-verify` option. This will skip the hooks, allowing
 you to push those changes.
+
+The `branch-history` check is run automatically through GitHub Actions. If a
+branch is known to have a large number of changes that are legitimate, the
+check can be overridden by setting  the label 'branch-history-override' on the
+pull request.

general/optparser.cpp

@@ -12,7 +12,9 @@
 #include "optparser.hpp"
 #include "../linalg/vector.hpp"
 #include "../general/communication.hpp"
+#include <cerrno>
 #include <cctype>
+#include <cstdlib>
 
 namespace mfem
 {
@@ -48,6 +50,19 @@ int isValidAsInt(char * s)
    return 1;
 }
 
+int isValidAsLongLong(char * s)
+{
+   if (s == NULL || *s == '\0')
+   {
+      return 0;
+   }
+
+   char *end = NULL;
+   errno = 0;
+   std::strtoll(s, &end, 10);
+   return errno == 0 && end != s && *end == '\0';
+}
+
 int isValidAsDouble(char * s)
 {
    // A valid floating point number for atof using the "C" locale is formed by
@@ -200,6 +215,10 @@ void OptionsParser::Parse()
                   isValid = isValidAsInt(argv[i]);
                   *(int *)(options[j].var_ptr) = atoi(argv[i++]);
                   break;
+               case LONG:
+                  isValid = isValidAsLongLong(argv[i]);
+                  *(long long *)(options[j].var_ptr) = std::strtoll(argv[i++], NULL, 10);
+                  break;
                case DOUBLE:
                   isValid = isValidAsDouble(argv[i]);
                   *(real_t *)(options[j].var_ptr) = atof(argv[i++]);
@@ -279,6 +298,10 @@ void OptionsParser::WriteValue(const Option &opt, std::ostream &os)
          os << *(int *)(opt.var_ptr);
          break;
 
+      case LONG:
+         os << *(long long *)(opt.var_ptr);
+         break;
+
       case DOUBLE:
          os << *(real_t *)(opt.var_ptr);
          break;
@@ -410,7 +433,7 @@ void OptionsParser::PrintHelp(ostream &os) const
    static const char *line_sep = "";
    static const char *types[] = { " <int>", " <double>", " <string>",
                                   " <string>", "", "", " '<int>...'",
-                                  " '<double>...'"
+                                  " '<double>...'", " <long long>"
                                 };
 
    os << indent << "-h" << seprtr << "--help" << descr_sep

general/optparser.hpp

@@ -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)

linalg/CMakeLists.txt

@@ -76,6 +76,12 @@ list(APPEND HDRS
   ttensor.hpp
   tensor.hpp
   vector.hpp
+  chpt/chpt.hpp
+  chpt/dynamic_checkpointing.hpp
+  chpt/file_checkpoint_storage.hpp
+  chpt/fixed_slot_checkpoint_storage.hpp
+  chpt/revolve_checkpointing.hpp
+  chpt/revolve_file_storage.hpp
   )
 
 if (MFEM_USE_MPI)

linalg/chpt/chpt.hpp

@@ -0,0 +1,21 @@
+// 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_CHPT_HPP
+#define MFEM_CHPT_HPP
+
+#include "dynamic_checkpointing.hpp"
+#include "file_checkpoint_storage.hpp"
+#include "fixed_slot_checkpoint_storage.hpp"
+#include "revolve_checkpointing.hpp"
+#include "revolve_file_storage.hpp"
+
+#endif // MFEM_CHPT_HPP

linalg/chpt/dynamic_checkpointing.hpp

@@ -0,0 +1,309 @@
+#ifndef MFEM_DYNAMIC_CHECKPOINTING_HPP
+#define MFEM_DYNAMIC_CHECKPOINTING_HPP
+
+#pragma once
+
+
+#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();
+         MFEM_VERIFY(!slots_[id],
+                     "InMemoryCheckpointStorage: free slot unexpectedly occupied.");
+
+         // Reuse the handle ID only; erased Snapshot objects are not reused.
+         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 (Wang–Moin–Iaccarino 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> info;
+      info.reserve(cps_.size());
+      for (const auto &kv : cps_)
+      {
+         info.push_back({kv.first, kv.second.level, storage_->IsValid(kv.second.h)});
+      }
+      return info;
+   }
+
+   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

linalg/chpt/file_checkpoint_storage.hpp

@@ -0,0 +1,304 @@
+#ifndef MFEM_FILE_CHECKPOINT_STORAGE_HPP
+#define MFEM_FILE_CHECKPOINT_STORAGE_HPP
+
+#pragma once
+#include "../vector.hpp"
+
+#include <cerrno>
+#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.HostRead();
+         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.HostWrite();
+         is.read(reinterpret_cast<char*>(data),
+                 (std::streamsize)(n * (std::int64_t)sizeof(mfem::real_t)));
+         v.Read();
+         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 || errno == ENOENT,
+                     "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

linalg/chpt/fixed_slot_checkpoint_storage.hpp

@@ -0,0 +1,408 @@
+#ifndef MFEM_FIXED_SLOT_CHECKPOINT_STORAGE_HPP
+#define MFEM_FIXED_SLOT_CHECKPOINT_STORAGE_HPP
+
+#pragma once
+
+#include "../vector.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 &snap) const
+   {
+      std::memcpy(&snap, 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.HostRead(), SlotBytes());
+   }
+
+   void Unpack(const void *src, mfem::Vector &v) const
+   {
+      v.SetSize(n_);
+      std::memcpy(v.HostWrite(), 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

linalg/chpt/revolve_checkpointing.hpp

@@ -0,0 +1,537 @@
+#ifndef MFEM_REVOLVE_CHECKPOINTING_HPP
+#define MFEM_REVOLVE_CHECKPOINTING_HPP
+
+#pragma once
+
+#include <cstdint>
+#include <cstring>
+#include <fstream>
+#include <string>
+#include <vector>
+#include <algorithm>
+
+namespace mfem
+{
+
+// -------------------------------
+// Fixed-slot storage: Memory
+// -------------------------------
+class RevolveFixedSlotMemoryStorage
+{
+private:
+   int max_slots_ = 0;
+   size_t slot_bytes_ = 0;
+   std::vector<uint8_t> data_;
+
+public:
+   RevolveFixedSlotMemoryStorage(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 RevolveFixedSlotFileStorage
+{
+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:
+   RevolveFixedSlotFileStorage(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");
+   }
+
+   ~RevolveFixedSlotFileStorage()
+   {
+      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;
+
+   // 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), 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

linalg/chpt/revolve_file_storage.hpp

@@ -0,0 +1,362 @@
+#ifndef MFEM_REVOLVE_SEPARATE_FILE_STORAGE_HPP
+#define MFEM_REVOLVE_SEPARATE_FILE_STORAGE_HPP
+
+#pragma once
+
+#include "../vector.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

linalg/linalg.hpp

@@ -40,6 +40,7 @@
 #include "filteredsolver.hpp"
 #include "ordering.hpp"
 #include "particlevector.hpp"
+#include "chpt/chpt.hpp"
 
 #ifdef MFEM_USE_AMGX
 #include "amgxsolver.hpp"

makefile

@@ -811,6 +811,7 @@ MINIAPPS_SUBDIRS = dpg/util hooke/operators hooke/preconditioners \
 FORMAT_FILES += $(foreach dir,$(TESTS_SUBDIRS),tests/$(dir)/*.?pp)
 FORMAT_FILES += $(foreach dir,$(UNIT_TESTS_SUBDIRS),tests/unit/$(dir)/*.?pp)
 FORMAT_FILES += $(foreach dir,$(MINIAPPS_SUBDIRS),miniapps/$(dir)/*.?pp)
+FORMAT_FILES += linalg/chpt/*.?pp miniapps/mtop/chpt/*.?pp
 FORMAT_FILES += config/cmake/config.hpp.in config/config.hpp.in mfem*.hpp
 FORMAT_EXCLUDE = general/tinyxml2.cpp tests/unit/catch.hpp
 FORMAT_LIST = $(filter-out $(FORMAT_EXCLUDE),$(wildcard $(FORMAT_FILES)))

miniapps/CMakeLists.txt

@@ -31,6 +31,7 @@ add_subdirectory(hdiv-linear-solver)
 add_subdirectory(hooke)
 add_subdirectory(meshing)
 add_subdirectory(mtop)
+add_subdirectory(mtop/chpt)
 add_subdirectory(multidomain)
 add_subdirectory(nurbs)
 add_subdirectory(parelag)

miniapps/mtop/chpt/CMakeLists.txt

@@ -0,0 +1,68 @@
+# 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_ex_scalar_chpt
+                 MAIN ex_scalar_field.cpp
+                 ${CHPT_COMMON_FILES}
+                 LIBRARIES mfem)
+
+add_mfem_miniapp(mtop_ex_adapt_chpt
+                 MAIN ex_scalar_adaptive.cpp
+                 ${CHPT_COMMON_FILES}
+                 LIBRARIES mfem)
+
+
+add_mfem_miniapp(mtop_ex_vector_chpt
+                 MAIN ex_vector_field.cpp
+                 ${CHPT_COMMON_FILES}
+                 LIBRARIES mfem)
+
+add_mfem_miniapp(mtop_ex_fixedslot
+                 MAIN ex_fixedslot_scalar_vector.cpp
+                 ${CHPT_COMMON_FILES}
+                 LIBRARIES mfem)
+
+# add_mfem_miniapp(mtop_ex_segmentio
+#                  MAIN ex_segmentio_scalar_vector_mfem.cpp
+#                  ${CHPT_COMMON_FILES}
+#                  LIBRARIES mfem)
+
+add_mfem_miniapp(mtop_ex_file
+                 MAIN ex_fileio_scalar_vector_mfem.cpp
+                 ${CHPT_COMMON_FILES}
+                 LIBRARIES mfem)
+
+
+add_mfem_miniapp(mtop_ex_revolve
+                 MAIN ex_revolve.cpp
+                 ${CHPT_COMMON_FILES}
+                 LIBRARIES mfem)
+
+add_mfem_miniapp(mtop_ex_dyn
+                 MAIN ex_dynamic_chpt.cpp
+                 ${CHPT_COMMON_FILES}
+                 LIBRARIES mfem)
+
+add_mfem_miniapp(mtop_ex_rev_fixedslot
+                 MAIN ex_revolve_chpt_fixedslot.cpp
+                 ${CHPT_COMMON_FILES}
+                 LIBRARIES mfem)
+

miniapps/mtop/chpt/ex_dynamic_chpt.cpp

@@ -0,0 +1,718 @@
+#include "mfem.hpp"
+
+#include <cmath>
+#include <iomanip>
+#include <string>
+
+using namespace mfem;
+
+
+/************************************************************************
+ * The forward state consists of vector and several additional values.
+ * The goal is to demonstrate 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 &snap) const
+   {
+      MFEM_VERIFY(src != nullptr, "Unpack: src is null.");
+
+      const unsigned char *b = static_cast<const unsigned char*>(src);
+
+      std::memcpy(&snap.time, b + 0*sizeof(mfem::real_t), sizeof(mfem::real_t));
+      std::memcpy(&snap.obj,  b + 1*sizeof(mfem::real_t), sizeof(mfem::real_t));
+      snap.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)
+      {
+         mfem::out<<"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)
+      {
+         mfem::out<<"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<<"\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 storage’s 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)
+      {
+         mfem::out<<"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)
+      {
+         mfem::out<<"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 that’s 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)
+      {
+         mfem::out<<"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)
+      {
+         mfem::out<<"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)
+      {
+         mfem::out<<"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)
+      {
+         mfem::out<<"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;
+}

miniapps/mtop/chpt/ex_fileio_scalar_vector_mfem.cpp

@@ -0,0 +1,313 @@
+#include "mfem.hpp"
+
+#include <cmath>
+#include <iomanip>
+
+using namespace mfem;
+
+static double LogisticStep(double u, double alpha, double dt)
+{
+   return u + dt * alpha * u * (1.0 - u);
+}
+
+static double LogisticJac(double u, double alpha, double dt)
+{
+   // d/du [u + dt*alpha*u*(1-u)] = 1 + dt*alpha*(1 - 2u)
+   return 1.0 + dt * alpha * (1.0 - 2.0*u);
+}
+
+#if MFEM_HAVE_FILESYSTEM
+static void PurgeDir(const std::string &dir)
+{
+   std::error_code ec;
+   mfem_fs::remove_all(mfem_fs::path(dir), ec);
+   mfem_fs::create_directories(mfem_fs::path(dir), ec);
+   MFEM_VERIFY(!ec, "Failed to purge/create directory.");
+}
+#endif
+
+int main(int argc, char *argv[])
+{
+   // ---------- shared options ----------
+   int s = 3;
+   double dt0 = 0.02;
+   double omega = 0.2;
+   double Tfinal = 1.0;
+   bool keep_files = false;
+   bool purge_dirs = true;
+
+   // scalar problem params
+   double alpha = 2.0;
+   double u0 = 0.2;
+   double target_s = 0.7;
+   double eps = 1e-7;
+
+   // vector problem params
+   int n = 16;
+   double target_v = 0.7;
+
+   std::string dir_scalar = "chk_scalar_files";
+   std::string dir_vector = "chk_vector_files";
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&s, "-s", "--checkpoints",
+                  "Checkpoint budget s (real checkpoints).");
+   args.AddOption(&dt0, "-dt0", "--dt0",
+                  "Base dt for dt(i)=dt0*(1+0.5*sin(omega*i)).");
+   args.AddOption(&omega, "-om", "--omega",
+                  "Omega for dt(i)=dt0*(1+0.5*sin(omega*i)).");
+   args.AddOption(&Tfinal, "-T", "--tfinal",
+                  "Stop when accumulated time reaches Tfinal.");
+   args.AddOption(&keep_files, "-k", "--keep-files", "-nk", "--no-keep-files",
+                  "Keep checkpoint files (debug).");
+   args.AddOption(&purge_dirs, "-p", "--purge-dirs", "-np", "--no-purge-dirs",
+                  "Purge checkpoint directories at start (recommended).");
+
+   args.AddOption(&alpha, "-a", "--alpha", "Logistic growth alpha.");
+   args.AddOption(&u0, "-u0", "--u0", "Scalar initial u0.");
+   args.AddOption(&target_s, "-ts", "--target-scalar", "Scalar target.");
+   args.AddOption(&eps, "-eps", "--fd-eps",
+                  "FD epsilon (scalar and vector directional).");
+
+   args.AddOption(&n, "-n", "--size", "Vector dimension.");
+   args.AddOption(&target_v, "-tv", "--target-vector",
+                  "Vector target value per component.");
+
+   args.AddOption(&dir_scalar, "-ds", "--dir-scalar",
+                  "Directory for scalar checkpoint files.");
+   args.AddOption(&dir_vector, "-dv", "--dir-vector",
+                  "Directory for vector checkpoint files.");
+
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(mfem::out);
+      return 1;
+   }
+   args.PrintOptions(mfem::out);
+
+   MFEM_VERIFY(s > 0, "Need s > 0.");
+   MFEM_VERIFY(dt0 > 0.0, "Need dt0 > 0.");
+   MFEM_VERIFY(Tfinal > 0.0, "Need Tfinal > 0.");
+
+#if MFEM_HAVE_FILESYSTEM
+   if (purge_dirs)
+   {
+      PurgeDir(dir_scalar);
+      PurgeDir(dir_vector);
+   }
+#else
+   MFEM_VERIFY(!purge_dirs, "purge_dirs requires <filesystem> support.");
+#endif
+
+   using StepS =
+      mfem::DynamicCheckpointing<double, mfem::FileCheckpointStorage<double>>::Step;
+   auto dt_func = [&](StepS i)
+   {
+      const double dt = dt0 * (1.0 + 0.5 * std::sin(omega * double(i)));
+      MFEM_VERIFY(dt > 0.0, "dt_func produced non-positive dt.");
+      return dt;
+   };
+
+   mfem::out << std::setprecision(15);
+
+   // ============================================================
+   // A) Scalar run (double) with FileCheckpointStorage
+   // ============================================================
+   {
+      mfem::FileCheckpointStorage<double> storage(dir_scalar, "ckpt_", ".bin", true,
+                                                  keep_files);
+      mfem::DynamicCheckpointing<double, mfem::FileCheckpointStorage<double>> ckpt(s,
+                                                                                   storage);
+
+      auto make_snapshot = [](const double &u) { return u; };
+      auto restore_snapshot = [](const double &snap, double &state) { state = snap; };
+
+      auto primal_step = [&](double &u, StepS i)
+      {
+         const double dt = dt_func(i);
+         u = LogisticStep(u, alpha, dt);
+      };
+
+      auto adjoint_step = [&](double &lambda, const double &u_i, StepS i)
+      {
+         const double dt = dt_func(i);
+         const double dF = LogisticJac(u_i, alpha, dt);
+         lambda *= dF;
+      };
+
+      double u = u0;
+      double t_phys = 0.0;
+      StepS i = 0;
+
+      while (t_phys < Tfinal)
+      {
+         ckpt.ForwardStep(i, u, primal_step, make_snapshot);
+         t_phys += dt_func(i);
+         ++i;
+      }
+      const StepS m = i;
+      const double u_m = u;
+      const double J = 0.5 * (u_m - target_s) * (u_m - target_s);
+
+      double lambda = (u_m - target_s); // terminal adjoint
+      double u_work = 0.0;
+
+      for (StepS j = m - 1; j >= 0; --j)
+      {
+         ckpt.BackwardStep(j, lambda, u_work,
+                           primal_step, adjoint_step,
+                           make_snapshot, restore_snapshot);
+         if (j == 0) { break; }
+      }
+
+      const double dJ_du0_adj = lambda;
+
+      // FD check
+      auto forward_only_J = [&](double u_init)
+      {
+         double uu = u_init;
+         double tt = 0.0;
+         StepS k = 0;
+         while (tt < Tfinal)
+         {
+            uu = LogisticStep(uu, alpha, dt_func(k));
+            tt += dt_func(k);
+            ++k;
+         }
+         const double r = uu - target_s;
+         return 0.5 * r * r;
+      };
+
+      const double Jp = forward_only_J(u0 + eps);
+      const double Jm = forward_only_J(u0 - eps);
+      const double dJ_du0_fd = (Jp - Jm) / (2.0 * eps);
+
+      mfem::out << "\n[FileIO][Scalar]\n";
+      mfem::out << "  m steps        = " << m << "\n";
+      mfem::out << "  u_m            = " << u_m << "\n";
+      mfem::out << "  J              = " << J << "\n";
+      mfem::out << "  dJ/du0 adjoint  = " << dJ_du0_adj << "\n";
+      mfem::out << "  dJ/du0 FD       = " << dJ_du0_fd << "\n";
+      mfem::out << "  abs err         = " << std::abs(dJ_du0_adj - dJ_du0_fd) << "\n";
+   }
+
+   // ============================================================
+   // B) Vector run (mfem::Vector) with FileCheckpointStorage
+   // ============================================================
+   {
+      mfem::FileCheckpointStorage<mfem::Vector> storage(dir_vector, "ckpt_", ".bin",
+                                                        true, keep_files);
+      mfem::DynamicCheckpointing<mfem::Vector, mfem::FileCheckpointStorage<mfem::Vector>>
+            ckpt(s, storage);
+
+      using StepV =
+         mfem::DynamicCheckpointing<mfem::Vector, mfem::FileCheckpointStorage<mfem::Vector>>::Step;
+
+      auto dt_func_v = [&](StepV i)
+      {
+         const double dt = dt0 * (1.0 + 0.5 * std::sin(omega * double(i)));
+         MFEM_VERIFY(dt > 0.0, "dt_func produced non-positive dt.");
+         return dt;
+      };
+
+      auto make_snapshot = [](const mfem::Vector &u) { return u; };
+      auto restore_snapshot = [](const mfem::Vector &snap, mfem::Vector &state) { state = snap; };
+
+      auto primal_step = [&](mfem::Vector &u, StepV i)
+      {
+         const double dt = dt_func_v(i);
+         for (int j = 0; j < u.Size(); ++j)
+         {
+            u[j] = LogisticStep(u[j], alpha, dt);
+         }
+      };
+
+      auto adjoint_step = [&](mfem::Vector &lambda, const mfem::Vector &u_i, StepV i)
+      {
+         const double dt = dt_func_v(i);
+         MFEM_ASSERT(lambda.Size() == u_i.Size(), "Size mismatch.");
+         for (int j = 0; j < lambda.Size(); ++j)
+         {
+            lambda[j] *= LogisticJac(u_i[j], alpha, dt);
+         }
+      };
+
+      mfem::Vector u0v(n), ut(n);
+      for (int j = 0; j < n; ++j) { u0v[j] = 0.2 + 0.05 * std::cos(0.7*(j+1)); }
+      ut = target_v;
+
+      mfem::Vector u = u0v;
+      double t_phys = 0.0;
+      StepV i = 0;
+
+      while (t_phys < Tfinal)
+      {
+         ckpt.ForwardStep(i, u, primal_step, make_snapshot);
+         t_phys += dt_func_v(i);
+         ++i;
+      }
+      const StepV m = i;
+      const mfem::Vector u_m = u;
+
+      mfem::Vector diff(u_m);
+      diff -= ut;
+      const double J = 0.5 * mfem::InnerProduct(diff, diff);
+
+      // backward
+      mfem::Vector lambda = diff; // terminal adjoint
+      mfem::Vector u_work(n);
+
+      for (StepV j = m - 1; j >= 0; --j)
+      {
+         ckpt.BackwardStep(j, lambda, u_work,
+                           primal_step, adjoint_step,
+                           make_snapshot, restore_snapshot);
+         if (j == 0) { break; }
+      }
+      const mfem::Vector &grad_u0 = lambda;
+
+      // directional FD check
+      mfem::Vector v(n);
+      for (int j = 0; j < n; ++j) { v[j] = std::sin(0.3*(j+1)) + 0.1; }
+
+      auto forward_only_J = [&](const mfem::Vector &u_init)
+      {
+         mfem::Vector uu = u_init;
+         double tt = 0.0;
+         StepV k = 0;
+         while (tt < Tfinal)
+         {
+            primal_step(uu, k);
+            tt += dt_func_v(k);
+            ++k;
+         }
+         mfem::Vector dd(uu);
+         dd -= ut;
+         return 0.5 * mfem::InnerProduct(dd, dd);
+      };
+
+      mfem::Vector u_plus(u0v), u_minus(u0v);
+      u_plus.Add(eps, v);
+      u_minus.Add(-eps, v);
+
+      const double Jp = forward_only_J(u_plus);
+      const double Jm = forward_only_J(u_minus);
+      const double dJ_dir_fd = (Jp - Jm) / (2.0 * eps);
+      const double dJ_dir_adj = mfem::InnerProduct(grad_u0, v);
+
+      mfem::out << "\n[FileIO][Vector]\n";
+      mfem::out << "  m steps             = " << m << "\n";
+      mfem::out << "  J                   = " << J << "\n";
+      mfem::out << "  v·grad adjoint      = " << dJ_dir_adj << "\n";
+      mfem::out << "  FD directional      = " << dJ_dir_fd << "\n";
+      mfem::out << "  abs err             = " << std::abs(dJ_dir_adj - dJ_dir_fd) <<
+                "\n";
+      mfem::out << "  ||u_m||_2           = " << u_m.Norml2() << "\n";
+      mfem::out << "  ||u_m - u_target||2 = " << diff.Norml2() << "\n";
+   }
+
+   mfem::out << "\nDone.\n";
+   return 0;
+}
+

miniapps/mtop/chpt/ex_fixedslot_scalar_vector.cpp

@@ -0,0 +1,385 @@
+#include "mfem.hpp"
+
+
+#include <cmath>
+#include <iomanip>
+#include <string>
+
+using namespace mfem;
+
+// Logistic ODE (scalar or componentwise):
+//   du/dt = alpha*u*(1-u)
+// Explicit Euler:
+//   u_{i+1} = u_i + dt_i * alpha*u_i*(1-u_i)
+//
+// Objective:
+//   Scalar: J = 0.5*(u_m - target)^2
+//   Vector: J = 0.5*||u_m - u_target||^2
+//
+// Discrete adjoint:
+//   dF/du = 1 + dt_i*alpha*(1 - 2*u_i)
+//   lambda_i = (dF/du at u_i) * lambda_{i+1}
+
+static inline double LogisticStep(const double u, const double alpha,
+                                  const double dt)
+{
+   return u + dt * alpha * u * (1.0 - u);
+}
+
+static inline double LogisticJac(const double u, const double alpha,
+                                 const double dt)
+{
+   return 1.0 + dt * alpha * (1.0 - 2.0*u);
+}
+
+// --------------------------
+// Scalar run (double)
+// --------------------------
+template <typename Storage>
+static void RunScalarFixedSlot(int s,
+                               Storage &storage,
+                               double alpha,
+                               double dt0,
+                               double omega,
+                               double Tfinal,
+                               double u0,
+                               double target,
+                               double eps)
+{
+   using CKPT = mfem::DynamicCheckpointing<double, Storage>;
+   using Step = typename CKPT::Step;
+
+   CKPT ckpt(s, storage);
+
+   auto dt_func = [&](Step i)
+   {
+      const double dt = dt0 * (1.0 + 0.5 * std::sin(omega * double(i)));
+      MFEM_VERIFY(dt > 0.0, "dt_func produced non-positive dt.");
+      return dt;
+   };
+
+   auto make_snapshot = [](const double &u) { return u; };
+   auto restore_snapshot = [](const double &snap, double &state) { state = snap; };
+
+   auto primal_step = [&](double &u, Step i)
+   {
+      const double dt = dt_func(i);
+      u = LogisticStep(u, alpha, dt);
+   };
+
+   auto adjoint_step = [&](double &lambda, const double &u_i, Step i)
+   {
+      const double dt = dt_func(i);
+      lambda *= LogisticJac(u_i, alpha, dt);
+   };
+
+   // Forward sweep (unknown m; stop on accumulated physical time)
+   double u = u0;
+   double t_phys = 0.0;
+   Step i = 0;
+
+   while (t_phys < Tfinal)
+   {
+      ckpt.ForwardStep(i, u, primal_step, make_snapshot);
+      t_phys += dt_func(i);
+      ++i;
+   }
+
+   const Step m = i;
+   const double u_m = u;
+   const double J = 0.5 * (u_m - target) * (u_m - target);
+
+   // Backward sweep
+   MFEM_VERIFY(m > 0, "Scalar run produced m=0 steps.");
+
+   double lambda = (u_m - target); // terminal condition dJ/du_m
+   double u_work = 0.0;
+
+   for (Step j = m - 1; j >= 0; --j)
+   {
+      ckpt.BackwardStep(j, lambda, u_work,
+                        primal_step, adjoint_step,
+                        make_snapshot, restore_snapshot);
+      if (j == 0) { break; }
+   }
+   const double dJ_du0_adj = lambda;
+
+   // FD check
+   auto forward_only_J = [&](double u_init)
+   {
+      double uu = u_init;
+      double tt = 0.0;
+      Step k = 0;
+      while (tt < Tfinal)
+      {
+         uu = LogisticStep(uu, alpha, dt_func(k));
+         tt += dt_func(k);
+         ++k;
+      }
+      const double r = uu - target;
+      return 0.5 * r * r;
+   };
+
+   const double Jp = forward_only_J(u0 + eps);
+   const double Jm = forward_only_J(u0 - eps);
+   const double dJ_du0_fd = (Jp - Jm) / (2.0 * eps);
+
+   mfem::out << "\n[FixedSlot][Scalar]\n";
+   mfem::out << "  m steps        = " << m << "\n";
+   mfem::out << "  t_phys         = " << t_phys << "\n";
+   mfem::out << "  u_m            = " << u_m << "\n";
+   mfem::out << "  J              = " << J << "\n";
+   mfem::out << "  dJ/du0 adjoint  = " << dJ_du0_adj << "\n";
+   mfem::out << "  dJ/du0 FD       = " << dJ_du0_fd << "\n";
+   mfem::out << "  abs err         = " << std::abs(dJ_du0_adj - dJ_du0_fd) << "\n";
+}
+
+// --------------------------
+// Vector run (mfem::Vector)
+// --------------------------
+template <typename Storage>
+static void RunVectorFixedSlot(int s,
+                               Storage &storage,
+                               int n,
+                               double alpha,
+                               double dt0,
+                               double omega,
+                               double Tfinal,
+                               double target_val,
+                               double eps)
+{
+   using CKPT = mfem::DynamicCheckpointing<mfem::Vector, Storage>;
+   using Step = typename CKPT::Step;
+
+   CKPT ckpt(s, storage);
+
+   auto dt_func = [&](Step i)
+   {
+      const double dt = dt0 * (1.0 + 0.5 * std::sin(omega * double(i)));
+      MFEM_VERIFY(dt > 0.0, "dt_func produced non-positive dt.");
+      return dt;
+   };
+
+   auto make_snapshot = [](const mfem::Vector &u) { return u; };
+   auto restore_snapshot = [](const mfem::Vector &snap, mfem::Vector &state) { state = snap; };
+
+   auto primal_step = [&](mfem::Vector &u, Step i)
+   {
+      const double dt = dt_func(i);
+      for (int j = 0; j < u.Size(); ++j)
+      {
+         u[j] = LogisticStep(u[j], alpha, dt);
+      }
+   };
+
+   auto adjoint_step = [&](mfem::Vector &lambda, const mfem::Vector &u_i, Step i)
+   {
+
+      const double dt = dt_func(i);
+      MFEM_ASSERT(lambda.Size() == u_i.Size(), "lambda and u_i size mismatch.");
+      for (int j = 0; j < lambda.Size(); ++j)
+      {
+         lambda[j] *= LogisticJac(u_i[j], alpha, dt);
+      }
+   };
+
+   // Initial condition and target
+   mfem::Vector u0(n), u_target(n);
+   for (int j = 0; j < n; ++j)
+   {
+      u0[j] = 0.2 + 0.05 * std::cos(0.7 * (j + 1));
+   }
+   u_target = target_val;
+
+   // Forward sweep
+   mfem::Vector u = u0;
+   double t_phys = 0.0;
+   Step i = 0;
+
+   while (t_phys < Tfinal)
+   {
+      ckpt.ForwardStep(i, u, primal_step, make_snapshot);
+      t_phys += dt_func(i);
+      ++i;
+   }
+
+   const Step m = i;
+   const mfem::Vector u_m = u;
+
+   mfem::Vector diff(u_m);
+   diff -= u_target;
+   const double J = 0.5 * mfem::InnerProduct(diff, diff);
+
+   // Backward sweep
+   MFEM_VERIFY(m > 0, "Vector run produced m=0 steps.");
+
+   mfem::Vector lambda = diff; // terminal adjoint = dJ/du_m
+   mfem::Vector u_work(n);
+
+   for (Step j = m - 1; j >= 0; --j)
+   {
+      ckpt.BackwardStep(j, lambda, u_work,
+                        primal_step, adjoint_step,
+                        make_snapshot, restore_snapshot);
+      if (j == 0) { break; }
+   }
+   const mfem::Vector &grad_u0 = lambda;
+
+   // Directional FD check
+   mfem::Vector v(n);
+   for (int j = 0; j < n; ++j) { v[j] = std::sin(0.3 * (j + 1)) + 0.1; }
+
+   auto forward_only_J = [&](const mfem::Vector &u_init)
+   {
+      mfem::Vector uu = u_init;
+      double tt = 0.0;
+      Step k = 0;
+      while (tt < Tfinal)
+      {
+         primal_step(uu, k);
+         tt += dt_func(k);
+         ++k;
+      }
+      mfem::Vector dd(uu);
+      dd -= u_target;
+      return 0.5 * mfem::InnerProduct(dd, dd);
+   };
+
+   mfem::Vector u_plus(u0), u_minus(u0);
+   u_plus.Add(eps, v);
+   u_minus.Add(-eps, v);
+
+   const double Jp = forward_only_J(u_plus);
+   const double Jm = forward_only_J(u_minus);
+   const double dJ_dir_fd = (Jp - Jm) / (2.0 * eps);
+   const double dJ_dir_adj = mfem::InnerProduct(grad_u0, v);
+
+   mfem::out << "\n[FixedSlot][Vector]\n";
+   mfem::out << "  n              = " << n << "\n";
+   mfem::out << "  m steps         = " << m << "\n";
+   mfem::out << "  t_phys          = " << t_phys << "\n";
+   mfem::out << "  J               = " << J << "\n";
+   mfem::out << "  ||u_m||_2        = " << u_m.Norml2() << "\n";
+   mfem::out << "  ||u_m-u_target|| = " << diff.Norml2() << "\n";
+   mfem::out << "  v·grad adjoint   = " << dJ_dir_adj << "\n";
+   mfem::out << "  FD directional   = " << dJ_dir_fd << "\n";
+   mfem::out << "  abs err          = " << std::abs(dJ_dir_adj - dJ_dir_fd) <<
+             "\n";
+}
+
+int main(int argc, char *argv[])
+{
+   // Backend selection:
+   //   0 = fixed-slot memory (single RAM block)
+   //   1 = fixed-slot file   (single file with fixed offsets)
+   int backend = 0;
+
+   // Common parameters
+   int s = 8;                 // number of REAL stored checkpoints
+   double alpha = 2.0;
+   double dt0 = 0.02;
+   double omega = 0.2;
+   double Tfinal = 1.0;
+   double eps = 1e-7;
+
+   // Scalar parameters
+   double u0 = 0.2;
+   double target_s = 0.7;
+
+   // Vector parameters
+   int n = 64;                // must stay fixed for fixed-slot vector packing
+   double target_v = 0.7;
+
+   // File backend parameters
+   std::string scalar_file = "scalar_fixedslots.bin";
+   std::string vector_file = "vector_fixedslots.bin";
+   bool truncate_files = true;
+   bool flush_on_store = true; // safer for demo correctness with iostreams
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&backend, "-b", "--backend",
+                  "Backend: 0=memory fixed-slots, 1=file fixed-slots (single file).");
+   args.AddOption(&s, "-s", "--checkpoints",
+                  "Checkpoint budget s (real checkpoints).");
+
+   args.AddOption(&alpha, "-a", "--alpha", "Logistic growth alpha.");
+   args.AddOption(&dt0, "-dt0", "--dt0",
+                  "Base dt for dt(i)=dt0*(1+0.5*sin(omega*i)).");
+   args.AddOption(&omega, "-om", "--omega",
+                  "Omega for dt(i)=dt0*(1+0.5*sin(omega*i)).");
+   args.AddOption(&Tfinal, "-T", "--tfinal",
+                  "Stop when accumulated time reaches Tfinal.");
+   args.AddOption(&eps, "-eps", "--fd-eps", "FD epsilon.");
+
+   args.AddOption(&u0, "-u0", "--u0", "Scalar initial condition u0.");
+   args.AddOption(&target_s, "-ts", "--target-scalar", "Scalar target.");
+
+   args.AddOption(&n, "-n", "--size", "Vector dimension (fixed).");
+   args.AddOption(&target_v, "-tv", "--target-vector",
+                  "Vector target value per component.");
+
+   args.AddOption(&scalar_file, "-sf", "--scalar-file",
+                  "File for scalar fixed-slot storage.");
+   args.AddOption(&vector_file, "-vf", "--vector-file",
+                  "File for vector fixed-slot storage.");
+   args.AddOption(&truncate_files, "-tr", "--truncate", "-ntr", "--no-truncate",
+                  "Truncate checkpoint files on startup.");
+   args.AddOption(&flush_on_store, "-fl", "--flush", "-nfl", "--no-flush",
+                  "Flush on each Store() (demo-safety; slower).");
+
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(mfem::out);
+      return 1;
+   }
+   args.PrintOptions(mfem::out);
+
+   MFEM_VERIFY(s > 0, "Need s > 0.");
+   MFEM_VERIFY(n > 0, "Need n > 0.");
+   MFEM_VERIFY(dt0 > 0.0, "Need dt0 > 0.");
+   MFEM_VERIFY(Tfinal > 0.0, "Need Tfinal > 0.");
+
+   mfem::out << std::setprecision(15);
+
+   if (backend == 0)
+   {
+      mfem::out << "\nUsing fixed-slot MEMORY backend (single RAM block)\n";
+
+      // Scalar: fixed-size POD => trivial packer
+      mfem::FixedSlotMemoryCheckpointStorage<double> stor_s(s);
+      RunScalarFixedSlot(s, stor_s, alpha, dt0, omega, Tfinal, u0, target_s, eps);
+
+      // Vector: fixed-size packing (n must remain constant)
+      mfem::FixedVectorPacker packer(n);
+      mfem::FixedSlotMemoryCheckpointStorage<mfem::Vector, mfem::FixedVectorPacker>
+      stor_v(s, packer);
+      RunVectorFixedSlot(s, stor_v, n, alpha, dt0, omega, Tfinal, target_v, eps);
+   }
+   else if (backend == 1)
+   {
+      mfem::out <<
+                "\nUsing fixed-slot FILE backend (single file with fixed offsets)\n";
+
+      // Scalar file
+      mfem::FixedSlotFileCheckpointStorage<double> stor_s(scalar_file, s,
+                                                          mfem::TrivialFixedPacker<double>(),
+                                                          truncate_files,
+                                                          flush_on_store);
+      RunScalarFixedSlot(s, stor_s, alpha, dt0, omega, Tfinal, u0, target_s, eps);
+
+      // Vector file (fixed-size packing with n)
+      mfem::FixedVectorPacker packer(n);
+      mfem::FixedSlotFileCheckpointStorage<mfem::Vector, mfem::FixedVectorPacker>
+      stor_v(vector_file, s, packer, truncate_files, flush_on_store);
+      RunVectorFixedSlot(s, stor_v, n, alpha, dt0, omega, Tfinal, target_v, eps);
+   }
+   else
+   {
+      MFEM_ABORT("Unknown backend. Use -b 0 (memory) or -b 1 (file).");
+   }
+
+   mfem::out << "\nDone.\n";
+   return 0;
+}
+

miniapps/mtop/chpt/ex_revolve.cpp

@@ -0,0 +1,183 @@
+#include "mfem.hpp"
+
+using namespace mfem;
+
+// ============================================================
+// 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)
+{
+   mfem::out << "\n--- Scalar REVOLVE example ---\n";
+   mfem::out << "Nsteps=" << Nsteps << ", Ncheckpoints=" << Ncheck << ", factor="
+             <<
+             factor << "\n";
+
+   const size_t snap_bytes = sizeof(double);
+   RevolveFixedSlotMemoryStorage storage(Ncheck, snap_bytes);
+   FixedStepRevolveCheckpointing<RevolveFixedSlotMemoryStorage> 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);
+
+   mfem::out << "uN = " << uN << "\n";
+   mfem::out << "lambda0 (computed) = " << lambda << "\n";
+   mfem::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)
+{
+   mfem::out << "\n--- mfem::Vector REVOLVE example ---\n";
+   mfem::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);
+   RevolveFixedSlotMemoryStorage storage(Ncheck, snap_bytes);
+   FixedStepRevolveCheckpointing<RevolveFixedSlotMemoryStorage> 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);
+   }
+
+   mfem::out << "||uN||_2      = " << uN.Norml2() << "\n";
+   mfem::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(mfem::out);
+      return 1;
+   }
+   args.PrintOptions(mfem::out);
+
+   RunScalarExample(Nsteps, Ncheck, factor);
+   RunVectorExample(Nsteps, Ncheck, dim, factor);
+
+   return 0;
+}

miniapps/mtop/chpt/ex_revolve_chpt_fixedslot.cpp

@@ -0,0 +1,174 @@
+#include "mfem.hpp"
+
+
+#include <cmath>
+#include <iomanip>
+#include <string>
+
+using namespace mfem;
+
+
+// the forward state consists of vector and
+// several additional values
+struct State
+{
+   mfem::real_t time = 0.0;
+   mfem::real_t obj  = 0.0;
+   mfem::Vector v;
+};
+
+/// class to manage the Snaphot Layout
+class StateSnapshotLayout
+{
+public:
+   // Constructor: n - size of the vector in the state
+   explicit StateSnapshotLayout(int n) : n_(n)
+   {
+      MFEM_VERIFY(n_ > 0, "StateSnapshotLayout: n must be > 0.");
+   }
+
+   // returns the total size of the snapshot in mfem::real_t
+   int N() const { return n_; }
+
+   // returns the total number of bytes necessary to store the state
+   int Bytes() const
+   {
+      return (2 + n_) * (int)sizeof(mfem::real_t);
+   }
+
+   // Pack State -> bytes (host bytes)
+   void Pack(const State &s, unsigned char *dst, int bytes) const
+   {
+      MFEM_VERIFY(dst != nullptr, "Pack: dst is null.");
+      MFEM_VERIFY(bytes == Bytes(), "Pack: snapshot byte size mismatch.");
+      MFEM_VERIFY(s.v.Size() == n_, "Pack: State.v size mismatch.");
+
+      // Ensure we read from host even if MFEM uses device memory.
+      const mfem::real_t *vh = s.v.HostRead();
+
+      std::memcpy(dst + 0*sizeof(mfem::real_t), &s.time, sizeof(mfem::real_t));
+      std::memcpy(dst + 1*sizeof(mfem::real_t), &s.obj,  sizeof(mfem::real_t));
+      std::memcpy(dst + 2*sizeof(mfem::real_t),
+                  vh, (std::size_t)n_ * sizeof(mfem::real_t));
+   }
+
+   // Unpack bytes -> State (host write)
+   void Unpack(const unsigned char *src, int bytes, State &state) const
+   {
+      MFEM_VERIFY(src != nullptr, "Unpack: src is null.");
+      MFEM_VERIFY(bytes == Bytes(), "Unpack: snapshot byte size mismatch.");
+
+      if (state.v.Size() != n_) { state.v.SetSize(n_); }
+      mfem::real_t *vh = state.v.HostWrite();
+
+      std::memcpy(&state.time, src + 0*sizeof(mfem::real_t), sizeof(mfem::real_t));
+      std::memcpy(&state.obj,  src + 1*sizeof(mfem::real_t), sizeof(mfem::real_t));
+      std::memcpy(vh,
+                  src + 2*sizeof(mfem::real_t),
+                  (std::size_t)n_ * sizeof(mfem::real_t));
+   }
+
+private:
+   int n_ = 0;
+};
+
+/* The REVOLVE manager (as provided earlier) uses this callback style:
+
+        * make_snapshot(const State&, uint8_t* out, size_t bytes)
+        * restore_snapshot(State&, const uint8_t* in, size_t bytes)
+        * primal_step(int step, State&)
+        * adjoint_step(int step, const State&, AdjointState&)
+
+        and a storage backend with:
+        * Save(slot, bytes)
+        * Load(slot, bytes)
+
+    REVOLVE manager will:
+        *request storing snapshots into checkpoint slots (takeshot)
+        *request restoring a checkpoint slot (restore)
+        *request recomputation forward (advance)
+        *then request the next adjoint step (firsturn / youturn)
+        all while using only (Ncheck) stored checkpoints.
+
+    Fixed-step REVOLVE manager:
+        * Snapshot is raw bytes packed/unpacked by your callbacks.
+        * Storage is a fixed indexed array (Save(slot) / Load(slot)),
+               because REVOLVE addresses checkpoints by slot index.
+
+    Best use when Nsteps is known in the begining of the simulation.
+*/
+
+int main(int argc, char *argv[])
+{
+   const int n = 100/* fixed State.v size */;
+   StateSnapshotLayout layout(n);
+   const std::size_t snapshot_bytes = (std::size_t)layout.Bytes();
+
+   const int Nsteps = 20 /* known number of time steps */;
+   const int Ncheck = 5 /* number of checkpoints (snaps) */;
+
+   // Fixed-slot file backend (single file)
+   RevolveFixedSlotFileStorage storage("revolve_ckpts.bin", Ncheck,
+                                       snapshot_bytes);
+   FixedStepRevolveCheckpointing<RevolveFixedSlotFileStorage>
+   cktp(Nsteps, Ncheck, snapshot_bytes, storage);
+
+   // Memory backend (single block) alternative:
+   // RevolveFixedSlotMemoryStorage storage(Ncheck, snapshot_bytes);
+   // FixedStepRevolveCheckpointing<RevolveFixedSlotMemoryStorage>
+   //                        cktp(Nsteps, Ncheck, snapshot_bytes, storage);
+
+
+   auto make_snapshot = [&](const State &s, uint8_t *outb, std::size_t bytes)
+   {
+      MFEM_VERIFY(bytes == snapshot_bytes, "make_snapshot: byte size mismatch");
+      layout.Pack(s, reinterpret_cast<unsigned char*>(outb), (int)bytes);
+   };
+
+   auto restore_snapshot = [&](State &s, const uint8_t *inb, std::size_t bytes)
+   {
+      MFEM_VERIFY(bytes == snapshot_bytes, "restore_snapshot: byte size mismatch");
+      layout.Unpack(reinterpret_cast<const unsigned char*>(inb), (int)bytes, s);
+   };
+
+   // set the work state
+   State u_work;
+   u_work.v.SetSize(n); u_work.v=0.0;
+
+   State u;
+   u.v.SetSize(n); u.v=0.0; u.obj=0.0; u.time=0.0;
+   // init u.time, u.obj, u.v ...
+
+   auto primal_step = [&](int step, State &u_inout)
+   {
+      (void)step;
+      // advance u_inout -> u_{step+1}
+   };
+
+   struct AdjointState { /* ... */ };
+   AdjointState lambda;
+
+   auto adjoint_step = [&](int step, const State &u_step,
+                           AdjointState &lambda_inout)
+   {
+      (void)step;
+      // update lambda_{step+1} -> lambda_step using u_step
+   };
+
+   // Forward: i = 0..Nsteps-1
+   for (int i = 0; i < Nsteps; ++i)
+   {
+      cktp.ForwardStep(i, u, primal_step, make_snapshot);
+   }
+
+   // Reverse: i = Nsteps-1..0
+   for (int i = Nsteps - 1; i >= 0; --i)
+   {
+      cktp.BackwardStep(i, lambda, u_work,
+                        primal_step, adjoint_step,
+                        make_snapshot, restore_snapshot);
+   }
+
+   mfem::out << "\nDone.\n";
+   return 0;
+};

miniapps/mtop/chpt/ex_scalar_adaptive.cpp

@@ -0,0 +1,212 @@
+#include "mfem.hpp"
+
+#include <cmath>
+#include <iomanip>
+
+using namespace mfem;
+
+// Scalar logistic ODE:
+//   du/dt = alpha*u*(1-u)
+// Explicit Euler with variable dt(i):
+//   u_{i+1} = u_i + dt_i*alpha*u_i*(1-u_i)
+// Objective:
+//   J = 0.5*(u_m - target)^2
+// Discrete adjoint (scalar):
+//   lambda_m = (u_m - target)
+//   lambda_i = (dF/du at u_i)^T * lambda_{i+1}
+//   dF/du = 1 + dt_i*alpha*(1 - 2*u_i)
+
+int main(int argc, char *argv[])
+{
+   int s = 3;
+   double alpha = 2.0;
+   double dt0 = 0.001;
+   double Tfinal = 1.0;
+   double u0 = 0.2;
+   double target = 0.7;
+   double eps = 1e-7;
+   double err = 1e-4;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&s,      "-s",   "--checkpoints",
+                  "Checkpoint budget s (real checkpoints).");
+   args.AddOption(&alpha,  "-a",   "--alpha",       "Logistic growth alpha.");
+   args.AddOption(&dt0,    "-dt0", "--dt0",         "Base for the time step dt.");
+   args.AddOption(&Tfinal, "-T",   "--tfinal",
+                  "Terminate when accumulated time reaches Tfinal.");
+   args.AddOption(&u0,     "-u0",  "--u0",          "Initial scalar state u0.");
+   args.AddOption(&target, "-ut",  "--target",
+                  "Target value in J=0.5*(u_m-target)^2.");
+   args.AddOption(&eps,    "-eps", "--fd-eps",      "Finite-difference epsilon.");
+   args.AddOption(&err,    "-err", "--time_err",
+                  "Allowed time integration error per time step.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(mfem::out);
+      return 1;
+   }
+   args.PrintOptions(mfem::out);
+
+   MFEM_VERIFY(s > 0, "Need s > 0.");
+   MFEM_VERIFY(dt0 > 0.0, "Need dt0 > 0.");
+   MFEM_VERIFY(Tfinal > 0.0, "Need Tfinal > 0.");
+
+   using Step = mfem::DynamicCheckpointing<double>::Step;
+
+   struct my_state
+   {
+      my_state(double a, double b, double c, double d=0.0, double tp_=0.0)
+      {
+         t=a; dt=b; u=c; up=d; tp=tp_;
+      }
+
+      //t,dt,u
+      double t;
+      double dt;
+      double u;
+      double tp;
+      double up;
+   };
+
+   // Checkpoint manager:
+   //   State   = double
+   //   Snapshot= double
+   mfem::DynamicCheckpointing<my_state> ckpt(s);
+
+   auto make_snapshot = [](const my_state &u) -> my_state { return u; };
+   auto restore_snapshot = [](const my_state &snap, my_state &state) { state = snap; };
+
+   double tmax=Tfinal;
+
+   auto primal_step = [&](my_state &su, Step i)
+   {
+      su.up=su.u;
+      su.tp=su.t;
+
+      double t=su.t;
+      double dt = su.dt;
+      bool flag=true;
+
+      if ((tmax-t)<su.dt)
+      {
+         dt=tmax-t;
+         flag=false;
+      }
+
+      double u=su.u;
+
+      double s0=alpha * u * (1.0 - u);
+      double ue=u+dt*s0;
+      double s1=alpha * ue * (1.0 - ue);
+      double uh=u + 0.5 * dt * (s0 + s1);
+      double ee = std::abs(uh-ue);
+
+      if ((ee < 0.5*err) && (flag))
+      {
+         dt=1.25*dt;
+         ue=u+dt*s0;
+         s1=alpha * ue * (1.0 - ue);
+         uh=u + 0.5 * dt * (s0 + s1);
+         ee = std::abs(uh-ue);
+      }
+
+      while (ee > err)
+      {
+         dt=0.5*dt;
+         ue=u+dt*s0;
+         s1=alpha * ue * (1.0 - ue);
+         uh=u + 0.5 * dt * (s0 + s1);
+         ee = std::abs(uh-ue);
+      }
+
+      //mfem::out<<" t="<<t+dt<<" dt="<<dt<<" err="<<ee<<std::endl;
+
+      su.t = t+dt;
+      su.dt = dt;
+      su.u = uh;
+   };
+
+   auto adjoint_step = [&](double &lambda, const my_state &u_i, Step i)
+   {
+      const double dt = u_i.dt;
+      // const double dF_du = 1.0 + dt * alpha * (1.0 - 2.0 * u_i.u);
+      // lambda = dF_du * lambda;
+
+      const double s0=alpha * (1.0 -2.0 * u_i.u) *lambda;
+      const double le= lambda+dt*s0;
+      const double s1= alpha * (1.0 -2.0 * u_i.up) * le;
+      lambda = lambda + 0.5 *dt *(s0+s1);
+   };
+
+   // ---------------- Forward sweep (unknown m) ----------------
+   my_state u(0.0, dt0, u0);
+   double t_phys = 0.0;
+
+   Step i = 0;
+   while (t_phys < Tfinal)
+   {
+      ckpt.ForwardStep(i, u, primal_step, make_snapshot);
+      t_phys = u.t;
+      ++i;
+   }
+
+   const Step m = i;
+   const double u_m = u.u;
+   const double J = 0.5 * (u_m - target) * (u_m - target);
+
+   mfem::out << std::setprecision(15);
+   mfem::out << "\n[Scalar] Forward finished:\n";
+   mfem::out << "  m (steps) = " << m << "\n";
+   mfem::out << "  t_phys    = " << t_phys << "\n";
+   mfem::out << "  u_m       = " << u_m << "\n";
+   mfem::out << "  J         = " << J << "\n\n";
+
+   mfem::out <<
+             "[Scalar] Checkpoint set after forward sweep (step, level, stored):\n";
+   for (const auto &cp : ckpt.GetCheckpointInfo())
+   {
+      mfem::out << "  step=" << cp.step
+                << ", level=" << cp.level
+                << ", stored=" << (cp.stored ? "yes" : "no")
+                << (cp.stored ? "" : " (placeholder)")
+                << "\n";
+   }
+   mfem::out << "\n";
+
+   MFEM_VERIFY(m > 0, "Forward produced m=0 steps; nothing to do.");
+
+   // ---------------- Backward sweep (adjoint) ----------------
+   double lambda = (u_m - target); // terminal adjoint = dJ/du_m
+   my_state u_work (0.0,0.0,0.0);            // scratch primal state u_i
+
+   for (Step j = m - 1; j >= 0; --j)
+   {
+      ckpt.BackwardStep(j, lambda, u_work,
+                        primal_step, adjoint_step,
+                        make_snapshot, restore_snapshot);
+      if (j == 0) { break; } // avoid signed underflow
+   }
+
+   const double dJ_du0_adjoint = lambda;
+
+   // analytic solution
+   const double sol=u0*exp(alpha*Tfinal)/(1.0-u0+u0*exp(alpha*Tfinal));
+   const double grd=exp(alpha*Tfinal)/std::pow(1.0-u0+u0*exp(alpha*Tfinal),2.0);
+
+   mfem::out << "true sol = "<<sol<<" \n";
+
+
+   const double dJ_du0_fd = (sol - target) * grd;
+
+   const double abs_err = std::abs(dJ_du0_adjoint - dJ_du0_fd);
+   const double rel_err = abs_err / (std::abs(dJ_du0_fd) + 1e-30);
+
+   mfem::out << "[Scalar] Gradient check (dJ/du0):\n";
+   mfem::out << "  adjoint = " << dJ_du0_adjoint << "\n";
+   mfem::out << "  true grad = " << dJ_du0_fd << "\n";
+   mfem::out << "  abs err = " << abs_err << "\n";
+   mfem::out << "  rel err = " << rel_err << "\n\n";
+
+   return 0;
+}

miniapps/mtop/chpt/ex_scalar_field.cpp

@@ -0,0 +1,168 @@
+#include "mfem.hpp"
+
+#include <cmath>
+#include <iomanip>
+
+using namespace mfem;
+
+// Scalar logistic ODE:
+//   du/dt = alpha*u*(1-u)
+// Explicit Euler with variable dt(i):
+//   u_{i+1} = u_i + dt_i*alpha*u_i*(1-u_i)
+// Objective:
+//   J = 0.5*(u_m - target)^2
+// Discrete adjoint (scalar):
+//   lambda_m = (u_m - target)
+//   lambda_i = (dF/du at u_i)^T * lambda_{i+1}
+//   dF/du = 1 + dt_i*alpha*(1 - 2*u_i)
+
+int main(int argc, char *argv[])
+{
+   int s = 3;
+   double alpha = 2.0;
+   double dt0 = 0.02;
+   double omega = 0.2;
+   double Tfinal = 1.0;
+   double u0 = 0.2;
+   double target = 0.7;
+   double eps = 1e-7;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&s,      "-s",   "--checkpoints",
+                  "Checkpoint budget s (real checkpoints).");
+   args.AddOption(&alpha,  "-a",   "--alpha",       "Logistic growth alpha.");
+   args.AddOption(&dt0,    "-dt0", "--dt0",
+                  "Base dt for dt(i)=dt0*(1+0.5*sin(omega*i)).");
+   args.AddOption(&omega,  "-om",  "--omega",
+                  "Omega in dt(i)=dt0*(1+0.5*sin(omega*i)).");
+   args.AddOption(&Tfinal, "-T",   "--tfinal",
+                  "Terminate when accumulated time reaches Tfinal.");
+   args.AddOption(&u0,     "-u0",  "--u0",          "Initial scalar state u0.");
+   args.AddOption(&target, "-ut",  "--target",
+                  "Target value in J=0.5*(u_m-target)^2.");
+   args.AddOption(&eps,    "-eps", "--fd-eps",      "Finite-difference epsilon.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(mfem::out);
+      return 1;
+   }
+   args.PrintOptions(mfem::out);
+
+   MFEM_VERIFY(s > 0, "Need s > 0.");
+   MFEM_VERIFY(dt0 > 0.0, "Need dt0 > 0.");
+   MFEM_VERIFY(Tfinal > 0.0, "Need Tfinal > 0.");
+
+   using Step = mfem::DynamicCheckpointing<double>::Step;
+
+   auto dt_func = [&](Step i)
+   {
+      const double dt = dt0 * (1.0 + 0.5 * std::sin(omega * double(i)));
+      MFEM_VERIFY(dt > 0.0, "dt_func produced non-positive dt.");
+      return dt;
+   };
+
+   // Checkpoint manager:
+   //   State   = double
+   //   Snapshot= double
+   mfem::DynamicCheckpointing<double> ckpt(s);
+
+   auto make_snapshot = [](const double &u) -> double { return u; };
+   auto restore_snapshot = [](const double &snap, double &state) { state = snap; };
+
+   auto primal_step = [&](double &u, Step i)
+   {
+      const double dt = dt_func(i);
+      u = u + dt * alpha * u * (1.0 - u);
+   };
+
+   auto adjoint_step = [&](double &lambda, const double &u_i, Step i)
+   {
+      const double dt = dt_func(i);
+      const double dF_du = 1.0 + dt * alpha * (1.0 - 2.0 * u_i);
+      lambda = dF_du * lambda;
+   };
+
+   // ---------------- Forward sweep (unknown m) ----------------
+   double u = u0;
+   double t_phys = 0.0;
+
+   Step i = 0;
+   while (t_phys < Tfinal)
+   {
+      ckpt.ForwardStep(i, u, primal_step, make_snapshot);
+      t_phys += dt_func(i);
+      ++i;
+   }
+
+   const Step m = i;
+   const double u_m = u;
+   const double J = 0.5 * (u_m - target) * (u_m - target);
+
+   mfem::out << std::setprecision(15);
+   mfem::out << "\n[Scalar] Forward finished:\n";
+   mfem::out << "  m (steps) = " << m << "\n";
+   mfem::out << "  t_phys    = " << t_phys << "\n";
+   mfem::out << "  u_m       = " << u_m << "\n";
+   mfem::out << "  J         = " << J << "\n\n";
+
+   mfem::out <<
+             "[Scalar] Checkpoint set after forward sweep (step, level, stored):\n";
+   for (const auto &cp : ckpt.GetCheckpointInfo())
+   {
+      mfem::out << "  step=" << cp.step
+                << ", level=" << cp.level
+                << ", stored=" << (cp.stored ? "yes" : "no")
+                << (cp.stored ? "" : " (placeholder)")
+                << "\n";
+   }
+   mfem::out << "\n";
+
+   MFEM_VERIFY(m > 0, "Forward produced m=0 steps; nothing to do.");
+
+   // ---------------- Backward sweep (adjoint) ----------------
+   double lambda = (u_m - target); // terminal adjoint = dJ/du_m
+   double u_work = 0.0;            // scratch primal state u_i
+
+   for (Step j = m - 1; j >= 0; --j)
+   {
+      ckpt.BackwardStep(j, lambda, u_work,
+                        primal_step, adjoint_step,
+                        make_snapshot, restore_snapshot);
+      if (j == 0) { break; } // avoid signed underflow
+   }
+
+   const double dJ_du0_adjoint = lambda;
+
+   // ---------------- Finite-difference gradient check ----------------
+   auto forward_only_J = [&](double u_init)
+   {
+      double uu = u_init;
+      double tt = 0.0;
+      Step k = 0;
+      while (tt < Tfinal)
+      {
+         primal_step(uu, k);
+         tt += dt_func(k);
+         ++k;
+      }
+      const double r = (uu - target);
+      return 0.5 * r * r;
+   };
+
+   const double Jp = forward_only_J(u0 + eps);
+   const double Jm = forward_only_J(u0 - eps);
+   const double dJ_du0_fd = (Jp - Jm) / (2.0 * eps);
+
+   const double abs_err = std::abs(dJ_du0_adjoint - dJ_du0_fd);
+   const double rel_err = abs_err / (std::abs(dJ_du0_fd) + 1e-30);
+
+   mfem::out << "[Scalar] Gradient check (dJ/du0):\n";
+   mfem::out << "  adjoint = " << dJ_du0_adjoint << "\n";
+   mfem::out << "  FD      = " << dJ_du0_fd << "\n";
+   mfem::out << "  abs err = " << abs_err << "\n";
+   mfem::out << "  rel err = " << rel_err << "\n\n";
+
+   return 0;
+}
+

miniapps/mtop/chpt/ex_vector_field.cpp

@@ -0,0 +1,199 @@
+#include "mfem.hpp"
+
+#include <cmath>
+#include <iomanip>
+
+using namespace mfem;
+
+// Vector logistic ODE (componentwise):
+//   du/dt = alpha*u*(1-u)   (applied to each component)
+// Explicit Euler:
+//   u_{i+1}[j] = u_i[j] + dt_i*alpha*u_i[j]*(1-u_i[j])
+// Objective:
+//   J = 0.5*||u_m - u_target||^2
+// Terminal adjoint:
+//   lambda_m = u_m - u_target
+// Adjoint step (componentwise):
+//   lambda_i[j] = (1 + dt_i*alpha*(1 - 2*u_i[j])) * lambda_{i+1}[j]
+
+int main(int argc, char *argv[])
+{
+   int s = 3;
+   int n = 16;
+   double alpha = 2.0;
+   double dt0 = 0.02;
+   double omega = 0.2;
+   double Tfinal = 1.0;
+   double target_val = 0.7;
+   double eps = 1e-7;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&s,         "-s",   "--checkpoints",
+                  "Checkpoint budget s (real checkpoints).");
+   args.AddOption(&n,         "-n",   "--size",        "Vector dimension n.");
+   args.AddOption(&alpha,     "-a",   "--alpha",       "Logistic growth alpha.");
+   args.AddOption(&dt0,       "-dt0", "--dt0",
+                  "Base dt for dt(i)=dt0*(1+0.5*sin(omega*i)).");
+   args.AddOption(&omega,     "-om",  "--omega",
+                  "Omega in dt(i)=dt0*(1+0.5*sin(omega*i)).");
+   args.AddOption(&Tfinal,    "-T",   "--tfinal",
+                  "Terminate when accumulated time reaches Tfinal.");
+   args.AddOption(&target_val,"-tv",  "--target",
+                  "Target value for each component.");
+   args.AddOption(&eps,       "-eps", "--fd-eps",
+                  "Finite-difference epsilon (directional).");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(mfem::out);
+      return 1;
+   }
+   args.PrintOptions(mfem::out);
+
+   MFEM_VERIFY(s > 0, "Need s > 0.");
+   MFEM_VERIFY(n > 0, "Need n > 0.");
+   MFEM_VERIFY(dt0 > 0.0, "Need dt0 > 0.");
+   MFEM_VERIFY(Tfinal > 0.0, "Need Tfinal > 0.");
+
+   using Step = mfem::DynamicCheckpointing<mfem::Vector>::Step;
+
+   auto dt_func = [&](Step i)
+   {
+      const double dt = dt0 * (1.0 + 0.5 * std::sin(omega * double(i)));
+      MFEM_VERIFY(dt > 0.0, "dt_func produced non-positive dt.");
+      return dt;
+   };
+
+   // Checkpoint manager:
+   //   State   = mfem::Vector
+   //   Snapshot= mfem::Vector
+   mfem::DynamicCheckpointing<mfem::Vector> ckpt(s);
+
+   auto make_snapshot = [](const mfem::Vector &u) -> mfem::Vector { return u; };
+   auto restore_snapshot = [](const mfem::Vector &snap, mfem::Vector &state) { state = snap; };
+
+   auto primal_step = [&](mfem::Vector &u, Step i)
+   {
+      const double dt = dt_func(i);
+      for (int j = 0; j < u.Size(); ++j)
+      {
+         const double uj = u[j];
+         u[j] = uj + dt * alpha * uj * (1.0 - uj);
+      }
+   };
+
+   auto adjoint_step = [&](mfem::Vector &lambda, const mfem::Vector &u_i, Step i)
+   {
+      const double dt = dt_func(i);
+      MFEM_ASSERT(lambda.Size() == u_i.Size(), "lambda and u_i size mismatch");
+      for (int j = 0; j < lambda.Size(); ++j)
+      {
+         const double dF_du = 1.0 + dt * alpha * (1.0 - 2.0 * u_i[j]);
+         lambda[j] *= dF_du;
+      }
+   };
+
+   // Initial condition and target
+   mfem::Vector u0(n), u_target(n);
+   for (int j = 0; j < n; ++j)
+   {
+      u0[j] = 0.2 + 0.05 * std::cos(0.7 * (j + 1));
+   }
+   u_target = target_val;
+
+   // ---------------- Forward sweep (unknown m) ----------------
+   mfem::Vector u = u0;
+   double t_phys = 0.0;
+
+   Step i = 0;
+   while (t_phys < Tfinal)
+   {
+      ckpt.ForwardStep(i, u, primal_step, make_snapshot);
+      t_phys += dt_func(i);
+      ++i;
+   }
+
+   const Step m = i;
+   const mfem::Vector u_m = u;
+
+   mfem::Vector diff(u_m);
+   diff -= u_target;
+   const double J = 0.5 * mfem::InnerProduct(diff, diff);
+
+   mfem::out << std::setprecision(15);
+   mfem::out << "\n[Vector] Forward finished:\n";
+   mfem::out << "  m (steps) = " << m << "\n";
+   mfem::out << "  t_phys    = " << t_phys << "\n";
+   mfem::out << "  J         = " << J << "\n";
+   mfem::out << "  ||u_m||   = " << u_m.Norml2() << "\n";
+   mfem::out << "  ||u_m-ut||= " << diff.Norml2() << "\n\n";
+
+   mfem::out <<
+             "[Vector] Checkpoint set after forward sweep (step, level, stored):\n";
+   for (const auto &cp : ckpt.GetCheckpointInfo())
+   {
+      mfem::out << "  step=" << cp.step
+                << ", level=" << cp.level
+                << ", stored=" << (cp.stored ? "yes" : "no")
+                << (cp.stored ? "" : " (placeholder)")
+                << "\n";
+   }
+   mfem::out << "\n";
+
+   MFEM_VERIFY(m > 0, "Forward produced m=0 steps; nothing to do.");
+
+   // ---------------- Backward sweep (adjoint) ----------------
+   mfem::Vector lambda = diff; // terminal = dJ/du_m
+   mfem::Vector u_work(n);     // scratch primal u_i
+
+   for (Step j = m - 1; j >= 0; --j)
+   {
+      ckpt.BackwardStep(j, lambda, u_work,
+                        primal_step, adjoint_step,
+                        make_snapshot, restore_snapshot);
+      if (j == 0) { break; }
+   }
+   const mfem::Vector &grad_u0 = lambda;
+
+   // ---------------- Directional FD check ----------------
+   mfem::Vector v(n);
+   for (int j = 0; j < n; ++j) { v[j] = std::sin(0.3 * (j + 1)) + 0.1; }
+
+   auto forward_only_J = [&](const mfem::Vector &u_init)
+   {
+      mfem::Vector uu = u_init;
+      double tt = 0.0;
+      Step k = 0;
+      while (tt < Tfinal)
+      {
+         primal_step(uu, k);
+         tt += dt_func(k);
+         ++k;
+      }
+      mfem::Vector dd(uu);
+      dd -= u_target;
+      return 0.5 * mfem::InnerProduct(dd, dd);
+   };
+
+   mfem::Vector u_plus(u0), u_minus(u0);
+   u_plus.Add(eps, v);
+   u_minus.Add(-eps, v);
+
+   const double Jp = forward_only_J(u_plus);
+   const double Jm = forward_only_J(u_minus);
+   const double dJ_dir_fd = (Jp - Jm) / (2.0 * eps);
+
+   const double dJ_dir_adj = mfem::InnerProduct(grad_u0, v);
+
+   const double abs_err = std::abs(dJ_dir_adj - dJ_dir_fd);
+   const double rel_err = abs_err / (std::abs(dJ_dir_fd) + 1e-30);
+
+   mfem::out << "[Vector] Directional derivative check:\n";
+   mfem::out << "  v·grad adjoint = " << dJ_dir_adj << "\n";
+   mfem::out << "  FD directional = " << dJ_dir_fd << "\n";
+   mfem::out << "  abs err        = " << abs_err << "\n";
+   mfem::out << "  rel err        = " << rel_err << "\n\n";
+
+   return 0;
+}
+

miniapps/mtop/chpt/makefile

@@ -0,0 +1,92 @@
+# 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.
+
+# Use the MFEM build directory
+MFEM_DIR ?= ../../..
+MFEM_BUILD_DIR ?= ../../..
+MFEM_INSTALL_DIR ?= ../../../mfem
+SRC = $(if $(MFEM_DIR:../../..=),$(MFEM_DIR)/miniapps/mtop/chpt/,)
+CONFIG_MK = $(or $(wildcard $(MFEM_BUILD_DIR)/config/config.mk),\
+   $(wildcard $(MFEM_INSTALL_DIR)/share/mfem/config.mk))
+
+MFEM_LIB_FILE = mfem_is_not_built
+-include $(CONFIG_MK)
+
+SEQ_MINIAPPS = mtop_ex_scalar_chpt mtop_ex_adapt_chpt mtop_ex_vector_chpt \
+   mtop_ex_fixedslot mtop_ex_file mtop_ex_revolve mtop_ex_dyn \
+   mtop_ex_rev_fixedslot
+PAR_MINIAPPS =
+MINIAPPS_CLEAN := $(SEQ_MINIAPPS) $(PAR_MINIAPPS)
+
+ifeq ($(MFEM_USE_MPI),NO)
+   MINIAPPS = $(SEQ_MINIAPPS)
+else
+   MINIAPPS = $(SEQ_MINIAPPS) $(PAR_MINIAPPS)
+endif
+
+.SUFFIXES:
+.SUFFIXES: .o .cpp .mk
+.PHONY: all clean clean-build clean-exec
+.PRECIOUS: %.o
+
+# Remove built-in rules
+%: %.cpp
+%.o: %.cpp
+
+all: $(MINIAPPS)
+
+# Rules for building the miniapps
+mtop_ex_scalar_chpt: $(SRC)ex_scalar_field.cpp $(MFEM_LIB_FILE) $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) $< -o $@ $(MFEM_LIBS)
+
+mtop_ex_adapt_chpt: $(SRC)ex_scalar_adaptive.cpp $(MFEM_LIB_FILE) $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) $< -o $@ $(MFEM_LIBS)
+
+mtop_ex_vector_chpt: $(SRC)ex_vector_field.cpp $(MFEM_LIB_FILE) $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) $< -o $@ $(MFEM_LIBS)
+
+mtop_ex_fixedslot: $(SRC)ex_fixedslot_scalar_vector.cpp $(MFEM_LIB_FILE) $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) $< -o $@ $(MFEM_LIBS)
+
+mtop_ex_file: $(SRC)ex_fileio_scalar_vector_mfem.cpp $(MFEM_LIB_FILE) $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) $< -o $@ $(MFEM_LIBS)
+
+mtop_ex_revolve: $(SRC)ex_revolve.cpp $(MFEM_LIB_FILE) $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) $< -o $@ $(MFEM_LIBS)
+
+mtop_ex_dyn: $(SRC)ex_dynamic_chpt.cpp $(MFEM_LIB_FILE) $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) $< -o $@ $(MFEM_LIBS)
+
+mtop_ex_rev_fixedslot: $(SRC)ex_revolve_chpt_fixedslot.cpp $(MFEM_LIB_FILE) $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) $< -o $@ $(MFEM_LIBS)
+
+MFEM_TESTS = MINIAPPS
+include $(MFEM_TEST_MK)
+
+TEST_NAME := MTop checkpointing miniapp
+%-test-seq: %
+	@$(call mfem-test,$<,, $(TEST_NAME))
+
+# Testing: "test" target and mfem-test* variables are defined in config/test.mk
+
+# Generate an error message if the MFEM library is not built and exit
+$(MFEM_LIB_FILE):
+	$(error The MFEM library is not built)
+
+clean: clean-build clean-exec
+
+clean-build:
+	rm -f *.o *~ $(MINIAPPS_CLEAN)
+	rm -rf *.dSYM *.TVD.*breakpoints
+
+clean-exec:
+	rm -rf chk_scalar_files chk_vector_files dyn_ckpts
+	rm -f dynamic_ckpts.bin revolve_ckpts.bin