Added support for MultiVector and cleaned up.

in some cases. This required changes in the internals of the class.

Added some new methods in class MultiVector.

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

@@ -142,10 +142,6 @@ 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
@@ -294,52 +290,6 @@ 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
@@ -355,14 +305,9 @@ jobs:
           metis-dir: ${{ env.METIS_TOP_DIR }}
           mfem-dir: ${{ env.MFEM_TOP_DIR }}
           precision: ${{ matrix.precision }}
-          config-options: ${{ matrix.config-opts }} ${{ env.CCACHE_CONFIG_OPTS }}
+          config-options: ${{ matrix.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

@@ -1,42 +0,0 @@
-# Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
-# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-# LICENSE and NOTICE for details. LLNL-CODE-806117.
-#
-# This file is part of the MFEM library. For more information and source code
-# availability visit https://mfem.org.
-#
-# MFEM is free software; you can redistribute it and/or modify it under the
-# terms of the BSD-3 license. We welcome feedback and contributions, see file
-# CONTRIBUTING.md for details.
----
-# 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,7 +13,6 @@ name: "Checks"
 
 permissions:
   actions: write
-  pull-requests: read
 
 on:
   push:
@@ -30,11 +29,6 @@ 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:
@@ -134,7 +128,10 @@ jobs:
 
   branch-history:
     if: |
-      github.ref != 'refs/heads/next' && github.ref != 'refs/heads/master'
+      github.ref != 'refs/heads/next' &&
+      github.ref != 'refs/heads/master' &&
+      (github.event_name == 'push' ||
+       github.event.pull_request.head.repo.full_name != github.repository)
     runs-on: ubuntu-latest
     steps:
       - name: checkout mfem
@@ -142,27 +139,7 @@ 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.
@@ -170,9 +147,3 @@ 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,8 +39,3 @@ 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.

fem/datacollection.hpp

@@ -25,21 +25,35 @@
 namespace mfem
 {
 
-/// Lightweight adaptor over an std::map from strings to pointer to T
-template<typename T>
-class NamedFieldsMap
+/// Lightweight adaptor over an std::map from type K to type to V
+template<typename K, typename V,
+         typename = typename std::enable_if<std::is_default_constructible<V>::value>::type>
+class GenericFieldMap
 {
+private:
+   static constexpr bool ValueIsPointer = std::is_pointer<V>::value;
+
 public:
-   typedef std::map<std::string, T*> MapType;
+   typedef std::map<K, V> MapType;
    typedef typename MapType::iterator iterator;
    typedef typename MapType::const_iterator const_iterator;
 
-   /// Register field @a field with name @a fname
-   /** Replace existing field associated with @a fname (and optionally
-       delete associated pointer if @a own_data is true) */
-   void Register(const std::string& fname, T* field, bool own_data)
+   /// Register field @a field with name @a key
+   /// Only enabled if the template parameter V is not a pointer
+   template<typename = std::enable_if<!ValueIsPointer, bool>>
+   void Register(const K& key, V field)
    {
-      T*& ref = field_map[fname];
+      field_map[key] = field;
+   }
+
+   /// Register field @a field with name @a key
+   /** Replace existing field associated with @a key (and optionally
+       delete associated pointer if @a own_data is true).
+       Only enabled if the template parameter V is a pointer*/
+   template<typename = std::enable_if<ValueIsPointer, bool>>
+   void Register(const K& key, V field, bool own_data)
+   {
+      V& ref = field_map[key];
       if (own_data)
       {
          delete ref; // if newly allocated -> ref is null -> OK
@@ -47,23 +61,40 @@ public:
       ref = field;
    }
 
-   /// Unregister association between field @a field and name @a fname
-   /** Optionally delete associated pointer if @a own_data is true */
-   void Deregister(const std::string& fname, bool own_data)
+   /// Unregister association between field @a field and name @a key
+   /// Only enabled if the template parameter V is not a pointer
+   template<typename = std::enable_if<!ValueIsPointer, bool>>
+   void Deregister(const K& key)
    {
-      iterator it = field_map.find(fname);
+      iterator it = field_map.find(key);
+      if ( it != field_map.end() )
+      {
+         field_map.erase(it);
+      }
+   }
+
+   /// Unregister association between field @a field and name @a key
+   /** Optionally delete associated pointer if @a own_data is true.
+       Only enabled if the template parameter V is a pointer */
+   template<typename = std::enable_if<ValueIsPointer, bool>>
+   void Deregister(const K& key, bool own_data)
+   {
+      iterator it = field_map.find(key);
       if ( it != field_map.end() )
       {
          if (own_data)
          {
             delete it->second;
+            it->second = nullptr;
          }
          field_map.erase(it);
       }
    }
 
    /// Clear all associations between names and fields
-   /** Delete associated pointers when @a own_data is true */
+   /** Delete associated pointers when @a own_data is true.
+       Only enabled if the template parameter V is a pointer */
+   template<typename = std::enable_if<ValueIsPointer, bool>>
    void DeleteData(bool own_data)
    {
       for (iterator it = field_map.begin(); it != field_map.end(); ++it)
@@ -76,22 +107,37 @@ public:
       }
    }
 
-   /// Predicate to check if a field is associated with name @a fname
-   bool Has(const std::string& fname) const
+   /// Predicate to check if a field is associated with name @a key
+   bool Has(const K& key) const
    {
-      return field_map.find(fname) != field_map.end();
+      return field_map.find(key) != field_map.end();
    }
 
-   /// Get a pointer to the field associated with name @a fname
-   /** @return Pointer to field associated with @a fname or NULL */
-   T* Get(const std::string& fname) const
+   /// Get a pointer to the field associated with name @a key
+   /** @return Field associated with @a key or NULL,
+       if value is pointer and key not found */
+   V Get(const K& key) const
    {
-      const_iterator it = field_map.find(fname);
-      return it != field_map.end() ? it->second : NULL;
+      const_iterator it = field_map.find(key);
+      if (it != field_map.end())
+      {
+         return it->second;
+      }
+      else
+      {
+         if constexpr (ValueIsPointer)
+         {
+            return nullptr;
+         }
+         else
+         {
+            return V(); // Return default-constructed value for non-pointer types
+         }
+      }
    }
 
    /// Returns a const reference to the underlying map
-   const MapType& GetMap() const { return field_map; }
+   const MapType &GetMap() const { return field_map; }
 
    /// Returns the number of registered fields
    int NumFields() const { return field_map.size(); }
@@ -106,21 +152,24 @@ public:
    /// Returns an end const iterator to the registered fields
    const_iterator end() const { return field_map.end(); }
 
-   /// Returns an iterator to the field @a fname
-   iterator find(const std::string& fname)
-   { return field_map.find(fname); }
+   /// Returns an iterator to the field @a key
+   iterator find(const K& key)
+   { return field_map.find(key); }
 
-   /// Returns a const iterator to the field @a fname
-   const_iterator find(const std::string& fname) const
-   { return field_map.find(fname); }
+   /// Returns a const iterator to the field @a key
+   const_iterator find(const K& key) const
+   { return field_map.find(key); }
 
-   /// Clears the map of registered fields without reclaiming memory
+   /// Clears the map of registered fields
    void clear() { field_map.clear(); }
 
 protected:
    MapType field_map;
 };
 
+/// Lightweight adaptor over an std::map from strings to pointer to T
+template<typename T>
+using NamedFieldsMap = GenericFieldMap<std::string, T*>;
 
 /** A class for collecting finite element data that is part of the same
     simulation. Currently, this class groups together grid functions (fields),

fem/dfem/util.hpp

@@ -671,6 +671,20 @@ public:
                     MPI_COMM_WORLD);
    }
 
+   Operator& GetGradient(const Vector &x0) const override
+   {
+      x = x0;
+      f.UseDevice(x.UseDevice());
+      xpev.UseDevice(x.UseDevice());
+
+      op.Mult(x, f);
+      const real_t xnorm_local = x.Norml2();
+      MPI_Allreduce(&xnorm_local, &xnorm, 1, MPITypeMap<real_t>::mpi_type, MPI_SUM,
+                    MPI_COMM_WORLD);
+
+      return const_cast<FDJacobian&>(*this);
+   }
+
    void Mult(const Vector &v, Vector &y) const override
    {
       // See [1] for choice of eps.
@@ -725,11 +739,11 @@ public:
 
 private:
    const Operator &op;
-   Vector x, f;
+   mutable Vector x, f;
    mutable Vector xpev;
    real_t lambda = 1.0e-6;
    real_t fixed_eps;
-   real_t xnorm;
+   mutable real_t xnorm;
 };
 
 /// @brief Find the index of a field descriptor in a vector of field descriptors.

fem/nonlinearform.hpp

@@ -224,6 +224,9 @@ public:
    /** @see GetGradient(const Vector &) */
    Operator &GetGradient(const Vector &x, bool finalize) const;
 
+   /// Suppress a warning about hiding overloaded virtual function.
+   using Operator::GetGradient;
+
    /// Update the NonlinearForm to propagate updates of the associated FE space.
    /** After calling this method, the essential boundary conditions need to be
        set again. */

general/optparser.cpp

@@ -12,9 +12,7 @@
 #include "optparser.hpp"
 #include "../linalg/vector.hpp"
 #include "../general/communication.hpp"
-#include <cerrno>
 #include <cctype>
-#include <cstdlib>
 
 namespace mfem
 {
@@ -50,19 +48,6 @@ 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
@@ -215,10 +200,6 @@ 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++]);
@@ -298,10 +279,6 @@ 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;
@@ -433,7 +410,7 @@ void OptionsParser::PrintHelp(ostream &os) const
    static const char *line_sep = "";
    static const char *types[] = { " <int>", " <double>", " <string>",
                                   " <string>", "", "", " '<int>...'",
-                                  " '<double>...'", " <long long>"
+                                  " '<double>...'"
                                 };
 
    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, LONG };
+   enum OptionType { INT, DOUBLE, STRING, STD_STRING, ENABLE, DISABLE, ARRAY, VECTOR };
 
 private:
    struct Option
@@ -98,14 +98,6 @@ 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

@@ -27,6 +27,7 @@ list(APPEND SRCS
   handle.cpp
   matrix.cpp
   mma.cpp
+  multivector.cpp
   ode.cpp
   operator.cpp
   ordering.cpp
@@ -63,6 +64,7 @@ list(APPEND HDRS
   linalg.hpp
   matrix.hpp
   mma.hpp
+  multivector.hpp
   ode.hpp
   operator.hpp
   ordering.hpp
@@ -76,12 +78,6 @@ 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

@@ -1,21 +0,0 @@
-// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
-// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-// LICENSE and NOTICE for details. LLNL-CODE-806117.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability visit https://mfem.org.
-//
-// MFEM is free software; you can redistribute it and/or modify it under the
-// terms of the BSD-3 license. We welcome feedback and contributions, see file
-// CONTRIBUTING.md for details.
-
-#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

@@ -1,309 +0,0 @@
-#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

@@ -1,304 +0,0 @@
-#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

@@ -1,408 +0,0 @@
-#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

@@ -1,537 +0,0 @@
-#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

@@ -1,362 +0,0 @@
-#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

@@ -15,6 +15,7 @@
 // Linear algebra header file
 
 #include "vector.hpp"
+#include "multivector.hpp"
 #include "operator.hpp"
 #include "matrix.hpp"
 #include "sparsemat.hpp"
@@ -40,7 +41,6 @@
 #include "filteredsolver.hpp"
 #include "ordering.hpp"
 #include "particlevector.hpp"
-#include "chpt/chpt.hpp"
 
 #ifdef MFEM_USE_AMGX
 #include "amgxsolver.hpp"

linalg/multivector.cpp

@@ -0,0 +1,60 @@
+// 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.
+
+#include "multivector.hpp"
+
+namespace mfem
+{
+
+MultiVector::MultiVector(const Array<int> &vector_sizes)
+{
+   SetSizes(vector_sizes);
+}
+
+MultiVector::MultiVector(const Array<int> &vector_sizes, MemoryType mt)
+{
+   SetSizes(vector_sizes, mt);
+}
+
+MultiVector::MultiVector(Vector &base, const Array<int> &vector_sizes)
+{
+   MakeRef(base, vector_sizes);
+}
+
+void MultiVector::SetSizes(const Array<int> &vector_sizes)
+{
+   blocks.resize(vector_sizes.Size());
+   for (int i = 0; i < vector_sizes.Size(); i++)
+   {
+      operator[](i).SetSize(vector_sizes[i]);
+   }
+}
+
+void MultiVector::SetSizes(const Array<int> &vector_sizes, MemoryType mt)
+{
+   blocks.resize(vector_sizes.Size());
+   for (int i = 0; i < vector_sizes.Size(); i++)
+   {
+      operator[](i).SetSize(vector_sizes[i], mt);
+   }
+}
+
+void MultiVector::MakeRef(Vector &base, const Array<int> &vector_sizes)
+{
+   blocks.resize(vector_sizes.Size());
+   for (int offset = 0, i = 0; i < vector_sizes.Size(); i++)
+   {
+      blocks[i].emplace<0>(base, offset, vector_sizes[i]);
+      offset += vector_sizes[i];
+   }
+}
+
+} // namespace mfem

linalg/multivector.hpp

@@ -0,0 +1,198 @@
+// 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_MULTIVECTOR_HPP
+#define MFEM_MULTIVECTOR_HPP
+
+#include "../general/array.hpp"
+#include "vector.hpp"
+#include <vector>
+#include <array>
+#include <variant>
+
+namespace mfem
+{
+
+/// Class representing an array of Vectors with generally different sizes.
+/** This class is similar to BlockVector with the following two main
+    differences:
+    - the data for the individual Vector blocks does not need to be part of one
+      big contiguous memory allocation;
+    - this class does not inherit from class Vector (as a consequence of the
+      first bullet).
+
+    Internally, each Vector block is represented as either:
+    - (default) a Vector object constructed and owned by this class; this
+      object, in turn, as any Vector object, can own its Memory allocation or
+      refer to a sub-Memory of another Memory object; or
+    - a pointer to an externally allocated Vector or classes derived from
+      Vector. */
+class MultiVector
+{
+private:
+   std::vector<std::variant<Vector,Vector*>> blocks;
+
+public:
+   /// Create an empty MultiVector with zero blocks.
+   MultiVector() = default;
+
+   /** @brief Create a MultiVector with @a num_blocks blocks. The individual
+       Vector blocks are default initialized, i.e. they all have size zero. */
+   MultiVector(int num_blocks)
+      : blocks(num_blocks) { }
+
+   /** @brief Construct a MultiVector with number of blocks and individual block
+       Vector sizes given by @a vector_sizes.
+
+       @note The memory of the individual Vector blocks is NOT initialized. */
+   MultiVector(const Array<int> &vector_sizes);
+
+   /** @brief Construct a MultiVector with number of blocks and individual block
+       Vector sizes given by @a vector_sizes. All Vector blocks use the
+       MemoryType @a mt.
+
+       @note The memory of the individual Vector blocks is NOT initialized. */
+   MultiVector(const Array<int> &vector_sizes, MemoryType mt);
+
+   /** @brief Construct a MultiVector referencing data within a given monolithic
+       Vector @a base.
+
+       With this constructor, the Memory flags of @a base and of the individual
+       Vector blocks may need to be explicitly synchronized when data is moved
+       between host and device. */
+   MultiVector(Vector &base, const Array<int> &vector_sizes);
+
+   /** @brief Construct a MultiVector referencing multiple Vectors given as
+       arguments.
+
+       The VectorTypes reference arguments are expected to be static_cast-able
+       to (Vector &) which is the case if the types are derived from Vector,
+       e.g. HypreParVector, GridFunction, etc.
+
+       With this constructor, operations on individual Vector blocks are
+       performed directly on the objects @a vs. In particular, there is no need
+       to synchronize the Memory flags of @a vs and the ones of the individual
+       Vector blocks when data is moved between host and device. */
+   template <typename... VectorTypes,
+             std::enable_if_t<
+                std::conjunction_v<
+                   std::is_convertible<VectorTypes&,Vector&>...>, bool> = true>
+   MultiVector(VectorTypes &...vs) { MakeRef(vs...); }
+
+   /// Return the number of Vectors in the MultiVector.
+   int NumBlocks() const { return blocks.size(); }
+
+   /** @brief Set the number of Vectors in the MultiVector. Existing Vector
+       blocks will remain unmodified. New Vector blocks will be default
+       initialized, i.e. they all have size zero. */
+   void SetNumBlocks(int num_blocks) { blocks.resize(num_blocks); }
+
+   /// Read-write access to the i-th Vector.
+   inline Vector &operator[](int i);
+
+   /// Read-only access to the i-th Vector.
+   inline const Vector &operator[](int i) const;
+
+   /** @brief Update the MultiVector according to the given @a vector_sizes.
+
+       This method can be used to add or remove blocks. The individual Vector
+       sizes are updated using the method Vector::SetSize(int). */
+   void SetSizes(const Array<int> &vector_sizes);
+
+   /** @brief Update the MultiVector according to the given @a vector_sizes and
+       MemoryType @a mt.
+
+       This method can be used to add or remove blocks. The individual Vector
+       sizes and MemoryType are updated using the method
+       Vector::SetSize(int, MemoryType). */
+   void SetSizes(const Array<int> &vector_sizes, MemoryType mt);
+
+   /** @brief Update the MultiVector to reference data within a given monolithic
+       Vector @a base.
+
+       After calling this method, the Memory flags of @a base and of the
+       individual Vector blocks may need to be explicitly synchronized when data
+       is moved between host and device.*/
+   void MakeRef(Vector &base, const Array<int> &vector_sizes);
+
+   /** @brief Update the @a i-th MultiVector block to reference data within the
+       given monolithic Vector @a base at the given @a offset and with the given
+       @a size.
+
+       After calling this method, the Memory flags of @a base and of the @a i-th
+       Vector block may need to be explicitly synchronized when data is moved
+       between host and device.*/
+   inline void MakeRef(int i, Vector &base, int offset, int size)
+   {
+      blocks[i].emplace<0>(base, offset, size);
+   }
+
+   /** @brief Update the MultiVector to reference multiple Vectors given as
+       arguments.
+
+       The VectorTypes reference arguments are expected to be static_cast-able
+       to (Vector &) which is the case if the types are derived from Vector,
+       e.g. HypreParVector, GridFunction, etc.
+
+       After calling this method, operations on individual Vector blocks are
+       performed directly on the objects @a vs. In particular, there is no need
+       to synchronize the Memory flags of @a vs and the ones of the individual
+       Vector blocks when data is moved between host and device. */
+   template <typename... VectorTypes,
+             std::enable_if_t<
+                std::conjunction_v<
+                   std::is_convertible<VectorTypes&,Vector&>...>, bool> = true>
+   inline void MakeRef(VectorTypes &...vs);
+
+   /** @brief Update the @a i-th MultiVector block to reference the given
+       Vector @a v.
+
+       After calling this method, operations on the @a i-th Vector block are
+       performed directly on the Vector @a v. In particular, there is no need
+       to synchronize the Memory flags of @a v and the ones of the @a i-th
+       Vector blocks when data is moved between host and device. */
+   inline void MakeRef(int i, Vector &v) { blocks[i] = &v; }
+};
+
+// Inline and template methods
+
+inline Vector &MultiVector::operator[](int i)
+{
+   auto &bi = blocks[i];
+   return (bi.index() == 0) ? std::get<0>(bi) : *std::get<1>(bi);
+}
+
+inline const Vector &MultiVector::operator[](int i) const
+{
+   auto &bi = blocks[i];
+   return (bi.index() == 0) ? std::get<0>(bi) : *std::get<1>(bi);
+}
+
+template <typename... VectorTypes,
+          std::enable_if_t<
+             std::conjunction_v<
+                std::is_convertible<VectorTypes&,Vector&>...>, bool>>
+inline void MultiVector::MakeRef(VectorTypes &...vs)
+{
+   blocks.resize(sizeof...(vs));
+   if constexpr (sizeof...(vs) > 0)
+   {
+      const std::array vs_p{&static_cast<Vector&>(vs)...};
+      for (std::size_t i = 0; i < sizeof...(vs); i++)
+      {
+         blocks[i] = vs_p[i];
+      }
+   }
+}
+
+} // namespace mfem
+
+#endif // MFEM_MULTIVECTOR_HPP

linalg/operator.cpp

@@ -111,6 +111,16 @@ void Operator::ArrayAddMultTranspose(const Array<const Vector *> &X,
    }
 }
 
+void Operator::Mult(const MultiVector &, MultiVector &)
+{
+   MFEM_ABORT("this method is not overriden for this class!");
+}
+
+Operator &Operator::GetGradient(const MultiVector &) const
+{
+   MFEM_ABORT("this method is not overriden for this class!");
+}
+
 void Operator::FormLinearSystem(const Array<int> &ess_tdof_list,
                                 Vector &x, Vector &b,
                                 Operator* &Aout, Vector &X, Vector &B,

linalg/operator.hpp

@@ -13,6 +13,7 @@
 #define MFEM_OPERATOR
 
 #include "vector.hpp"
+#include "multivector.hpp"
 
 namespace mfem
 {
@@ -129,6 +130,16 @@ public:
    virtual void ArrayAddMultTranspose(const Array<const Vector *> &X,
                                       Array<Vector *> &Y, const real_t a = 1.0) const;
 
+   /** @brief Operator application, y = A(x), where the input @a x and the
+       output @a y are MultiVector objects, i.e. they generally use
+       non-contiguous memory representation.
+
+       The typical use case for this method are block operators like
+       DifferentiableOperator.
+
+       The base class implementation for the method is to generate an error. */
+   virtual void Mult(const MultiVector &x, MultiVector &y);
+
    /** @brief Evaluate the gradient operator at the point @a x. The default
        behavior in class Operator is to generate an error. */
    virtual Operator &GetGradient(const Vector &x) const
@@ -137,6 +148,16 @@ public:
       return const_cast<Operator &>(*this);
    }
 
+   /** @brief Evaluate the gradient operator at the point @a x. The input @a x
+       is provided as a MultiVector, i.e. it generally uses non-contiguous
+       memory representation.
+
+       The typical use case for this method are block operators like
+       DifferentiableOperator.
+
+       The base class implementation for the method is to generate an error. */
+   virtual Operator &GetGradient(const MultiVector &x) const;
+
    /** @brief Computes the diagonal entries into @a diag. Typically, this
        operation only makes sense for linear Operator%s. In some cases, only an
        approximation of the diagonal is computed. */

makefile

@@ -811,7 +811,6 @@ 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

@@ -22,6 +22,7 @@ add_subdirectory(common)
 add_subdirectory(contact)
 add_subdirectory(dfem)
 add_subdirectory(diag-smoothers)
+add_subdirectory(multiapp)
 add_subdirectory(dpg)
 add_subdirectory(electromagnetics)
 add_subdirectory(fluids/navier)
@@ -31,7 +32,6 @@ 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/hdiv-linear-solver/change_basis.hpp

@@ -80,6 +80,8 @@ public:
    // limitations
    void MultRT_2D(const Vector &x, Vector &y, Mode mode) const;
    void MultRT_3D(const Vector &x, Vector &y, Mode mode) const;
+   // suppress warning about hiding overloaded virtual function:
+   using Operator::Mult;
 };
 
 } // namespace mfem

miniapps/mtop/chpt/CMakeLists.txt

@@ -1,68 +0,0 @@
-# Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
-# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-# LICENSE and NOTICE for details. LLNL-CODE-806117.
-#
-# This file is part of the MFEM library. For more information and source code
-# availability visit https://mfem.org.
-#
-# MFEM is free software; you can redistribute it and/or modify it under the
-# terms of the BSD-3 license. We welcome feedback and contributions, see file
-# CONTRIBUTING.md for details.
-
-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

@@ -1,718 +0,0 @@
-#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

@@ -1,313 +0,0 @@
-#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

@@ -1,385 +0,0 @@
-#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

@@ -1,183 +0,0 @@
-#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

@@ -1,174 +0,0 @@
-#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

@@ -1,212 +0,0 @@
-#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

@@ -1,168 +0,0 @@
-#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

@@ -1,199 +0,0 @@
-#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

@@ -1,92 +0,0 @@
-# Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
-# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-# LICENSE and NOTICE for details. LLNL-CODE-806117.
-#
-# This file is part of the MFEM library. For more information and source code
-# availability visit https://mfem.org.
-#
-# MFEM is free software; you can redistribute it and/or modify it under the
-# terms of the BSD-3 license. We welcome feedback and contributions, see file
-# CONTRIBUTING.md for details.
-
-# 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

miniapps/multiapp/CMakeLists.txt

@@ -0,0 +1,52 @@
+# 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.
+
+set(MESH_FILES
+  backward-facing-step.msh
+  channel-cylinder.msh
+)
+
+# Add a target to copy the mesh files from the source directory; used by sample
+# runs.
+set(SRC_MESH_FILES)
+foreach(MESH_FILE ${MESH_FILES})
+  list(APPEND SRC_MESH_FILES ${CMAKE_CURRENT_SOURCE_DIR}/${MESH_FILE})
+endforeach()
+add_custom_command(OUTPUT data_is_copied
+  COMMAND ${CMAKE_COMMAND} -E copy_if_different ${SRC_MESH_FILES} .
+  COMMAND ${CMAKE_COMMAND} -E touch data_is_copied
+  COMMENT "Copying multiapp miniapps data files ...")
+add_custom_target(copy_miniapps_multiapp_data DEPENDS data_is_copied)
+
+
+list(APPEND MULTIAPP_COMMON_SOURCES
+	multiapp.cpp)
+
+list(APPEND MULTIAPP_COMMON_HEADERS
+	multiapp.hpp)
+
+set(MULTIAPP_COMMON_FILES
+	EXTRA_SOURCES ${MULTIAPP_COMMON_SOURCES}
+	EXTRA_HEADERS ${MULTIAPP_COMMON_HEADERS})
+
+
+# Parallel apps.
+if (MFEM_USE_MPI)
+  add_mfem_miniapp(coupled-diffusion
+    MAIN coupled-diffusion.cpp
+    ${MFEM_MINIAPPS_COMMON_HEADERS} ${MULTIAPP_COMMON_FILES}
+    LIBRARIES mfem-common)
+  add_dependencies(coupled-diffusion copy_miniapps_multiapp_data)
+
+# Add parallel tests.
+  # if (MFEM_ENABLE_TESTING)
+  # endif()
+endif()

miniapps/multiapp/coupled-diffusion.cpp

@@ -0,0 +1,895 @@
+#include "mfem.hpp"
+#include "multiapp.hpp"
+#include <fstream>
+using namespace mfem;
+using namespace std;
+
+struct CaseContext
+{
+   int ser_ref = 1;         // Serial mesh refinement
+   int order = 3;           // Finite element order
+   bool visualization = true;// Visualization on/off
+   int grad_mode = 1;       // Gradient mode for the coupled operator - 0: finite difference,
+                            //                                          1: back/forward propagation
+   bool coupled = true;     // Coupled (true) vs. uncoupled (false) solves
+   int nl_iter = 50;        // Maximum number of nonlinear iterations
+   int lin_iter = 2000;     // Maximum number of linear iterations
+
+#if defined(MFEM_USE_DOUBLE)
+   real_t tol_nsolve = 1e-4;
+   real_t tol_lsolve = 1e-6;
+#elif defined(MFEM_USE_SINGLE)
+   real_t tol_nsolve = 1e-3;
+   real_t tol_lsolve = 1e-3;
+#else
+#error "Only single and double precision are supported!"
+   real_t tol_nsolve = 0;
+   real_t tol_lsolve = 0;
+#endif
+} ctx;
+
+void SetSolverParameters(IterativeSolver *solver, real_t rtol, real_t atol , int max_it,
+                         int print_level, bool iterative_mode);
+
+
+/// A functional diffusion coefficient (i.e., k(T))
+class FunctionalCoefficient : public Coefficient
+{
+public:
+   enum Mode { FUNC = 0, GRAD = 1};
+
+protected:
+   ParGridFunction *T_gf = nullptr;
+   real_t kref = 1.0;
+   real_t a0 = 0.0, a1 = 0.0, a2 = 0.0;
+   int findex = 0;
+   Mode mode = Mode::FUNC; // otherwise, grad
+
+public:
+   FunctionalCoefficient(ParGridFunction *T_gf, real_t kref): 
+                         T_gf(T_gf), kref(kref) { }
+
+   FunctionalCoefficient(ParGridFunction *T_gf, real_t kref, real_t a0):
+                         T_gf(T_gf), kref(kref), a0(a0) { findex = 1; }
+
+   FunctionalCoefficient(ParGridFunction *T_gf, real_t kref,
+                         real_t a0, real_t a1, real_t a2): T_gf(T_gf),
+                         kref(kref), a0(a0), a1(a1), a2(a2) { findex = 2; }
+
+    real_t Exponential(real_t x, bool eval_f) const
+    { 
+      real_t f = kref*exp(a0*x);
+      return (eval_f ? f : a0*f);
+    }
+    real_t Polynomial(real_t x, bool eval_f) const
+    { 
+      return (eval_f ? kref*(a0 + a1*x + a2*x*x) : kref*(a1 + 2*a2*x));
+    }
+
+    void SetMode(Mode mode) { this->mode = mode; }
+    Mode GetMode() const { return mode; }
+    
+    void UpdateGridFunction(ParGridFunction *gf) { T_gf = gf; }
+
+    real_t Eval(real_t x, bool eval_f) const
+    {
+        switch (findex)
+        {
+        case 1:
+         return Exponential(x, eval_f);
+        case 2:
+         return Polynomial(x, eval_f);
+        default:
+         return kref;
+        }
+    }
+
+    real_t Eval(ElementTransformation &Tr,
+              const IntegrationPoint &ip) override
+    {
+        real_t T = T_gf ? T_gf->GetValue(Tr, ip) : 0.0;
+        bool eval_f = (mode == Mode::FUNC);
+        return Eval(T, eval_f);
+    }
+};
+
+
+/// A coefficient defined by the product of grid functions, e.g. k(T) = prod_i x_i
+class GridFunctionProductCoefficient : public Coefficient
+{
+protected:
+   std::vector<ParGridFunction*> &x;
+
+public:
+   GridFunctionProductCoefficient(std::vector<ParGridFunction*> &x) : x(x) { }
+
+   real_t Eval(ElementTransformation &Tr, const IntegrationPoint &ip) override
+   {
+      real_t prod = 1.0;
+      for(size_t i = 0; i < x.size(); i++)
+      {
+         real_t val = x[i]->GetValue(Tr, ip);
+         prod *= val;
+      }
+      return prod;
+   }
+};
+
+class CoefficientIntegrator : public NonlinearFormIntegrator
+{
+protected:
+   FunctionalCoefficient *func = nullptr;
+   Vector shape;
+
+public:
+   CoefficientIntegrator(FunctionalCoefficient *func) : func(func) { }
+
+
+   void SetCoefficient(FunctionalCoefficient *f) { func = f; }
+
+   void AssembleElementVector(const FiniteElement &el,
+                              ElementTransformation &Tr,
+                              const Vector &elfun, Vector &elvect)
+   {
+      int dof = el.GetDof();
+      shape.SetSize(dof);
+      elvect.SetSize(dof);
+      elvect = 0.0;
+
+      const IntegrationRule *ir = &el.GetNodes();
+
+      for (int i = 0; i < ir->GetNPoints(); i++)
+      {
+         const IntegrationPoint &ip = ir->IntPoint(i);
+         el.CalcShape(ip, shape);
+         Tr.SetIntPoint(&ip);
+         real_t x = elfun * shape; // Evaluate the function at the integration point
+         real_t fval = func->Eval(x, true);
+         for (int j = 0; j < dof; j++)
+         {
+            elvect(j) += fval * shape(j);
+         }
+      }
+   }
+
+   void AssembleElementGrad(const FiniteElement &el, ElementTransformation &Tr,
+                            const Vector &elfun, DenseMatrix &elmat)
+   {
+      int dof = el.GetDof();
+      shape.SetSize(dof);
+      elmat.SetSize(dof);
+      elmat = 0.0;
+
+      const IntegrationRule *ir = &el.GetNodes();
+
+      for (int i = 0; i < ir->GetNPoints(); i++)
+      {
+         const IntegrationPoint &ip = ir->IntPoint(i);
+         el.CalcShape(ip, shape);
+         Tr.SetIntPoint(&ip);
+         real_t x = elfun * shape; // Evaluate the function at the integration point
+         real_t dfdx = func->Eval(x, false); // Evaluate the derivative of the function at the integration point
+         for (int j = 0; j < dof; j++)
+         {
+            elmat(j,j) += dfdx * shape(j); // Diagonal contribution to the Jacobian
+         }
+      }
+   }
+};
+
+
+class NonlinearDiffusionIntegrator : public NonlinearFormIntegrator
+{
+protected:
+   Coefficient *k;
+   Coefficient *dk;
+
+   Vector u, vec, shape;
+   DenseMatrix dshape, dshapedxt, adjJ;
+public:
+   NonlinearDiffusionIntegrator(Coefficient *kappa, Coefficient *dkappa) :
+                                k(kappa), dk(dkappa) { }
+
+   virtual void AssembleElementVector(const FiniteElement &el,
+                                      ElementTransformation &Tr,
+                                      const Vector &elfun, Vector &elvec)
+   {
+      int dim = el.GetDim();
+      int dof = el.GetDof();
+      real_t w;
+
+      elvec.SetSize(dof);
+      elvec = 0.0;
+
+      const IntegrationRule *ir = IntRule ? IntRule : &IntRules.Get(el.GetGeomType(), 2*el.GetOrder());
+      u.SetSize(dim);
+      vec.SetSize(dim);
+      dshape.SetSize(dof, dim);
+      adjJ.SetSize(dim, dim);
+
+      for (int i = 0; i < ir->GetNPoints(); i++)
+      {
+         const IntegrationPoint &ip = ir->IntPoint(i);
+         el.CalcDShape(ip, dshape);
+
+         Tr.SetIntPoint(&ip);
+         CalcAdjugate(Tr.Jacobian(), adjJ);
+         w = ip.weight / Tr.Weight();
+
+         dshape.MultTranspose(elfun, u);
+         adjJ.MultTranspose(u, vec);
+         if(k)
+         {
+            w *= k->Eval(Tr, ip);
+         }
+
+         vec *= w;
+         adjJ.Mult(vec, u);
+         dshape.AddMult(u, elvec);
+      }
+   }
+
+   void AssembleElementGrad(const FiniteElement &el, ElementTransformation &Tr,
+                            const Vector &elfun, DenseMatrix &elmat)
+   {
+      int dim = el.GetDim();
+      int dof = el.GetDof();
+      real_t w, k0 = 0.0, dk0 = 0.0;
+
+      elmat.SetSize(dof);
+      elmat = 0.0;
+
+      const IntegrationRule *ir = IntRule ? IntRule : &IntRules.Get(el.GetGeomType(), 2*el.GetOrder());
+      u.SetSize(dim);
+      shape.SetSize(dof);
+      vec.SetSize(dof);
+      dshape.SetSize(dof, dim);
+      dshapedxt.SetSize(dof, dim);
+
+      // f = grad(psi) * k(u) * grad(T)
+      // df/dT = grad(psi) ( k(u0) * grad(T) + k'(u0) * grad(u0) * T )
+      for (int i = 0; i < ir->GetNPoints(); i++)
+      {
+         const IntegrationPoint &ip = ir->IntPoint(i);
+         el.CalcShape(ip, shape);
+         el.CalcDShape(ip, dshape);
+
+         Tr.SetIntPoint(&ip);
+         w = ip.weight / Tr.Weight();
+
+         Mult(dshape, Tr.AdjugateJacobian(), dshapedxt);
+
+         k0  = k ? k->Eval(Tr, ip) : 0.0;
+         dk0 = dk ? dk->Eval(Tr, ip) : 0.0;
+
+         if(k0 != 0.0) // grad(psi) * k(u0) * grad(T)
+         {
+            real_t kdT = w*k0;
+            AddMult_a_AAt(kdT, dshapedxt, elmat);
+         }
+
+         if(dk0 != 0.0) // grad(psi) * (k'(T0) * grad(T0)) * T
+         {
+            dk0 = w*dk->Eval(Tr, ip);
+            dshapedxt.MultTranspose(elfun, u); // grad(T0) in physical space
+            u *= dk0; // k'(T0) * grad(T0)
+            dshapedxt.Mult(u, vec); // grad(psi) * k'(T0) * grad(T0)
+            AddMultVWt(vec, shape, elmat); // grad(psi) * k'(T0) * grad(T0) * T
+         }
+      }
+   }
+};
+
+/// An application that takes an input field T, and computes an output field k(T)
+// represented by the FunctionalCoefficient class.
+class DiffusionCoefficient : public GraphNode
+{
+public:
+   using Mode = FunctionalCoefficient::Mode;
+
+protected:
+   ParFiniteElementSpace &fes;
+   mutable ParGridFunction T, k;
+   mutable FunctionalCoefficient *kc;
+   // mutable Vector tdof, kdof, dk_dof, dT_dof;
+   mutable Mode mode = Mode::FUNC;
+
+   mutable ParNonlinearForm Nform;
+   mutable Operator *J = nullptr; // Jacobian for the nonlinear form
+
+   CoefficientIntegrator *coeff_integrator = nullptr;
+
+public:
+   DiffusionCoefficient(ParFiniteElementSpace &fes) :
+                        GraphNode(fes.GetTrueVSize()), fes(fes), T(&fes), k(&fes),
+                        kc(new FunctionalCoefficient(&T, 1.0, 5.0e-2)),
+                        Nform(&fes),
+                        coeff_integrator(new CoefficientIntegrator(kc))
+   {
+      k = 0.0;
+      T = 0.0;
+      k.ProjectCoefficient(*kc);
+
+      // Testing with the nonlinear form framework to compute k(T) and dk/dT
+      Nform.AddDomainIntegrator(coeff_integrator); // Transfer ownership
+      Nform.SetGradientType(Operator::Type::Hypre_ParCSR);
+      Nform.Setup();
+
+      SetInputOffsets(Array<int>({0, fes.GetTrueVSize()}));
+      SetOutputOffsets(Array<int>({0, fes.GetTrueVSize()}));
+   }
+
+   void SetMode(Mode mode) { this->mode = mode; }
+
+   FunctionalCoefficient* GetCoefficient() { return kc; }
+
+   void SetCoefficient(FunctionalCoefficient *fc)
+   {
+      if(kc) delete kc;
+      kc = fc;
+      kc->SetMode(mode);
+      kc->UpdateGridFunction(&T);
+      coeff_integrator->SetCoefficient(kc);
+   }
+
+   void Mult(const Vector &x, Vector &y) const override
+   {
+      BlockVector xb(x.GetData(), InputOffsets());
+      BlockVector yb(y.GetData(), OutputOffsets());
+
+      MultiVector xmv(1), ymv(1);
+      xmv.MakeRef(0, xb.GetBlock(0));
+      ymv.MakeRef(0, yb.GetBlock(0));
+
+      const_cast<DiffusionCoefficient*>(this)->Mult(xmv, ymv);
+   }
+
+   void Mult(const MultiVector &x, MultiVector &y) override
+   {
+      const Vector &tdof = x[0];
+      Vector &kdof = y[0];
+
+      Nform.Mult(tdof, kdof);
+      if(exec_mode == GraphNode::GRADIENT_MODE)
+      {
+         J = &Nform.GetGradient(tdof); // Store jacobian for JVP
+      }
+      else
+      {
+         J = nullptr; // Clear the Jacobian if not in gradient mode
+      }
+   }
+
+   // TODO: Possibly delete and only support MultiVector version of GradientMult
+   void GradientMult(const Vector &x, const Vector &dx, Vector &dy) const override
+   {
+      BlockVector xb(x.GetData(), InputOffsets());
+      BlockVector dxb(dx.GetData(), InputOffsets());
+      BlockVector dyb(dy.GetData(), OutputOffsets());
+
+      MultiVector xmv(1), dxmv(1), dymv(1);
+      xmv.MakeRef(0, xb.GetBlock(0));
+      dxmv.MakeRef(0, dxb.GetBlock(0));
+      dymv.MakeRef(0, dyb.GetBlock(0));
+      const_cast<DiffusionCoefficient*>(this)->GradientMult(xmv, dxmv, dymv);
+   }
+
+   void GradientMult(const MultiVector &x, const MultiVector &dx, MultiVector &dy) const override
+   {
+      const Vector &tdof = x[0];
+      const Vector &xadj = dx[0];
+      Vector &yadj = dy[0];
+
+      if(J)
+      {
+         J->Mult(xadj, yadj);
+      }
+      else
+      {
+         J = &Nform.GetGradient(tdof); // Store jacobian for JVP
+         J->Mult(xadj, yadj);
+      }
+   }
+
+   ~DiffusionCoefficient() override
+   {
+      if(kc) delete kc;
+   }
+};
+
+/// An application that takes n input fields x_i, and computes an output 
+/// field prod(x) := y = prod_i x_i.
+/// Also provides the derivative dy/dx_i = prod_{j!=i} x_j * dx_i/dx for i = 0,...,n-1.
+class ProductGridFunctions : public GraphNode
+{
+protected:
+
+   ParFiniteElementSpace &fes;
+   mutable std::vector<ParGridFunction*> x_gf;
+   mutable Vector dfdx;
+   mutable ParGridFunction y_gf;
+   mutable GridFunctionProductCoefficient prod_coeff;
+
+public:
+   ProductGridFunctions(ParFiniteElementSpace &fes, int n) :
+                     //   GraphNode(fes.GetTrueVSize()),
+                       GraphNode(fes.GetTrueVSize(), fes.GetTrueVSize() * n),
+                       fes(fes), x_gf(n),
+                       y_gf(&fes), prod_coeff(x_gf)
+   {
+      Array<int> offsets(n+1);
+      offsets[0] = 0;
+      for (int i = 0; i < n; i++)
+      {
+         x_gf[i]  = new ParGridFunction(&fes);
+         *x_gf[i] = 0.0;
+         offsets[i+1] = offsets[i] + fes.GetTrueVSize();
+      }
+      y_gf = 0.0;
+      y_gf.ProjectCoefficient(prod_coeff);
+
+      SetInputOffsets(offsets);
+      SetOutputOffsets(Array<int>({0, fes.GetTrueVSize()}));
+   }
+
+   void Mult(const Vector &x, Vector &y) const override
+   {
+      BlockVector xb(x.GetData(), InputOffsets());
+      BlockVector yb(y.GetData(), OutputOffsets());
+
+      MultiVector xmv(x_gf.size()), ymv(1);
+      for (size_t i = 0; i < x_gf.size(); i++)
+      {
+         xmv.MakeRef(i, xb.GetBlock(i));
+      }
+      ymv.MakeRef(0, yb.GetBlock(0));
+      const_cast<ProductGridFunctions*>(this)->Mult(xmv, ymv);
+   }
+
+   void Mult(const MultiVector &x, MultiVector &y) override
+   {
+      for (size_t i = 0; i < x_gf.size(); i++)
+      {
+         const Vector &x_dof = x[i];
+         x_gf[i]->SetFromTrueDofs(x_dof);
+      }
+
+      Field *out_field = OutputField(0);
+      Vector &y_dof = y[0];
+      y_gf.ProjectCoefficient(prod_coeff);
+      y_gf.GetTrueDofs(y_dof);
+   }
+
+   // TODO: Possibly delete and only support MultiVector version of GradientMult
+   void GradientMult(const Vector &x, const Vector &dx, Vector &dy) const override
+   {
+      BlockVector xb(x.GetData(), InputOffsets());
+      BlockVector dxb(dx.GetData(), InputOffsets());
+      BlockVector dyb(dy.GetData(), OutputOffsets());
+
+      MultiVector xmv(x_gf.size()), dxmv(x_gf.size()), dymv(1);
+      for (size_t i = 0; i < x_gf.size(); i++)
+      {
+         xmv.MakeRef(i, xb.GetBlock(i));
+         dxmv.MakeRef(i, dxb.GetBlock(i));
+      }
+      dymv.MakeRef(0, dyb.GetBlock(0));
+      const_cast<ProductGridFunctions*>(this)->GradientMult(xmv, dxmv, dymv);
+   }
+
+   void GradientMult(const MultiVector &x, const MultiVector &dx, MultiVector &dy) const override
+   {
+      // Jacobian vector product for y = prod_i x_i is:
+      // dy/dx = sum_i (prod_{j!=i} x_j * dx_i/dx)
+      for (size_t i = 0; i < x_gf.size(); i++)
+      {
+         const Vector &x_dof = x[i];
+         x_gf[i]->SetFromTrueDofs(x_dof); // Set all x_i
+      }
+
+      Vector &jvp = dy[0];
+      jvp = 0.0;
+      for (size_t i = 0; i < x_gf.size(); i++)
+      {
+         const Vector &x_dof = x[i];
+         const Vector &dx_dof = dx[i]; // Get dx_i/dx
+
+         x_gf[i]->SetFromTrueDofs(dx_dof); // Set x_i = dx_i/dx for i-th term in the product
+         y_gf.ProjectCoefficient(prod_coeff); // Recompute product with x_i replaced by dx_i/dx
+         y_gf.GetTrueDofs(dfdx); // Get prod_{j!=i} x_j * dx_i/dx for i-th term
+         jvp += dfdx; // Accumulate contribution from i-th term
+         x_gf[i]->SetFromTrueDofs(x_dof); // reset to original value for next iteration
+      }
+   }
+
+   ~ProductGridFunctions() override
+   {
+      for (size_t i = 0; i < x_gf.size(); i++)
+      {
+         if(x_gf[i]) delete x_gf[i];
+      }
+   }
+};
+
+
+/// An application that represents the nonlinear diffusion operator: f(T) = -Div(k(u) grad(T)) 
+/// with input field T and k, and output field f(T).
+class DiffusionOperator : public GraphNode
+{
+public:
+
+   // Mesh and finite element space
+   ParMesh &mesh;
+   ParFiniteElementSpace &fes;
+
+   /// Essential dof array.
+   Array<int> ess_tdofs;
+
+   /// Grid functions for the temperature and heat flux
+   mutable ParGridFunction T, k, dk;
+   mutable GridFunctionCoefficient k_gfc, dk_gfc;
+   mutable ParNonlinearForm Nform;
+   mutable ParLinearForm bform;
+   mutable Vector b;
+
+   ConstantCoefficient zero_coeff, one_coeff;
+
+   mutable FunctionalCoefficient *kc = nullptr;
+   mutable HypreParMatrix *dfdk_mat = nullptr, *dfdT_mat = nullptr;
+
+public:
+
+   DiffusionOperator(ParFiniteElementSpace &fes_) :
+                     // GraphNode(fes_.GetTrueVSize()),
+                     GraphNode(fes_.GetTrueVSize(),2*fes_.GetTrueVSize()),
+                     mesh(*fes_.GetParMesh()), fes(fes_),
+                     T(&fes), k(&fes), dk(&fes),
+                     k_gfc(&k), dk_gfc(&dk),
+                     Nform(&fes), bform(&fes),
+                     zero_coeff(0.0), one_coeff(1.0)
+   {
+      fes.GetBoundaryTrueDofs(ess_tdofs);
+      T = 0.0;
+      k = 0.0;
+      dk = 0.0;
+
+      bform.AddDomainIntegrator(new DomainLFIntegrator(one_coeff));
+      Nform.AddDomainIntegrator(new NonlinearDiffusionIntegrator(&k_gfc, &dk_gfc));
+      Nform.SetGradientType(Operator::Type::Hypre_ParCSR);
+
+      b.SetSize(fes.GetTrueVSize()); b = 0.0;
+      Assemble();
+
+      SetInputOffsets(Array<int>({0, fes.GetTrueVSize(), 2*fes.GetTrueVSize()}));
+      SetOutputOffsets(Array<int>({0, fes.GetTrueVSize()}));
+   }
+
+   void SetCoefficient(FunctionalCoefficient *fc) { kc = fc; }
+
+   void Assemble()
+   {
+      AssembleLinearForms();
+      AssembleBilinearForms();
+      AssembleNonlinearForms();
+   }
+
+   void AssembleBilinearForms()
+   {}
+
+   void AssembleNonlinearForms()
+   {
+      Nform.SetEssentialTrueDofs(ess_tdofs);
+      Nform.Setup();
+   }
+
+   void AssembleLinearForms()
+   {
+      bform.Assemble();
+      bform.ParallelAssemble(b);
+   }
+
+   void Mult(const Vector &x, Vector &y) const override
+   {
+      BlockVector xb(x.GetData(), InputOffsets());
+      BlockVector yb(y.GetData(), OutputOffsets());
+
+      MultiVector xmv(2), ymv(1);
+      xmv.MakeRef(0, xb.GetBlock(0));
+      xmv.MakeRef(1, xb.GetBlock(1));
+      ymv.MakeRef(0, yb.GetBlock(0));
+
+      const_cast<DiffusionOperator*>(this)->Mult(xmv, ymv);
+   }
+
+   void Mult(const MultiVector &x, MultiVector &y) override
+   {
+      const Vector &tdofs = x[0];
+      const Vector &kdofs = x[1];
+      Vector &fdofs = y[0];
+
+      k.SetFromTrueDofs(kdofs); // update for use in k_gfc
+
+      if(exec_mode == GraphNode::GRADIENT_MODE)
+      {
+         if(dfdT_mat) delete dfdT_mat;
+         if(dfdk_mat) delete dfdk_mat;
+
+         dk = 0.0;
+         k.SetFromTrueDofs(kdofs);
+         Operator* grad = &Nform.GetGradient(tdofs);
+         dfdT_mat = new HypreParMatrix(dynamic_cast<const HypreParMatrix&>(*grad)); // deep copy
+
+         dk = 1.0;
+         k  = 0.0;
+         grad = &Nform.GetGradient(tdofs);
+         dfdk_mat = new HypreParMatrix(dynamic_cast<const HypreParMatrix&>(*grad)); // deep copy
+      }
+      else
+      {
+         if(dfdT_mat) { delete dfdT_mat; dfdT_mat = nullptr; }
+         if(dfdk_mat) { delete dfdk_mat; dfdk_mat = nullptr; }
+      }
+
+      Nform.Mult(tdofs, fdofs);
+      fdofs.SetSubVector(ess_tdofs, 0.0);
+   }
+
+   // Exact block jacobian [df/dT, df/dk]
+   Operator& GetGradient(const Vector &x) const override
+   {
+      MFEM_ABORT("GetGradient not implemented for DiffusionOperator");
+   }
+
+   // TODO: Possibly delete and only support MultiVector version of GradientMult
+   void GradientMult(const Vector &x, const Vector &dx, Vector &dy) const override
+   {
+      BlockVector xb(x.GetData(), InputOffsets());
+      BlockVector dxb(dx.GetData(), InputOffsets());
+      BlockVector dyb(dy.GetData(), OutputOffsets());
+
+      Vector &Tadj = dxb.GetBlock(0);
+      Vector &kadj = dxb.GetBlock(1);
+      Vector &yadj = dyb.GetBlock(0);
+
+      Vector &tdofs = xb.GetBlock(0);
+      Vector &kdofs = xb.GetBlock(1);
+
+      dfdT_mat->Mult(Tadj, yadj);
+      dfdk_mat->AddMult(kadj, yadj);
+   }
+
+   void GradientMult(const MultiVector &x, const MultiVector &dx, MultiVector &dy) const override
+   {
+      const Vector &Tadj = dx[0];
+      const Vector &kadj = dx[1];
+      Vector &yadj = dy[0];
+
+      const Vector &tdofs = x[0];
+      const Vector &kdofs = x[1];
+
+      dfdT_mat->Mult(Tadj, yadj);
+      dfdk_mat->AddMult(kadj, yadj);
+   }
+
+   /// @brief Destroy the DiffusionOperator object
+   ~DiffusionOperator() override
+   {
+      if(dfdT_mat) delete dfdT_mat;
+      if(dfdk_mat) delete dfdk_mat;
+   }
+};
+
+
+int main(int argc, char *argv[])
+{
+   Mpi::Init();
+   Hypre::Init();
+
+   using GradMode = DAGraph::GradMode;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&ctx.order, "-o", "--order",
+                  "Finite element order (polynomial degree).");
+   args.AddOption(&ctx.visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&ctx.ser_ref, "-rs", "--serial-refine",
+                  "Number of times to refine the mesh in serial.");
+
+   args.AddOption(&ctx.grad_mode, "-gm", "--grad-mode",
+                  "Gradient mode for the coupled operator (0: exact, 1: finite difference, 2: algorithmic differentiation)");
+   args.AddOption(&ctx.coupled, "-cp", "--coupled", "-ucp", "--uncoupled",
+                  "Coupled (true) vs. uncoupled (false) solves.");
+   args.ParseCheck();
+
+
+   int order = ctx.order;
+   std::string mesh_file = "../../data/star.mesh";
+   Mesh *serial_mesh = new Mesh(mesh_file);
+   int dim = serial_mesh->Dimension();
+
+   for (int i = 0; i < ctx.ser_ref; ++i) { serial_mesh->UniformRefinement(); }
+   serial_mesh->SetCurvature(order, false, dim, Ordering::byNODES);
+
+   ParMesh pmesh(MPI_COMM_WORLD, *serial_mesh);
+   delete serial_mesh;
+   pmesh.UniformRefinement();
+
+   // Finite element spaces
+   H1_FECollection fec(order, dim);
+   ParFiniteElementSpace fes(&pmesh, &fec);
+
+   // Build all operator nodes
+   DiffusionCoefficient diff_coeff_1(fes);
+   diff_coeff_1.SetName("k(T1)");
+   diff_coeff_1.SetCoefficient(new FunctionalCoefficient(nullptr, 1.0, 3.5e-2));
+   // diff_coeff_1.SetCoefficient(new FunctionalCoefficient(nullptr, 1.0, 1.0, 0.1, 0.0));
+
+   DiffusionCoefficient diff_coeff_2(fes);
+   diff_coeff_2.SetName("k(T2)");
+   diff_coeff_2.SetCoefficient(new FunctionalCoefficient(nullptr, 1.0, 1.0, 2.0, 0.0));
+   // diff_coeff_2.SetCoefficient(new FunctionalCoefficient(nullptr, 1.5, 2.5e-2));
+
+   ProductGridFunctions prod_coeff(fes, 2);
+   prod_coeff.SetName("k(T1,T2)");
+
+   DiffusionOperator diff_op1(fes);
+   diff_op1.SetName("Div(k(T1,T2) grad(T1))");
+   diff_op1.SetCoefficient(diff_coeff_1.GetCoefficient());
+
+   DiffusionOperator diff_op2(fes);
+   diff_op2.SetName("Div(k(T1,T2) grad(T2))");
+   diff_op2.SetCoefficient(diff_coeff_2.GetCoefficient());
+
+
+   // Build the DAG in any order, and then sort it to ensure the correct execution order
+   DAGraph dag(5);
+   dag.AddOperator(&diff_coeff_1);
+   dag.AddOperator(&diff_op1, fes.GetTrueVSize());
+   dag.AddOperator(&diff_op2, fes.GetTrueVSize());
+   dag.AddOperator(&diff_coeff_2);
+   dag.AddOperator(&prod_coeff);
+
+   Vector k1vec(fes.GetTrueVSize()); k1vec = 0.0;
+   Vector k2vec(fes.GetTrueVSize()); k2vec = 0.0;
+   Vector kpvec(fes.GetTrueVSize()); kpvec = 0.0;
+
+   Vector k1adj(fes.GetTrueVSize()); k1adj = 0.0;
+   Vector k2adj(fes.GetTrueVSize()); k2adj = 0.0;
+   // Vector kpadj(fes.GetTrueVSize()); kpadj = 0.0;
+
+   // Input fields get data from 'x' in DAGraph::Mult(x, y)
+   Field T1_field(nullptr, nullptr);
+   Field T2_field(nullptr, nullptr);
+
+   // Write space for data and adjoint only needed
+   // for the intermediate fields k1, k2, and k_prod
+   Field k1_field(&k1vec, &k1adj);
+   Field k2_field(&k2vec, &k2adj);
+   Field kp_field(&kpvec, &kpvec); // can use same space for data & adjoint
+
+   // Output fields get data from 'y' in DAGraph::Mult(x, y)
+   Field f1_field(nullptr, nullptr);
+   Field f2_field(nullptr, nullptr);
+
+
+   // Add input and output to the DAG
+   int sz = fes.GetTrueVSize();
+   dag.AddInput(&T1_field, sz);
+   dag.AddInput(&T2_field, sz);
+   dag.AddOutput(&f1_field, sz);
+   dag.AddOutput(&f2_field, sz);
+
+   // Form connections between the nodes in the DAG
+   diff_coeff_1.AddInput(&T1_field);
+   diff_coeff_1.AddOutput(&k1_field);
+
+   diff_coeff_2.AddInput(&T2_field);
+   diff_coeff_2.AddOutput(&k2_field);
+
+   prod_coeff.AddInputs(&k1_field, &k2_field);
+   prod_coeff.AddOutput(&kp_field);
+
+   diff_op1.AddInput(&T1_field);
+   diff_op1.AddOutput(&f1_field);
+
+   diff_op2.AddInput(&T2_field);
+   diff_op2.AddOutput(&f2_field);
+
+   if(ctx.coupled)
+   {
+      diff_op1.AddInput(&kp_field); // kp_field
+      diff_op2.AddInput(prod_coeff.OutputField(0)); // Can also use kp_field directly
+   }
+   else
+   {
+      diff_op1.AddInput(&k1_field); // Can also use diff_coeff_1.OutputField(0)
+      diff_op2.AddInput(&k2_field); // Can also use diff_coeff_2.OutputField(0)
+   }
+
+   // Assemble DAG: topological sort, validate nodes, etc.
+   dag.Assemble();
+
+   std::string output_prefix = ctx.coupled ? "Coupled_Diffusion" : "Uncoupled_Diffusion";
+
+   if(Mpi::Root())
+   {
+      std::ofstream fout(output_prefix+"-dag.txt");
+      fout << "{\n";
+      dag.Save(fout);
+      fout << "}\n";
+      fout << std::flush;
+      fout.close();
+   }
+
+   // Set initial guess and boundary conditions for T1 and T2
+   Array<int> ess_tdofs;
+   fes.GetBoundaryTrueDofs(ess_tdofs);
+
+   int T1_idx = 0;
+   int T2_idx = 1;
+
+   BlockVector xb(dag.InputOffsets());
+   BlockVector yb(dag.OutputOffsets());
+
+   xb.GetBlock(T1_idx).Randomize();
+   xb.GetBlock(T2_idx).Randomize();
+   xb.GetBlock(T1_idx).SetSubVector(ess_tdofs, 0.0);
+   xb.GetBlock(T2_idx).SetSubVector(ess_tdofs, 0.0);
+
+   // Build the nonlinear solver and linear solver for the DAG
+   NewtonSolver newton_solver(pmesh.GetComm());
+   GMRESSolver linear_solver(pmesh.GetComm());
+   linear_solver.SetKDim(500);
+   SetSolverParameters(&newton_solver, ctx.tol_nsolve, 0.0, ctx.nl_iter, 1, true);
+   SetSolverParameters(&linear_solver, ctx.tol_lsolve, 0.0, ctx.lin_iter, 1, false);
+
+   newton_solver.SetPreconditioner(linear_solver);
+   linear_solver.SetPrintLevel(1);
+
+   // Set the gradient mode for the DAG and solve the coupled system
+   GradMode gm = static_cast<GradMode>(ctx.grad_mode);
+   dag.SetGradientMode(gm);
+   newton_solver.SetOperator(dag);
+   newton_solver.Mult(xb, yb);
+
+   ParaViewDataCollection *pv = nullptr;
+   if (ctx.visualization)
+   {
+      std::string pv_prefix;
+      switch (ctx.grad_mode)
+      {
+         case 0: pv_prefix = "FD"; break;
+         case 1: pv_prefix = "MF"; break;
+         default: pv_prefix = "Unknown"; break;
+      }
+
+      pv = new ParaViewDataCollection(output_prefix+"-"+pv_prefix, &pmesh);
+      pv->SetLevelsOfDetail(order);
+      pv->SetDataFormat(VTKFormat::BINARY);
+      pv->SetHighOrderOutput(true);
+
+      ParGridFunction T1_gf(&fes);
+      ParGridFunction T2_gf(&fes);
+      T1_gf.SetFromTrueDofs(yb.GetBlock(T1_idx));
+      T2_gf.SetFromTrueDofs(yb.GetBlock(T2_idx));
+
+      pv->RegisterField("T1", &T1_gf);
+      pv->RegisterField("T2", &T2_gf);
+      pv->Save();
+      delete pv;
+   }
+
+   std::cout << "Finished solving the coupled diffusion problem." << std::endl;
+   return 0;
+}
+
+void SetSolverParameters(IterativeSolver *solver, real_t rtol, real_t atol,
+                         int max_it, int print_level, bool iterative_mode)
+{
+    solver->SetRelTol(rtol);
+    solver->SetAbsTol(atol);
+    solver->SetMaxIter(max_it);
+    solver->SetPrintLevel(print_level);
+    solver->iterative_mode = iterative_mode;
+}

miniapps/multiapp/multiapp.cpp

@@ -0,0 +1,873 @@
+// 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.
+
+#include "multiapp.hpp"
+
+namespace mfem
+{
+
+DAGraph::~DAGraph()
+{
+    for(int i=0; i < nnodes; i++)
+    {
+        if(node_owned[i] && nodes[i]) delete nodes[i];
+    }
+    if(grad) delete grad;
+}
+
+void DAGraph::Assemble()
+{
+    // Sort graph nodes topologically to ensure correct execution order
+    // Ordering is not unique, hence, id->index maps are needed
+    TopologicalSort();
+
+    // Collect all fields from the nodes into the field map
+    CollectFieldMaps();
+
+    // Compute depth of the graph nodes
+    ComputeDepth();
+
+    // Validate each node
+    for (auto &node : nodes)
+    {
+        ValidateNode(*node);
+    }
+
+    // Update width and height of the DAG from offsets
+    // Check that the input and output offsets are consistent
+    ValidateOffsets();
+    width  = input_offsets.Last();
+    height = output_offsets.Last();
+
+    // Delete any existing gradient operator as node ordering may have changed
+    if (grad) delete grad;
+
+    assembled = true;
+}
+
+void DAGraph::ValidateOffsets()
+{
+    // Check that the input and output offsets are consistent
+    // with the number of inputs and outputs
+    if(InputFields().Size() > 1)
+    {
+        MFEM_ASSERT(input_offsets.Size() == InputFields().Size() + 1,
+                    "Input offsets size inconsistent with number of input fields");
+    }
+    else
+    {
+        input_offsets = Array<int>({0, nodes[0]->Width()});
+    }
+
+    if(OutputFields().Size() > 1)
+    {
+        MFEM_ASSERT(output_offsets.Size() == OutputFields().Size() + 1,
+                    "Output offsets size inconsistent with number of output fields");
+    }
+    else
+    {
+        output_offsets = Array<int>({0, nodes.Last()->Height()});
+    }
+}
+
+void DAGraph::ValidateNode(GraphNode &node)
+{
+    // Validate that the node's input and output fields are consistent with the graph's field map
+    auto inputs = node.InputFields();
+    auto outputs = node.OutputFields();
+
+    // Check offsets match width and height of the node
+    MFEM_ASSERT(node.InputOffsets().Last() == node.Width(),
+                "Node ID: " << node.ID() << " input offsets do not match node width.");
+    MFEM_ASSERT(node.OutputOffsets().Last() == node.Height(),
+                "Node ID: " << node.ID() << " output offsets do not match node height.");
+
+    // Check number of input and output fields match the offsets
+    MFEM_ASSERT(node.InputOffsets().Size() == inputs.Size() + 1,
+                "Node input offsets size inconsistent with number of input fields");
+    MFEM_ASSERT(node.OutputOffsets().Size() == outputs.Size() + 1,
+                "Node output offsets size inconsistent with number of output fields");
+
+    // Check that all input and output fields are registered in the graph's field map
+    for(auto input_field : inputs)
+    {
+        MFEM_ASSERT(fid_to_index.Has(input_field->ID()),
+                    "Input field ID " << input_field->ID() << " not found in graph's field map");
+    }
+    for(auto output_field : outputs)
+    {
+        MFEM_ASSERT(fid_to_index.Has(output_field->ID()),
+                    "Output field ID " << output_field->ID() << " not found in graph's field map");
+    }
+}
+
+void DAGraph::TopologicalSort()
+{
+    Array<int> sorted_indices;
+    sorted_indices.Reserve(nnodes);
+
+    Array<bool> visited(nnodes);
+    visited = false; // Initialize all nodes as unvisited
+
+    // Perform a depth-first search to sort the nodes topologically
+    std::function<void(int)> DepthFirstSearch = [&](int node_index)
+    {
+        if(visited[node_index]) return;
+        visited[node_index] = true;
+        auto node = nodes[node_index];
+        // Visit all nodes that this node depends on
+        for(auto input_field : node->InputFields())
+        {
+            for(int j=0; j < nnodes; j++)
+            {
+                auto other_node = nodes[j];
+                if(other_node == node) continue;
+                for(auto output_field : other_node->OutputFields())
+                {
+                    if(input_field->ID() == output_field->ID()) // Compare by unique ID
+                    {
+                        DepthFirstSearch(j);
+                    }
+                }
+            }
+        }
+        sorted_indices.push_back(node_index);
+    };
+
+    for(int i=0; i < nnodes; i++)
+    {
+        DepthFirstSearch(i);
+    }
+
+    nodes.Permute(sorted_indices);
+    node_owned.Permute(sorted_indices);
+
+    // Update the node indices after sorting
+    for(int i=0; i < nnodes; i++)
+    {
+        nodes[i]->SetNodeIndex(i);
+    }
+
+    sorted = true;
+}
+
+void DAGraph::ComputeDepth()
+{
+    // Compute depth of ordered nodes
+    node_depth.SetSize(nnodes);
+    node_depth = 0;
+    for(int i=0; i < nnodes; i++)
+    {
+        int max_depth = 0;
+        auto node = nodes[i];
+        for(auto input_field : node->InputFields())
+        {
+            for(int j=0; j < i; j++)
+            {
+                auto other_node = nodes[j];
+                if(other_node == node) continue;
+                for(auto output_field : other_node->OutputFields())
+                {
+                    if(input_field->ID() == output_field->ID()) // Compare by unique ID
+                    {
+                        max_depth = std::max(max_depth, node_depth[j] + 1);
+                    }
+                }
+            }
+        }
+        node_depth[i] = max_depth;
+    }
+}
+
+void DAGraph::CollectFieldMaps()
+{
+    MFEM_ASSERT(sorted, "DAGraph must be topologically sorted before collecting fields");
+
+    fid_to_index.clear();
+    fid_to_field.clear();
+
+    int nfields = 0;
+    for (auto f : InputFields())
+    {
+        fid_to_index.Register(f->ID(), nfields++);
+        fid_to_field.Register(f->ID(), f);
+    }
+
+    for (auto &node : nodes)
+    {
+        for (auto f : node->OutputFields())
+        {
+            if (!fid_to_index.Has(f->ID()))
+            {
+                fid_to_index.Register(f->ID(), nfields++);
+            }
+            if (!fid_to_field.Has(f->ID()))
+            {
+                fid_to_field.Register(f->ID(), f);
+            }
+        }
+    }
+    // TODO: Possibly add all intermediate fields from nodes to the graph's FieldCollection
+}
+
+void DAGraph::Mult(const Vector &x, Vector &y) const
+{
+    MFEM_ASSERT(width == x.Size(), "Input vector size (" << x.Size()
+                << ") must match matrix width (" << width << ")");
+
+    MFEM_ASSERT(height == y.Size(), "Output vector size (" << y.Size()
+                << ") must match matrix height (" << height << ")");
+
+    auto inputs  = InputFields();
+    auto outputs = OutputFields();
+
+    BlockVector xb(x.GetData(), input_offsets);
+    BlockVector yb(y.GetData(), output_offsets);
+    MultiVector xmv(inputs.Size()), ymv(outputs.Size());
+
+    // Set the data pointers of the input and output fields
+    // of the graph to point to the corresponding blocks of
+    // the input and output vectors
+    for(int i=0; i < inputs.Size(); i++)
+    {
+        xmv.MakeRef(i, xb.GetBlock(i));
+    }
+
+    for(int i=0; i < outputs.Size(); i++)
+    {
+        ymv.MakeRef(i, yb.GetBlock(i));
+    }
+
+    const_cast<DAGraph*>(this)->Mult(xmv, ymv);
+}
+
+void DAGraph::Mult(const MultiVector &x, MultiVector &y)
+{
+    auto inputs  = InputFields();
+    auto outputs = OutputFields();
+
+    MFEM_ASSERT(inputs.Size() == x.NumBlocks(), "Number of input blocks (" << x.NumBlocks()
+                << ") must match number of input fields (" << inputs.Size() << ")");
+
+    MFEM_ASSERT(outputs.Size() == y.NumBlocks(), "Number of output blocks (" << y.NumBlocks()
+                << ") must match number of output fields (" << outputs.Size() << ")");
+
+    for(int i=0; i < inputs.Size(); i++)
+    {
+        inputs[i]->SetData(const_cast<Vector*>(&x[i]));
+    }
+    for (int i=0; i < outputs.Size(); i++)
+    {
+        outputs[i]->SetData(&y[i]);
+    }
+
+    auto index_map = GetFieldIdToIndexMap();
+    auto fld_map = GetFieldIdToFieldMap();
+    int nfields = index_map.NumFields();
+    MultiVector ymv(nfields); // TODO: Should this be a member function?
+
+    // Assemble the multivector from the individual fields based on their IDs
+    // This multivector contains all input, output, and intermediate fields in the graph
+    for (auto const& [id, idx] : index_map)
+    {
+        if (fld_map.Has(id))
+        {
+            auto field = fld_map.Get(id);
+            ymv.MakeRef(idx, *field->Data());
+        }
+        else
+        {
+            MFEM_ABORT("Field ID " << id << " not found in field map");
+        }
+    }
+    
+    Execute(x, ymv);
+
+    for(auto &f : inputs)
+    {
+        f->SetData(nullptr);
+    }
+    for(auto &f : outputs)
+    {
+        f->SetData(nullptr);
+    }
+}
+
+void DAGraph::Execute(const MultiVector &x, MultiVector &y) const
+{
+    MFEM_ASSERT(assembled, "DAGraph must be assembled before calling Execute()");
+
+    MFEM_ASSERT(x.NumBlocks() == InputFields().Size(),
+                "Number of input blocks (" << x.NumBlocks()
+                << ") must match number of input fields (" << InputFields().Size() << ")");
+    
+    auto index_map = GetFieldIdToIndexMap();
+    
+    MFEM_ASSERT(y.NumBlocks() == index_map.NumFields(),
+                "Number of output blocks (" << y.NumBlocks()
+                << ") must match number of fields (" << index_map.NumFields() << ")");
+
+    auto inputs  = InputFields();
+    for(int i=0; i < inputs.Size(); i++)
+    {
+        int idx = index_map.Get(inputs[i]->ID());
+        if(&y[idx] != &x[i]) // copy data, if address is different
+        {
+            y[idx] = x[i];
+        }
+    }
+
+    if(input_type == InputType::VECTOR)
+    {
+        x_node.SetSize(MaxWidth());
+        y_node.SetSize(MaxHeight());
+
+        for (auto node : nodes)
+        {
+            x_node.SetSize(node->Width());
+            y_node.SetSize(node->Height());
+
+            // Assemble input fields into a single vector for the node
+            auto node_inputs = node->InputFields();
+            auto ioffsets = node->InputOffsets();
+            for (int i=0; i < node_inputs.Size(); i++)
+            {
+                auto in_field = node_inputs[i];
+                int idx = index_map.Get(in_field->ID());
+                x_node.SetVector(y[idx],ioffsets[i]);
+            }
+
+            node->Mult(x_node, y_node);
+
+            // Disassemble output vector back
+            auto node_outputs = node->OutputFields();
+            BlockVector ynb(y_node.GetData(), node->OutputOffsets());
+            for (int i=0; i < node_outputs.Size(); i++)
+            {
+                auto out_field = node_outputs[i];
+                int idx = index_map.Get(out_field->ID());
+                y[idx] = ynb.GetBlock(i);
+            }
+        }
+    }
+    else if(input_type == InputType::MULTIVECTOR)
+    {
+        for (auto node : nodes)
+        {
+            auto node_inputs = node->InputFields();
+            auto node_outputs = node->OutputFields();
+            xmv_node.SetNumBlocks(node_inputs.Size());
+            ymv_node.SetNumBlocks(node_outputs.Size());
+
+            for (int i=0; i < node_inputs.Size(); i++)
+            {
+                int idx = index_map.Get(node_inputs[i]->ID());
+                xmv_node.MakeRef(i, y[idx]);
+            }
+            for (int i=0; i < node_outputs.Size(); i++)
+            {
+                int idx = index_map.Get(node_outputs[i]->ID());
+                ymv_node.MakeRef(i, y[idx]);
+            }
+            node->Mult(xmv_node, ymv_node);
+        }
+    }
+    else if(input_type == InputType::NONE)
+    {
+        Vector x_unused, y_unused;
+        for (auto node : nodes)
+        {
+            node->Mult(x_unused, y_unused);
+        }
+    }
+    else
+    {
+        MFEM_ABORT("DAGraph::Execute() not implemented for input type: "
+                    << static_cast<int>(input_type));
+    }
+}
+
+Operator& DAGraph::GetGradient(const Vector &x) const
+{
+    // TODO: Should/could be removed
+    if(grad_mode == GradMode::FINITE_DIFF)
+    {
+        if(!grad)
+        {
+            grad = new future::FDJacobian(*this, x, 1e-6);
+        }
+        else
+        {
+            grad->GetGradient(x); // Update the FDJacobian with new point x
+        }
+        return *grad;
+    }
+
+    MFEM_ASSERT(static_cast<int>(grad_mode) < static_cast<int>(GradMode::NONE),
+                "DAGraph::GetGradient() called with invalid grad_mode: "
+                << static_cast<int>(grad_mode));
+
+    if(!grad)
+    {
+        grad = new GraphGradient(const_cast<DAGraph&>(*this));
+    }
+
+    if(grad_mode == GradMode::ASSEMBLED)
+    {
+        return grad->GetGradient(x); // Assemble the Jacobian matrix
+    }
+    else // GradMode::MATRIX_FREE
+    {
+        dynamic_cast<GraphGradient*>(grad)->Update(x); // Update the GraphGradient with new point x
+    }
+
+    return *grad;
+}
+
+GraphGradient::GraphGradient(DAGraph &dag) : Operator(dag.Height(), dag.Width()),
+                                             graph(&dag)
+{
+    MFEM_ASSERT(graph->IsAssembled(), "GraphGradient requires an assembled DAGraph.");
+    MFEM_ASSERT(graph->IsSorted(), "GraphGradient requires a topologically sorted DAGraph.");
+
+    auto index_map = graph->GetFieldIdToIndexMap();
+    auto field_map = graph->GetFieldIdToFieldMap();
+
+    MFEM_ASSERT(index_map.NumFields() == field_map.NumFields(),
+                "Mismatch in number of fields between index_map and field_map");
+
+    x_work.DeleteAll(); // Clear any existing pointers
+    x_work.SetSize(index_map.NumFields());
+    x_work = nullptr; // Initialize all pointers to nullptr
+    xlin.SetNumBlocks(index_map.NumFields());
+
+    for (auto const& [id, idx] : index_map)
+    {
+        MFEM_ASSERT(idx >= 0 && idx < x_work.Size(), "Index out of bounds for field ID: " << id);
+        MFEM_ASSERT(field_map.Has(id), "Field ID not found in field_map: " << id);
+
+        if(x_work[idx] == nullptr)
+        {
+            x_work[idx] = new Vector(); // Allocate a new Vector for this field
+        }
+        xlin.MakeRef(idx, *x_work[idx]); // Make xlin refer to the allocated Vector
+    }
+
+}
+
+void GraphGradient::Update(const Vector &x)
+{
+    MFEM_ASSERT(graph != nullptr, "GraphGradient operator requires a non-null DAGraph pointer.");
+
+    auto set_exec_mode = [&](DAGraph::ExecutionMode mode)
+    {
+        for (auto &node : graph->Nodes())
+        {
+            node->SetExecutionMode(mode);
+        }
+    };
+
+    auto inputs = graph->InputFields();
+    BlockVector xb(x.GetData(), graph->InputOffsets());
+    MultiVector xmv(inputs.Size());
+    for(int i=0; i < inputs.Size(); i++)
+    {
+        xmv.MakeRef(i, xb.GetBlock(i));
+    }
+
+    set_exec_mode(DAGraph::ExecutionMode::GRADIENT_MODE);
+    graph->Execute(xmv, xlin); // Forward pass to populate fields for gradient computations
+    set_exec_mode(DAGraph::ExecutionMode::DEFAULT_MODE); // Reset execution mode for forward pass
+}
+
+void GraphGradient::Mult(const Vector &x, Vector &y) const
+{
+    MFEM_ASSERT(x.Size() == graph->Width(), "Input vector size (" << x.Size()
+                << ") must match graph width (" << graph->Width() << ")");
+
+    MFEM_ASSERT(y.Size() == graph->Height(), "Output vector size (" << y.Size()
+                << ") must match graph height (" << graph->Height() << ")");
+
+    auto in_offsets  = graph->InputOffsets();
+    auto out_offsets = graph->OutputOffsets();
+
+    auto inputs  = graph->InputFields();
+    auto outputs = graph->OutputFields();
+
+    BlockVector xb(x.GetData(), in_offsets);
+    BlockVector yb(y.GetData(), out_offsets);
+    MultiVector xmv(inputs.Size()), ymv(outputs.Size());
+
+    for(int i=0; i < inputs.Size(); i++)
+    {
+        xmv.MakeRef(i, xb.GetBlock(i));
+    }
+
+    for(int i=0; i < outputs.Size(); i++)
+    {
+        ymv.MakeRef(i, yb.GetBlock(i));
+    }
+
+    const_cast<GraphGradient*>(this)->Mult(xmv, ymv); // Forward mode: compute JVP, y = J(z) * x
+}
+
+void GraphGradient::Mult(const MultiVector &x, MultiVector &y)
+{
+    auto inputs  = graph->InputFields();
+    auto outputs = graph->OutputFields();
+
+    MFEM_ASSERT(inputs.Size() == x.NumBlocks(), "Number of input blocks (" << x.NumBlocks()
+                << ") must match number of input fields (" << inputs.Size() << ")");
+
+    MFEM_ASSERT(outputs.Size() == y.NumBlocks(), "Number of output blocks (" << y.NumBlocks()
+                << ") must match number of output fields (" << outputs.Size() << ")");
+
+    for(int i=0; i < inputs.Size(); i++)
+    {
+        inputs[i]->SetAdjoint(const_cast<Vector*>(&x[i]));
+    }
+    for (int i=0; i < outputs.Size(); i++)
+    {
+        outputs[i]->SetAdjoint(&y[i]);
+    }
+
+    auto index_map = graph->GetFieldIdToIndexMap();
+    auto fld_map = graph->GetFieldIdToFieldMap();
+    int nfields = index_map.NumFields();
+    MultiVector ymv(nfields); // TODO: Should this be a member function?
+
+    // Assemble the multivector from the individual fields based on their IDs
+    // This multivector contains all input, output, and intermediate fields in the graph
+    for (auto const& [id, idx] : index_map)
+    {
+        if (fld_map.Has(id))
+        {
+            auto field = fld_map.Get(id);
+            ymv.MakeRef(idx, *field->Adjoint());
+        }
+        else
+        {
+            MFEM_ABORT("Field ID " << id << " not found in field map");
+        }
+    }
+
+    Forward(x, ymv); // Forward mode: compute JVP, y = J(z) * x
+
+    for (auto &f : inputs)
+    {
+        f->SetAdjoint(nullptr);
+    }
+    for (auto &f : outputs)
+    {
+        f->SetAdjoint(nullptr);
+    }
+}
+
+void GraphGradient::MultTranspose(const Vector &x, Vector &y) const
+{
+    MFEM_ASSERT(x.Size() == graph->Height(), "Input vector size (" << x.Size()
+                << ") must match graph height (" << graph->Height() << ")");
+    MFEM_ASSERT(y.Size() == graph->Width(), "Output vector size (" << y.Size()
+                << ") must match graph width (" << graph->Width() << ")");
+
+    auto in_offsets  = graph->InputOffsets();
+    auto out_offsets = graph->OutputOffsets();
+
+    auto inputs  = graph->InputFields();
+    auto outputs = graph->OutputFields();
+
+    BlockVector xb(x.GetData(), out_offsets);
+    BlockVector yb(y.GetData(), in_offsets);
+    MultiVector xmv(outputs.Size()), ymv(inputs.Size());
+
+    for(int i=0; i < inputs.Size(); i++)
+    {
+        xmv.MakeRef(i, xb.GetBlock(i));
+    }
+
+    for(int i=0; i < outputs.Size(); i++)
+    {
+        ymv.MakeRef(i, yb.GetBlock(i));
+    }
+    const_cast<GraphGradient*>(this)->MultTranspose(xmv, ymv); // Reverse mode: compute VJP, y = J(z)^T * x
+}
+
+void GraphGradient::MultTranspose(const MultiVector &x, MultiVector &y)
+{
+    auto inputs  = graph->InputFields();
+    auto outputs = graph->OutputFields();
+
+    MFEM_ASSERT(outputs.Size() == x.NumBlocks(), "Number of input blocks (" << x.NumBlocks()
+                << ") must match number of output fields (" << outputs.Size() << ")");
+
+    MFEM_ASSERT(inputs.Size() == y.NumBlocks(), "Number of output blocks (" << y.NumBlocks()
+                << ") must match number of input fields (" << inputs.Size() << ")");
+
+    for(int i=0; i < outputs.Size(); i++)
+    {
+        outputs[i]->SetAdjoint(const_cast<Vector*>(&x[i]));
+    }
+    for (int i=0; i < inputs.Size(); i++)
+    {
+        inputs[i]->SetAdjoint(&y[i]);
+    }
+
+    auto index_map = graph->GetFieldIdToIndexMap();
+    auto fld_map = graph->GetFieldIdToFieldMap();
+    int nfields = index_map.NumFields();
+    MultiVector ymv(nfields); // TODO: Should this be a member function?
+
+    for(auto const& [id, idx] : index_map)
+    {
+        if (fld_map.Has(id))
+        {
+            auto field = fld_map.Get(id);
+            ymv.MakeRef(idx, *field->Adjoint());
+        }
+        else
+        {
+            MFEM_ABORT("Field ID " << id << " not found in field map");
+        }
+    }
+
+    Reverse(x, ymv); // Reverse mode: compute VJP, y = J(z)^T * x
+
+    for (auto &f : outputs)
+    {
+        f->SetAdjoint(nullptr);
+    }
+    for (auto &f : inputs)
+    {
+        f->SetAdjoint(nullptr);
+    }
+}
+
+void GraphGradient::Forward(const MultiVector &x, MultiVector &y) const
+{
+    MFEM_ASSERT(x.NumBlocks() == graph->InputFields().Size(),
+                "Number of input blocks (" << x.NumBlocks()
+                << ") must match number of input fields (" << graph->InputFields().Size() << ")");
+    
+    auto in_type = graph->GetInputType();
+    auto index_map  = graph->GetFieldIdToIndexMap();
+    auto field_map = graph->GetFieldIdToFieldMap();
+
+    MFEM_ASSERT(y.NumBlocks() == index_map.NumFields(),
+                "Number of output blocks (" << y.NumBlocks()
+                << ") must match number of fields (" << index_map.NumFields() << ")");
+
+    auto inputs  = graph->InputFields();
+    for(int i=0; i < inputs.Size(); i++)
+    {
+        int idx = index_map.Get(inputs[i]->ID());
+        if(&y[idx] != &x[i]) // copy data, if address is different
+        {
+            y[idx] = x[i];
+        }
+    }
+
+    if(in_type == InputType::VECTOR)
+    {
+        x0.SetSize(graph->MaxWidth());
+        dx.SetSize(graph->MaxWidth());
+        dy.SetSize(graph->MaxHeight());
+
+        auto nodes = graph->Nodes();
+        for (auto node : nodes)
+        {
+            x0.SetSize(node->Width());
+            dx.SetSize(node->Width());
+            dy.SetSize(node->Height());
+
+            // Assemble input fields into a single vector for the node
+            auto node_inputs = node->InputFields();
+            auto ioffsets = node->InputOffsets();
+            for(int i=0; i < node_inputs.Size(); i++)
+            {
+                auto in_field = node_inputs[i];
+                MFEM_ASSERT(index_map.Has(in_field->ID()), "Input field ID not found in index_map");
+                int idx = index_map.Get(in_field->ID());
+                x0.SetVector(xlin[idx], ioffsets[i]);
+                dx.SetVector(y[idx], ioffsets[i]);
+            }
+
+            node->GradientMult(x0, dx, dy); // Compute JVP for the node
+
+            // Disassemble output vector back
+            auto node_outputs = node->OutputFields();
+            BlockVector ynb(dy.GetData(), node->OutputOffsets());
+            for(int i=0; i < node_outputs.Size(); i++)
+            {
+                auto out_field = node_outputs[i];
+                MFEM_ASSERT(index_map.Has(out_field->ID()), "Output field ID not found in index_map");
+                int idx = index_map.Get(out_field->ID());
+                y[idx] = ynb.GetBlock(i);
+            }
+        }
+    }
+    else if(in_type == InputType::MULTIVECTOR)
+    {
+        auto nodes = graph->Nodes();
+        for (auto node : nodes)
+        {
+            auto node_inputs = node->InputFields();
+            auto node_outputs = node->OutputFields();
+            x0_mv.SetNumBlocks(node_inputs.Size());
+            dx_mv.SetNumBlocks(node_inputs.Size());
+            dy_mv.SetNumBlocks(node_outputs.Size());
+
+            for(int i=0; i < node_inputs.Size(); i++)
+            {
+                int idx = index_map.Get(node_inputs[i]->ID());
+                x0_mv.MakeRef(i, xlin[idx]);
+                dx_mv.MakeRef(i, y[idx]);
+            }
+            for(int i=0; i < node_outputs.Size(); i++)
+            {
+                int idx = index_map.Get(node_outputs[i]->ID());
+                dy_mv.MakeRef(i, y[idx]);
+            }
+            node->GradientMult(x0_mv, dx_mv, dy_mv); // Compute JVP for the node
+        }
+    }
+    else if(in_type == InputType::NONE)
+    {
+        Vector x_unused, dx_unused, dy_unused;
+        auto nodes = graph->Nodes();
+        for (auto node : nodes)
+        {
+            node->GradientMult(x_unused, dx_unused, dy_unused);
+        }
+    }
+    else
+    {
+        MFEM_ABORT("GraphGradient::Forward() not implemented for input type: "
+                    << static_cast<int>(in_type));
+    }
+}
+
+void GraphGradient::Reverse(const MultiVector &x, MultiVector &y) const
+{
+    MFEM_ASSERT(x.NumBlocks() == graph->OutputFields().Size(),
+                "Number of input blocks (" << x.NumBlocks()
+                << ") must match number of output fields (" << graph->OutputFields().Size() << ")");
+
+    auto in_type = graph->GetInputType();
+    auto index_map  = graph->GetFieldIdToIndexMap();
+    auto field_map = graph->GetFieldIdToFieldMap();
+    int nnodes  = graph->Size();
+
+    MFEM_ASSERT(y.NumBlocks() == index_map.NumFields(),
+                "Number of output blocks (" << y.NumBlocks()
+                << ") must match number of fields (" << index_map.NumFields() << ")");
+
+    auto outputs  = graph->OutputFields();
+    for(int i=0; i < outputs.Size(); i++)
+    {
+        int idx = index_map.Get(outputs[i]->ID());
+        if(&y[idx] != &x[i]) // copy data, if address is different
+        {
+            y[idx] = x[i];
+        }
+    }
+
+    if(in_type == InputType::VECTOR)
+    {
+        x0.SetSize(graph->MaxWidth());
+        dx.SetSize(graph->MaxHeight());
+        dy.SetSize(graph->MaxWidth());
+
+        for (int i=nnodes-1; i >= 0; i--)
+        {
+            auto node = graph->GetNode(i);
+            x0.SetSize(node->Width());
+            dx.SetSize(node->Height());
+            dy.SetSize(node->Width());
+
+            auto node_inputs = node->InputFields();
+            auto ioffsets = node->InputOffsets();
+            for(int i=0; i < node_inputs.Size(); i++)
+            {
+                auto in_field = node_inputs[i];
+                MFEM_ASSERT(index_map.Has(in_field->ID()), "Input field ID not found in index_map");
+                int idx = index_map.Get(in_field->ID());
+                x0.SetVector(xlin[idx], ioffsets[i]);
+            }
+
+            auto node_outputs = node->OutputFields();
+            auto ooffsets = node->OutputOffsets();
+            for(int i=0; i < node_outputs.Size(); i++)
+            {
+                auto out_field = node_outputs[i];
+                MFEM_ASSERT(index_map.Has(out_field->ID()), "Output field ID not found in index_map");
+                int idx = index_map.Get(out_field->ID());
+                dx.SetVector(y[idx], ooffsets[i]);
+            }
+
+            node->GradientMultTranspose(x0, dx, dy); // Compute JVP for the node
+
+            BlockVector dynb(dy.GetData(), node->InputOffsets());
+            for(int i=0; i < node_inputs.Size(); i++)
+            {
+                int idx = index_map.Get(node_inputs[i]->ID());
+                y[idx] = dynb.GetBlock(i);
+            }
+        }
+    }
+    else if(in_type == InputType::MULTIVECTOR)
+    {
+        for (int i=nnodes-1; i >= 0; i--)
+        {
+            auto node = graph->GetNode(i);
+            auto node_inputs = node->InputFields();
+            auto node_outputs = node->OutputFields();
+            x0_mv.SetNumBlocks(node_inputs.Size());
+            dx_mv.SetNumBlocks(node_outputs.Size());
+            dy_mv.SetNumBlocks(node_inputs.Size());
+
+            for(int i=0; i < node_inputs.Size(); i++)
+            {
+                int idx = index_map.Get(node_inputs[i]->ID());
+                x0_mv.MakeRef(i, xlin[idx]);
+                dy_mv.MakeRef(i, y[idx]);
+            }
+            for(int i=0; i < node_outputs.Size(); i++)
+            {
+                int idx = index_map.Get(node_outputs[i]->ID());
+                dx_mv.MakeRef(i, y[idx]);
+            }
+            node->GradientMultTranspose(x0_mv, dx_mv, dy_mv); // Compute JVP for the node
+        }
+    }
+    else if(in_type == InputType::NONE)
+    {
+        Vector x_unused, dx_unused, dy_unused;
+        for (int i=nnodes-1; i >= 0; i--)
+        {
+            auto node = graph->GetNode(i);
+            node->GradientMultTranspose(x_unused, dx_unused, dy_unused); // Compute VJP for the node
+        }
+    }
+    else
+    {
+        MFEM_ABORT("GraphGradient::Reverse() not implemented for input type: "
+                    << static_cast<int>(in_type));
+    }
+}
+
+Operator& GraphGradient::GetGradient(const Vector &x) const
+{
+    // Used to build Jacobian matrix
+    MFEM_ABORT("GraphGradient::GetGradient() not implemented");
+}
+
+
+} // namespace mfem

miniapps/multiapp/multiapp.hpp

@@ -0,0 +1,838 @@
+// 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_MULTIAPP_HPP
+#define MFEM_MULTIAPP_HPP
+
+#include "mfem.hpp"
+
+namespace mfem
+{
+
+/// Forward declarations needed below
+class Field;
+class FieldCollection;
+class GraphNode;
+class DAGraph;
+class GraphGradient;
+
+
+/// @brief Base class for storing data (Vector) and distinguishing
+/// fields variables
+class Field
+{
+public:
+    enum Type ///< Not used for now, but could be used to distinguish between input/output fields
+    {
+        INPUT , ///< Input field
+        OUTPUT, ///< Output field
+        DEFAULT ///< Any field
+    };
+
+    friend class GraphNode;
+
+private:
+    Type type = Type::DEFAULT;
+    inline static int next_id = 0;
+
+protected:
+    Vector *data = nullptr;
+    Vector *adjoint = nullptr; // For storing derivative info
+    int id = -1; // initialized to invalid id
+
+    std::string name; // Optional name for the field
+    Operator *oper  = nullptr; // Operator that outputs this field
+
+    int GetValidID(int id_, int lb=0, int ub = std::numeric_limits<int>::max())
+    {
+        return (id_ >= lb && id_ <= ub) ? id_ : next_id++;
+    }
+
+public:
+
+    ///@brief Constructor for a Field of type Type with optional ID
+    Field(Vector *field, Vector *adjoint, Type type, int id_ = -1) :
+          type(type), data(field), adjoint(adjoint), id(GetValidID(id_)),
+          name("Field_" + std::to_string(id)) { }
+
+    ///@brief Constructor for a Field of Default type with optional ID
+    Field(Vector *field, Vector *adjoint, int id_ = -1) : 
+          Field(field, adjoint, Type::DEFAULT, id_) { }
+
+    ///@brief Constructor for an input field
+    Field(Vector *field, int id_ = -1) :
+          Field(field, nullptr, Type::DEFAULT, id_) { }
+
+    ///@brief Constructor for a Field of type Type
+    Field(Vector *field, Type type, int id_ = -1) :
+          Field(field, nullptr, type, id_) { }
+
+    ///@brief Get the stored internally stored data pointer
+    Vector* Data() const { return data; }
+    Vector* Adjoint() const { return adjoint; }
+    Operator* GetOperator() const { return oper; }
+
+    ///@brief Set the internally stored data pointer
+    virtual void SetData(Vector *field) { data = field; }
+    virtual void SetAdjoint(Vector *adj) { adjoint = adj; }
+    virtual void SetOperator(Operator *op) { oper = op; }
+
+    virtual void GetData(Vector &field) const { field = *data; }
+    virtual void GetAdjoint(Vector &adj) const { adj = *adjoint; }
+
+    std::string Name() const { return name; }
+    void SetName(const std::string &n) { name = n; }
+    int ID() const { return id; }
+
+    void SetID(int i)
+    {
+        MFEM_ASSERT(i >= 0, "ID must be non-negative.");
+        id = i;
+    }
+
+    bool IsInput() const {return (type == Type::INPUT);}
+    bool IsOutput() const {return (type == Type::OUTPUT);}
+    bool IsDefault() const {return (type == Type::DEFAULT);}
+
+    virtual ~Field() = default;
+
+protected:
+
+    ///@brief Set the type of the field (prevents changing type of input/output fields)
+    void SetType(Type t)
+    {
+        type = t;
+    }
+};
+
+/// @brief A collection of Fields, each identified by a name
+class FieldCollection
+{
+public:
+    using FieldMap = GenericFieldMap<std::string, Field*>;
+    using IndexMap = GenericFieldMap<std::string, int>;
+
+private:
+    std::string name; /// Name of the collection
+    Operator *oper = nullptr; /// Operator associated with this collection (not owned)
+    FieldMap fields; /// Map from field name to Field pointer
+    IndexMap index_map; /// Map from field name to index in input/output vectors
+
+    Array<Field*> input_fields;  // Input fields for this node
+    Array<Field*> output_fields; // Output fields for this node
+
+public:
+
+    FieldCollection() = default;
+
+    /// @brief Constructor with collection name and optional associated operator
+    FieldCollection(std::string collection_name, Operator *op = nullptr):
+                    name(collection_name), oper(op) {}
+
+    /// @brief Constructor with associated operator and default collection name
+    FieldCollection(Operator *op) : name("FieldCollection"), oper(op) {}
+
+    /// @brief Get the number of fields in the collection
+    int Size() const { return fields.NumFields(); }
+
+    /// @brief Set the name of the collection
+    void SetName(const std::string &collection_name) { name = collection_name;}
+
+    /// @brief Get the name of the collection
+    std::string Name() const { return name; }
+
+    /// @brief Set the operator associated with this collection
+    void SetOperator(Operator *op){ oper = op; }
+
+    /// @brief Get the operator associated with this collection
+    const Operator* GetOperator() const { return oper; }
+
+    /// @brief Get the field associated with the given name, or nullptr if not found
+    Field* GetField(const std::string &field_name) const
+    {
+        return fields.Get(field_name);
+    }
+
+    /// @brief Add a field to the collection with a given name and ownership flag
+    void AddField(const std::string &field_name, Field *field, bool own = false)
+    {
+        if(fields.Has(field_name))
+        {
+            MFEM_WARNING("FieldCollection::AddField: Field with name "
+                         << field_name << " already exists. Replacing existing field.");
+        }
+        fields.Register(field_name, field, own);
+    }
+
+    void AddInput(const std::string &field_name,
+                  Field *field, bool own = false)
+    {
+        bool has_field = fields.Has(field_name);
+        bool has_index = index_map.Has(field_name);
+        if(has_field && has_index)
+        {
+            int i = index_map.Get(field_name);
+            input_fields[i] = field;
+        }
+        else
+        {
+            input_fields.push_back(field);
+            index_map.Register(field_name, input_fields.Size() - 1);
+        }
+        AddField(field_name, field, own);
+    }
+
+    void AddOutput(const std::string &field_name,
+                   Field *field, bool own = false)
+    {
+        bool has_field = fields.Has(field_name);
+        bool has_index = index_map.Has(field_name);
+        if(has_field && has_index)
+        {
+            int i = index_map.Get(field_name);
+            output_fields[i] = field;
+        }
+        else
+        {
+            output_fields.push_back(field);
+            index_map.Register(field_name, output_fields.Size() - 1);
+        }
+
+        AddField(field_name, field, own);
+        if(field->GetOperator() == nullptr)
+        {
+            field->SetOperator(oper);
+        }
+    }
+
+    Array<Field*>& InputFields() { return input_fields; }
+    Array<Field*>& OutputFields() { return output_fields; }
+
+    Field* InputField(int i) const { return input_fields[i]; }
+    Field *InputField(const std::string &field_name) const
+    {
+        bool has_index = index_map.Has(field_name);
+        if(!has_index)
+        {
+            MFEM_WARNING("FieldCollection::InputField: Field with name "
+                         << field_name << " does not exist in the collection.");
+            return nullptr;
+        }
+
+        int index = index_map.Get(field_name);
+        MFEM_VERIFY(index >= 0 && index < input_fields.Size(),
+                    "FieldCollection::InputField: Invalid index for field name: "
+                    << field_name << ".");
+        return input_fields[index];
+    }
+
+    Field* OutputField(int i) const { return output_fields[i]; }
+    Field *OutputField(const std::string &field_name) const
+    {
+        bool has_index = index_map.Has(field_name);
+        if(!has_index)
+        {
+            MFEM_WARNING("FieldCollection::OutputField: Field with name "
+                         << field_name << " does not exist in the collection.");
+            return nullptr;
+        }
+
+        int index = index_map.Get(field_name);
+        MFEM_VERIFY(index >= 0 && index < output_fields.Size(),
+                    "FieldCollection::OutputField: Invalid index for field name: "
+                    << field_name << ".");
+        return output_fields[index];
+    }
+
+    FieldMap &Fields() { return fields; }
+    FieldMap Fields() const { return fields; }
+
+    virtual void Save (std::ostream &out) const
+    {
+        out << "\"Fields\":\n";
+        out << "{\n";
+        for (auto f = fields.begin(); f != fields.end(); ++f)
+        {
+            std::string f_name = f->first;
+            Field *f_obj = f->second;
+            // out << "  " << f_name << ": ID " << f_obj->ID() << ",\n";
+            // out << f_obj->ID() << ": " << f_name << ",\n";
+            out << '\"' << f_obj->ID() << "\": \"" << f_name << "\"";
+            if(f != std::prev(fields.end())) out << ",";
+            out << "\n";
+        }
+        out << "},\n";
+
+        out << "\"Inputs\":\n";
+        out << "{\n";
+        for (int i = 0; i < input_fields.Size(); ++i)
+        {
+            Field *f_obj = input_fields[i];
+            out << '\"' << f_obj->ID() << "\": \"" << f_obj->Name() << "\"";
+            if(i != input_fields.Size() - 1) out << ",";
+            out << "\n";
+        }
+        out << "},\n";
+
+        out << "\"Outputs\":\n";
+        out << "{\n";
+        for (int i = 0; i < output_fields.Size(); ++i)
+        {
+            Field *f_obj = output_fields[i];
+            out << '\"' << f_obj->ID() << "\": \"" << f_obj->Name() << "\"";
+            if(i != output_fields.Size() - 1) out << ",";
+            out << "\n";
+        }
+        out << "}\n";
+    }
+
+    Field* HasField(const Field &field) const
+    {
+        for (auto f = fields.begin(); f != fields.end(); ++f)
+        {
+            if(f->second == &field)
+            {
+                return f->second;
+            }
+        }
+        return nullptr;
+    }
+
+    Field* HasField(const std::string &field_name) const
+    {
+        return fields.Get(field_name);
+    }
+
+    Field* HasField(const int id) const
+    {
+        for (auto f = fields.begin(); f != fields.end(); ++f)
+        {
+            if(f->second->ID() == id)
+            {
+                return f->second;
+            }
+        }
+        return nullptr;
+    }
+
+    ~FieldCollection(){}
+
+};
+
+
+class GraphNode : public Operator
+{
+public:
+    enum ExecutionMode
+    {
+        GRADIENT_MODE, ///< Node is being executed as part of a gradient evaluation
+        DEFAULT_MODE   ///< Node is being executed as default, operator evaluation
+    };
+
+private:
+    inline static int next_id = 0;
+
+protected:
+    int id = -1;
+    int node_index = -1;
+    mutable ExecutionMode exec_mode = DEFAULT_MODE;
+
+    std::string name;
+    mutable FieldCollection field_collection; ///< Collection of fields associated with this node
+
+    // Offsets to be used for operation on BlockVector
+    Array<int> input_offsets;  ///< Offsets for input fields
+    Array<int> output_offsets; ///< Offsets for output fields
+
+    int GetValidID(int id_, int lb=0, int ub = std::numeric_limits<int>::max())
+    {
+        return (id_ >= lb && id_ <= ub) ? id_ : next_id++;
+    }
+
+public:
+
+    GraphNode(int h, int w) : Operator(h,w), id(GetValidID(-1)),
+                              name("Node_" + std::to_string(id)),
+                              field_collection(this) { }
+
+    GraphNode(int s = 0) : GraphNode(s, s) { }
+
+    void SetNodeIndex(int index){ node_index = index; }
+    int GetNodeIndex() const { return node_index; }
+
+    void SetExecutionMode(ExecutionMode mode) { exec_mode = mode; }
+    ExecutionMode GetExecutionMode() const { return exec_mode; }
+
+    void SetName(const std::string &name_) { name = name_; }
+    std::string Name() const { return name; }
+
+    void SetID(int id_) { id = id_; }
+    int ID() const { return id; }
+
+    FieldCollection::FieldMap& Fields() { return field_collection.Fields(); }
+    Field* Fields(const std::string &f) { return field_collection.GetField(f); }
+
+    FieldCollection::FieldMap Fields() const { return field_collection.Fields(); }
+    Field* Fields(const std::string &f) const { return field_collection.GetField(f); }
+
+    Array<Field*>& InputFields() const { return field_collection.InputFields(); }
+    Array<Field*>& OutputFields() const { return field_collection.OutputFields(); }
+    Field* InputField(int i) const { return field_collection.InputField(i); }
+    Field* OutputField(int i) const { return field_collection.OutputField(i); }
+
+
+    virtual void AddInput(const std::string &field_name,
+                          Field *field, bool own = false)
+    { field_collection.AddInput(field_name, field, own); }
+
+    virtual void AddInput(Field *field, bool own = false)
+    { AddInput(field->Name(), field, own); }
+
+    template<bool OwnInputs = false,
+             typename... Args,
+             bool AreFields = std::conjunction<std::is_base_of<Field, std::remove_pointer_t<Args>> ...>::value,
+             typename std::enable_if<AreFields, bool>::type = true >
+    void AddInputs(Args... args)
+    {
+        ((AddInput(std::forward<Args>(args), OwnInputs)), ...);
+    }
+
+    virtual void AddOutput(const std::string &field_name,
+                           Field *field, bool own = false)
+    { field_collection.AddOutput(field_name, field, own); }
+
+    virtual void AddOutput(Field *field, bool own = false)
+    { AddOutput(field->Name(), field, own); }
+
+    template<bool OwnOutputs = false,
+             typename... Args,
+             bool AreFields = std::conjunction<std::is_base_of<Field, std::remove_pointer_t<Args>> ...>::value,
+             typename std::enable_if<AreFields, bool>::type = true >
+    void AddOutputs(Args... args)
+    {
+        ((AddOutput(std::forward<Args>(args), OwnOutputs)), ...);
+    }
+
+    virtual void Save (std::ostream &out) const
+    {
+        out << "\"Node-" << id << "\" : " << std::endl;
+        out << "{\n";
+        out << "\"Name\": \"" << name << "\",\n";
+        field_collection.Save(out);
+        out << "}";
+    }
+
+    virtual void Mult(const Vector &x, Vector &y) const override
+    {
+        MFEM_ABORT("GraphNode::Mult() not implemented");
+    }
+
+    virtual void Mult(const MultiVector &x, MultiVector &y) override
+    {
+        MFEM_ABORT("GraphNode::Mult(MultiVector) not implemented");
+    }
+
+    using Operator::GetGradient;
+
+    // TODO: Possibly remove this and only support MultiVector version of GradientMult
+    virtual void GradientMult(const Vector &x, const Vector &dx, Vector &dy) const
+    {
+        MFEM_ABORT("GraphNode::GradientMult() not implemented");
+        GetGradient(x).Mult(dx, dy);
+    }
+
+    virtual void GradientMult(const MultiVector &x, const MultiVector &dx, MultiVector &dy) const
+    {
+        MFEM_ABORT("GraphNode::GradientMult() not implemented");
+        GetGradient(x).Mult(dx, dy);
+    }
+
+    // TODO: Possibly remove this and only support MultiVector version of GradientMultTranspose
+    virtual void GradientMultTranspose(const Vector &x, const Vector &dx, Vector &dy) const
+    {
+        MFEM_ABORT("GraphNode::GradientMultTranspose() not implemented");
+        GetGradient(x).MultTranspose(dx, dy);
+    }
+
+    virtual void GradientMultTranspose(const MultiVector &x, const MultiVector &dx, MultiVector &dy) const
+    {
+        MFEM_ABORT("GraphNode::GradientMultTranspose() not implemented");
+        // GetGradient(x).MultTranspose(dx, dy); // Not yet implemented
+    }
+
+    /// @brief Return the input offsets for block starts.
+    Array<int>& InputOffsets() { return input_offsets; }
+
+    /// @brief Read only access to the input offsets for block starts.
+    const Array<int>& InputOffsets() const { return input_offsets; }
+
+    void SetInputOffsets(const Array<int> &offsets) { input_offsets = offsets; }
+
+    /// @brief Return the output offsets for block starts.
+    Array<int>& OutputOffsets() { return output_offsets; }
+
+    /// @brief Read only access to the output offsets for block starts.
+    const Array<int>& OutputOffsets() const { return output_offsets; }
+
+    void SetOutputOffsets(const Array<int> &offsets) { output_offsets = offsets; }
+
+    virtual ~GraphNode() = default;
+};
+
+
+/**
+   @brief An abstract, type-erased class to define the interface for 
+   operators, not inherited from @a GraphNode. It performs SFINAE 
+   checks for stored operator's member functions and override the Mult
+   to call the stored object's functions.
+ */
+template <typename OpType>
+class AbstractOperator : public GraphNode
+{
+protected:
+    /// Define a template class 'check' to test for the existence of member functions
+    template <typename C>
+    class CheckMember{
+        private:
+
+        /// @brief A type trait to check if the erased class has the function Mult
+        /// with the needed signatures.
+        template<class T>
+        using Mult = decltype(std::declval<T&>().Mult(std::declval<const Vector&>(),
+                                                      std::declval<Vector&>()));
+
+        template<class T>
+        using MultPtr = decltype(std::declval<T&>().Mult(std::declval<const int>(),
+                                                         std::declval<const real_t*>(),
+                                                         std::declval<const int>(),
+                                                         std::declval<real_t*>()));
+        // ---------------------------------------------------------------------
+        
+        template <typename T, template<typename> typename Func, typename R>
+        static constexpr auto Check(T*) -> typename std::is_same< Func<T>, R>::type;
+
+        template <typename, template<typename> typename, typename >
+        static constexpr std::false_type Check(...);
+
+        // --- Check for the existence of the member functions
+        typedef decltype(Check<C,Mult,void>(0)) Has_Mult;
+        typedef decltype(Check<C,MultPtr,void>(0)) Has_MultPtr;
+    public:
+        static constexpr bool HasMult  = Has_Mult::value;
+        static constexpr bool HasMultPtr  = Has_MultPtr::value;
+    };
+
+    OpType *op;  ///< Pointer to the operator
+
+public:
+
+    constexpr bool HasExecute(){return CheckMember<OpType>::HasStep;}
+    constexpr bool HasMult(){return CheckMember<OpType>::HasMult;}
+
+
+    /// @brief Constructor for the type-erased AbstractOperator class
+    AbstractOperator(OpType *op_, int h, int w) : GraphNode(h,w), op(op_)
+    { }
+
+    /// @brief Constructor for the type-erased AbstractOperator class.
+    AbstractOperator(OpType *op_, int s = 0) : AbstractOperator(op_,s,s) {}
+
+    /**
+       @brief Perform Mult operation with the stored operator, if it exists.
+     */
+    void Mult(const Vector &x, Vector &y) const override
+    {
+        if constexpr (CheckMember<OpType>::HasMult)
+        {
+            op->Mult(x,y);
+        }
+        else if constexpr (CheckMember<OpType>::HasMultPtr)
+        {
+            op->Mult(x.Size(), x.GetData(), y.Size(), y.GetData());
+        }
+        else
+        {
+            MFEM_ABORT("The AbstractOperator does not have the function, "
+                       "Mult(const Vector&, Vector&) or "
+                       "Mult(int, double*, int, double*).");
+        }
+    }
+};
+
+
+/**
+   @brief A class to store and coupled multiple operators together.
+ */
+class DAGraph : public GraphNode
+{
+public:
+
+    using IntToIntMap = GenericFieldMap<int, int>;
+    using IntToFieldMap = GenericFieldMap<int, Field*>;
+
+    enum class GradMode
+    {
+        FINITE_DIFF = 0,  ///< Finite difference Jacobian
+        MATRIX_FREE = 1,  ///< Matrix-free Jacobian
+        ASSEMBLED = 2,    ///< Assembled Jacobian
+        NONE = 3          ///< Not implemented
+    };
+
+    enum InputType
+    {
+        VECTOR,      ///< Asemble the input blockvector from individual fields
+        MULTIVECTOR, ///< Asemble the multivector from individual fields
+        NONE         ///< No input
+    };
+
+protected:
+    Array<GraphNode*> nodes; ///< Vector of individual operators
+    Array<bool> node_owned; ///< Whether the operators are owned
+    Array<int> node_depth; ///< Depth of each operator in the graph
+
+    int max_width  = 0;     ///< Largest operator width
+    int max_height = 0;     ///< Largest operator height
+    int nnodes     = 0;     ///< The number of nodes
+    bool sorted    = false; ///< True if the nodes are topologically sorted
+    bool assembled = false; ///< True if the graph is assembled
+
+    GradMode grad_mode = GradMode::MATRIX_FREE; ///< Gradient mode for the graph
+    mutable Operator *grad = nullptr; ///< Gradient operator
+
+    InputType input_type = InputType::MULTIVECTOR; ///< Input type for the graph
+    mutable Vector x_node, y_node; ///< Temporary vectors for evaluating nodes
+    mutable MultiVector xmv_node, ymv_node; ///< Temporary multivectors for evaluating nodes
+
+    IntToFieldMap fid_to_field; ///< Map from Field ID to Field pointer
+    IntToIntMap fid_to_index; ///< Map from ID to index in an array; needed since ordering is not unique
+
+    friend class GraphGradient;
+
+public:
+    /**
+       @brief Construct a new CoupledOperator object.
+       @param nop Total number of operators to couple
+     */
+    DAGraph(const int nop) : GraphNode()
+    {
+        nodes.Reserve(nop);
+        node_owned.Reserve(nop);
+    }
+
+    /**
+       @brief Construct a new CoupledOperator object for an 
+       abstract non/mfem operator.
+     */
+    template <class OpType>
+    DAGraph(const OpType &op) : DAGraph(1)
+    {
+        AddOperator(op);
+    }
+
+    /**
+       @brief Add an operator to the list of coupled operator and
+       return pointer to it. Not owned unless it's not derived from GraphNode.
+     */
+    template <class OpType>
+    GraphNode* AddOperator(OpType *op_, int h, int w)
+    {
+        // Add operator to list of operators
+        if constexpr(std::is_base_of<GraphNode, OpType>::value)
+        {
+            nodes.push_back(op_);
+            node_owned.Append(false);
+        } 
+        else
+        {
+            nodes.push_back(new AbstractOperator<OpType>(op_,h,w));
+            node_owned.Append(true);
+        }
+        nnodes++;
+
+        // Update size of the coupled operator and the block offsets
+        GraphNode* op = nodes.Last();
+        op->SetNodeIndex(nnodes-1); // Set the index of the operator
+
+        int ht = op->Height();
+        int wt = op->Width();
+
+        max_width = std::max(max_width, wt);
+        max_height = std::max(max_height, ht);
+        sorted = false;
+
+        return op;
+    }
+
+    /// @brief Add an operator to the list of coupled operator and return pointer to it.
+    template <class OpType>
+    GraphNode* AddOperator(OpType *op_, int s = 0) { return AddOperator(op_,s,s);}
+
+    /// @brief Get the number of coupled operators
+    int Size(){return nnodes;}
+
+    /// @brief Get the size of the largest operator
+    int MaxWidth() const {return max_width;}
+    int MaxHeight() const {return max_height;}
+
+    IntToIntMap &GetFieldIdToIndexMap() { return fid_to_index; }
+    IntToIntMap GetFieldIdToIndexMap() const { return fid_to_index; }
+
+    IntToFieldMap &GetFieldIdToFieldMap() { return fid_to_field; }
+    IntToFieldMap GetFieldIdToFieldMap() const { return fid_to_field; }
+
+    /// @brief Get the operator at index @a i
+    GraphNode* GetNode(const int i)
+    {
+        MFEM_ASSERT(i >= 0 && i < nnodes,
+               "index [" << i << "] is out of range [0," << nnodes << ")");
+        return nodes[i];
+    }
+
+    Array<GraphNode*>& Nodes() { return nodes; }
+
+    /// @brief Specify whether the operator at index @a i is owned.
+    void OwnNode(const int i, bool own = true)
+    {
+        MFEM_ASSERT(i >= 0 && i < nnodes,
+               "index [" << i << "] is out of range [0," << nnodes << ")");
+        node_owned[i] = own;
+    }
+
+    void Assemble();
+    bool IsAssembled() const { return assembled; }
+
+    void TopologicalSort();
+    bool IsSorted() const { return sorted; }
+
+    void ComputeDepth();
+
+    void ValidateOffsets();
+
+    void ValidateNode(GraphNode &node);
+
+    void CollectFieldMaps();
+
+    using GraphNode::AddInput;
+    void AddInput(Field *field, int sz, bool own = false)
+    {
+        if(input_offsets.Size() == 0)
+        {   // First entry
+            input_offsets.Append(0);
+        }
+        input_offsets.Append(input_offsets.Last() + sz);
+        AddInput(field, own);
+    }
+
+    using GraphNode::AddOutput;
+    void AddOutput(Field *field, int sz, bool own = false)
+    {
+        if(output_offsets.Size() == 0)
+        {   // First entry
+            output_offsets.Append(0);
+        }
+        output_offsets.Append(output_offsets.Last() + sz);
+        AddOutput(field, own);
+    }
+
+    /// @brief Set the gradient mode for the coupled operator
+    void SetGradientMode(GradMode mode)
+    {
+        if(mode != grad_mode)
+        {
+            if(grad) { delete grad; grad = nullptr; }
+            grad_mode = mode;
+        }
+    }
+
+    void SetInputType(InputType type) { input_type = type; }
+    InputType GetInputType() const { return input_type; }
+
+    /**
+       @brief Apply the operator to the vector @a x 
+       and return the result in @a y.
+     */
+    virtual void Mult(const Vector &x, Vector &y) const override;
+
+    virtual void Mult(const MultiVector &x, MultiVector &y) override;
+
+    virtual void Execute(const MultiVector &x, MultiVector &y) const;
+
+    virtual void Save (std::ostream &out) const
+    {
+        out << "\"DAGraph\":\n";
+        out << "{\n";
+        // out << "\"nodes\" : " << nnodes << ",\n";
+        out << "\"Nodes\":\n";
+        out << "{\n";
+        for (int i = 0; i < nodes.Size(); i++)
+        {
+            nodes[i]->Save(out);
+            if(i != nodes.Size()-1) out << ",";
+            out << "\n";
+        }
+        out << "},\n"; // End of Nodes
+        field_collection.Save(out);
+        out << "}\n";
+    }
+
+    Operator& GetGradient(const Vector &x) const override;
+
+    /// @brief Destroy the Coupled Application object
+    ~DAGraph();
+};
+
+
+
+class GraphGradient : public Operator
+{
+public:
+    using InputType = DAGraph::InputType;
+
+protected:
+    mutable DAGraph *graph = nullptr; ///< Pointer to the DAGraph for which this is the gradient operator
+    Array<Vector*> x_work; ///< Array to store linearization point (intermediate fields)
+    mutable MultiVector xlin;
+    mutable Vector x0, dx, dy;
+    mutable MultiVector x0_mv, dx_mv, dy_mv;
+
+public:
+    GraphGradient(DAGraph &dag);
+
+    void Update(const Vector &x);
+
+    void Mult(const Vector &x, Vector &y) const override;
+
+    void Mult(const MultiVector &x, MultiVector &y) override;
+
+    void MultTranspose(const Vector &x, Vector &y) const override;
+
+    void MultTranspose(const MultiVector &x, MultiVector &y);
+
+    Operator &GetGradient(const Vector &x) const override;
+
+    void Forward(const MultiVector &x, MultiVector &y) const;
+
+    void Reverse(const MultiVector &x, MultiVector &y) const;
+
+    ~GraphGradient()
+    {
+        for (auto &v : x_work)
+        {
+            if(v) { delete v; v = nullptr; }
+        }
+        x_work.DeleteAll();
+    }
+};
+
+} //mfem namespace
+
+#endif