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."

CHANGELOG

@@ -68,10 +68,6 @@ Linear and nonlinear solvers
   PRefinement multigrid methods for problems posed on trace spaces (see e.g. the
   DPG miniapps).
 
-- Added interface to MUMPS direct solver for complex-valued problems.
-  Its usage is demonstrated in ex25p. See http://mumps.enseeiht.fr/ for more details.
-  Supported versions >= 5.1.1.
-
 GPU computing
 -------------
 - Added device assembly support for 3D H(curl) VectorFEDomainLFIntegrator.

CMakeLists.txt

@@ -76,9 +76,7 @@ set(XSDK_ENABLE_Fortran OFF)
 # Check if we need to enable C or Fortran.
 if (MFEM_USE_CONDUIT OR
     MFEM_USE_SIDRE OR
-    MFEM_USE_PETSC OR
-    MFEM_USE_MUMPS OR
-    MFEM_USE_COMPLEX_MUMPS)
+    MFEM_USE_PETSC)
    # This seems to be needed by:
    #  * find_package(BLAS REQUIRED) and
    #  * find_package(HDF5 REQUIRED) needed, in turn, by:
@@ -86,7 +84,7 @@ if (MFEM_USE_CONDUIT OR
    #  * find_package(PETSc REQUIRED)
    set(XSDK_ENABLE_C ON)
 endif()
-if (MFEM_USE_STRUMPACK OR MFEM_USE_MUMPS OR MFEM_USE_COMPLEX_MUMPS)
+if (MFEM_USE_STRUMPACK OR MFEM_USE_MUMPS)
   # Just needed to find the MPI_Fortran libraries to link with
   set(XSDK_ENABLE_Fortran ON)
 endif()
@@ -344,7 +342,7 @@ if (MFEM_USE_MPI)
     endif()
   endif()
 else()
-  set(PKGS_NEED_MPI SUPERLU MUMPS COMPLEX_MUMPS PETSC SLEPC STRUMPACK PUMI)
+  set(PKGS_NEED_MPI SUPERLU MUMPS PETSC SLEPC STRUMPACK PUMI)
   foreach(PKG IN LISTS PKGS_NEED_MPI)
     if (MFEM_USE_${PKG})
       message(STATUS "Disabling package ${PKG} - requires MPI")
@@ -417,7 +415,7 @@ if (MFEM_USE_SUPERLU)
 endif()
 
 # MUMPS can only be enabled in parallel
-if (MFEM_USE_MUMPS OR MFEM_USE_COMPLEX_MUMPS)
+if (MFEM_USE_MUMPS)
   if (MFEM_USE_MPI)
     find_package(MUMPS REQUIRED mumps_common pord)
     set(MFEM_MUMPS_VERSION ${MUMPS_VERSION})

config/cmake/MFEMConfig.cmake.in

@@ -34,7 +34,6 @@ set(MFEM_USE_SUNDIALS @MFEM_USE_SUNDIALS@)
 set(MFEM_USE_SUITESPARSE @MFEM_USE_SUITESPARSE@)
 set(MFEM_USE_SUPERLU @MFEM_USE_SUPERLU@)
 set(MFEM_USE_MUMPS @MFEM_USE_MUMPS@)
-set(MFEM_USE_COMPLEX_MUMPS @MFEM_USE_COMPLEX_MUMPS@)
 set(MFEM_USE_STRUMPACK @MFEM_USE_STRUMPACK@)
 set(MFEM_USE_CUDSS @MFEM_USE_CUDSS@)
 set(MFEM_USE_GINKGO @MFEM_USE_GINKGO@)

config/cmake/config.hpp.in

@@ -103,7 +103,6 @@
 
 // Enable MFEM functionality based on the MUMPS library.
 #cmakedefine MFEM_USE_MUMPS
-#cmakedefine MFEM_USE_COMPLEX_MUMPS
 #cmakedefine MFEM_MUMPS_VERSION @MFEM_MUMPS_VERSION@
 
 // Enable MFEM functionality based on the STRUMPACK library.

config/cmake/modules/FindMUMPS.cmake

@@ -17,31 +17,13 @@
 
 include(MfemCmakeUtilities)
 
-# Decide headers/libs by MFEM precision
+# Toggle which precision of MUMPS to use depending on the precision of MFEM.
 if (MFEM_USE_DOUBLE)
-  set(_rmumps_header dmumps_c.h)
-  set(_rmumps_lib    dmumps)
-  set(_cmumps_header zmumps_c.h)
-  set(_cmumps_lib    zmumps)
-elseif (MFEM_USE_SINGLE)
-  set(_rmumps_header smumps_c.h)
-  set(_rmumps_lib    smumps)
-  set(_cmumps_header cmumps_c.h)
-  set(_cmumps_lib    cmumps)
-endif()
-
-# Choose which header/lib mfem_find_package should use as the "primary" one.
-# If both enabled, prefer the real one as primary (either is fine).
-if (MFEM_USE_MUMPS)
-  set(_mumps_header ${_rmumps_header})
-  set(_mumps_lib    ${_rmumps_lib})
-elseif (MFEM_USE_COMPLEX_MUMPS)
-  set(_mumps_header ${_cmumps_header})
-  set(_mumps_lib    ${_cmumps_lib})
-else()
-  # Should not happen in practice because FindMUMPS is only called when enabled,
-  set(_mumps_header ${_rmumps_header})
-  set(_mumps_lib    ${_rmumps_lib})
+  set(_mumps_header dmumps_c.h)
+  set(_mumps_lib dmumps)
+elseif(MFEM_USE_SINGLE)
+  set(_mumps_header smumps_c.h)
+  set(_mumps_lib smumps)
 endif()
 
 mfem_find_package(MUMPS MUMPS MUMPS_DIR
@@ -49,35 +31,8 @@ mfem_find_package(MUMPS MUMPS MUMPS_DIR
   "Paths to headers required by MUMPS."
   "Libraries required by MUMPS."
   ADD_COMPONENT mumps_common "include" ${_mumps_header} "lib" mumps_common
-  ADD_COMPONENT pord        "include" ${_mumps_header} "lib" pord)
+  ADD_COMPONENT pord "include" ${_mumps_header} "lib" pord)
 
-# If BOTH real and complex are enabled, ensure BOTH solver libs are linked.
-if (MUMPS_FOUND AND MFEM_USE_MUMPS AND MFEM_USE_COMPLEX_MUMPS)
-  # Find the "other" solver library and append it.
-  find_library(_mfem_other_mumps_solver
-    NAMES ${_cmumps_lib}   
-    HINTS ${MUMPS_DIR}
-    PATH_SUFFIXES lib lib64
-    NO_DEFAULT_PATH)
-
-  if (NOT _mfem_other_mumps_solver)
-    # Fall back to system search
-    find_library(_mfem_other_mumps_solver NAMES ${_cmumps_lib})
-  endif()
-
-  if (NOT _mfem_other_mumps_solver)
-    message(FATAL_ERROR
-      "MFEM_USE_MUMPS=ON and MFEM_USE_COMPLEX_MUMPS=ON, but could not find "
-      "the complex solver library '${_cmumps_lib}' in MUMPS_DIR='${MUMPS_DIR}'.")
-  endif()
-
-  # Put solver libs first (important for static link order)
-  # MUMPS_LIBRARIES contains the primary solver already + common + pord.
-  # We prepend the other solver.
-  list(INSERT MUMPS_LIBRARIES 0 ${_mfem_other_mumps_solver})
-endif()
-
-# Version detection 
 if (MUMPS_FOUND AND (NOT MUMPS_VERSION))
   try_run(MUMPS_VERSION_RUN_RESULT MUMPS_VERSION_COMPILE_RESULT
           ${CMAKE_CURRENT_BINARY_DIR}/config

config/cmake/modules/MfemCmakeUtilities.cmake

@@ -876,12 +876,12 @@ function(mfem_export_mk_files)
       MFEM_USE_ZLIB MFEM_USE_LIBUNWIND MFEM_USE_LAPACK MFEM_THREAD_SAFE
       MFEM_USE_LEGACY_OPENMP MFEM_USE_OPENMP MFEM_USE_MEMALLOC MFEM_USE_SUNDIALS
       MFEM_USE_SUITESPARSE MFEM_USE_SUPERLU MFEM_USE_SUPERLU5 MFEM_USE_MUMPS
-      MFEM_USE_COMPLEX_MUMPS MFEM_USE_STRUMPACK MFEM_USE_GINKGO MFEM_USE_AMGX
-      MFEM_USE_MAGMA MFEM_USE_GNUTLS MFEM_USE_NETCDF MFEM_USE_PETSC
-      MFEM_USE_SLEPC MFEM_USE_MPFR MFEM_USE_SIDRE MFEM_USE_FMS MFEM_USE_CONDUIT
-      MFEM_USE_PUMI MFEM_USE_HIOP MFEM_USE_GSLIB MFEM_USE_CUDA MFEM_USE_HIP
-      MFEM_USE_RAJA MFEM_USE_OCCA MFEM_USE_CEED MFEM_USE_CALIPER MFEM_USE_UMPIRE
-      MFEM_USE_SIMD MFEM_USE_ADIOS2 MFEM_USE_MKL_CPARDISO MFEM_USE_MKL_PARDISO
+      MFEM_USE_STRUMPACK MFEM_USE_GINKGO MFEM_USE_AMGX MFEM_USE_MAGMA
+      MFEM_USE_GNUTLS MFEM_USE_NETCDF MFEM_USE_PETSC MFEM_USE_SLEPC
+      MFEM_USE_MPFR MFEM_USE_SIDRE MFEM_USE_FMS MFEM_USE_CONDUIT MFEM_USE_PUMI
+      MFEM_USE_HIOP MFEM_USE_GSLIB MFEM_USE_CUDA MFEM_USE_HIP MFEM_USE_RAJA
+      MFEM_USE_OCCA MFEM_USE_CEED MFEM_USE_CALIPER MFEM_USE_UMPIRE MFEM_USE_SIMD
+      MFEM_USE_ADIOS2 MFEM_USE_MKL_CPARDISO MFEM_USE_MKL_PARDISO
       MFEM_USE_ADFORWARD MFEM_USE_CODIPACK MFEM_USE_BENCHMARK MFEM_USE_PARELAG
       MFEM_USE_TRIBOL MFEM_USE_MOONOLITH MFEM_USE_ALGOIM MFEM_USE_ENZYME
       MFEM_USE_HDF5)

config/config.hpp

@@ -140,9 +140,6 @@ constexpr real_t operator""_r(unsigned long long v)
 #ifdef MFEM_USE_MUMPS
 #error Building with MUMPS (MFEM_USE_MUMPS=YES) requires MPI (MFEM_USE_MPI=YES)
 #endif
-#ifdef MFEM_USE_COMPLEX_MUMPS
-#error Building with COMPLEX_MUMPS (MFEM_USE_COMPLEX_MUMPS=YES) requires MPI (MFEM_USE_MPI=YES)
-#endif
 #ifdef MFEM_USE_STRUMPACK
 #error Building with STRUMPACK (MFEM_USE_STRUMPACK=YES) requires MPI (MFEM_USE_MPI=YES)
 #endif

config/config.hpp.in

@@ -103,7 +103,6 @@
 
 // Enable MFEM functionality based on the MUMPS library.
 // #define MFEM_USE_MUMPS
-// #define MFEM_USE_COMPLEX_MUMPS
 // #define MFEM_MUMPS_VERSION @MFEM_MUMPS_VERSION@
 
 // Enable MFEM functionality based on the STRUMPACK library.

config/defaults.cmake

@@ -37,7 +37,6 @@ option(MFEM_USE_SUITESPARSE "Enable SuiteSparse usage" OFF)
 option(MFEM_USE_SUPERLU "Enable SuperLU_DIST usage" OFF)
 option(MFEM_USE_SUPERLU5 "Use the old SuperLU_DIST 5.1 version" OFF)
 option(MFEM_USE_MUMPS "Enable MUMPS usage" OFF)
-option(MFEM_USE_COMPLEX_MUMPS "Enable COMPLEX_MUMPS usage" OFF)
 option(MFEM_USE_STRUMPACK "Enable STRUMPACK usage" OFF)
 option(MFEM_USE_CUDSS "Enable cuDSS usage" OFF)
 option(MFEM_USE_GINKGO "Enable Ginkgo usage" OFF)
@@ -154,7 +153,7 @@ set(SuperLUDist_REQUIRED_PACKAGES "MPI" "ParMETIS" "METIS"
 set(MUMPS_DIR "${MFEM_DIR}/../MUMPS_5.5.0" CACHE PATH
     "Path to the MUMPS library.")
 # MUMPS may also depend on "OpenMP", depending on how it was compiled.
-set(MUMPS_REQUIRED_PACKAGES "MPI" "MPI_Fortran" "METIS"
+set(MUMPS_REQUIRED_PACKAGES "MPI" "MPI_Fortran" "ParMETIS" "METIS"
     "ScaLAPACK" "LAPACK" "BLAS" CACHE STRING
     "Additional packages required by MUMPS.")
 # If the MPI package does not find all required Fortran libraries:

config/defaults.mk

@@ -152,7 +152,6 @@ MFEM_USE_SUITESPARSE   = NO
 MFEM_USE_SUPERLU       = NO
 MFEM_USE_SUPERLU5      = NO
 MFEM_USE_MUMPS         = NO
-MFEM_USE_COMPLEX_MUMPS = NO
 MFEM_USE_STRUMPACK     = NO
 MFEM_USE_CUDSS         = NO
 MFEM_USE_GINKGO        = NO
@@ -250,7 +249,7 @@ ifeq (YES,$(MFEM_USE_HIP))
 endif
 
 # METIS library configuration
-ifeq ($(MFEM_USE_SUPERLU)$(MFEM_USE_STRUMPACK)$(MFEM_USE_MUMPS)$(MFEM_USE_COMPLEX_MUMPS),NONONONO)
+ifeq ($(MFEM_USE_SUPERLU)$(MFEM_USE_STRUMPACK)$(MFEM_USE_MUMPS),NONONO)
    ifeq ($(MFEM_USE_METIS_5),NO)
      METIS_DIR = @MFEM_DIR@/../metis-4.0
      METIS_OPT =
@@ -354,19 +353,13 @@ MPI_FORTRAN_LIB = -lmpifort
 # MUMPS library configuration
 MUMPS_DIR = @MFEM_DIR@/../MUMPS_5.5.0
 MUMPS_OPT = -I$(MUMPS_DIR)/include
-COMPLEX_MUMPS_OPT = $(MUMPS_OPT)
-
-MUMPS_COMMON_PATH = $(XLINKER)-rpath,$(MUMPS_DIR)/lib -L$(MUMPS_DIR)/lib
-MUMPS_COMMON_LIB = -lmumps_common -lpord $(SCALAPACK_LIB) $(LAPACK_LIB) $(MPI_FORTRAN_LIB)
+MUMPS_LIB = $(XLINKER)-rpath,$(MUMPS_DIR)/lib -L$(MUMPS_DIR)/lib
 ifeq ($(MFEM_USE_SINGLE),YES)
-  MUMPS_SOLVER_LIB         = -lsmumps
-  COMPLEX_MUMPS_SOLVER_LIB = -lcmumps
+   MUMPS_LIB += -lsmumps
 else
-  MUMPS_SOLVER_LIB         = -ldmumps
-  COMPLEX_MUMPS_SOLVER_LIB = -lzmumps
+   MUMPS_LIB += -ldmumps
 endif
-MUMPS_LIB         = $(MUMPS_COMMON_PATH) $(MUMPS_SOLVER_LIB) $(MUMPS_COMMON_LIB)
-COMPLEX_MUMPS_LIB = $(MUMPS_COMMON_PATH) $(COMPLEX_MUMPS_SOLVER_LIB) $(MUMPS_COMMON_LIB)
+MUMPS_LIB += -lmumps_common -lpord $(SCALAPACK_LIB) $(LAPACK_LIB) $(MPI_FORTRAN_LIB)
 
 # STRUMPACK library configuration
 STRUMPACK_DIR = @MFEM_DIR@/../STRUMPACK-build

config/get_mumps_version.cpp

@@ -10,18 +10,10 @@
 // CONTRIBUTING.md for details.
 
 #ifdef MFEM_USE_SINGLE
-#ifdef MFEM_USE_COMPLEX_MUMPS
-#include "cmumps_c.h"
-#else
 #include "smumps_c.h"
-#endif
-#else
-#ifdef MFEM_USE_COMPLEX_MUMPS
-#include "zmumps_c.h"
 #else
 #include "dmumps_c.h"
 #endif
-#endif
 #include <string>
 #include <iostream>
 #include <algorithm>

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.

config/makefile

@@ -42,7 +42,7 @@ GHV_FLAGS = $(MFEM_CXXFLAGS) $(subst @MFEM_DIR@,$(if $(MFEM_DIR),$(MFEM_DIR),..)
 SMX = $(if $(MFEM_USE_PUMI:NO=),MFEM_USE_SIMMETRIX)
 SMX_PATH = $(PUMI_DIR)/include/gmi_sim.h
 SMX_FILE = $(subst @MFEM_DIR@,$(if $(MFEM_DIR),$(MFEM_DIR),..),$(SMX_PATH))
-MUMPS = $(if $(filter YES,$(MFEM_USE_MUMPS) $(MFEM_USE_COMPLEX_MUMPS)),YES,)
+MUMPS = $(MFEM_USE_MUMPS:NO=)
 GMV_CXX ?= $(MFEM_CXX)
 GMV = get_mumps_version
 GMV_FLAGS = $(MFEM_CXXFLAGS) $(subst @MFEM_DIR@,$(if $(MFEM_DIR),$(MFEM_DIR),..),$(MUMPS_OPT))

examples/ex25p.cpp

@@ -200,7 +200,7 @@ int main(int argc, char *argv[])
    args.AddOption(&slu_solver, "-slu", "--superlu", "-no-slu",
                   "--no-superlu", "Use the SuperLU Solver.");
 #endif
-#if defined(MFEM_USE_MUMPS) || defined(MFEM_USE_COMPLEX_MUMPS)
+#ifdef MFEM_USE_MUMPS
    args.AddOption(&mumps_solver, "-mumps", "--mumps-solver", "-no-mumps",
                   "--no-mumps-solver", "Use the MUMPS Solver.");
 #endif
@@ -502,24 +502,16 @@ int main(int argc, char *argv[])
       delete A;
    }
 #endif
-#if defined(MFEM_USE_MUMPS) || defined(MFEM_USE_COMPLEX_MUMPS)
+#ifdef MFEM_USE_MUMPS
    if (!pa && mumps_solver)
    {
-      ComplexHypreParMatrix *Ahc = Ah.As<ComplexHypreParMatrix>();
-#ifdef MFEM_USE_COMPLEX_MUMPS
-      ComplexMUMPSSolver mumps(MPI_COMM_WORLD);
-      mumps.SetPrintLevel(0);
-      mumps.SetOperator(*Ahc);
-      mumps.Mult(B, X);
-#else
-      HypreParMatrix *A = Ahc->GetSystemMatrix();
+      HypreParMatrix *A = Ah.As<ComplexHypreParMatrix>()->GetSystemMatrix();
       MUMPSSolver mumps(A->GetComm());
       mumps.SetPrintLevel(0);
       mumps.SetMatrixSymType(MUMPSSolver::MatType::UNSYMMETRIC);
       mumps.SetOperator(*A);
       mumps.Mult(B, X);
       delete A;
-#endif
    }
 #endif
    // 16a. Set up the parallel Bilinear form a(.,.) for the preconditioner

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. */

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

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"
@@ -65,7 +66,7 @@
 #include "hypre_parcsr.hpp"
 #include "hypre.hpp"
 
-#if defined(MFEM_USE_MUMPS) || defined(MFEM_USE_COMPLEX_MUMPS)
+#ifdef MFEM_USE_MUMPS
 #include "mumps.hpp"
 #endif
 

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/mumps.cpp

@@ -12,9 +12,11 @@
 #include "../config/config.hpp"
 #include "../general/communication.hpp"
 
-#if defined(MFEM_USE_MUMPS) || defined(MFEM_USE_COMPLEX_MUMPS)
+#ifdef MFEM_USE_MUMPS
+#ifdef MFEM_USE_MPI
+
 #include "mumps.hpp"
-#include <unordered_map>
+
 #include <algorithm>
 
 #if MFEM_MUMPS_VERSION >= 530
@@ -33,14 +35,9 @@
 #define MUMPS_INFO(I) info[(I) -1]
 #define MUMPS_INFOG(I) infog[(I) -1]
 
-#endif // MFEM_USE_MUMPS || MFEM_USE_COMPLEX_MUMPS
-
-
 namespace mfem
 {
 
-#ifdef MFEM_USE_MUMPS
-
 MUMPSSolver::MUMPSSolver(MPI_Comm comm_)
 {
    Init(comm_);
@@ -670,760 +667,7 @@ void MUMPSSolver::RedistributeSol(const int *rmap, const real_t *x,
 }
 #endif
 
-#endif // MFEM_USE_MUMPS
-
-
-#ifdef MFEM_USE_COMPLEX_MUMPS
-
-ComplexMUMPSSolver::ComplexMUMPSSolver(MPI_Comm comm_)
-{
-   Init(comm_);
-}
-
-ComplexMUMPSSolver::ComplexMUMPSSolver(const Operator &op)
-{
-   auto APtr = dynamic_cast<const ComplexHypreParMatrix *>(&op);
-   MFEM_VERIFY(APtr,
-               "ComplexMUMPSSolver requires a ComplexHypreParMatrix operator");
-   SetOperator(op);
-}
-
-void ComplexMUMPSSolver::Init(MPI_Comm comm_)
-{
-   comm = comm_;
-   MPI_Comm_size(comm, &numProcs);
-   MPI_Comm_rank(comm, &myid);
-
-   print_level = 2;
-   row_start = 0;
-
-   id = nullptr;
-
-#if MFEM_MUMPS_VERSION >= 530
-   irhs_loc = nullptr;
-   isol_loc = nullptr;
-   rhs_loc  = nullptr;
-   sol_loc  = nullptr;
-#else
-   global_num_rows = 0;
-   recv_counts = nullptr;
-   displs = nullptr;
-   rhs_glob = nullptr;
-   rhs_glob_r = nullptr;
-   rhs_glob_i = nullptr;
-#endif
-}
-
-ComplexMUMPSSolver::~ComplexMUMPSSolver()
-{
-#if MFEM_MUMPS_VERSION >= 530
-   delete [] irhs_loc;
-   delete [] isol_loc;
-   delete [] rhs_loc;
-   delete [] sol_loc;
-#else
-   delete [] recv_counts;
-   delete [] displs;
-   delete [] rhs_glob;
-   delete [] rhs_glob_r;
-   delete [] rhs_glob_i;
-#endif
-
-   if (id)
-   {
-      id->job = -2;
-      mumps_call();
-      delete id;
-      id = nullptr;
-   }
-}
-
-void ComplexMUMPSSolver::SetOperator(const Operator &op)
-{
-   auto APtr = dynamic_cast<const ComplexHypreParMatrix *>(&op);
-   MFEM_VERIFY(APtr,
-               "ComplexMUMPSSolver requires a ComplexHypreParMatrix operator");
-
-   height = op.Height();
-   width  = op.Width();
-
-   const HypreParMatrix *Ar = (APtr->hasRealPart()) ? &APtr->real() : nullptr;
-   const HypreParMatrix *Ai = (APtr->hasImagPart()) ? &APtr->imag() : nullptr;
-
-   MFEM_VERIFY(Ar || Ai, "ComplexMUMPSSolver: both real and imag parts are null.");
-
-   // Pick communicator from the non-null part
-   MPI_Comm op_comm = (Ar ? Ar->GetComm() : Ai->GetComm());
-
-   // Comm setup/check
-   if (comm == MPI_COMM_NULL) { Init(op_comm); }
-   else
-   {
-      int cmp = MPI_UNEQUAL;
-      MPI_Comm_compare(comm, op_comm, &cmp);
-      MFEM_VERIFY(cmp != MPI_UNEQUAL, "MPI Comm mismatch");
-   }
-
-   // HostRead only if non-null
-   if (Ar) { Ar->HostRead(); }
-   if (Ai) { Ai->HostRead(); }
-
-   // hypre parcsr pointers
-   hypre_ParCSRMatrix *parcsr_op_r = nullptr;
-   hypre_ParCSRMatrix *parcsr_op_i = nullptr;
-
-   if (Ar) { parcsr_op_r = (hypre_ParCSRMatrix*) const_cast<HypreParMatrix&>(*Ar); }
-   if (Ai) { parcsr_op_i = (hypre_ParCSRMatrix*) const_cast<HypreParMatrix&>(*Ai); }
-
-   // Merge diag+offd for whichever exists
-   hypre_CSRMatrix *csr_op_r = nullptr;
-   hypre_CSRMatrix *csr_op_i = nullptr;
-
-   if (parcsr_op_r) { csr_op_r = hypre_MergeDiagAndOffd(parcsr_op_r); }
-   if (parcsr_op_i) { csr_op_i = hypre_MergeDiagAndOffd(parcsr_op_i); }
-
-#if MFEM_HYPRE_VERSION >= 21600
-   if (csr_op_r) { hypre_CSRMatrixBigJtoJ(csr_op_r); }
-   if (csr_op_i) { hypre_CSRMatrixBigJtoJ(csr_op_i); }
-#endif
-
-   // Determine local/global sizes and row_start from an existing part
-   const int n_loc = internal::to_int((csr_op_r ? csr_op_r->num_rows :
-                                       csr_op_i->num_rows));
-   row_start = internal::to_int((parcsr_op_r ? parcsr_op_r->first_row_index
-                                 : parcsr_op_i->first_row_index));
-   const int global_n = internal::to_int((parcsr_op_r ?
-                                          parcsr_op_r->global_num_rows
-                                          : parcsr_op_i->global_num_rows));
-
-   // Use nullptr checks
-   const int *Ir = csr_op_r ? csr_op_r->i : nullptr;
-   const int *Jr = csr_op_r ? csr_op_r->j : nullptr;
-   const real_t *Vr = csr_op_r ? (const real_t*)csr_op_r->data : nullptr;
-
-   const int *Ii = csr_op_i ? csr_op_i->i : nullptr;
-   const int *Ji = csr_op_i ? csr_op_i->j : nullptr;
-   const real_t *Vi = csr_op_i ? (const real_t*)csr_op_i->data : nullptr;
-
-   // Build union COO
-   std::vector<int> Icoo, Jcoo;
-   std::vector<mumps_complex_t> Zcoo;
-
-   size_t nnz_r = csr_op_r ? (size_t)csr_op_r->num_nonzeros : 0;
-   size_t nnz_i = csr_op_i ? (size_t)csr_op_i->num_nonzeros : 0;
-   Icoo.reserve(nnz_r + nnz_i);
-   Jcoo.reserve(nnz_r + nnz_i);
-   Zcoo.reserve(nnz_r + nnz_i);
-
-   BuildUnionCOO(n_loc, row_start, Ir, Jr, Vr, Ii, Ji, Vi, Icoo, Jcoo, Zcoo);
-
-   const int nnz = (int)Icoo.size();
-   int *I = new int[nnz];
-   int *J = new int[nnz];
-   mumps_complex_t *A = new mumps_complex_t[nnz];
-
-   std::copy(Icoo.begin(), Icoo.end(), I);
-   std::copy(Jcoo.begin(), Jcoo.end(), J);
-   std::copy(Zcoo.begin(), Zcoo.end(), A);
-
-   // New ComplexMUMPS object or reuse an existing one
-   if (!id || !reorder_reuse)
-   {
-      if (id)
-      {
-         id->job = -2;
-         mumps_call();
-         delete id;
-         id = nullptr;
-      }
-
-#ifdef MFEM_USE_SINGLE
-      id = new CMUMPS_STRUC_C();
-#else
-      id = new ZMUMPS_STRUC_C();
-#endif
-
-      id->sym = 0; // general complex
-      id->par = 1;
-      id->comm_fortran = (MUMPS_INT)MPI_Comm_c2f(comm);
-
-      // Init
-      id->job = -1;
-      mumps_call();
-
-      // Set parameters
-      SetParameters();
-
-      // Attach matrix
-      id->n       = global_n;
-      id->nnz_loc = nnz;
-      id->irn_loc = I;
-      id->jcn_loc = J;
-      id->a_loc   = A;
-
-      // Analysis (ordering + symbolic)
-      id->job = 1;
-      mumps_call();
-   }
-   else
-   {
-      // Reuse symbolic factorization / ordering
-      MFEM_VERIFY(id->n == global_n,
-                  "ReorderingReuse requires same global size (id->n mismatch)");
-
-      // Update matrix pointers (pattern is assumed compatible)
-      id->nnz_loc = nnz;
-      id->irn_loc = I;
-      id->jcn_loc = J;
-      id->a_loc   = A;
-   }
-
-   // Factorization
-   id->job = 2;
-   {
-      const int mem_relax_lim = 200;
-      while (true)
-      {
-         mumps_call();
-         if (id->MUMPS_INFOG(1) < 0)
-         {
-            if (id->MUMPS_INFOG(1) == -8 || id->MUMPS_INFOG(1) == -9)
-            {
-               id->MUMPS_ICNTL(14) += 20;
-               MFEM_VERIFY(id->MUMPS_ICNTL(14) <= mem_relax_lim,
-                           "Memory relaxation limit reached for ComplexMUMPSSolver factorization");
-               if (myid == 0 && print_level > 0)
-               {
-                  out << "Re-running ComplexMUMPSSolver factorization with memory relaxation "
-                      << id->MUMPS_ICNTL(14) << '\n';
-               }
-            }
-            else
-            {
-               MFEM_ABORT("Error during ComplexMUMPSSolver numerical factorization");
-            }
-         }
-         else { break; }
-      }
-   }
-
-   // Done with input storage
-   if (csr_op_r) { hypre_CSRMatrixDestroy(csr_op_r);}
-   if (csr_op_i) { hypre_CSRMatrixDestroy(csr_op_i);}
-   delete [] I;
-   delete [] J;
-   delete [] A;
-
-   // Post-factorization RHS/SOL setup
-   id->nrhs = -1;
-
-#if MFEM_MUMPS_VERSION >= 530
-   // Distributed RHS/SOL sizes
-   id->nloc_rhs = n_loc;
-   id->lrhs_loc = n_loc;
-   id->lsol_loc = id->MUMPS_INFO(23);
-
-   delete [] irhs_loc;
-   irhs_loc = new int[id->lrhs_loc];
-   for (int i = 0; i < n_loc; i++)
-   {
-      irhs_loc[i] = row_start + i + 1;
-   }
-   id->irhs_loc = irhs_loc;
-
-   delete [] isol_loc;
-   isol_loc = new int[id->lsol_loc];
-   id->isol_loc = isol_loc;
-
-   row_starts.SetSize(numProcs);
-   MPI_Allgather(&row_start, 1, MPI_INT, row_starts, 1, MPI_INT, comm);
-
-   // Reset cached buffers
-   delete [] rhs_loc; rhs_loc = nullptr;
-   delete [] sol_loc; sol_loc = nullptr;
-   rhs1_buf.clear();
-
-#else
-   // Centralized RHS/SOL on root
-   id->lrhs = id->n;
-
-   global_num_rows = id->n;
-
-   if (myid == 0)
-   {
-      delete [] recv_counts;
-      delete [] displs;
-      recv_counts = new int[numProcs];
-      displs = new int[numProcs];
-
-      delete [] rhs_glob;   rhs_glob = nullptr;
-      delete [] rhs_glob_r; rhs_glob_r = nullptr;
-      delete [] rhs_glob_i; rhs_glob_i = nullptr;
-   }
-
-   MPI_Gather(&n_loc, 1, MPI_INT, recv_counts, 1, MPI_INT, 0, comm);
-
-   if (myid == 0)
-   {
-      displs[0] = 0;
-      int s = 0;
-      for (int k = 0; k < numProcs - 1; k++)
-      {
-         s += recv_counts[k];
-         displs[k+1] = s;
-      }
-   }
-#endif
-}
-
-void ComplexMUMPSSolver::InitRhsSol(int nrhs) const
-{
-#if MFEM_MUMPS_VERSION >= 530
-
-   MFEM_VERIFY(id, "InitRhsSol called before SetOperator");
-
-   if (id->nrhs != nrhs)
-   {
-      delete [] rhs_loc;
-      delete [] sol_loc;
-
-      rhs_loc = new mumps_complex_t[(size_t)nrhs * (size_t)id->lrhs_loc];
-      sol_loc = new mumps_complex_t[(size_t)nrhs * (size_t)id->lsol_loc];
-
-      id->rhs_loc = rhs_loc;
-      id->sol_loc = sol_loc;
-   }
-   id->nrhs = nrhs;
-
-#else
-   MFEM_VERIFY(id, "InitRhsSol called before SetOperator");
-
-   id->nrhs = nrhs;
-   id->lrhs = id->n;
-
-   if (myid == 0)
-   {
-      const size_t N = (size_t)nrhs * (size_t)global_num_rows;
-
-      delete [] rhs_glob;
-      delete [] rhs_glob_r;
-      delete [] rhs_glob_i;
-
-      rhs_glob   = new mumps_complex_t[N];
-      rhs_glob_r = new real_t[N];
-      rhs_glob_i = new real_t[N];
-
-      id->rhs = rhs_glob;
-   }
-#endif
-}
-
-void ComplexMUMPSSolver::Mult(const Vector &x, Vector &y) const
-{
-   Array<const Vector *> X(1);
-   Array<Vector *> Y(1);
-   X[0] = &x;
-   Y[0] = &y;
-   ArrayMult(X, Y);
-}
-
-void ComplexMUMPSSolver::ArrayMult(const Array<const Vector *> &X,
-                                   Array<Vector *> &Y) const
-{
-   MFEM_ASSERT(X.Size() == Y.Size(),
-               "Number of columns mismatch in ComplexMUMPSSolver::Mult!");
-   MFEM_VERIFY(id, "ComplexMUMPSSolver::ArrayMult called before SetOperator");
-
-   InitRhsSol(X.Size());
-
-#if MFEM_MUMPS_VERSION >= 530
-   MFEM_VERIFY(irhs_loc && isol_loc, "RHS/SOL maps not initialized");
-   MFEM_VERIFY(rhs_loc && sol_loc, "RHS/SOL buffers not initialized");
-   const int n_loc = id->lrhs_loc;
-   const int nrhs  = id->nrhs;
-
-   // Pack all RHS
-   int xisign = (conv == ComplexOperator::BLOCK_SYMMETRIC) ? -1 : 1;
-   for (int i = 0; i < nrhs; i++)
-   {
-      MFEM_ASSERT(X[i], "Missing Vector in Mult!");
-      X[i]->HostRead();
-      MFEM_VERIFY(X[i]->Size() == 2*n_loc, "RHS size mismatch");
-
-      const real_t *xdata = X[i]->GetData();
-      const real_t *xr = xdata;
-      const real_t *xi = xdata + n_loc;
-
-      mumps_complex_t *dst = rhs_loc + i * n_loc;
-      for (int j = 0; j < n_loc; j++)
-      {
-         dst[j].r = xr[j];
-         dst[j].i = xisign * xi[j];
-      }
-   }
-
-   id->rhs_loc  = rhs_loc;
-   id->sol_loc  = sol_loc;
-   id->irhs_loc = irhs_loc;
-   id->isol_loc = isol_loc;
-
-   // MUMPS solve
-   id->job = 3;
-   mumps_call();
-
-   const int lsol = id->lsol_loc;
-
-   // Redistribute each solution column into Y
-   for (int i = 0; i < nrhs; i++)
-   {
-      MFEM_ASSERT(Y[i], "Missing output Vector in Mult!");
-      Y[i]->HostWrite();
-      MFEM_VERIFY(Y[i]->Size() == 2*n_loc, "Output size mismatch");
-
-      const mumps_complex_t *xcol = sol_loc + i * lsol;
-      RedistributeSol(isol_loc, xcol, Y[i]->GetData(), n_loc, lsol);
-   }
-
-#else // MFEM_MUMPS_VERSION < 530
-
-   const int nrhs = id->nrhs;
-
-   MFEM_VERIFY(X.Size() > 0 && X[0], "Missing RHS");
-   const int n_loc = X[0]->Size()/2;
-
-   for (int i = 0; i < nrhs; i++)
-   {
-      MFEM_ASSERT(X[i], "Missing Vector in Mult!");
-      X[i]->HostRead();
-      MFEM_VERIFY(X[i]->Size() == 2*n_loc, "RHS size mismatch");
-   }
-
-   // Gather each RHS column (real+imag separately) into root staging
-   for (int i = 0; i < nrhs; i++)
-   {
-      const real_t *xdata = X[i]->GetData();
-
-      MPI_Gatherv(xdata, n_loc, MPITypeMap<real_t>::mpi_type,
-                  rhs_glob_r + i * global_num_rows,
-                  recv_counts, displs, MPITypeMap<real_t>::mpi_type,
-                  0, comm);
-
-      MPI_Gatherv(xdata + n_loc, n_loc, MPITypeMap<real_t>::mpi_type,
-                  rhs_glob_i + i * global_num_rows,
-                  recv_counts, displs, MPITypeMap<real_t>::mpi_type,
-                  0, comm);
-   }
-
-   // Pack into MUMPS complex RHS on root: id->rhs is in-place
-   if (myid == 0)
-   {
-      for (int i = 0; i < nrhs; i++)
-      {
-         mumps_complex_t *dst = rhs_glob + i * global_num_rows;
-         const real_t *rr = rhs_glob_r + i * global_num_rows;
-         const real_t *ri = rhs_glob_i + i * global_num_rows;
-
-         for (int j = 0; j < global_num_rows; j++)
-         {
-            dst[j].r = rr[j];
-            dst[j].i = ri[j];
-         }
-      }
-      id->rhs = rhs_glob;
-   }
-
-   // Solve
-   id->job = 3;
-   mumps_call();
-
-   // Unpack to real/imag
-   if (myid == 0)
-   {
-      for (int i = 0; i < nrhs; i++)
-      {
-         const mumps_complex_t *src = rhs_glob + i * global_num_rows;
-         real_t *rr = rhs_glob_r + i * global_num_rows;
-         real_t *ri = rhs_glob_i + i * global_num_rows;
-
-         for (int j = 0; j < global_num_rows; j++)
-         {
-            rr[j] = src[j].r;
-            ri[j] = src[j].i;
-         }
-      }
-   }
-
-   // Scatter each RHS solution
-   for (int i = 0; i < nrhs; i++)
-   {
-      MFEM_ASSERT(Y[i], "Missing Vector in Mult!");
-      Y[i]->HostWrite();
-      MFEM_VERIFY(Y[i]->Size() == 2*n_loc, "Output size mismatch");
-
-      real_t *ydata = Y[i]->GetData();
-
-      MPI_Scatterv(rhs_glob_r + i * global_num_rows,
-                   recv_counts, displs, MPITypeMap<real_t>::mpi_type,
-                   ydata, n_loc, MPITypeMap<real_t>::mpi_type,
-                   0, comm);
-
-      MPI_Scatterv(rhs_glob_i + i * global_num_rows,
-                   recv_counts, displs, MPITypeMap<real_t>::mpi_type,
-                   ydata + n_loc, n_loc, MPITypeMap<real_t>::mpi_type,
-                   0, comm);
-   }
-
-#endif
-}
-
-void ComplexMUMPSSolver::MultTranspose(const Vector &x, Vector &y) const
-{
-   MFEM_VERIFY(id, "MultTranspose called before SetOperator");
-
-   // Transpose solve
-   id->MUMPS_ICNTL(9) = 0;
-   Mult(x, y);
-   id->MUMPS_ICNTL(9) = 1;
-}
-
-void ComplexMUMPSSolver::ArrayMultTranspose(const Array<const Vector *> &X,
-                                            Array<Vector *> &Y) const
-{
-   MFEM_VERIFY(id, "ArrayMultTranspose called before SetOperator");
-
-   // Transpose solve
-   id->MUMPS_ICNTL(9) = 0;
-   ArrayMult(X, Y);
-   id->MUMPS_ICNTL(9) = 1;
-}
-
-void ComplexMUMPSSolver::SetParameters()
-{
-   // Output stream for error messages
-   id->MUMPS_ICNTL(1) = 6;
-   // Output stream for diagnostic printing local to each proc
-   id->MUMPS_ICNTL(2) = 0;
-   // Output stream for global info
-   id->MUMPS_ICNTL(3) = 6;
-   // Level of error printing
-   id->MUMPS_ICNTL(4) = print_level;
-
-   // Input matrix format (assembled)
-   id->MUMPS_ICNTL(5) = 0;
-   // Use A or A^T
-   id->MUMPS_ICNTL(9) = 1;
-   // Iterative refinement (disabled)
-   id->MUMPS_ICNTL(10) = 0;
-   // Error analysis-statistics (disabled)
-   id->MUMPS_ICNTL(11) = 0;
-   // Use of ScaLAPACK (disabled)
-   id->MUMPS_ICNTL(13) = 0;
-   // Workspace relaxation (% increase)
-   id->MUMPS_ICNTL(14) = 20;
-   // OpenMP threads (default)
-   id->MUMPS_ICNTL(16) = 0;
-   // Matrix input format (distributed)
-   id->MUMPS_ICNTL(18) = 3;
-   // Schur complement (none)
-   id->MUMPS_ICNTL(19) = 0;
-
-#if MFEM_MUMPS_VERSION >= 530
-   // Distributed RHS
-   id->MUMPS_ICNTL(20) = 10;
-   // Distributed Sol
-   id->MUMPS_ICNTL(21) = 1;
-#else
-   // Centralized RHS
-   id->MUMPS_ICNTL(20) = 0;
-   // Centralized Sol
-   id->MUMPS_ICNTL(21) = 0;
-#endif
-
-   // Out-of-core (disabled)
-   id->MUMPS_ICNTL(22) = 0;
-   // Max size of working memory (default)
-   id->MUMPS_ICNTL(23) = 0;
-
-   switch (reorder_method)
-   {
-      case ReorderingStrategy::AUTOMATIC:
-         id->MUMPS_ICNTL(28) = 0;
-         id->MUMPS_ICNTL(7) = 7;
-         id->MUMPS_ICNTL(29) = 0;
-         break;
-      case ReorderingStrategy::AMD:
-         id->MUMPS_ICNTL(28) = 1;
-         id->MUMPS_ICNTL(7) = 0;
-         break;
-      case ReorderingStrategy::AMF:
-         id->MUMPS_ICNTL(28) = 1;
-         id->MUMPS_ICNTL(7) = 2;
-         break;
-      case ReorderingStrategy::PORD:
-         id->MUMPS_ICNTL(28) = 1;
-         id->MUMPS_ICNTL(7) = 4;
-         break;
-      case ReorderingStrategy::METIS:
-         id->MUMPS_ICNTL(28) = 1;
-         id->MUMPS_ICNTL(7) = 5;
-         break;
-      case ReorderingStrategy::PARMETIS:
-         id->MUMPS_ICNTL(28) = 2;
-         id->MUMPS_ICNTL(29) = 2;
-         break;
-      case ReorderingStrategy::SCOTCH:
-         id->MUMPS_ICNTL(28) = 1;
-         id->MUMPS_ICNTL(7) = 3;
-         break;
-      case ReorderingStrategy::PTSCOTCH:
-         id->MUMPS_ICNTL(28) = 2;
-         id->MUMPS_ICNTL(29) = 1;
-         break;
-      default:
-         break; // This should be unreachable
-   }
-
-}
-
-void ComplexMUMPSSolver::BuildUnionCOO(const int n_loc,
-                                       const int row_start_,
-                                       const int *Ir, const int *Jr, const real_t *Vr,
-                                       const int *Ii, const int *Ji, const real_t *Vi,
-                                       std::vector<int> &Icoo,
-                                       std::vector<int> &Jcoo,
-                                       std::vector<mumps_complex_t> &Zcoo) const
-{
-   for (int r = 0; r < n_loc; ++r)
-   {
-      std::unordered_map<int, std::pair<real_t, real_t>> row;
-
-      const int rr0 = Ir ? Ir[r]   : 0;
-      const int rr1 = Ir ? Ir[r+1] : 0;
-      const int ii0 = Ii ? Ii[r]   : 0;
-      const int ii1 = Ii ? Ii[r+1] : 0;
-
-      row.reserve((rr1 - rr0) + (ii1 - ii0));
-
-      if (Ir)
-      {
-         for (int p = rr0; p < rr1; ++p) { row[Jr[p]].first += Vr[p]; }
-      }
-      if (Ii)
-      {
-         for (int p = ii0; p < ii1; ++p) { row[Ji[p]].second += Vi[p]; }
-      }
-
-      for (const auto &kv : row)
-      {
-         Icoo.push_back(row_start_ + r + 1);
-         Jcoo.push_back(kv.first + 1);
-         Zcoo.push_back(mumps_complex_t{kv.second.first, kv.second.second});
-      }
-   }
-}
-
-#if MFEM_MUMPS_VERSION >= 530
-int ComplexMUMPSSolver::GetRowRank(int i, const Array<int> &row_starts_) const
-{
-   if (row_starts_.Size() == 1) { return 0; }
-   auto up = std::upper_bound(row_starts_.begin(), row_starts_.end(), i);
-   return (int)std::distance(row_starts_.begin(), up) - 1;
-}
-
-void ComplexMUMPSSolver::RedistributeSol(const int *row_map,
-                                         const mumps_complex_t *x,
-                                         real_t *y_ri,
-                                         int n_loc,
-                                         int lsol_loc) const
-{
-   int *send_count = new int[numProcs]();
-   for (int i = 0; i < lsol_loc; i++)
-   {
-      const int j = row_map[i] - 1;
-      const int row_rank = GetRowRank(j, row_starts);
-      if (myid == row_rank) { continue; }
-      send_count[row_rank]++;
-   }
-
-   int *recv_count = new int[numProcs];
-   MPI_Alltoall(send_count, 1, MPI_INT, recv_count, 1, MPI_INT, comm);
-
-   int *send_displ = new int[numProcs]; send_displ[0] = 0;
-   int *recv_displ = new int[numProcs]; recv_displ[0] = 0;
-
-   int sbuff_size = send_count[numProcs-1];
-   int rbuff_size = recv_count[numProcs-1];
-   for (int k = 0; k < numProcs - 1; k++)
-   {
-      send_displ[k+1] = send_displ[k] + send_count[k];
-      recv_displ[k+1] = recv_displ[k] + recv_count[k];
-      sbuff_size += send_count[k];
-      rbuff_size += recv_count[k];
-   }
-
-   int   *sendbuf_index = new int[sbuff_size];
-   real_t *sendbuf_r    = new real_t[sbuff_size];
-   real_t *sendbuf_i    = new real_t[sbuff_size];
-   int   *soffs         = new int[numProcs]();
-
-   for (int i = 0; i < lsol_loc; i++)
-   {
-      const int j = row_map[i] - 1;
-      const int row_rank = GetRowRank(j, row_starts);
-
-      const real_t xr = (real_t)x[i].r;
-      const real_t xi = (real_t)x[i].i;
-
-      if (myid == row_rank)
-      {
-         const int local_index = j - row_start;
-         y_ri[local_index]        = xr;
-         y_ri[local_index+n_loc]  = xi;
-      }
-      else
-      {
-         const int k = send_displ[row_rank] + soffs[row_rank];
-         sendbuf_index[k] = j;
-         sendbuf_r[k] = xr;
-         sendbuf_i[k] = xi;
-         soffs[row_rank]++;
-      }
-   }
-
-   int    *recvbuf_index = new int[rbuff_size];
-   real_t *recvbuf_r     = new real_t[rbuff_size];
-   real_t *recvbuf_i     = new real_t[rbuff_size];
-
-   MPI_Alltoallv(sendbuf_index, send_count, send_displ, MPI_INT,
-                 recvbuf_index, recv_count, recv_displ, MPI_INT, comm);
-
-   MPI_Alltoallv(sendbuf_r, send_count, send_displ, MPITypeMap<real_t>::mpi_type,
-                 recvbuf_r, recv_count, recv_displ, MPITypeMap<real_t>::mpi_type, comm);
-
-   MPI_Alltoallv(sendbuf_i, send_count, send_displ, MPITypeMap<real_t>::mpi_type,
-                 recvbuf_i, recv_count, recv_displ, MPITypeMap<real_t>::mpi_type, comm);
-
-   for (int i = 0; i < rbuff_size; i++)
-   {
-      const int local_index = recvbuf_index[i] - row_start;
-      y_ri[local_index]       = recvbuf_r[i];
-      y_ri[local_index+n_loc] = recvbuf_i[i];
-   }
-
-   delete [] recvbuf_i;
-   delete [] recvbuf_r;
-   delete [] recvbuf_index;
-   delete [] soffs;
-   delete [] sendbuf_i;
-   delete [] sendbuf_r;
-   delete [] sendbuf_index;
-   delete [] recv_displ;
-   delete [] send_displ;
-   delete [] recv_count;
-   delete [] send_count;
-}
-#endif // MFEM_MUMPS_VERSION >= 530
-
-#endif // MFEM_USE_COMPLEX_MUMPS
-
 } // namespace mfem
+
+#endif // MFEM_USE_MPI
+#endif // MFEM_USE_MUMPS

linalg/mumps.hpp

@@ -14,34 +14,22 @@
 
 #include "../config/config.hpp"
 
-#if defined(MFEM_USE_MUMPS) || defined(MFEM_USE_COMPLEX_MUMPS)
+#ifdef MFEM_USE_MUMPS
+#ifdef MFEM_USE_MPI
+
 #include "operator.hpp"
 #include "hypre.hpp"
 #include <mpi.h>
-#endif
 
-#ifdef MFEM_USE_MUMPS
 #ifdef MFEM_USE_SINGLE
 #include "smumps_c.h"
 #else
 #include "dmumps_c.h"
 #endif
-#endif
-
-#ifdef MFEM_USE_COMPLEX_MUMPS
-#include "complex_operator.hpp"
-#ifdef MFEM_USE_SINGLE
-#include "cmumps_c.h"
-#else
-#include "zmumps_c.h"
-#endif
-#include <vector>
-#endif
 
 namespace mfem
 {
 
-#ifdef MFEM_USE_MUMPS
 /**
  * @brief MUMPS: A Parallel Sparse Direct Solver
  *
@@ -64,7 +52,7 @@ public:
    /// Specify the reordering strategy for the MUMPS solver
    enum ReorderingStrategy
    {
-      /// Let MUMPS automatically decide the reordering strategy
+      /// Let MUMPS automatically decide the reording strategy
       AUTOMATIC = 0,
       /// Approximate Minimum Degree with auto quasi-dense row detection is used
       AMD,
@@ -110,7 +98,7 @@ public:
    void Mult(const Vector &x, Vector &y) const;
 
    /**
-    * @brief Solve $ Y_i = Op^{-1} X_i $
+    * @brief Solve $ Y_i = Op^{-T} X_i $
     *
     * @param X Array of RHS vectors
     * @param Y Array of Solution vectors
@@ -141,8 +129,8 @@ public:
     * - 0:  No output printed
     * - 1:  Only errors printed
     * - 2:  Errors, warnings, and main stats printed
-    * - 3:  Errors, warnings, main stats, and terse diagnostics printed
-    * - 4:  Errors, warnings, main stats, diagnostics, and input/output printed
+    * - 3:  Errors, warning, main stats, and terse diagnostics printed
+    * - 4:  Errors, warning, main stats, diagnostics, and input/output printed
     *
     * @param print_lvl Print level, default is 2
     *
@@ -270,247 +258,8 @@ private:
 #endif
 }; // mfem::MUMPSSolver class
 
-#endif // MFEM_USE_MUMPS
-
-#ifdef MFEM_USE_COMPLEX_MUMPS
-/**
- * @brief Complex MUMPS: Parallel sparse direct solver for ComplexHypreParMatrix
- *
- * Notes:
- *  - Expects Operator to be a ComplexHypreParMatrix.
- *  - Complex vectors are assumed packed as [Re; Im] in a real Vector.
- *  - SetOperator(): analysis + factorization
- *  - Mult()       : solve
- */
-class ComplexMUMPSSolver : public Solver
-{
-public:
-   /// Specify the reordering strategy
-   enum ReorderingStrategy
-   {
-      /// Let MUMPS automatically decide the reordering strategy
-      AUTOMATIC = 0,
-      /// Approximate Minimum Degree with auto quasi-dense row detection is used
-      AMD,
-      /// Approximate Minimum Fill method will be used
-      AMF,
-      /// The PORD library will be used
-      PORD,
-      /// The METIS library will be used
-      METIS,
-      /// The ParMETIS library will be used
-      PARMETIS,
-      /// The Scotch library will be used
-      SCOTCH,
-      /// The PTScotch library will be used
-      PTSCOTCH
-   };
-
-   /**
-    * @brief Constructor with MPI_Comm parameter.
-    */
-   ComplexMUMPSSolver(MPI_Comm comm_);
-   /**
-    * @brief Constructor with a ComplexHypreParMatrix Operator.
-    */
-   ComplexMUMPSSolver(const Operator &op);
-
-   /**
-    * @brief Set the Operator and perform factorization
-    *
-    * @a op needs to be of type ComplexHypreParMatrix.
-    *
-    * @param op Operator used in factorization and solve
-    */
-   void SetOperator(const Operator &op);
-
-   /**
-    * @brief Solve $ y = Op^{-1} x $
-    *
-    * @param x RHS vector
-    * @param y Solution vector
-    */
-   void Mult(const Vector &x, Vector &y) const;
-   /**
-    * @brief Solve $ Y_i = Op^{-1} X_i $
-    *
-    * @param X Array of RHS vectors
-    * @param Y Array of Solution vectors
-    */
-   void ArrayMult(const Array<const Vector *> &X, Array<Vector *> &Y) const;
-   /**
-    * @brief Transpose Solve $ y = Op^{-T} x $
-    * @note This is not a Hermitian/conjugate-transpose solve.
-    *
-    * @param x RHS vector
-    * @param y Solution vector
-    */
-   void MultTranspose(const Vector &x, Vector &y) const;
-
-   /**
-    * @brief Transpose Solve $ Y_i = Op^{-T} X_i $
-    * @note This is not a Hermitian/conjugate-transpose solve.
-    *
-    * @param X Array of RHS vectors
-    * @param Y Array of Solution vectors
-    */
-   void ArrayMultTranspose(const Array<const Vector *> &X,
-                           Array<Vector *> &Y) const;
-
-   /**
-    * @brief Set the error print level for MUMPS
-    *
-    * Supported values are:
-    * - 0:  No output printed
-    * - 1:  Only errors printed
-    * - 2:  Errors, warnings, and main stats printed
-    * - 3:  Errors, warnings, main stats, and terse diagnostics printed
-    * - 4:  Errors, warnings, main stats, diagnostics, and input/output printed
-    *
-    * @param print_lvl Print level, default is 2
-    *
-    * @note This method has to be called before SetOperator
-    */
-   void SetPrintLevel(int print_lvl) { print_level = print_lvl;}
-
-   /**
-    * @brief Set the reordering strategy
-    *
-    * Supported reorderings are: ComplexMUMPSSolver::AUTOMATIC,
-    * ComplexMUMPSSolver::AMD, ComplexMUMPSSolver::AMF,
-    * ComplexMUMPSSolver::PORD, ComplexMUMPSSolver::METIS,
-    * ComplexMUMPSSolver::PARMETIS, ComplexMUMPSSolver::SCOTCH,
-    * and ComplexMUMPSSolver::PTSCOTCH
-    *
-    * @param method Reordering method
-    *
-    * @note This method has to be called before SetOperator
-    */
-   void SetReorderingStrategy(ReorderingStrategy method) { reorder_method = method; }
-
-   /**
-    * @brief Set the flag controlling reuse of the symbolic factorization
-    * for multiple operators
-    *
-    * @param reuse Flag to reuse symbolic factorization
-    *
-    * @note This method has to be called before repeated calls to SetOperator
-    */
-   void SetReorderingReuse(bool reuse) { reorder_reuse = reuse; }
-
-   ~ComplexMUMPSSolver();
-
-private:
-   // MPI communicator
-   MPI_Comm comm = MPI_COMM_NULL;
-
-   // Number of procs
-   int numProcs;
-
-   // MPI rank
-   int myid;
-
-   // Parameter controlling the printing level
-   int print_level = 0;
-
-   // Parameter controlling the reordering strategy
-   ReorderingStrategy reorder_method = ReorderingStrategy::AUTOMATIC;
-
-   // Parameter controlling whether or not to reuse the symbolic factorization
-   // for multiple calls to SetOperator
-   bool reorder_reuse = false;
-
-   // Local row offsets
-   int row_start;
-
-   // ComplexMUMPS object
-#ifdef MFEM_USE_SINGLE
-   CMUMPS_STRUC_C *id = nullptr;
-   using mumps_complex_t = mumps_complex;
-#else
-   ZMUMPS_STRUC_C *id = nullptr;
-   using mumps_complex_t = mumps_double_complex;
-#endif
-
-   /// Method for initialization
-   void Init(MPI_Comm comm_);
-
-   /// Method for setting ComplexMUMPS internal parameters
-   void SetParameters();
-
-   /// Method for configuring storage for distributed/centralized
-   /// RHS and solution
-   void InitRhsSol(int nrhs) const;
-
-   /// Method for calling the single/double ComplexMUMPS solver
-   inline void mumps_call() const
-   {
-#ifdef MFEM_USE_SINGLE
-      cmumps_c(id);
-#else
-      zmumps_c(id);
-#endif
-   }
-
-   /// Method for building the COO format of the combined complex operator
-   /// from the real and imaginary parts. This is particularly useful when
-   /// real and imaginary parts have different sparsity patterns.
-   void BuildUnionCOO(const int n_loc,
-                      const int row_start,
-                      const int *Ir, const int *Jr, const real_t *Vr,
-                      const int *Ii, const int *Ji, const real_t *Vi,
-                      std::vector<int> &Icoo,
-                      std::vector<int> &Jcoo,
-                      std::vector<mumps_complex_t> &Zcoo) const;
-
-#if MFEM_MUMPS_VERSION >= 530
-   // Row offsets on all procs
-   Array<int> row_starts;
-
-   // Local RHS row indices
-   int *irhs_loc = nullptr;
-
-   // Local solution row map returned by MUMPS
-   int *isol_loc = nullptr;
-
-   // Cached buffers
-   mutable mumps_complex_t *rhs_loc = nullptr;
-   mutable mumps_complex_t *sol_loc = nullptr;
-
-   // RHS buffers
-   mutable std::vector<mumps_complex_t> rhs1_buf;
-
-   // These two methods are needed to distribute the local solution
-   // vectors returned by MUMPS to the original MFEM parallel partition
-   int GetRowRank(int i, const Array<int> &row_starts_) const;
-
-   void RedistributeSol(const int *row_map,
-                        const mumps_complex_t *x,
-                        real_t *y_ri,
-                        int n_loc,
-                        int lsol_loc) const;
-
-   ComplexOperator::Convention conv = ComplexOperator::Convention::HERMITIAN;
-
-#else
-   // Root-gather path
-   int global_num_rows;
-
-   // Arrays needed for MPI_Gatherv and MPI_Scatterv
-   int *recv_counts = nullptr;
-   int *displs = nullptr;
-
-   // Complex RHS/solution on root
-   mutable mumps_complex_t *rhs_glob = nullptr;
-
-   // Cached real/imag staging on root
-   mutable real_t *rhs_glob_r = nullptr;
-   mutable real_t *rhs_glob_i = nullptr;
-#endif
-};
-
-#endif // MFEM_USE_COMPLEX_MUMPS
-
 } // namespace mfem
 
+#endif // MFEM_USE_MPI
+#endif // MFEM_USE_MUMPS
 #endif // MFEM_MUMPS

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

@@ -229,7 +229,7 @@ CXXFLAGS ?= $(OPTIM_FLAGS)
 # MPI configuration
 ifneq ($(MFEM_USE_MPI),YES)
    HOST_CXX = $(CXX)
-   PKGS_NEED_MPI = SUPERLU MUMPS COMPLEX_MUMPS STRUMPACK PETSC PUMI SLEPC MKL_CPARDISO
+   PKGS_NEED_MPI = SUPERLU MUMPS STRUMPACK PETSC PUMI SLEPC MKL_CPARDISO
    $(foreach mpidep,$(PKGS_NEED_MPI),$(if $(MFEM_USE_$(mpidep):NO=),\
      $(warning *** [MPI is OFF] setting MFEM_USE_$(mpidep) = NO)\
      $(eval override MFEM_USE_$(mpidep)=NO),))
@@ -304,7 +304,7 @@ ifeq ($(MFEM_USE_LEGACY_OPENMP),YES)
 endif
 
 # List of MFEM dependencies, that require the *_LIB variable to be non-empty
-MFEM_REQ_LIB_DEPS = SUPERLU MUMPS COMPLEX_MUMPS METIS FMS CONDUIT SIDRE LAPACK SUNDIALS\
+MFEM_REQ_LIB_DEPS = SUPERLU MUMPS METIS FMS CONDUIT SIDRE LAPACK SUNDIALS\
  SUITESPARSE STRUMPACK GINKGO GNUTLS HDF5 NETCDF SLEPC PETSC MPFR PUMI HIOP\
  GSLIB OCCA CEED RAJA UMPIRE MKL_CPARDISO MKL_PARDISO AMGX MAGMA CALIPER PARELAG\
  TRIBOL BENCHMARK MOONOLITH ALGOIM CUDSS
@@ -374,7 +374,7 @@ MFEM_DEFINES = MFEM_VERSION MFEM_VERSION_STRING MFEM_GIT_STRING MFEM_USE_MPI\
  MFEM_USE_PUMI MFEM_USE_HIOP MFEM_USE_GSLIB MFEM_USE_CUDA MFEM_USE_HIP\
  MFEM_USE_OCCA MFEM_USE_MOONOLITH MFEM_USE_CEED MFEM_USE_RAJA MFEM_USE_UMPIRE\
  MFEM_USE_SIMD MFEM_USE_ADIOS2 MFEM_USE_MKL_CPARDISO MFEM_USE_MKL_PARDISO MFEM_USE_AMGX\
- MFEM_USE_MAGMA MFEM_USE_MUMPS MFEM_USE_COMPLEX_MUMPS MFEM_USE_ADFORWARD MFEM_USE_CODIPACK MFEM_USE_CALIPER\
+ MFEM_USE_MAGMA MFEM_USE_MUMPS MFEM_USE_ADFORWARD MFEM_USE_CODIPACK MFEM_USE_CALIPER\
  MFEM_USE_BENCHMARK MFEM_USE_PARELAG MFEM_USE_TRIBOL MFEM_USE_ALGOIM MFEM_USE_ENZYME\
  MFEM_SOURCE_DIR MFEM_INSTALL_DIR MFEM_SHARED_BUILD MFEM_USE_DOUBLE MFEM_USE_SINGLE\
  MFEM_USE_CUDSS MFEM_CUDSS_COMM_LIB MFEM_CUDSS_THREADING_LIB
@@ -412,7 +412,7 @@ MFEM_INSTALL_DIR = $(abspath $(MFEM_PREFIX))
 # If we have 'config' target, export variables used by config/makefile
 ifneq (,$(filter config,$(MAKECMDGOALS)))
    export $(MFEM_DEFINES) MFEM_DEFINES $(MFEM_CONFIG_VARS) MFEM_CONFIG_VARS
-   export VERBOSE HYPRE_OPT PUMI_DIR MUMPS_OPT COMPLEX_MUMPS_OPT GSLIB_OPT CUDSS_OPT
+   export VERBOSE HYPRE_OPT PUMI_DIR MUMPS_OPT GSLIB_OPT CUDSS_OPT
 endif
 
 # If we have 'install' target, export variables used by config/makefile
@@ -742,7 +742,6 @@ status info:
 	$(info MFEM_USE_SUPERLU       = $(MFEM_USE_SUPERLU))
 	$(info MFEM_USE_SUPERLU5      = $(MFEM_USE_SUPERLU5))
 	$(info MFEM_USE_MUMPS         = $(MFEM_USE_MUMPS))
-	$(info MFEM_USE_COMPLEX_MUMPS = $(MFEM_USE_COMPLEX_MUMPS))
 	$(info MFEM_USE_STRUMPACK     = $(MFEM_USE_STRUMPACK))
 	$(info MFEM_USE_CUDSS         = $(MFEM_USE_CUDSS))
 	$(info MFEM_USE_GINKGO        = $(MFEM_USE_GINKGO))

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)

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/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

tests/unit/linalg/test_direct_solvers.cpp

@@ -188,6 +188,8 @@ TEST_CASE("Serial Direct Solvers", "[GPU]")
 
 TEST_CASE("Parallel Direct Solvers", "[Parallel], [GPU]")
 {
+   int rank;
+   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    const int ne = 4;
    for (int dim = 1; dim < 4; ++dim)
    {
@@ -355,113 +357,3 @@ TEST_CASE("Parallel Direct Solvers", "[Parallel], [GPU]")
 }
 
 #endif
-
-
-#ifdef MFEM_USE_COMPLEX_MUMPS
-
-TEST_CASE("ComplexMUMPS Solver", "[Parallel], [GPU]")
-{
-   const int ne = 4;
-   const int order = 3;
-   auto dim = GENERATE(1, 2, 3);
-   auto conv = GENERATE(ComplexOperator::HERMITIAN,
-                        ComplexOperator::BLOCK_SYMMETRIC);
-   CAPTURE(dim, conv);
-
-   Mesh mesh;
-   if (dim == 1)
-   {
-      mesh = Mesh::MakeCartesian1D(ne, 1.0);
-   }
-   else if (dim == 2)
-   {
-      mesh = Mesh::MakeCartesian2D(
-                ne, ne, Element::QUADRILATERAL, 1, 1.0, 1.0);
-   }
-   else
-   {
-      mesh = Mesh::MakeCartesian3D(
-                ne, ne, ne, Element::HEXAHEDRON, 1.0, 1.0, 1.0);
-   }
-
-   ParMesh pmesh(MPI_COMM_WORLD, mesh);
-   mesh.Clear();
-   H1_FECollection fec(order, dim);
-   ParFiniteElementSpace fespace(&pmesh, &fec);
-   Array<int> ess_tdof_list, ess_bdr;
-   if (pmesh.bdr_attributes.Size())
-   {
-      ess_bdr.SetSize(pmesh.bdr_attributes.Max());
-      ess_bdr = 1;
-      fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
-   }
-
-   ConstantCoefficient one(1.0);
-   ConstantCoefficient negone(-1.0);
-   ConstantCoefficient two(2.0);
-
-   ComplexLinearForm b(&fespace, conv);
-   b.AddDomainIntegrator(new DomainLFIntegrator(one), new DomainLFIntegrator(two));
-   b.Assemble();
-
-   ParSesquilinearForm a_r(&fespace, conv);
-   a_r.AddDomainIntegrator(new DiffusionIntegrator(one), nullptr);
-   a_r.AddDomainIntegrator(new MassIntegrator(one), nullptr);
-   a_r.Assemble();
-
-   ParSesquilinearForm a_i(&fespace, conv);
-   a_i.AddDomainIntegrator(nullptr, new DiffusionIntegrator(one));
-   a_i.AddDomainIntegrator(nullptr, new MassIntegrator(one));
-   a_i.Assemble();
-
-   ParSesquilinearForm a_c(&fespace, conv);
-   a_c.AddDomainIntegrator(new DiffusionIntegrator(one), new MassIntegrator(two));
-   a_c.AddDomainIntegrator(new MassIntegrator(negone),nullptr);
-   a_c.Assemble();
-
-   ParComplexGridFunction x_c(&fespace);
-   ParComplexGridFunction x_r(&fespace);
-   ParComplexGridFunction x_i(&fespace);
-   x_c = 0.0; x_r = 0.0; x_i = 0.0;
-
-   OperatorPtr Ac, Ar, Ai;
-   Vector Bc, Br, Bi, Xr, Xc, Xi;
-   a_c.FormLinearSystem(ess_tdof_list, x_c, b, Ac, Xc, Bc);
-   a_r.FormLinearSystem(ess_tdof_list, x_r, b, Ar, Xr, Br);
-   a_i.FormLinearSystem(ess_tdof_list, x_i, b, Ai, Xi, Bi);
-
-   ComplexHypreParMatrix *Ahc = Ac.As<ComplexHypreParMatrix>();
-   ComplexHypreParMatrix *Ahr = Ar.As<ComplexHypreParMatrix>();
-   ComplexHypreParMatrix *Ahi = Ai.As<ComplexHypreParMatrix>();
-
-   ComplexMUMPSSolver cmumps(MPI_COMM_WORLD);
-   cmumps.SetPrintLevel(0);
-
-   cmumps.SetOperator(*Ahc);
-   cmumps.Mult(Bc, Xc);
-
-   cmumps.SetOperator(*Ahr);
-   cmumps.Mult(Br, Xr);
-
-   cmumps.SetOperator(*Ahi);
-   cmumps.Mult(Bi, Xi);
-
-   Vector Yc(Xc.Size()), Yr(Xr.Size()), Yi(Xi.Size());
-   Ahc->Mult(Xc, Yc); Ahr->Mult(Xr, Yr); Ahi->Mult(Xi, Yi);
-
-   if (conv == ComplexOperator::BLOCK_SYMMETRIC)
-   {
-      int n = Bc.Size()/2;
-      real_t * dataBc = Bc.GetData(); Vector Bcimag(&dataBc[n],n);
-      real_t * dataBr = Br.GetData(); Vector Brimag(&dataBr[n],n);
-      real_t * dataBi = Bi.GetData(); Vector Biimag(&dataBi[n],n);
-      Bcimag *= -1; Brimag *= -1; Biimag *= -1;
-   }
-
-   Yc -= Bc; REQUIRE(Yc.Norml2() < 1.e-12);
-   Yr -= Br; REQUIRE(Yr.Norml2() < 1.e-12);
-   Yi -= Bi; REQUIRE(Yi.Norml2() < 1.e-12);
-
-} // Test case "ComplexMUMPS Solver"
-
-#endif // MFEM_USE_COMPLEX_MUMPS