Add GPU scan support for Array<T>::PartialSum

Remove explicit template instantiations

.github/workflows/stale.yml

@@ -29,3 +29,47 @@ jobs:
         operations-per-run: 500
         exempt-issue-labels: "bug,WIP,ready-for-review,in-review,in-next"
         exempt-pr-labels: "bug,WIP,ready-for-review,in-review,in-next"
+
+  # Stale action for PRs with "in-review" label.
+  stale-in-review-pr:
+
+    runs-on: ubuntu-latest
+    permissions:
+      issues: write
+      pull-requests: write
+      actions: write
+
+    steps:
+    - uses: actions/stale@v9
+      with:
+        repo-token: ${{ secrets.GITHUB_TOKEN }}
+        stale-pr-message: ':warning: This PR has been automatically marked as stale because it has not had any activity in the last 150 days. *If no activity occurs in the next 30 days, it will be automatically closed.* Thank you for your contributions.'
+        only-pr-labels: "in-review"
+        days-before-pr-stale: 150
+        days-before-pr-close: 30
+        days-before-issue-stale: -1
+        days-before-issue-close: -1
+        stale-pr-label: 'stale'
+        operations-per-run: 500
+
+  # Stale action for PRs with "WIP" label.
+  stale-wip-pr:
+
+    runs-on: ubuntu-latest
+    permissions:
+      issues: write
+      pull-requests: write
+      actions: write
+
+    steps:
+    - uses: actions/stale@v9
+      with:
+        repo-token: ${{ secrets.GITHUB_TOKEN }}
+        stale-pr-message: ':warning: This PR has been automatically marked as stale because it has not had any activity in the last 300 days. *If no activity occurs in the next 30 days, it will be automatically closed.* Thank you for your contributions.'
+        only-pr-labels: "WIP"
+        days-before-pr-stale: 300
+        days-before-pr-close: 30
+        days-before-issue-stale: -1
+        days-before-issue-close: -1
+        stale-pr-label: 'stale'
+        operations-per-run: 500

.gitignore

@@ -208,10 +208,13 @@ miniapps/electromagnetics/volta
 miniapps/electromagnetics/tesla
 miniapps/electromagnetics/maxwell
 miniapps/electromagnetics/joule
+miniapps/electromagnetics/lorentz
 miniapps/electromagnetics/Volta-AMR*
 miniapps/electromagnetics/Tesla-AMR*
 miniapps/electromagnetics/Maxwell-Parallel*
 miniapps/electromagnetics/Joule_[0-9]*
+miniapps/electromagnetics/Lorentz_[0-9]*
+miniapps/electromagnetics/Lorentz.dat
 
 miniapps/gslib/field-diff
 miniapps/gslib/field-interp

CHANGELOG

@@ -91,6 +91,10 @@ New and updated examples and miniapps
 - Added a new miniapp (tools/gridfunction-bounds) to compute piecewise linear
   bounds on a given high-order grid function.
 
+- Added a new miniapp (electromagnetics/lorentz) which computes the trajectory
+  of a charged particle, subject to Lorentz forces, in electrostatic and/or
+  magnetostatic fields as computed by the volta or tesla miniapps.
+
 API changes:
 -----------
 - mfem::internal::tensor and mfem::internal::dual have been moved to

CMakeLists.txt

@@ -278,6 +278,11 @@ if (MFEM_USE_OPENMP OR MFEM_USE_LEGACY_OPENMP)
   endif()
 endif()
 
+# Umpire (must be included before hypre, so hypre can use it if needed)
+if (MFEM_USE_UMPIRE)
+   find_package(UMPIRE REQUIRED)
+endif()
+
 # MPI -> hypre; PETSc (optional)
 if (MFEM_USE_MPI)
   find_package(MPI REQUIRED)
@@ -502,11 +507,6 @@ if (MFEM_USE_RAJA)
    find_package(RAJA REQUIRED)
 endif()
 
-# UMPIRE
-if (MFEM_USE_UMPIRE)
-   find_package(UMPIRE REQUIRED)
-endif()
-
 # GOOGLE-BENCHMARK
 if (MFEM_USE_BENCHMARK)
    find_package(Benchmark REQUIRED)

config/cmake/modules/FindHYPRE.cmake

@@ -38,51 +38,91 @@ if (HYPRE_FOUND OR TARGET HYPRE)
 endif()
 
 if (HYPRE_FETCH OR FETCH_TPLS)
+  # Collect all HYPRE_ENABLE variables and pass them to hypre, assuming they are BOOL.
+  set(HYPRE_CMAKE_OPTIONS "")
+  get_cmake_property(all_vars VARIABLES)
+  foreach(var ${all_vars})
+    if(var MATCHES "^HYPRE_ENABLE")
+      list(APPEND HYPRE_CMAKE_OPTIONS "-D${var}:BOOL=${${var}}")
+    endif()
+  endforeach()
+
   set(HYPRE_FETCH_VERSION 2.33.0)
+  set(HYPRE_FETCH_TAG "v${HYPRE_FETCH_VERSION}" CACHE STRING "Tag, branch, or commit for HYPRE")
   add_library(HYPRE STATIC IMPORTED)
   # set options and associated dependencies
-  set(CMAKE_OPTIONS)
-  list(APPEND CMAKE_OPTIONS -DCMAKE_BUILD_TYPE:STRING=${CMAKE_BUILD_TYPE})
+  list(APPEND HYPRE_CMAKE_OPTIONS -DCMAKE_BUILD_TYPE:STRING=${CMAKE_BUILD_TYPE})
   if (MFEM_USE_CUDA)
-    list(APPEND CMAKE_OPTIONS -DHYPRE_WITH_CUDA:BOOL=ON)
+    list(APPEND HYPRE_CMAKE_OPTIONS -DHYPRE_ENABLE_CUDA:BOOL=ON -DCMAKE_CUDA_ARCHITECTURES:STRING=${CMAKE_CUDA_ARCHITECTURES})
     find_package(CUDAToolkit REQUIRED)
     target_link_libraries(HYPRE INTERFACE CUDA::cusparse CUDA::curand CUDA::cublas)
   elseif (MFEM_USE_HIP)
-    list(APPEND CMAKE_OPTIONS -DHYPRE_WITH_HIP:BOOL=ON)
+    list(APPEND HYPRE_CMAKE_OPTIONS -DHYPRE_ENABLE_HIP:BOOL=ON)
     find_package(rocsparse REQUIRED)
     find_package(rocrand REQUIRED)
     target_link_libraries(HYPRE INTERFACE rocsparse rocrand)
   endif()
+  if (MFEM_USE_CUDA OR MFEM_USE_HIP)
+    if (MFEM_USE_UMPIRE)
+      if (EXISTS ${umpire_DIR})
+        list(APPEND HYPRE_CMAKE_OPTIONS -DHYPRE_ENABLE_UMPIRE:BOOL=ON -Dumpire_DIR:PATH=${umpire_DIR})
+      else()
+        message(FATAL_ERROR "MFEM_USE_UMPIRE=ON, however umpire_DIR isn't visible to HYPRE")
+      endif()
+    else()
+      list(APPEND HYPRE_CMAKE_OPTIONS -DHYPRE_ENABLE_UMPIRE:BOOL=OFF)
+      message(WARNING
+"================================================================================
+ Umpire is disabled while building HYPRE with GPU support.
+ This is not recommended for performance reasons!
+ Consider enabling Umpire with -DMFEM_USE_UMPIRE=ON and providing -DUMPIRE_DIR.
+================================================================================")
+    endif()
+  endif()
   if (MFEM_USE_SINGLE)
-    list(APPEND CMAKE_OPTIONS -DHYPRE_ENABLE_SINGLE:BOOL=ON)
+    list(APPEND HYPRE_CMAKE_OPTIONS -DHYPRE_ENABLE_SINGLE:BOOL=ON)
   endif()
   # define external project and create future include directory so it is present
   # to pass CMake checks at end of MFEM configuration step
-  message(STATUS "Will fetch HYPRE ${HYPRE_FETCH_VERSION} to be built with ${CMAKE_OPTIONS}")
-  set(PREFIX ${CMAKE_BINARY_DIR}/fetch/hypre)
+  message(STATUS "Will fetch HYPRE ${HYPRE_FETCH_TAG} to be built with ${HYPRE_CMAKE_OPTIONS}")
+  set(HYPRE_INSTALL ${CMAKE_BINARY_DIR}/fetch/hypre)
   include(ExternalProject)
   ExternalProject_Add(hypre
     GIT_REPOSITORY https://github.com/hypre-space/hypre.git
-    GIT_TAG v${HYPRE_FETCH_VERSION}
+    GIT_TAG ${HYPRE_FETCH_TAG}
     GIT_SHALLOW TRUE
+    GIT_PROGRESS TRUE
     UPDATE_DISCONNECTED TRUE
     SOURCE_SUBDIR src
-    PREFIX ${PREFIX}
-    CMAKE_CACHE_ARGS -DCMAKE_INSTALL_PREFIX:PATH=${PREFIX} -DCMAKE_INSTALL_LIBDIR:PATH=lib ${CMAKE_OPTIONS})
-  file(MAKE_DIRECTORY ${PREFIX}/include)
+    PREFIX ${HYPRE_INSTALL}
+    BUILD_COMMAND ${CMAKE_COMMAND} --build . -- -j${CMAKE_BUILD_PARALLEL_LEVEL}
+    CMAKE_CACHE_ARGS -DCMAKE_INSTALL_PREFIX:PATH=${HYPRE_INSTALL} -DCMAKE_INSTALL_LIBDIR:PATH=lib ${HYPRE_CMAKE_OPTIONS})
+  file(MAKE_DIRECTORY ${HYPRE_INSTALL}/include)
   # set imported library target properties
   add_dependencies(HYPRE hypre)
   set_target_properties(HYPRE PROPERTIES
-    IMPORTED_LOCATION ${PREFIX}/lib/libHYPRE.a
-    INTERFACE_INCLUDE_DIRECTORIES ${PREFIX}/include)
+    IMPORTED_LOCATION ${HYPRE_INSTALL}/lib/libHYPRE.a
+    INTERFACE_INCLUDE_DIRECTORIES ${HYPRE_INSTALL}/include)
   # convert HYPRE version to integer
-  string(REGEX MATCHALL "[0-9]+" HYPRE_SPLIT_VERSION ${HYPRE_FETCH_VERSION})
-  list(GET HYPRE_SPLIT_VERSION 0 HYPRE_MAJOR_VERSION)
-  list(GET HYPRE_SPLIT_VERSION 1 HYPRE_MINOR_VERSION)
-  list(GET HYPRE_SPLIT_VERSION 2 HYPRE_PATCH_VERSION)
-  math(EXPR HYPRE_VERSION "10000*${HYPRE_MAJOR_VERSION} + 100*${HYPRE_MINOR_VERSION} + ${HYPRE_PATCH_VERSION}")
-  # set cache variables that would otherwise be set after mfem_find_package call
-  set(HYPRE_VERSION ${HYPRE_VERSION} CACHE STRING "HYPRE version." FORCE)
+  if (HYPRE_FETCH_TAG MATCHES "^v?([0-9]+)\\.([0-9]+)\\.([0-9]+)$")
+    # Exact release tag X.Y.Z
+    string(REGEX MATCHALL "[0-9]+" HYPRE_SPLIT_VERSION "${HYPRE_FETCH_TAG}")
+  elseif (HYPRE_FETCH_VERSION MATCHES "([0-9]+)\\.([0-9]+)(\\.([0-9]+))?")
+    string(REGEX MATCHALL "[0-9]+" HYPRE_SPLIT_VERSION "${HYPRE_FETCH_VERSION}")
+  else (NOT DEFINED HYPRE_VERSION)
+    message(FATAL_ERROR "Unable to find HYPRE release version. Please provide it via -DHYPRE_VERSION")
+  endif()
+  if (HYPRE_SPLIT_VERSION AND NOT DEFINED HYPRE_VERSION)
+    list(GET HYPRE_SPLIT_VERSION 0 HYPRE_MAJOR_VERSION)
+    list(GET HYPRE_SPLIT_VERSION 1 HYPRE_MINOR_VERSION)
+    if (HYPRE_SPLIT_VERSION GREATER 2)
+      list(GET HYPRE_SPLIT_VERSION 2 HYPRE_PATCH_VERSION)
+    else()
+      set(HYPRE_PATCH_VERSION 0)
+    endif()
+    math(EXPR HYPRE_VERSION "10000*${HYPRE_MAJOR_VERSION} + 100*${HYPRE_MINOR_VERSION} + ${HYPRE_PATCH_VERSION}")
+    set(HYPRE_VERSION ${HYPRE_VERSION} CACHE STRING "HYPRE version." FORCE)
+  endif()
   return()
 endif()
 

config/cmake/modules/MfemCmakeUtilities.cmake

@@ -932,12 +932,14 @@ function(mfem_export_mk_files)
   endif()
   set(MFEM_BUILD_TAG "${CMAKE_SYSTEM}")
   set(MFEM_PREFIX "${CMAKE_INSTALL_PREFIX}")
-  # For the next 4 variable, these are the values for the build-tree version of
+  # For the next 4 variables, these are the values for the build-tree version of
   # 'config.mk'
   set(MFEM_INC_DIR "${PROJECT_BINARY_DIR}")
   set(MFEM_LIB_DIR "${PROJECT_BINARY_DIR}")
   set(MFEM_TEST_MK "${PROJECT_SOURCE_DIR}/config/test.mk")
   set(MFEM_CONFIG_EXTRA "MFEM_BUILD_DIR ?= ${PROJECT_BINARY_DIR}")
+  # TODO: CUDA/HIP support:
+  set(MFEM_XLINKER "${CMAKE_CXX_LINKER_WRAPPER_FLAG}")
   set(MFEM_MPIEXEC ${MPIEXEC})
   if (NOT MFEM_MPIEXEC)
     set(MFEM_MPIEXEC "mpirun")

config/config.mk.in

@@ -88,6 +88,7 @@ MFEM_BUILD_TAG = @MFEM_BUILD_TAG@
 MFEM_PREFIX    = @MFEM_PREFIX@
 MFEM_INC_DIR   = @MFEM_INC_DIR@
 MFEM_LIB_DIR   = @MFEM_LIB_DIR@
+MFEM_XLINKER   = @MFEM_XLINKER@
 
 # Location of test.mk
 MFEM_TEST_MK = @MFEM_TEST_MK@

config/defaults.mk

@@ -57,7 +57,7 @@ CUDA_DIR = $(or $(CUDA_HOME),$(patsubst %/,%,$(dir \
 CLANG_CUDA_FLAGS = -xcuda --cuda-path=$(CUDA_DIR) --cuda-gpu-arch=$(CUDA_ARCH)
 # flags for nvcc
 NVCC_FLAGS = -x=cu --expt-extended-lambda --expt-relaxed-constexpr \
- -arch=$(CUDA_ARCH)
+ -arch=$(CUDA_ARCH) -isystem "$(CUDA_DIR)/include"
 # Prefixes for passing flags to the host compiler and linker when using
 # CUDA_CXX=nvcc
 CUDA_XCOMPILER = -Xcompiler=

data/mfem.mesh

@@ -0,0 +1,593 @@
+MFEM mesh v1.0
+
+# Created by: Pointwise
+
+# MFEM Geometry Types:
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+# PRISM       = 6
+
+dimension
+2
+
+elements
+160
+1 3 1 164 163 0
+1 3 164 165 162 163
+1 3 2 166 164 1
+1 3 166 132 165 164
+1 3 3 167 166 2
+1 3 167 131 132 166
+1 3 4 168 167 3
+1 3 168 130 131 167
+1 3 5 169 168 4
+1 3 169 129 130 168
+1 3 6 170 169 5
+1 3 170 128 129 169
+1 3 171 172 170 6
+1 3 172 127 128 170
+1 3 124 125 172 171
+1 3 125 126 127 172
+1 3 162 165 173 161
+1 3 165 132 133 173
+1 3 161 173 174 160
+1 3 173 133 134 174
+1 3 160 174 175 159
+1 3 174 134 135 175
+1 3 6 7 176 171
+1 3 7 8 177 176
+1 3 171 176 123 124
+1 3 176 177 122 123
+1 3 159 175 178 158
+1 3 175 135 136 178
+1 3 158 178 179 157
+1 3 178 136 137 179
+1 3 157 179 180 156
+1 3 179 137 138 180
+1 3 122 177 181 121
+1 3 177 8 182 181
+1 3 8 9 183 182
+1 3 9 10 184 183
+1 3 10 11 185 184
+1 3 11 12 186 185
+1 3 12 13 187 186
+1 3 13 14 15 187
+1 3 121 181 119 120
+1 3 181 182 118 119
+1 3 182 183 117 118
+1 3 183 184 188 117
+1 3 184 185 109 188
+1 3 185 186 108 109
+1 3 186 187 189 108
+1 3 187 15 16 189
+1 3 109 110 190 188
+1 3 110 111 191 190
+1 3 111 112 113 191
+1 3 188 190 116 117
+1 3 190 191 115 116
+1 3 191 113 114 115
+1 3 189 192 107 108
+1 3 192 193 106 107
+1 3 193 194 105 106
+1 3 194 195 104 105
+1 3 195 196 103 104
+1 3 16 17 192 189
+1 3 17 18 193 192
+1 3 18 19 194 193
+1 3 19 20 195 194
+1 3 20 21 196 195
+1 3 97 98 197 96
+1 3 98 99 198 197
+1 3 99 100 199 198
+1 3 100 101 200 199
+1 3 101 102 201 200
+1 3 102 103 202 201
+1 3 103 196 203 202
+1 3 196 21 22 203
+1 3 96 197 204 95
+1 3 197 198 39 204
+1 3 198 199 38 39
+1 3 199 200 205 38
+1 3 200 201 32 205
+1 3 201 202 31 32
+1 3 202 203 206 31
+1 3 203 22 23 206
+1 3 32 33 207 205
+1 3 33 34 35 207
+1 3 205 207 37 38
+1 3 207 35 36 37
+1 3 39 40 208 204
+1 3 40 41 209 208
+1 3 41 42 210 209
+1 3 42 43 211 210
+1 3 43 44 212 211
+1 3 204 208 94 95
+1 3 208 209 93 94
+1 3 209 210 92 93
+1 3 210 211 91 92
+1 3 211 212 90 91
+1 3 90 212 213 89
+1 3 212 44 214 213
+1 3 44 45 215 214
+1 3 45 46 216 215
+1 3 46 47 217 216
+1 3 47 48 218 217
+1 3 48 49 219 218
+1 3 49 50 51 219
+1 3 89 213 87 88
+1 3 213 214 86 87
+1 3 214 215 85 86
+1 3 215 216 84 85
+1 3 216 217 83 84
+1 3 217 218 82 83
+1 3 218 219 220 82
+1 3 219 51 52 220
+1 3 53 221 220 52
+1 3 221 81 82 220
+1 3 54 222 221 53
+1 3 222 80 81 221
+1 3 55 223 222 54
+1 3 223 79 80 222
+1 3 26 27 224 25
+1 3 27 28 29 224
+1 3 25 224 225 24
+1 3 224 29 30 225
+1 3 24 225 206 23
+1 3 225 30 31 206
+1 3 154 155 226 153
+1 3 155 156 180 226
+1 3 153 226 227 152
+1 3 226 180 138 227
+1 3 152 227 228 151
+1 3 227 138 139 228
+1 3 151 228 229 150
+1 3 228 139 140 229
+1 3 150 229 230 149
+1 3 229 140 141 230
+1 3 149 230 231 148
+1 3 230 141 142 231
+1 3 148 231 232 147
+1 3 231 142 143 232
+1 3 147 232 145 146
+1 3 232 143 144 145
+1 3 56 233 223 55
+1 3 233 78 79 223
+1 3 57 234 233 56
+1 3 234 77 78 233
+1 3 58 235 234 57
+1 3 235 76 77 234
+1 3 61 236 59 60
+1 3 236 235 58 59
+1 3 62 237 236 61
+1 3 237 76 235 236
+1 3 63 238 237 62
+1 3 238 75 76 237
+1 3 64 239 238 63
+1 3 239 74 75 238
+1 3 65 240 239 64
+1 3 240 73 74 239
+1 3 66 241 240 65
+1 3 241 72 73 240
+1 3 67 242 241 66
+1 3 242 71 72 241
+1 3 68 69 242 67
+1 3 69 70 71 242
+
+boundary
+164
+3 1 0 1
+3 1 1 2
+3 1 2 3
+3 1 3 4
+3 1 4 5
+3 1 5 6
+3 1 6 7
+3 1 7 8
+3 1 8 9
+3 1 9 10
+3 1 10 11
+3 1 11 12
+3 1 12 13
+3 1 13 14
+3 1 16 17
+3 1 17 18
+3 1 18 19
+3 1 19 20
+3 1 20 21
+3 1 21 22
+3 1 22 23
+3 1 23 24
+3 1 24 25
+3 1 25 26
+3 1 26 27
+3 1 27 28
+3 1 28 29
+3 1 29 30
+3 1 30 31
+3 1 31 32
+3 1 32 33
+3 1 33 34
+3 1 34 35
+3 1 35 36
+3 1 36 37
+3 1 37 38
+3 1 38 39
+3 1 39 40
+3 1 40 41
+3 1 41 42
+3 1 42 43
+3 1 43 44
+3 1 49 50
+3 1 48 49
+3 1 47 48
+3 1 46 47
+3 1 45 46
+3 1 44 45
+3 1 52 53
+3 1 53 54
+3 1 54 55
+3 1 57 58
+3 1 56 57
+3 1 55 56
+3 1 60 61
+3 1 61 62
+3 1 62 63
+3 1 63 64
+3 1 64 65
+3 1 65 66
+3 1 66 67
+3 1 67 68
+3 1 75 76
+3 1 74 75
+3 1 73 74
+3 1 72 73
+3 1 71 72
+3 1 70 71
+3 1 76 77
+3 1 77 78
+3 1 78 79
+3 1 81 82
+3 1 80 81
+3 1 79 80
+3 1 82 83
+3 1 83 84
+3 1 84 85
+3 1 85 86
+3 1 86 87
+3 1 87 88
+3 1 94 95
+3 1 93 94
+3 1 92 93
+3 1 91 92
+3 1 90 91
+3 1 96 97
+3 1 95 96
+3 1 97 98
+3 1 98 99
+3 1 99 100
+3 1 100 101
+3 1 101 102
+3 1 102 103
+3 1 107 108
+3 1 106 107
+3 1 105 106
+3 1 104 105
+3 1 103 104
+3 1 108 109
+3 1 109 110
+3 1 110 111
+3 1 111 112
+3 1 112 113
+3 1 113 114
+3 1 114 115
+3 1 115 116
+3 1 116 117
+3 1 119 120
+3 1 118 119
+3 1 117 118
+3 1 131 132
+3 1 130 131
+3 1 129 130
+3 1 128 129
+3 1 127 128
+3 1 126 127
+3 1 132 133
+3 1 133 134
+3 1 134 135
+3 1 137 138
+3 1 136 137
+3 1 135 136
+3 1 138 139
+3 1 139 140
+3 1 140 141
+3 1 141 142
+3 1 142 143
+3 1 143 144
+3 1 147 148
+3 1 146 147
+3 1 153 154
+3 1 152 153
+3 1 151 152
+3 1 150 151
+3 1 149 150
+3 1 148 149
+3 1 156 157
+3 1 157 158
+3 1 158 159
+3 1 161 162
+3 1 160 161
+3 1 159 160
+2 1 69 70
+2 1 68 69
+3 1 88 89
+3 1 89 90
+3 1 121 122
+3 1 120 121
+3 1 123 124
+3 1 122 123
+3 1 125 126
+3 1 124 125
+1 1 144 145
+1 1 145 146
+3 1 15 16
+3 1 14 15
+3 1 50 51
+3 1 51 52
+3 1 59 60
+3 1 58 59
+3 1 154 155
+3 1 155 156
+3 1 163 0
+3 1 162 163
+
+vertices
+243
+2
+4 4
+4 3.5
+4 3
+4 2.5
+4 2
+4 1.5
+4 1
+4.5 1
+5 1
+5 1.5
+5 2
+5 2.5
+5 3
+5 3.5
+5 4
+5.500 4
+6 4
+6.500 4
+7 4
+7.5 4
+8 4
+8.5 4
+9 4
+9.5 4
+10 4
+10.5 4
+11 4
+11 3.5
+11 3
+10.5 3
+10 3
+9.5 3
+9.5 2.5
+10 2.5
+10.5 2.5
+10.5 2
+10.5 1.5
+10 1.5
+9.5 1.5
+9.5 1
+10 1
+10.5 1
+11 1
+11.5 1
+12 1
+12 1.5
+12 2
+12 2.5
+12 3
+12 3.5
+12 4
+12.5 4
+13 4
+13.333 3.75
+13.666 3.5
+14.000 3.25
+14.333 3.5
+14.666 3.75
+15.000 4
+15.500 4
+16.000 4
+16.000 3.5
+16.000 3
+16.000 2.5
+16.000 2
+16.000 1.5
+16.000 1
+16.000 0.5
+16.000 0
+15.500 0
+15.000 0
+15.000 0.5000000000000002
+15.000 1
+15.000 1.5
+15.000 2
+15.000 2.5
+15.000 3
+14.666 2.75
+14.333 2.5
+14.000 2.25
+13.666 2.5
+13.333 2.75
+13 3
+13 2.5
+13 2
+13 1.5
+13 1
+13 0.500
+13 0
+12.5 0
+12 0
+11.5 0
+11 0
+10.5 0
+10 0
+9.5 0
+9 0
+8.5 0
+8.5 0.5
+8.5 1
+8.5 1.5
+8.5 2
+8.5 2.5
+8.5 3
+8 3
+7.5 3
+7 3
+6.500 3
+6 3
+6 2.5
+6.5 2.5
+7 2.5
+7.5 2.5
+7.5 2
+7.5 1.5
+7.000 1.5
+6.5 1.5
+6 1.5
+6 1
+6 0.5
+6 0
+5.5 0
+5 0
+4.5 0
+4 0
+3.5 0
+3 0
+3 0.500
+3 1
+3 1.5
+3 2
+3 2.5
+3 3
+2.666 2.75
+2.333 2.5
+2.000 2.25
+1.666 2.5
+1.333 2.75
+1.000 3
+1.000 2.5
+1.000 2
+1.000 1.5
+1.000 1
+1.000 0.5000
+1.000 0
+0.5000 0
+0.0000 0
+0.0000 0.5
+0.0000 1
+0.0000 1.5
+0.0000 2
+0.0000 2.5
+0.0000 3
+0.0000 3.5
+0.0000 4
+0.5000 4
+1.000 4
+1.333 3.75
+1.666 3.5
+2.000 3.25
+2.333 3.5
+2.666 3.75
+3 4
+3.5 4
+3.5 3.5
+3 3.5
+3.5 3
+3.5 2.5
+3.5 2
+3.5 1.5
+3.5 1
+4 0.5
+3.5 0.5
+2.666 3.25
+2.333 3
+2.000 2.75
+4.5 0.5
+5 0.5
+1.666 3
+1.333 3.25
+1.000 3.5
+5.5 0.5
+5.500 1
+5.500 1.5
+5.500 2
+5.500 2.5
+5.500 3
+5.500 3.5
+6 2
+6 3.5
+6.5 2
+7 2
+6.5 3.5
+7 3.5
+7.5 3.5
+8 3.5
+8.5 3.5
+9 0.5
+9 1
+9 1.5
+9 2
+9 2.5
+9 3
+9 3.5
+9.5 0.5
+9.5 2
+9.5 3.5
+10 2
+10 0.5
+10.5 0.5
+11 0.5
+11.5 0.5
+12 0.5
+12.5 0.500
+12.5 1
+12.5 1.5
+12.5 2
+12.5 2.5
+12.5 3
+12.5 3.5
+13 3.5
+13.333 3.250
+13.666 3
+14.000 2.75
+10.5 3.5
+10 3.5
+0.500 3.5
+0.500 3
+0.500 2.5
+0.500 2
+0.500 1.5
+0.500 1
+0.500 0.5
+14.333 3
+14.666 3.25
+15.000 3.5
+15.500 3.5
+15.500 3
+15.500 2.5
+15.500 2
+15.500 1.5
+15.500 1
+15.500 0.5

doc/CodeDocumentation.dox

@@ -202,6 +202,7 @@ namespace mfem {
  * - <a class="el" href="tesla_8cpp_source.html">Tesla</a>: simple magnetostatics simulation code
  * - <a class="el" href="maxwell_8cpp_source.html">Maxwell</a>: simple transient full-wave electromagnetics simulation code
  * - <a class="el" href="joule_8cpp_source.html">Joule</a>: transient magnetics and Joule heating miniapp
+ * - <a class="el" href="lorentz_8cpp_source.html">Lorentz</a>: simple particle tracking code based on the Lorentz force
  * - <a class="el" href="classmfem_1_1navier_1_1NavierSolver.html">Navier</a>: solve the transient incompressible Navier-Stokes equations
  * - <a class="el" href="mobius-strip_8cpp_source.html">Mobius Strip</a>: generate various Mobius strip-like meshes
  * - <a class="el" href="klein-bottle_8cpp_source.html">Klein Bottle</a>: generate three types of Klein bottle surfaces

fem/fe_coll.cpp

@@ -2456,8 +2456,7 @@ RT_FECollection::RT_FECollection(const int order, const int dim,
       const char *cb_name = BasisType::Name(cb_type); // this may abort
       MFEM_ABORT("unknown closed BasisType: " << cb_name);
    }
-   if (Quadrature1D::CheckOpen(op_type) == Quadrature1D::Invalid &&
-       ob_type != BasisType::IntegratedGLL)
+   if (Quadrature1D::CheckOpen(op_type) == Quadrature1D::Invalid)
    {
       const char *ob_name = BasisType::Name(ob_type); // this may abort
       MFEM_ABORT("unknown open BasisType: " << ob_name);
@@ -2784,8 +2783,7 @@ ND_FECollection::ND_FECollection(const int p, const int dim,
    int cp_type = BasisType::GetQuadrature1D(cb_type);
 
    // Error checking
-   if (Quadrature1D::CheckOpen(op_type) == Quadrature1D::Invalid &&
-       ob_type != BasisType::IntegratedGLL)
+   if (Quadrature1D::CheckOpen(op_type) == Quadrature1D::Invalid)
    {
       const char *ob_name = BasisType::Name(ob_type);
       MFEM_ABORT("Invalid open basis point type: " << ob_name);

fem/intrules.cpp

@@ -947,6 +947,7 @@ int Quadrature1D::CheckOpen(int type)
       case OpenUniform:
       case ClosedUniform:
       case OpenHalfUniform:
+      case ClosedGL:
          return type; // all types can work as open
       default:
          return Invalid;

fem/nonlinearform.cpp

@@ -436,7 +436,7 @@ void NonlinearForm::Mult(const Vector &x, Vector &y) const
    // In parallel, the result is in 'py' which is an alias for 'aux2'.
 }
 
-Operator &NonlinearForm::GetGradient(const Vector &x) const
+Operator &NonlinearForm::GetGradient(const Vector &x, bool finalize) const
 {
    if (ext)
    {
@@ -644,6 +644,8 @@ Operator &NonlinearForm::GetGradient(const Vector &x) const
       }
    }
 
+   if (!finalize) { return *Grad; }
+
    if (!Grad->Finalized())
    {
       Grad->Finalize(skip_zeros);
@@ -1203,7 +1205,14 @@ const BlockVector &BlockNonlinearForm::Prolongate(const BlockVector &bx) const
       aux1.Update(block_offsets);
       for (int s = 0; s < fes.Size(); s++)
       {
-         P[s]->Mult(bx.GetBlock(s), aux1.GetBlock(s));
+         if (P[s])
+         {
+            P[s]->Mult(bx.GetBlock(s), aux1.GetBlock(s));
+         }
+         else
+         {
+            aux1.GetBlock(s) = bx.GetBlock(s);
+         }
       }
       return aux1;
    }
@@ -1232,11 +1241,16 @@ void BlockNonlinearForm::Mult(const Vector &x, Vector &y) const
       {
          cP[s]->MultTranspose(pby.GetBlock(s), by.GetBlock(s));
       }
+      else if (needs_prolongation)
+      {
+         by.GetBlock(s) = pby.GetBlock(s);
+      }
       by.GetBlock(s).SetSubVector(*ess_tdofs[s], 0.0);
    }
 }
 
-void BlockNonlinearForm::ComputeGradientBlocked(const BlockVector &bx) const
+void BlockNonlinearForm::ComputeGradientBlocked(const BlockVector &bx,
+                                                bool finalize) const
 {
    const int skip_zeros = 0;
    Array<Array<int> *> vdofs(fes.Size());
@@ -1490,7 +1504,7 @@ void BlockNonlinearForm::ComputeGradientBlocked(const BlockVector &bx) const
       }
    }
 
-   if (!Grads(0,0)->Finalized())
+   if (finalize && !Grads(0,0)->Finalized())
    {
       for (int i=0; i<fes.Size(); ++i)
       {
@@ -1529,7 +1543,23 @@ Operator &BlockNonlinearForm::GetGradient(const Vector &x) const
          for (int s2 = 0; s2 < fes.Size(); ++s2)
          {
             delete cGrads(s1, s2);
-            cGrads(s1, s2) = RAP(*cP[s1], *Grads(s1, s2), *cP[s2]);
+            if (cP[s1] && cP[s2])
+            {
+               cGrads(s1, s2) = RAP(*cP[s1], *Grads(s1, s2), *cP[s2]);
+            }
+            else if (cP[s1])
+            {
+               cGrads(s1, s2) = TransposeMult(*cP[s1], *Grads(s1, s2));
+            }
+            else if (cP[s2])
+            {
+               cGrads(s1, s2) = mfem::Mult(*Grads(s1, s2), *cP[s2]);
+            }
+            else
+            {
+               cGrads(s1, s2) = NULL;
+               continue;
+            }
             mGrads(s1, s2) = cGrads(s1, s2);
          }
       }

fem/nonlinearform.hpp

@@ -217,7 +217,12 @@ public:
        In general, @a x may have non-homogeneous essential boundary values.
 
        The state @a x must be a true-dof vector. */
-   Operator &GetGradient(const Vector &x) const override;
+   Operator &GetGradient(const Vector &x) const override { return GetGradient(x, true); }
+
+   /** @brief Compute the gradient Operator of the NonlinearForm corresponding
+       to the state @a x with optional finalization and elimintaion. */
+   /** @see GetGradient(const Vector &) */
+   Operator &GetGradient(const Vector &x, bool finalize) const;
 
    /// Update the NonlinearForm to propagate updates of the associated FE space.
    /** After calling this method, the essential boundary conditions need to be
@@ -308,7 +313,7 @@ protected:
    void MultBlocked(const BlockVector &bx, BlockVector &by) const;
 
    /// Specialized version of GetGradient() for BlockVector
-   void ComputeGradientBlocked(const BlockVector &bx) const;
+   void ComputeGradientBlocked(const BlockVector &bx, bool finalize = true) const;
 
 public:
    /// Construct an empty BlockNonlinearForm. Initialize with SetSpaces().

fem/pbilinearform.cpp

@@ -151,6 +151,15 @@ void ParBilinearForm::ParallelRAP(SparseMatrix &loc_A, OperatorHandle &A,
    }
 }
 
+HypreParMatrix *ParBilinearForm::ParallelAssembleInternalMatrix()
+{
+   if (p_mat.Ptr() == NULL)
+   {
+      ParallelAssemble(p_mat, mat);
+   }
+   return p_mat.As<HypreParMatrix>();
+}
+
 void ParBilinearForm::ParallelAssemble(OperatorHandle &A, SparseMatrix *A_local)
 {
    A.Clear();
@@ -333,6 +342,15 @@ void ParBilinearForm
    A.EliminateRowsCols(dof_list, X, B);
 }
 
+void ParBilinearForm::ParallelEliminateEssentialBC(
+   const Array<int> &bdr_attr_is_ess, const HypreParVector &X, HypreParVector &B)
+{
+   Array<int> dof_list;
+   pfes->GetEssentialTrueDofs(bdr_attr_is_ess, dof_list);
+
+   p_mat.As<HypreParMatrix>()->EliminateRowsCols(dof_list, X, B);
+}
+
 HypreParMatrix *ParBilinearForm::
 ParallelEliminateEssentialBC(const Array<int> &bdr_attr_is_ess,
                              HypreParMatrix &A) const
@@ -344,6 +362,26 @@ ParallelEliminateEssentialBC(const Array<int> &bdr_attr_is_ess,
    return A.EliminateRowsCols(dof_list);
 }
 
+void ParBilinearForm::ParallelEliminateEssentialBC(const Array<int>
+                                                   &bdr_attr_is_ess)
+{
+   Array<int> tdofs_list;
+   pfes->GetEssentialTrueDofs(bdr_attr_is_ess, tdofs_list);
+
+   ParallelEliminateTDofs(tdofs_list);
+}
+
+void ParBilinearForm::ParallelEliminateTDofs(const Array<int> &tdofs_list)
+{
+   p_mat_e.EliminateRowsCols(p_mat, tdofs_list);
+}
+
+void ParBilinearForm::ParallelEliminateTDofsInRHS(
+   const Array<int> &tdofs_list, const Vector &x, Vector &b)
+{
+   p_mat.EliminateBC(p_mat_e, tdofs_list, x, b);
+}
+
 void ParBilinearForm::TrueAddMult(const Vector &x, Vector &y, const real_t a)
 const
 {
@@ -485,7 +523,7 @@ void ParBilinearForm::FormLinearSystem(
       HypreParVector true_X(pfes), true_B(pfes);
       P.MultTranspose(b, true_B);
       R.Mult(x, true_X);
-      p_mat.EliminateBC(p_mat_e, ess_tdof_list, true_X, true_B);
+      ParallelEliminateTDofsInRHS(ess_tdof_list, true_X, true_B);
       R.MultTranspose(true_B, b);
       hybridization->ReduceRHS(true_B, B);
       X.SetSize(B.Size());
@@ -498,17 +536,11 @@ void ParBilinearForm::FormLinearSystem(
       B.SetSize(X.Size());
       P.MultTranspose(b, B);
       R.Mult(x, X);
-      p_mat.EliminateBC(p_mat_e, ess_tdof_list, X, B);
+      ParallelEliminateTDofsInRHS(ess_tdof_list, X, B);
       if (!copy_interior) { X.SetSubVectorComplement(ess_tdof_list, 0.0); }
    }
 }
 
-void ParBilinearForm::EliminateVDofsInRHS(
-   const Array<int> &vdofs, const Vector &x, Vector &b)
-{
-   p_mat.EliminateBC(p_mat_e, vdofs, x, b);
-}
-
 void ParBilinearForm::FormSystemMatrix(const Array<int> &ess_tdof_list,
                                        OperatorHandle &A)
 {
@@ -553,7 +585,7 @@ void ParBilinearForm::FormSystemMatrix(const Array<int> &ess_tdof_list,
          mat = NULL;
          delete mat_e;
          mat_e = NULL;
-         p_mat_e.EliminateRowsCols(p_mat, ess_tdof_list);
+         ParallelEliminateTDofs(ess_tdof_list);
       }
       if (hybridization)
       {
@@ -615,36 +647,180 @@ void ParBilinearForm::Update(FiniteElementSpace *nfes)
    p_mat_e.Clear();
 }
 
-
-HypreParMatrix *ParMixedBilinearForm::ParallelAssemble()
+void ParMixedBilinearForm::pAllocMat()
 {
-   // construct the block-diagonal matrix A
-   HypreParMatrix *A =
-      new HypreParMatrix(trial_pfes->GetComm(),
-                         test_pfes->GlobalVSize(),
-                         trial_pfes->GlobalVSize(),
-                         test_pfes->GetDofOffsets(),
-                         trial_pfes->GetDofOffsets(),
-                         mat);
+   const int trial_nbr_size = trial_pfes->GetFaceNbrVSize();
+   const int test_nbr_size = test_pfes->GetFaceNbrVSize();
 
-   HypreParMatrix *rap = RAP(test_pfes->Dof_TrueDof_Matrix(), A,
-                             trial_pfes->Dof_TrueDof_Matrix());
-
-   delete A;
-
-   return rap;
+   if (keep_nbr_block)
+   {
+      mat = new SparseMatrix(height + test_nbr_size, width + trial_nbr_size);
+   }
+   else
+   {
+      mat = new SparseMatrix(height, width + trial_nbr_size);
+   }
 }
 
-void ParMixedBilinearForm::ParallelAssemble(OperatorHandle &A)
+void ParMixedBilinearForm::AssembleSharedFaces(int skip_zeros)
 {
-   // construct the rectangular block-diagonal matrix dA
-   OperatorHandle dA(A.Type());
-   dA.MakeRectangularBlockDiag(trial_pfes->GetComm(),
-                               test_pfes->GlobalVSize(),
-                               trial_pfes->GlobalVSize(),
-                               test_pfes->GetDofOffsets(),
-                               trial_pfes->GetDofOffsets(),
-                               mat);
+   ParMesh *pmesh = trial_pfes->GetParMesh();
+   FaceElementTransformations *T;
+   Array<int> tr_vdofs1, tr_vdofs2, tr_vdofs_all;
+   Array<int> te_vdofs1, te_vdofs2, te_vdofs_all;
+   DenseMatrix elemmat;
+
+   int nfaces = pmesh->GetNSharedFaces();
+   for (int i = 0; i < nfaces; i++)
+   {
+      T = pmesh->GetSharedFaceTransformations(i);
+      int Elem2NbrNo = T->Elem2No - pmesh->GetNE();
+      trial_pfes->GetElementVDofs(T->Elem1No, tr_vdofs1);
+      test_pfes->GetElementVDofs(T->Elem1No, te_vdofs1);
+      trial_pfes->GetFaceNbrElementVDofs(Elem2NbrNo, tr_vdofs2);
+      test_pfes->GetFaceNbrElementVDofs(Elem2NbrNo, te_vdofs2);
+
+      tr_vdofs1.Copy(tr_vdofs_all);
+      for (int j = 0; j < tr_vdofs2.Size(); j++)
+      {
+         if (tr_vdofs2[j] >= 0)
+         {
+            tr_vdofs2[j] += width;
+         }
+         else
+         {
+            tr_vdofs2[j] -= width;
+         }
+      }
+      tr_vdofs_all.Append(tr_vdofs2);
+
+      if (keep_nbr_block)
+      {
+         te_vdofs1.Copy(te_vdofs_all);
+         for (int j = 0; j < te_vdofs2.Size(); j++)
+         {
+            if (te_vdofs2[j] >= 0)
+            {
+               te_vdofs2[j] += height;
+            }
+            else
+            {
+               te_vdofs2[j] -= height;
+            }
+         }
+         te_vdofs_all.Append(te_vdofs2);
+      }
+
+      for (int k = 0; k < interior_face_integs.Size(); k++)
+      {
+         interior_face_integs[k]->
+         AssembleFaceMatrix(*trial_pfes->GetFE(T->Elem1No),
+                            *test_pfes->GetFE(T->Elem1No),
+                            *trial_pfes->GetFaceNbrFE(Elem2NbrNo),
+                            *test_pfes->GetFaceNbrFE(Elem2NbrNo),
+                            *T, elemmat);
+         if (keep_nbr_block)
+         {
+            mat->AddSubMatrix(te_vdofs_all, tr_vdofs_all, elemmat, skip_zeros);
+         }
+         else
+         {
+            mat->AddSubMatrix(te_vdofs1, tr_vdofs_all, elemmat, skip_zeros);
+         }
+      }
+   }
+}
+
+void ParMixedBilinearForm::Assemble(int skip_zeros)
+{
+   if (interior_face_integs.Size())
+   {
+      trial_pfes->ExchangeFaceNbrData();
+      test_pfes->ExchangeFaceNbrData();
+      if (!ext && mat == NULL)
+      {
+         pAllocMat();
+      }
+   }
+
+   MixedBilinearForm::Assemble(skip_zeros);
+
+   if (!ext && interior_face_integs.Size() > 0)
+   {
+      AssembleSharedFaces(skip_zeros);
+   }
+}
+
+HypreParMatrix *ParMixedBilinearForm::ParallelAssembleInternalMatrix()
+{
+   if (p_mat.Ptr() == NULL)
+   {
+      ParallelAssemble(p_mat, mat);
+   }
+   return p_mat.As<HypreParMatrix>();
+}
+
+HypreParMatrix *ParMixedBilinearForm::ParallelAssemble(SparseMatrix *m)
+{
+   OperatorHandle Mh(Operator::Hypre_ParCSR);
+   ParallelAssemble(Mh, m);
+   Mh.SetOperatorOwner(false);
+   return Mh.As<HypreParMatrix>();
+}
+
+void ParMixedBilinearForm::ParallelAssemble(OperatorHandle &A,
+                                            SparseMatrix *A_local)
+{
+   A.Clear();
+
+   if (A_local == NULL) { return; }
+   MFEM_VERIFY(A_local->Finalized(), "the local matrix must be finalized");
+
+   OperatorHandle dA(A.Type()), hdA;
+
+   if (interior_face_integs.Size() == 0)
+   {
+      // construct the rectangular block-diagonal matrix dA
+      dA.MakeRectangularBlockDiag(trial_pfes->GetComm(),
+                                  test_pfes->GlobalVSize(),
+                                  trial_pfes->GlobalVSize(),
+                                  test_pfes->GetDofOffsets(),
+                                  trial_pfes->GetDofOffsets(),
+                                  A_local);
+   }
+   else
+   {
+      // handle the case when 'a' contains off-diagonal
+      const int lvrows = test_pfes->GetVSize();
+      const int lvcols = trial_pfes->GetVSize();
+      const HYPRE_BigInt *face_nbr_glob_lcol = trial_pfes->GetFaceNbrGlobalDofMap();
+      const HYPRE_BigInt lcol_offset = trial_pfes->GetMyDofOffset();
+
+      Array<HYPRE_BigInt> glob_J(A_local->NumNonZeroElems());
+      const int *J = A_local->GetJ();
+      for (int i = 0; i < glob_J.Size(); i++)
+      {
+         if (J[i] < lvcols)
+         {
+            glob_J[i] = J[i] + lcol_offset;
+         }
+         else
+         {
+            glob_J[i] = face_nbr_glob_lcol[J[i] - lvcols];
+         }
+      }
+
+      // TODO - construct dA directly in the A format
+      hdA.Reset(
+         new HypreParMatrix(trial_pfes->GetComm(), lvrows, test_pfes->GlobalVSize(),
+                            trial_pfes->GlobalVSize(), A_local->GetI(), glob_J,
+                            A_local->GetData(), test_pfes->GetDofOffsets(),
+                            trial_pfes->GetDofOffsets()));
+      // - hdA owns the new HypreParMatrix
+      // - the above constructor copies all input arrays
+      glob_J.DeleteAll();
+      dA.ConvertFrom(hdA);
+   }
 
    OperatorHandle P_test(A.Type()), P_trial(A.Type());
 
@@ -670,6 +846,44 @@ void ParMixedBilinearForm::TrueAddMult(const Vector &x, Vector &y,
    test_pfes->Dof_TrueDof_Matrix()->MultTranspose(a, Yaux, 1.0, y);
 }
 
+void ParMixedBilinearForm::ParallelEliminateTrialEssentialBC(
+   const Array<int> &bdr_attr_is_ess)
+{
+   Array<int> trial_tdof_list;
+   trial_pfes->GetEssentialTrueDofs(bdr_attr_is_ess, trial_tdof_list);
+
+   ParallelEliminateTrialTDofs(trial_tdof_list);
+}
+
+void ParMixedBilinearForm::ParallelEliminateTrialTDofs(
+   const Array<int> &trial_tdof_list)
+{
+   HypreParMatrix *temp = p_mat.As<HypreParMatrix>()->EliminateCols(
+                             trial_tdof_list);
+   p_mat_e.Reset(temp, true);
+}
+
+void ParMixedBilinearForm::ParallelEliminateTrialTDofsInRHS(
+   const Array<int> &trial_tdof_list, const Vector &x, Vector &b)
+{
+   p_mat_e.As<HypreParMatrix>()->Mult(-1.0, x, 1.0, b);
+}
+
+void ParMixedBilinearForm::ParallelEliminateTestEssentialBC(
+   const Array<int> &bdr_attr_is_ess)
+{
+   Array<int> test_tdof_list;
+   test_pfes->GetEssentialTrueDofs(bdr_attr_is_ess, test_tdof_list);
+
+   ParallelEliminateTestTDofs(test_tdof_list);
+}
+
+void ParMixedBilinearForm::ParallelEliminateTestTDofs(
+   const Array<int> &test_tdof_list)
+{
+   p_mat.As<HypreParMatrix>()->EliminateRows(test_tdof_list);
+}
+
 void ParMixedBilinearForm::FormRectangularSystemMatrix(
    const Array<int>
    &trial_tdof_list,
@@ -690,10 +904,8 @@ void ParMixedBilinearForm::FormRectangularSystemMatrix(
       mat = NULL;
       delete mat_e;
       mat_e = NULL;
-      HypreParMatrix *temp =
-         p_mat.As<HypreParMatrix>()->EliminateCols(trial_tdof_list);
-      p_mat.As<HypreParMatrix>()->EliminateRows(test_tdof_list);
-      p_mat_e.Reset(temp, true);
+      ParallelEliminateTrialTDofs(trial_tdof_list);
+      ParallelEliminateTestTDofs(test_tdof_list);
    }
 
    A = p_mat;
@@ -723,7 +935,7 @@ void ParMixedBilinearForm::FormRectangularLinearSystem(
    test_P->MultTranspose(b, B);
    trial_R->Mult(x, X);
 
-   p_mat_e.As<HypreParMatrix>()->Mult(-1.0, X, 1.0, B);
+   ParallelEliminateTrialTDofsInRHS(trial_tdof_list, X, B);
    B.SetSubVector(test_tdof_list, 0.0);
 }
 

fem/pbilinearform.hpp

@@ -73,7 +73,7 @@ public:
    /** When set to true and the ParBilinearForm has interior face integrators,
        the local SparseMatrix will include the rows (in addition to the columns)
        corresponding to face-neighbor dofs. The default behavior is to disregard
-       those rows. Must be called before the first Assemble call. */
+       those rows. Must be called before the first Assemble() call. */
    void KeepNbrBlock(bool knb = true) { keep_nbr_block = knb; }
 
    /** @brief Set the operator type id for the parallel matrix/operator when
@@ -101,6 +101,14 @@ public:
        diagonal for this case. */
    void AssembleDiagonal(Vector &diag) const override;
 
+   /// Returns the matrix assembled on the true dofs, i.e. P^t A P.
+   /** The returned matrix is the internal one, owned by the form. It is not
+       reassembled if it has been already constructed. If FormSystemMatrix()
+       has been called before, it is the system matrix with eliminated
+       essential DOFs, otherwise the parallel matrix is assembled here without
+       the elimination process. */
+   HypreParMatrix *ParallelAssembleInternalMatrix();
+
    /// Returns the matrix assembled on the true dofs, i.e. P^t A P.
    /** The returned matrix has to be deleted by the caller. */
    HypreParMatrix *ParallelAssemble() { return ParallelAssemble(mat); }
@@ -146,6 +154,13 @@ public:
                                      const HypreParVector &X,
                                      HypreParVector &B) const;
 
+   /// Eliminate essential boundary DOFs from the parallel system matrix.
+   /** The array @a bdr_attr_is_ess marks boundary attributes that constitute
+       the essential part of the boundary. */
+   void ParallelEliminateEssentialBC(const Array<int> &bdr_attr_is_ess,
+                                     const HypreParVector &X,
+                                     HypreParVector &B);
+
    /// Eliminate essential boundary DOFs from a parallel assembled matrix @a A.
    /** The array @a bdr_attr_is_ess marks boundary attributes that constitute
        the essential part of the boundary. The eliminated part is stored in a
@@ -157,6 +172,12 @@ public:
    HypreParMatrix *ParallelEliminateEssentialBC(const Array<int> &bdr_attr_is_ess,
                                                 HypreParMatrix &A) const;
 
+   /// Eliminate essential boundary DOFs from the parallel system matrix.
+   /** The array @a bdr_attr_is_ess marks boundary attributes that constitute
+       the essential part of the boundary. This method relies on
+       ParallelEliminateTDofs(const Array<int> &), see it for details. */
+   void ParallelEliminateEssentialBC(const Array<int> &bdr_attr_is_ess);
+
    /// Eliminate essential true DOFs from a parallel assembled matrix @a A.
    /** Given a list of essential true dofs and the parallel assembled matrix
        @a A, eliminate the true dofs from the matrix, storing the eliminated
@@ -169,6 +190,28 @@ public:
                                           HypreParMatrix &A) const
    { return A.EliminateRowsCols(tdofs_list); }
 
+   /// Eliminate essential true DOFs from the parallel system matrix.
+   /** Given a list of essential true dofs, eliminate the true dofs from
+       the parallel assembled system matrix, storing the eliminated part
+       internally. This method works in conjunction with
+       ParallelEliminateTDofsInRHS() and allows elimination of boundary
+       conditions in multiple right-hand sides. */
+   void ParallelEliminateTDofs(const Array<int> &tdofs_list);
+
+   /** @brief Use the stored eliminated part of the parallel system matrix for
+       elimination of boundary conditions in the r.h.s. */
+   /** Given a list of essential true dofs, eliminate the true dofs from the
+       right-hand side @a b using the solution vector @a x and the previously
+       stored eliminated part of the parallel assembled system matrix produced
+       by ParallelEliminateTDofs(const Array<int> &). */
+   void ParallelEliminateTDofsInRHS(const Array<int> &tdofs, const Vector &x,
+                                    Vector &b);
+
+   /// @deprecated Use ParallelEliminateTDofsInRHS() instead.
+   MFEM_DEPRECATED void EliminateVDofsInRHS(const Array<int> &vdofs,
+                                            const Vector &x, Vector &b)
+   { ParallelEliminateTDofsInRHS(vdofs, x, b); }
+
    /** @brief Compute @a y += @a a (P^t A P) @a x, where @a x and @a y are
        vectors on the true dofs. */
    void TrueAddMult(const Vector &x, Vector &y, const real_t a = 1.0) const;
@@ -238,8 +281,6 @@ public:
 
    void Update(FiniteElementSpace *nfes = NULL) override;
 
-   void EliminateVDofsInRHS(const Array<int> &vdofs, const Vector &x, Vector &b);
-
    virtual ~ParBilinearForm() { }
 };
 
@@ -257,6 +298,13 @@ protected:
    /// Matrix and eliminated matrix
    OperatorHandle p_mat, p_mat_e;
 
+   bool keep_nbr_block;
+
+   // Allocate mat - called when (mat == NULL && fbfi.Size() > 0)
+   void pAllocMat();
+
+   void AssembleSharedFaces(int skip_zeros = 1);
+
 private:
    /// Copy construction is not supported; body is undefined.
    ParMixedBilinearForm(const ParMixedBilinearForm &);
@@ -276,6 +324,7 @@ public:
    {
       trial_pfes = trial_fes;
       test_pfes  = test_fes;
+      keep_nbr_block = false;
    }
 
    /** @brief Create a ParMixedBilinearForm on the given FiniteElementSpace%s
@@ -295,15 +344,89 @@ public:
    {
       trial_pfes = trial_fes;
       test_pfes  = test_fes;
+      keep_nbr_block = false;
    }
 
+   /** When set to true and the ParMixedBilinearForm has interior face
+       integrators, the local SparseMatrix will include the rows (in addition
+       to the columns) corresponding to face-neighbor dofs. The default
+       behavior is to disregard those rows. Must be called before the first
+       Assemble() call. */
+   void KeepNbrBlock(bool knb = true) { keep_nbr_block = knb; }
+
+   /// Assemble the local matrix
+   void Assemble(int skip_zeros = 1);
+
    /// Returns the matrix assembled on the true dofs, i.e. P_test^t A P_trial.
-   HypreParMatrix *ParallelAssemble();
+   /** The returned matrix is the internal one, owned by the form. It is not
+       reassembled if it has been already constructed. If
+       FormRectangularSystemMatrix() has been called before, it is the system
+       matrix with eliminated essential DOFs, otherwise the parallel matrix is
+       assembled here without the elimination process. */
+   HypreParMatrix *ParallelAssembleInternalMatrix();
+
+   /// Returns the matrix assembled on the true dofs, i.e. P_test^t A P_trial.
+   /** The returned matrix has to be deleted by the caller. */
+   HypreParMatrix *ParallelAssemble() { return ParallelAssemble(mat); }
+
+   /** @brief Returns the eliminated matrix assembled on the true dofs, i.e.
+       P_test^t A_local P_trial. */
+   /** The returned matrix has to be deleted by the caller. */
+   HypreParMatrix *ParallelAssembleElim() { return ParallelAssemble(mat_e); }
+
+   /** @brief Return the matrix @a m assembled on the true dofs, i.e. P_test^t
+       A_local P_trial. */
+   /** The returned matrix has to be deleted by the caller. */
+   HypreParMatrix *ParallelAssemble(SparseMatrix *m);
 
    /** @brief Returns the matrix assembled on the true dofs, i.e.
        @a A = P_test^t A_local P_trial, in the format (type id) specified by
        @a A. */
-   void ParallelAssemble(OperatorHandle &A);
+   void ParallelAssemble(OperatorHandle &A) { ParallelAssemble(A, mat); }
+
+   /** Returns the eliminated matrix assembled on the true dofs, i.e.
+       @a A_elim = P^t A_elim_local P in the format (type id) specified by @a A.
+    */
+   void ParallelAssembleElim(OperatorHandle &A_elim)
+   { ParallelAssemble(A_elim, mat_e); }
+
+   /** Returns the matrix @a A_local assembled on the true dofs, i.e.
+       @a A = P_test^t A_local P_trial in the format (type id) specified by
+       @a A. */
+   void ParallelAssemble(OperatorHandle &A, SparseMatrix *A_local);
+
+   /// Eliminate essential boundary trial DOFs from the parallel system matrix.
+   /** The array @a bdr_attr_is_ess marks boundary attributes that constitute
+       the essential part of the boundary. This method relies on
+       ParallelEliminateTrialTDofs(const Array<int> &), see it for details. */
+   void ParallelEliminateTrialEssentialBC(const Array<int> &bdr_attr_is_ess);
+
+   /// Eliminate essential trial true DOFs from the parallel system matrix.
+   /** Given a list of essential trial true dofs, eliminate the trial true dofs
+       from the parallel assembled system matrix, storing the eliminated part
+       internally. This method works in conjunction with
+       ParallelEliminateTrialTDofsInRHS() and allows elimination of boundary
+       conditions in multiple right-hand sides. */
+   void ParallelEliminateTrialTDofs(const Array<int> &trial_tdof_list);
+
+   /** @brief Use the stored eliminated part of the parallel system matrix for
+       elimination of boundary conditions in the r.h.s. */
+   /** Given a list of essential trial true dofs, eliminate the trial true dofs
+       from the right-hand side @a B using the solution vector @a X and the
+       previously stored eliminated part of the parallel assembled system
+       matrix produced by ParallelEliminateTrialTDofs(const Array<int> &). */
+   void ParallelEliminateTrialTDofsInRHS(const Array<int> &trial_tdof_list,
+                                         const Vector &X, Vector &B);
+
+   /// Eliminate essential boundary test DOFs from the parallel system matrix.
+   /** The array @a bdr_attr_is_ess marks boundary attributes that constitute
+       the essential part of the boundary. */
+   void ParallelEliminateTestEssentialBC(const Array<int> &bdr_attr_is_ess);
+
+   /// Eliminate essential test true DOFs from the parallel system matrix.
+   /** Given a list of essential test true dofs, eliminate the test true dofs
+       from the parallel assembled system matrix. */
+   void ParallelEliminateTestTDofs(const Array<int> &test_tdof_list);
 
    using MixedBilinearForm::FormRectangularSystemMatrix;
    using MixedBilinearForm::FormRectangularLinearSystem;

fem/pnonlinearform.cpp

@@ -105,6 +105,59 @@ const SparseMatrix &ParNonlinearForm::GetLocalGradient(const Vector &x) const
    return *Grad;
 }
 
+void ParNonlinearForm::GradientSharedFaces(const Vector &x,
+                                           int skip_zeros) const
+{
+   ParFiniteElementSpace *pfes = ParFESpace();
+   ParMesh *pmesh = pfes->GetParMesh();
+   FaceElementTransformations *T;
+   Array<int> vdofs1, vdofs2, vdofs_all;
+   DenseMatrix elemmat;
+   Vector el_x, nbr_x, face_x;
+   const Vector &px = Prolongate(x);
+
+   ParGridFunction pgf(pfes, const_cast<Vector&>(px), 0);
+   pgf.ExchangeFaceNbrData();
+
+   int nfaces = pmesh->GetNSharedFaces();
+   for (int i = 0; i < nfaces; i++)
+   {
+      T = pmesh->GetSharedFaceTransformations(i);
+      int Elem2NbrNo = T->Elem2No - pmesh->GetNE();
+
+      pfes->GetElementVDofs(T->Elem1No, vdofs1);
+      pfes->GetFaceNbrElementVDofs(Elem2NbrNo, vdofs2);
+      face_x.SetSize(vdofs1.Size() + vdofs2.Size());
+
+      el_x.MakeRef(face_x, 0, vdofs1.Size());
+      pgf.GetSubVector(vdofs1, el_x);
+
+      nbr_x.MakeRef(face_x, vdofs1.Size(), vdofs2.Size());
+      pgf.FaceNbrData().GetSubVector(vdofs2, nbr_x);
+
+      vdofs1.Copy(vdofs_all);
+      for (int j = 0; j < vdofs2.Size(); j++)
+      {
+         if (vdofs2[j] >= 0)
+         {
+            vdofs2[j] += height;
+         }
+         else
+         {
+            vdofs2[j] -= height;
+         }
+      }
+      vdofs_all.Append(vdofs2);
+      for (int k = 0; k < fnfi.Size(); k++)
+      {
+         fnfi[k]->AssembleFaceGrad(*pfes->GetFE(T->Elem1No),
+                                   *pfes->GetFaceNbrFE(Elem2NbrNo),
+                                   *T, face_x, elemmat);
+         Grad->AddSubMatrix(vdofs1, vdofs_all, elemmat, skip_zeros);
+      }
+   }
+}
+
 Operator &ParNonlinearForm::GetGradient(const Vector &x) const
 {
    if (NonlinearForm::ext) { return NonlinearForm::GetGradient(x); }
@@ -112,19 +165,61 @@ Operator &ParNonlinearForm::GetGradient(const Vector &x) const
    ParFiniteElementSpace *pfes = ParFESpace();
 
    pGrad.Clear();
+   OperatorHandle dA(pGrad.Type()), Ph(pGrad.Type()), hdA;
 
-   NonlinearForm::GetGradient(x); // (re)assemble Grad, no b.c.
-
-   OperatorHandle dA(pGrad.Type()), Ph(pGrad.Type());
-
-   if (fnfi.Size() == 0)
+   if (fnfi.Size())
    {
-      dA.MakeSquareBlockDiag(pfes->GetComm(), pfes->GlobalVSize(),
-                             pfes->GetDofOffsets(), Grad);
+      const int skip_zeros = 0;
+
+      pfes->ExchangeFaceNbrData();
+      if (Grad == NULL)
+      {
+         int nbr_size = pfes->GetFaceNbrVSize();
+         Grad = new SparseMatrix(pfes->GetVSize(), pfes->GetVSize() + nbr_size);
+      }
+
+      NonlinearForm::GetGradient(x, false); // (re)assemble Grad, no b.c.
+
+      GradientSharedFaces(x, skip_zeros);
+
+      Grad->Finalize(skip_zeros);
+
+      // handle the case when 'a' contains off-diagonal
+      int lvsize = pfes->GetVSize();
+      const HYPRE_BigInt *face_nbr_glob_ldof = pfes->GetFaceNbrGlobalDofMap();
+      HYPRE_BigInt ldof_offset = pfes->GetMyDofOffset();
+
+      Array<HYPRE_BigInt> glob_J(Grad->NumNonZeroElems());
+      int *J = Grad->GetJ();
+      for (int i = 0; i < glob_J.Size(); i++)
+      {
+         if (J[i] < lvsize)
+         {
+            glob_J[i] = J[i] + ldof_offset;
+         }
+         else
+         {
+            glob_J[i] = face_nbr_glob_ldof[J[i] - lvsize];
+         }
+      }
+
+      // TODO - construct dA directly in the A format
+      hdA.Reset(
+         new HypreParMatrix(pfes->GetComm(), lvsize, pfes->GlobalVSize(),
+                            pfes->GlobalVSize(), Grad->GetI(), glob_J,
+                            Grad->GetData(), pfes->GetDofOffsets(),
+                            pfes->GetDofOffsets()));
+      // - hdA owns the new HypreParMatrix
+      // - the above constructor copies all input arrays
+      glob_J.DeleteAll();
+      dA.ConvertFrom(hdA);
    }
    else
    {
-      MFEM_ABORT("TODO: assemble contributions from shared face terms");
+      NonlinearForm::GetGradient(x); // (re)assemble Grad, no b.c.
+
+      dA.MakeSquareBlockDiag(pfes->GetComm(), pfes->GlobalVSize(),
+                             pfes->GetDofOffsets(), Grad);
    }
 
    // RAP the local gradient dA.
@@ -271,7 +366,70 @@ void ParBlockNonlinearForm::Mult(const Vector &x, Vector &y) const
 
    if (fnfi.Size() > 0)
    {
-      MFEM_ABORT("TODO: assemble contributions from shared face terms");
+      // Terms over shared interior faces in parallel.
+      ParMesh *pmesh = ParFESpace(0)->GetParMesh();
+      FaceElementTransformations *tr;
+
+      Array<Array<int> *>vdofs(fes.Size());
+      Array<Array<int> *>vdofs2(fes.Size());
+      Array<Vector *> el_x(fes.Size());
+      Array<const Vector *> el_x_const(fes.Size());
+      Array<Vector *> el_y(fes.Size());
+      Array<const FiniteElement *> fe(fes.Size());
+      Array<const FiniteElement *> fe2(fes.Size());
+      Array<ParGridFunction *> pgfs(fes.Size());
+      for (int s=0; s<fes.Size(); ++s)
+      {
+         el_x_const[s] = el_x[s] = new Vector();
+         el_y[s] = new Vector();
+         vdofs[s] = new Array<int>;
+         vdofs2[s] = new Array<int>;
+         pgfs[s] = new ParGridFunction(const_cast<ParFiniteElementSpace*>(ParFESpace(s)),
+                                       xs.GetBlock(s));
+         pgfs[s]->ExchangeFaceNbrData();
+      }
+
+      const int n_shared_faces = pmesh->GetNSharedFaces();
+      for (int i = 0; i < n_shared_faces; i++)
+      {
+         tr = pmesh->GetSharedFaceTransformations(i, true);
+         int Elem2NbrNo = tr->Elem2No - pmesh->GetNE();
+
+         for (int s=0; s<fes.Size(); ++s)
+         {
+            const ParFiniteElementSpace *pfes = ParFESpace(s);
+            fe[s] = pfes->GetFE(tr->Elem1No);
+            fe2[s] = pfes->GetFaceNbrFE(Elem2NbrNo);
+
+            pfes->GetElementVDofs(tr->Elem1No, *(vdofs[s]));
+            pfes->GetFaceNbrElementVDofs(Elem2NbrNo, *(vdofs2[s]));
+
+            el_x[s]->SetSize(vdofs[s]->Size() + vdofs2[s]->Size());
+            xs.GetBlock(s).GetSubVector(*(vdofs[s]), el_x[s]->GetData());
+            pgfs[s]->FaceNbrData().GetSubVector(*(vdofs2[s]),
+                                                el_x[s]->GetData() + vdofs[s]->Size());
+         }
+
+         for (int k = 0; k < fnfi.Size(); ++k)
+         {
+            fnfi[k]->AssembleFaceVector(fe, fe2, *tr, el_x_const, el_y);
+
+            for (int s=0; s<fes.Size(); ++s)
+            {
+               if (el_y[s]->Size() == 0) { continue; }
+               ys.GetBlock(s).AddElementVector(*(vdofs[s]), *el_y[s]);
+            }
+         }
+      }
+
+      for (int s=0; s<fes.Size(); ++s)
+      {
+         delete pgfs[s];
+         delete vdofs2[s];
+         delete vdofs[s];
+         delete el_y[s];
+         delete el_x[s];
+      }
    }
 
    for (int s=0; s<fes.Size(); ++s)
@@ -328,6 +486,106 @@ void ParBlockNonlinearForm::SetGradientType(Operator::Type tid)
    }
 }
 
+void ParBlockNonlinearForm::GradientSharedFaces(const BlockVector &xs,
+                                                int skip_zeros) const
+{
+   // Terms over shared interior faces in parallel.
+   ParMesh *pmesh = ParFESpace(0)->GetParMesh();
+   FaceElementTransformations *tr;
+
+   Array<Array<int> *>vdofs(fes.Size());
+   Array<Array<int> *>vdofs2(fes.Size());
+   Array<Array<int> *>vdofs_all(fes.Size());
+   Array<Vector *> el_x(fes.Size());
+   Array<const Vector *> el_x_const(fes.Size());
+   Array2D<DenseMatrix *> elmats(fes.Size(), fes.Size());
+   Array<const FiniteElement *> fe(fes.Size());
+   Array<const FiniteElement *> fe2(fes.Size());
+   Array<ParGridFunction *> pgfs(fes.Size());
+
+   for (int s1=0; s1<fes.Size(); ++s1)
+   {
+      el_x_const[s1] = el_x[s1] = new Vector();
+      vdofs[s1] = new Array<int>;
+      vdofs2[s1] = new Array<int>;
+      vdofs_all[s1] = new Array<int>;
+      pgfs[s1] = new ParGridFunction(
+         const_cast<ParFiniteElementSpace*>(ParFESpace(s1)),
+         const_cast<Vector&>(xs.GetBlock(s1)));
+      pgfs[s1]->ExchangeFaceNbrData();
+      for (int s2=0; s2<fes.Size(); ++s2)
+      {
+         elmats(s1,s2) = new DenseMatrix();
+      }
+   }
+
+   const int n_shared_faces = pmesh->GetNSharedFaces();
+   for (int i = 0; i < n_shared_faces; i++)
+   {
+      tr = pmesh->GetSharedFaceTransformations(i, true);
+      int Elem2NbrNo = tr->Elem2No - pmesh->GetNE();
+
+      for (int s=0; s<fes.Size(); ++s)
+      {
+         const ParFiniteElementSpace *pfes = ParFESpace(s);
+         fe[s] = pfes->GetFE(tr->Elem1No);
+         fe2[s] = pfes->GetFaceNbrFE(Elem2NbrNo);
+
+         pfes->GetElementVDofs(tr->Elem1No, *(vdofs[s]));
+         pfes->GetFaceNbrElementVDofs(Elem2NbrNo, *(vdofs2[s]));
+
+         el_x[s]->SetSize(vdofs[s]->Size() + vdofs2[s]->Size());
+         xs.GetBlock(s).GetSubVector(*(vdofs[s]), el_x[s]->GetData());
+         pgfs[s]->FaceNbrData().GetSubVector(*(vdofs2[s]),
+                                             el_x[s]->GetData() + vdofs[s]->Size());
+
+         vdofs[s]->Copy(*vdofs_all[s]);
+
+         const int lvsize = pfes->GetVSize();
+         for (int j = 0; j < vdofs2[s]->Size(); j++)
+         {
+            if ((*vdofs2[s])[j] >= 0)
+            {
+               (*vdofs2[s])[j] += lvsize;
+            }
+            else
+            {
+               (*vdofs2[s])[j] -= lvsize;
+            }
+         }
+         vdofs_all[s]->Append(*(vdofs2[s]));
+      }
+
+      for (int k = 0; k < fnfi.Size(); ++k)
+      {
+         fnfi[k]->AssembleFaceGrad(fe, fe2, *tr, el_x_const, elmats);
+
+         for (int s1=0; s1<fes.Size(); ++s1)
+         {
+            for (int s2=0; s2<fes.Size(); ++s2)
+            {
+               if (elmats(s1,s2)->Height() == 0) { continue; }
+               Grads(s1,s2)->AddSubMatrix(*vdofs[s1], *vdofs_all[s2],
+                                          *elmats(s1,s2), skip_zeros);
+            }
+         }
+      }
+   }
+
+   for (int s1=0; s1<fes.Size(); ++s1)
+   {
+      delete pgfs[s1];
+      delete vdofs_all[s1];
+      delete vdofs2[s1];
+      delete vdofs[s1];
+      delete el_x[s1];
+      for (int s2=0; s2<fes.Size(); ++s2)
+      {
+         delete elmats(s1,s2);
+      }
+   }
+}
+
 BlockOperator & ParBlockNonlinearForm::GetGradient(const Vector &x) const
 {
    if (pBlockGrad == NULL)
@@ -347,49 +605,155 @@ BlockOperator & ParBlockNonlinearForm::GetGradient(const Vector &x) const
       }
    }
 
-   GetLocalGradient(x); // gradients are stored in 'Grads'
+   // xs_true is not modified, so const_cast is okay
+   xs_true.Update(const_cast<Vector &>(x), block_trueOffsets);
+   xs.Update(block_offsets);
+
+   for (int s=0; s<fes.Size(); ++s)
+   {
+      fes[s]->GetProlongationMatrix()->Mult(
+         xs_true.GetBlock(s), xs.GetBlock(s));
+   }
 
    if (fnfi.Size() > 0)
    {
-      MFEM_ABORT("TODO: assemble contributions from shared face terms");
-   }
+      const int skip_zeros = 0;
 
-   for (int s1=0; s1<fes.Size(); ++s1)
-   {
-      for (int s2=0; s2<fes.Size(); ++s2)
+      for (int s=0; s<fes.Size(); ++s)
       {
-         OperatorHandle dA(phBlockGrad(s1,s2)->Type()),
-                        Ph(phBlockGrad(s1,s2)->Type()),
-                        Rh(phBlockGrad(s1,s2)->Type());
+         const_cast<ParFiniteElementSpace*>(pfes[s])->ExchangeFaceNbrData();
+      }
 
-         if (s1 == s2)
+      for (int s1=0; s1<fes.Size(); ++s1)
+      {
+         for (int s2=0; s2<fes.Size(); ++s2)
          {
-            dA.MakeSquareBlockDiag(pfes[s1]->GetComm(), pfes[s1]->GlobalVSize(),
-                                   pfes[s1]->GetDofOffsets(), Grads(s1,s1));
-            Ph.ConvertFrom(pfes[s1]->Dof_TrueDof_Matrix());
-            phBlockGrad(s1,s1)->MakePtAP(dA, Ph);
-
-            OperatorHandle Ae;
-            Ae.EliminateRowsCols(*phBlockGrad(s1,s1), *ess_tdofs[s1]);
+            if (Grads(s1,s2) == NULL)
+            {
+               int nbr_size = pfes[s2]->GetFaceNbrVSize();
+               Grads(s1,s2) = new SparseMatrix(pfes[s1]->GetVSize(),
+                                               pfes[s2]->GetVSize() + nbr_size);
+            }
          }
-         else
+      }
+
+      // (re)assemble Grad without b.c. into 'Grads'
+      BlockNonlinearForm::ComputeGradientBlocked(xs, false);
+
+      GradientSharedFaces(xs, skip_zeros);
+
+      // finalize the gradients
+      for (int s1=0; s1<fes.Size(); ++s1)
+         for (int s2=0; s2<fes.Size(); ++s2)
          {
-            dA.MakeRectangularBlockDiag(pfes[s1]->GetComm(),
-                                        pfes[s1]->GlobalVSize(),
-                                        pfes[s2]->GlobalVSize(),
-                                        pfes[s1]->GetDofOffsets(),
-                                        pfes[s2]->GetDofOffsets(),
-                                        Grads(s1,s2));
-            Rh.ConvertFrom(pfes[s1]->Dof_TrueDof_Matrix());
-            Ph.ConvertFrom(pfes[s2]->Dof_TrueDof_Matrix());
-
-            phBlockGrad(s1,s2)->MakeRAP(Rh, dA, Ph);
-
-            phBlockGrad(s1,s2)->EliminateRows(*ess_tdofs[s1]);
-            phBlockGrad(s1,s2)->EliminateCols(*ess_tdofs[s2]);
+            Grads(s1,s2)->Finalize(skip_zeros);
          }
 
-         pBlockGrad->SetBlock(s1, s2, phBlockGrad(s1,s2)->Ptr());
+      for (int s1=0; s1<fes.Size(); ++s1)
+      {
+         for (int s2=0; s2<fes.Size(); ++s2)
+         {
+            OperatorHandle hdA;
+            OperatorHandle dA(phBlockGrad(s1,s2)->Type()),
+                           Ph(phBlockGrad(s1,s2)->Type()),
+                           Rh(phBlockGrad(s1,s2)->Type());
+
+            // handle the case when 'a' contains off-diagonal
+            int lvsize = pfes[s2]->GetVSize();
+            const HYPRE_BigInt *face_nbr_glob_ldof =
+               const_cast<ParFiniteElementSpace*>(pfes[s2])->GetFaceNbrGlobalDofMap();
+            HYPRE_BigInt ldof_offset = pfes[s2]->GetMyDofOffset();
+
+            Array<HYPRE_BigInt> glob_J(Grads(s1,s2)->NumNonZeroElems());
+            int *J = Grads(s1,s2)->GetJ();
+            for (int i = 0; i < glob_J.Size(); i++)
+            {
+               if (J[i] < lvsize)
+               {
+                  glob_J[i] = J[i] + ldof_offset;
+               }
+               else
+               {
+                  glob_J[i] = face_nbr_glob_ldof[J[i] - lvsize];
+               }
+            }
+
+            // TODO - construct dA directly in the A format
+            hdA.Reset(
+               new HypreParMatrix(pfes[s2]->GetComm(), pfes[s1]->GetVSize(),
+                                  pfes[s1]->GlobalVSize(), pfes[s2]->GlobalVSize(),
+                                  Grads(s1,s2)->GetI(), glob_J, Grads(s1,s2)->GetData(),
+                                  pfes[s1]->GetDofOffsets(), pfes[s2]->GetDofOffsets()));
+            // - hdA owns the new HypreParMatrix
+            // - the above constructor copies all input arrays
+            glob_J.DeleteAll();
+            dA.ConvertFrom(hdA);
+
+            if (s1 == s2)
+            {
+               Ph.ConvertFrom(pfes[s1]->Dof_TrueDof_Matrix());
+               phBlockGrad(s1,s1)->MakePtAP(dA, Ph);
+
+               OperatorHandle Ae;
+               Ae.EliminateRowsCols(*phBlockGrad(s1,s1), *ess_tdofs[s1]);
+            }
+            else
+            {
+               Rh.ConvertFrom(pfes[s1]->Dof_TrueDof_Matrix());
+               Ph.ConvertFrom(pfes[s2]->Dof_TrueDof_Matrix());
+
+               phBlockGrad(s1,s2)->MakeRAP(Rh, dA, Ph);
+
+               phBlockGrad(s1,s2)->EliminateRows(*ess_tdofs[s1]);
+               phBlockGrad(s1,s2)->EliminateCols(*ess_tdofs[s2]);
+            }
+
+            pBlockGrad->SetBlock(s1, s2, phBlockGrad(s1,s2)->Ptr());
+         }
+      }
+   }
+   else
+   {
+      // (re)assemble Grad without b.c. into 'Grads'
+      BlockNonlinearForm::ComputeGradientBlocked(xs);
+
+      for (int s1=0; s1<fes.Size(); ++s1)
+      {
+         for (int s2=0; s2<fes.Size(); ++s2)
+         {
+            OperatorHandle dA(phBlockGrad(s1,s2)->Type()),
+                           Ph(phBlockGrad(s1,s2)->Type()),
+                           Rh(phBlockGrad(s1,s2)->Type());
+
+            if (s1 == s2)
+            {
+               dA.MakeSquareBlockDiag(pfes[s1]->GetComm(), pfes[s1]->GlobalVSize(),
+                                      pfes[s1]->GetDofOffsets(), Grads(s1,s1));
+               Ph.ConvertFrom(pfes[s1]->Dof_TrueDof_Matrix());
+               phBlockGrad(s1,s1)->MakePtAP(dA, Ph);
+
+               OperatorHandle Ae;
+               Ae.EliminateRowsCols(*phBlockGrad(s1,s1), *ess_tdofs[s1]);
+            }
+            else
+            {
+               dA.MakeRectangularBlockDiag(pfes[s1]->GetComm(),
+                                           pfes[s1]->GlobalVSize(),
+                                           pfes[s2]->GlobalVSize(),
+                                           pfes[s1]->GetDofOffsets(),
+                                           pfes[s2]->GetDofOffsets(),
+                                           Grads(s1,s2));
+               Rh.ConvertFrom(pfes[s1]->Dof_TrueDof_Matrix());
+               Ph.ConvertFrom(pfes[s2]->Dof_TrueDof_Matrix());
+
+               phBlockGrad(s1,s2)->MakeRAP(Rh, dA, Ph);
+
+               phBlockGrad(s1,s2)->EliminateRows(*ess_tdofs[s1]);
+               phBlockGrad(s1,s2)->EliminateCols(*ess_tdofs[s2]);
+            }
+
+            pBlockGrad->SetBlock(s1, s2, phBlockGrad(s1,s2)->Ptr());
+         }
       }
    }
 

fem/pnonlinearform.hpp

@@ -29,6 +29,8 @@ protected:
    mutable ParGridFunction X, Y;
    mutable OperatorHandle pGrad;
 
+   void GradientSharedFaces(const Vector &x, int skip_zeros = 1) const;
+
 public:
    ParNonlinearForm(ParFiniteElementSpace *pf);
 
@@ -81,6 +83,8 @@ protected:
    mutable Array2D<OperatorHandle *> phBlockGrad;
    mutable BlockOperator *pBlockGrad;
 
+   void GradientSharedFaces(const BlockVector &xs, int skip_zeros) const;
+
 public:
    /// Computes the energy of the system
    real_t GetEnergy(const Vector &x) const override;

general/array.cpp

@@ -1,214 +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.
-
-
-// Abstract array data type
-
-#include "array.hpp"
-#include "../general/forall.hpp"
-#include <fstream>
-#include <type_traits>
-
-namespace mfem
-{
-
-template <class T>
-void Array<T>::Print(std::ostream &os, int width) const
-{
-   for (int i = 0; i < size; i++)
-   {
-      os << data[i];
-      if ( !((i+1) % width) || i+1 == size )
-      {
-         os << '\n';
-      }
-      else
-      {
-         os << " ";
-      }
-   }
-}
-
-template <class T>
-void Array<T>::Save(std::ostream &os, int fmt) const
-{
-   if (fmt == 0)
-   {
-      os << size << '\n';
-   }
-   for (int i = 0; i < size; i++)
-   {
-      os << operator[](i) << '\n';
-   }
-}
-
-template <class T>
-void Array<T>::Load(std::istream &in, int fmt)
-{
-   if (fmt == 0)
-   {
-      int new_size;
-      in >> new_size;
-      SetSize(new_size);
-   }
-   for (int i = 0; i < size; i++)
-   {
-      in >> operator[](i);
-   }
-}
-
-template <class T>
-T Array<T>::Max() const
-{
-   MFEM_ASSERT(size > 0, "Array is empty with size " << size);
-
-   T max = operator[](0);
-   for (int i = 1; i < size; i++)
-   {
-      if (max < operator[](i))
-      {
-         max = operator[](i);
-      }
-   }
-
-   return max;
-}
-
-template <class T>
-T Array<T>::Min() const
-{
-   MFEM_ASSERT(size > 0, "Array is empty with size " << size);
-
-   T min = operator[](0);
-   for (int i = 1; i < size; i++)
-   {
-      if (operator[](i) < min)
-      {
-         min = operator[](i);
-      }
-   }
-
-   return min;
-}
-
-// Partial Sum
-template <class T>
-void Array<T>::PartialSum()
-{
-   T sum = static_cast<T>(0);
-   for (int i = 0; i < size; i++)
-   {
-      sum+=operator[](i);
-      operator[](i) = sum;
-   }
-}
-
-template <class T>
-void Array<T>::Abs()
-{
-   static_assert(std::is_arithmetic<T>::value, "Use with arithmetic types!");
-   const bool useDevice = UseDevice();
-   const int N = size;
-   auto y = ReadWrite(useDevice);
-   mfem::forall_switch(useDevice, N, [=] MFEM_HOST_DEVICE (int i)
-   {
-      y[i] = std::abs(y[i]);
-   });
-}
-
-// Sum
-template <class T>
-T Array<T>::Sum() const
-{
-   T sum = static_cast<T>(0);
-   for (int i = 0; i < size; i++)
-   {
-      sum+=operator[](i);
-   }
-
-   return sum;
-}
-
-template <class T>
-int Array<T>::IsSorted() const
-{
-   T val_prev = operator[](0), val;
-   for (int i = 1; i < size; i++)
-   {
-      val=operator[](i);
-      if (val < val_prev)
-      {
-         return 0;
-      }
-      val_prev = val;
-   }
-
-   return 1;
-}
-
-template <class T>
-bool Array<T>::IsConstant() const
-{
-   if (size < 2) { return true; }
-   const T v0 = data[0];
-   for (int i = 1; i < size; i++)
-   {
-      if (data[i] != v0)
-      {
-         return false;
-      }
-   }
-
-   return true;
-}
-
-template <class T>
-void Array2D<T>::Load(const char *filename, int fmt)
-{
-   std::ifstream in;
-   in.open(filename, std::ifstream::in);
-   MFEM_VERIFY(in.is_open(), "File " << filename << " does not exist.");
-   Load(in, fmt);
-   in.close();
-}
-
-template <class T>
-void Array2D<T>::Print(std::ostream &os, int width_)
-{
-   int height = this->NumRows();
-   int width  = this->NumCols();
-
-   for (int i = 0; i < height; i++)
-   {
-      os << "[row " << i << "]\n";
-      for (int j = 0; j < width; j++)
-      {
-         os << (*this)(i,j);
-         if ( (j+1) == width_ || (j+1) % width_ == 0 )
-         {
-            os << '\n';
-         }
-         else
-         {
-            os << ' ';
-         }
-      }
-   }
-}
-
-template class Array<char>;
-template class Array<int>;
-template class Array<long long>;
-template class Array<real_t>;
-template class Array2D<int>;
-template class Array2D<real_t>;
-
-} // namespace mfem

general/array.hpp

@@ -16,9 +16,13 @@
 #include "mem_manager.hpp"
 #include "device.hpp"
 #include "error.hpp"
+#include "forall.hpp"
 #include "globals.hpp"
+#include "reducers.hpp"
+#include "scan.hpp"
 
 #include <iostream>
+#include <fstream>
 #include <cstdlib>
 #include <cstring>
 #include <algorithm>
@@ -135,6 +139,8 @@ public:
    /// Return the device flag of the Memory object used by the Array
    bool UseDevice() const { return data.UseDevice(); }
 
+   void UseDevice(bool use_dev) { data.UseDevice(use_dev); }
+
    /// Return true if the data will be deleted by the Array
    inline bool OwnsData() const { return data.OwnsHostPtr(); }
 
@@ -275,11 +281,11 @@ public:
 
    /** @brief Find the maximal element in the array, using the comparison
        operator `<` for class T. */
-   T Max() const;
+   inline T Max() const;
 
    /** @brief Find the minimal element in the array, using the comparison
        operator `<` for class T. */
-   T Min() const;
+   inline T Min() const;
 
    /// Sorts the array in ascending order. This requires operator< to be defined for T.
    void Sort() { std::sort((T*)data, data + size); }
@@ -297,22 +303,22 @@ public:
    }
 
    /// Return 1 if the array is sorted from lowest to highest.  Otherwise return 0.
-   int IsSorted() const;
+   inline int IsSorted() const;
 
    /// Does the Array have Size zero.
    bool IsEmpty() const { return Size() == 0; }
 
    /// Return true if all entries of the array are the same.
-   bool IsConstant() const;
+   inline bool IsConstant() const;
 
    /// Fill the entries of the array with the cumulative sum of the entries.
-   void PartialSum();
+   inline void PartialSum();
 
    /// Replace each entry of the array with its absolute value.
-   void Abs();
+   inline void Abs();
 
    /// Return the sum of all the array entries using the '+'' operator for class 'T'.
-   T Sum() const;
+   inline T Sum() const;
 
    /// Set all entries of the array to the provided constant.
    inline void operator=(const T &a);
@@ -797,8 +803,14 @@ template <typename T> template <typename CT>
 inline Array<T> &Array<T>::operator=(const Array<CT> &src)
 {
    SetSize(src.Size());
-   for (int i = 0; i < size; i++) { (*this)[i] = T(src[i]); }
-   return *this;
+
+   const bool use_dev = UseDevice() || src.UseDevice();
+   const auto x = src.Read(use_dev);
+   auto y = Write(use_dev);
+   mfem::forall_switch(use_dev, size, [=] MFEM_HOST_DEVICE (int i)
+   {
+      y[i] = x[i];
+   });
 }
 
 template <class T>
@@ -1014,19 +1026,24 @@ template <class T>
 inline void Array<T>::GetSubArray(int offset, int sa_size, Array<T> &sa) const
 {
    sa.SetSize(sa_size);
-   for (int i = 0; i < sa_size; i++)
+   const bool use_dev = UseDevice() || sa.UseDevice();
+   const auto x = Read(use_dev);
+   auto y = sa.Write(use_dev);
+   mfem::forall_switch(use_dev, sa_size, [=] MFEM_HOST_DEVICE (int i)
    {
-      sa[i] = (*this)[offset+i];
-   }
+      y[i] = x[offset + i];
+   });
 }
 
 template <class T>
 inline void Array<T>::operator=(const T &a)
 {
-   for (int i = 0; i < size; i++)
+   const bool use_dev = UseDevice();
+   auto x = Write(use_dev);
+   mfem::forall_switch(use_dev, size, [=] MFEM_HOST_DEVICE (int i)
    {
-      data[i] = a;
-   }
+      x[i] = a;
+   });
 }
 
 template <class T>
@@ -1035,6 +1052,153 @@ inline void Array<T>::Assign(const T *p)
    data.CopyFromHost(p, Size());
 }
 
+template <class T>
+inline void Array<T>::Print(std::ostream &os, int width) const
+{
+   for (int i = 0; i < size; i++)
+   {
+      os << data[i];
+      if ( !((i+1) % width) || i+1 == size )
+      {
+         os << '\n';
+      }
+      else
+      {
+         os << " ";
+      }
+   }
+}
+
+template <class T>
+inline void Array<T>::Save(std::ostream &os, int fmt) const
+{
+   if (fmt == 0)
+   {
+      os << size << '\n';
+   }
+   for (int i = 0; i < size; i++)
+   {
+      os << operator[](i) << '\n';
+   }
+}
+
+template <class T>
+void Array<T>::Load(std::istream &in, int fmt)
+{
+   if (fmt == 0)
+   {
+      int new_size;
+      in >> new_size;
+      SetSize(new_size);
+   }
+   for (int i = 0; i < size; i++)
+   {
+      in >> operator[](i);
+   }
+}
+
+template <class T>
+inline T Array<T>::Max() const
+{
+   MFEM_ASSERT(size > 0, "Array is empty with size " << size);
+
+   T max = operator[](0);
+   for (int i = 1; i < size; i++)
+   {
+      if (max < operator[](i))
+      {
+         max = operator[](i);
+      }
+   }
+
+   return max;
+}
+
+template <class T>
+inline T Array<T>::Min() const
+{
+   MFEM_ASSERT(size > 0, "Array is empty with size " << size);
+
+   T min = operator[](0);
+   for (int i = 1; i < size; i++)
+   {
+      if (operator[](i) < min)
+      {
+         min = operator[](i);
+      }
+   }
+
+   return min;
+}
+
+// Partial Sum
+template <class T>
+inline void Array<T>::PartialSum()
+{
+   auto data_ptr = ReadWrite(UseDevice());
+   InclusiveScan(UseDevice(), data_ptr, data_ptr, size);
+}
+
+template <class T>
+inline void Array<T>::Abs()
+{
+   static_assert(std::is_arithmetic<T>::value, "Use with arithmetic types!");
+   const bool useDevice = UseDevice();
+   const int N = size;
+   auto y = ReadWrite(useDevice);
+   mfem::forall_switch(useDevice, N, [=] MFEM_HOST_DEVICE (int i)
+   {
+      y[i] = std::abs(y[i]);
+   });
+}
+
+// Sum
+template <class T>
+inline T Array<T>::Sum() const
+{
+   T sum = static_cast<T>(0);
+   if (size > 0)
+   {
+      const auto m_data = Read(UseDevice());
+      reduce(size, sum, [=] MFEM_HOST_DEVICE(int i, T &r) { r += m_data[i]; },
+      /*  */ SumReducer<T> {}, UseDevice());
+   }
+   return sum;
+}
+
+template <class T>
+inline int Array<T>::IsSorted() const
+{
+   T val_prev = operator[](0), val;
+   for (int i = 1; i < size; i++)
+   {
+      val=operator[](i);
+      if (val < val_prev)
+      {
+         return 0;
+      }
+      val_prev = val;
+   }
+
+   return 1;
+}
+
+template <class T>
+inline bool Array<T>::IsConstant() const
+{
+   if (size < 2) { return true; }
+   const T v0 = data[0];
+   for (int i = 1; i < size; i++)
+   {
+      if (data[i] != v0)
+      {
+         return false;
+      }
+   }
+
+   return true;
+}
+
 
 template <class T>
 inline const T &Array2D<T>::operator()(int i, int j) const
@@ -1074,6 +1238,40 @@ inline T *Array2D<T>::operator[](int i)
    return &array1d[i*N];
 }
 
+template <class T>
+void Array2D<T>::Load(const char *filename, int fmt)
+{
+   std::ifstream in;
+   in.open(filename, std::ifstream::in);
+   MFEM_VERIFY(in.is_open(), "File " << filename << " does not exist.");
+   Load(in, fmt);
+   in.close();
+}
+
+template <class T>
+void Array2D<T>::Print(std::ostream &os, int width_)
+{
+   int height = this->NumRows();
+   int width  = this->NumCols();
+
+   for (int i = 0; i < height; i++)
+   {
+      os << "[row " << i << "]\n";
+      for (int j = 0; j < width; j++)
+      {
+         os << (*this)(i,j);
+         if ( (j+1) == width_ || (j+1) % width_ == 0 )
+         {
+            os << '\n';
+         }
+         else
+         {
+            os << ' ';
+         }
+      }
+   }
+}
+
 
 template <class T>
 inline void Swap(Array2D<T> &a, Array2D<T> &b)

general/reducers.hpp

@@ -12,7 +12,6 @@
 #ifndef MFEM_REDUCERS_HPP
 #define MFEM_REDUCERS_HPP
 
-#include "array.hpp"
 #include "forall.hpp"
 
 #include <cmath>
@@ -514,6 +513,33 @@ template<class B, class R> struct reduction_kernel
       }
    }
 };
+
+template <class T>
+class ReductionWorkspace
+{
+   Memory<T> workspace;
+
+   static ReductionWorkspace &Instance()
+   {
+      static ReductionWorkspace instance;
+      return instance;
+   }
+
+   ~ReductionWorkspace() { workspace.Delete(); }
+
+public:
+   static T *Get(int num_blocks)
+   {
+      ReductionWorkspace &instance = Instance();
+      if (instance.workspace.Capacity() < num_blocks)
+      {
+         instance.workspace.Delete();
+         instance.workspace.New(num_blocks, MemoryType::HOST_PINNED);
+      }
+      return instance.workspace;
+   }
+};
+
 }
 
 /**
@@ -529,8 +555,7 @@ template<class B, class R> struct reduction_kernel
  @tparam T value_type to operate on
  */
 template <class T, class B, class R>
-void reduce(int N, T &res, B &&body, const R &reducer, bool use_dev,
-            Array<T> &workspace)
+void reduce(int N, T &res, B &&body, const R &reducer, bool use_dev)
 {
    if (N == 0)
    {
@@ -567,13 +592,7 @@ void reduce(int N, T &res, B &&body, const R &reducer, bool use_dev,
 
       red_type red{nullptr, std::forward<B>(body), reducer, N, items_per_thread};
       // allocate res to fit block_size entries
-      auto mt = workspace.GetMemory().GetMemoryType();
-      if (mt != MemoryType::HOST_PINNED && mt != MemoryType::MANAGED)
-      {
-         mt = MemoryType::HOST_PINNED;
-      }
-      workspace.SetSize(nblocks, mt);
-      auto work = workspace.HostWrite();
+      auto work = internal::ReductionWorkspace<T>::Get(nblocks);
       red.work = work;
       forall_2D(nblocks, block_size, 1, std::move(red));
       // wait for results

general/scan.hpp

@@ -28,8 +28,37 @@
 
 namespace mfem
 {
-/// Equivalent to InclusiveScan(use_dev, d_in, d_out, num_items, workspace,
-/// std::plus<>{})
+
+namespace internal
+{
+class ScanWorkspace
+{
+   Memory<std::byte> workspace;
+   static ScanWorkspace &Instance()
+   {
+      static ScanWorkspace instance;
+      return instance;
+   }
+   ~ScanWorkspace() { workspace.Delete(); }
+public:
+   static std::byte *Get(int num_bytes)
+   {
+      ScanWorkspace &instance = Instance();
+      if (Size() < num_bytes)
+      {
+         instance.workspace.Delete();
+         instance.workspace.New(num_bytes);
+      }
+      return instance.workspace.Write(MemoryClass::DEVICE, Size());
+   }
+   static int Size()
+   {
+      return Instance().workspace.Capacity();
+   }
+};
+}
+
+/// Equivalent to InclusiveScan(use_dev, d_in, d_out, num_items, std::plus<>{})
 template <class InputIt, class OutputIt>
 void InclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items)
 {
@@ -37,12 +66,12 @@ void InclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items)
 #if defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP)
    if (use_dev && mfem::Device::Allows(Backend::CUDA_MASK | Backend::HIP_MASK))
    {
-      static Array<std::byte> workspace;
-      size_t bytes = workspace.Size();
-      if (bytes)
+      using internal::ScanWorkspace;
+      size_t bytes = ScanWorkspace::Size();
+      if (bytes > 0)
       {
          auto err = MFEM_CUB_NAMESPACE::DeviceScan::InclusiveSum(
-                       workspace.Write(), bytes, d_in, d_out, num_items);
+                       ScanWorkspace::Get(bytes), bytes, d_in, d_out, num_items);
 #if defined(MFEM_USE_CUDA)
          if (err == cudaSuccess)
          {
@@ -57,11 +86,12 @@ void InclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items)
       }
       // try allocating a larger buffer
       bytes = 0;
+      // get size of buffer
       MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::InclusiveSum(
                         nullptr, bytes, d_in, d_out, num_items));
-      workspace.SetSize(bytes);
+      // resize buffer (in ScanWorkspace::Get) and try again
       MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::InclusiveSum(
-                        workspace.Write(), bytes, d_in, d_out, num_items));
+                        ScanWorkspace::Get(bytes), bytes, d_in, d_out, num_items));
       return;
    }
 #endif
@@ -101,12 +131,13 @@ void InclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
 #if defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP)
    if (use_dev && mfem::Device::Allows(Backend::CUDA_MASK | Backend::HIP_MASK))
    {
-      static Array<std::byte> workspace;
-      size_t bytes = workspace.Size();
-      if (bytes)
+      using internal::ScanWorkspace;
+      size_t bytes = ScanWorkspace::Size();
+      if (bytes > 0)
       {
          auto err = MFEM_CUB_NAMESPACE::DeviceScan::InclusiveScan(
-                       workspace.Write(), bytes, d_in, d_out, scan_op, num_items);
+                       ScanWorkspace::Get(bytes), bytes, d_in, d_out, scan_op,
+                       num_items);
 #if defined(MFEM_USE_CUDA)
          if (err == cudaSuccess)
          {
@@ -123,9 +154,9 @@ void InclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
       bytes = 0;
       MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::InclusiveScan(
                         nullptr, bytes, d_in, d_out, scan_op, num_items));
-      workspace.SetSize(bytes);
       MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::InclusiveScan(
-                        workspace.Write(), bytes, d_in, d_out, scan_op, num_items));
+                        ScanWorkspace::Get(bytes), bytes, d_in, d_out, scan_op,
+                        num_items));
       return;
    }
 #endif
@@ -164,13 +195,13 @@ void ExclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
 #if defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP)
    if (use_dev && mfem::Device::Allows(Backend::CUDA_MASK | Backend::HIP_MASK))
    {
-      static Array<std::byte> workspace;
-      size_t bytes = workspace.Size();
+      using internal::ScanWorkspace;
+      size_t bytes = ScanWorkspace::Size();
       if (bytes)
       {
          auto err = MFEM_CUB_NAMESPACE::DeviceScan::ExclusiveScan(
-                       workspace.Write(), bytes, d_in, d_out, scan_op, init_value,
-                       num_items);
+                       ScanWorkspace::Get(bytes), bytes, d_in, d_out, scan_op,
+                       init_value, num_items);
 #if defined(MFEM_USE_CUDA)
          if (err == cudaSuccess)
          {
@@ -187,10 +218,9 @@ void ExclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
       bytes = 0;
       MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::ExclusiveScan(
                         nullptr, bytes, d_in, d_out, scan_op, init_value, num_items));
-      workspace.SetSize(bytes);
       MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::ExclusiveScan(
-                        workspace.Write(), bytes, d_in, d_out, scan_op, init_value,
-                        num_items));
+                        ScanWorkspace::Get(bytes), bytes, d_in, d_out, scan_op,
+                        init_value, num_items));
       return;
    }
 #endif
@@ -213,7 +243,7 @@ void ExclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
 }
 
 /// Equivalent to ExclusiveScan(use_dev, d_in, d_out, num_items, init_value,
-/// workspace, std::plus<>{})
+/// std::plus<>{})
 template <class InputIt, class OutputIt, class T>
 void ExclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
                    T init_value)

linalg/batched/magma.cpp

@@ -167,7 +167,7 @@ void MagmaBatchedLinAlg::Invert(DenseTensor &A) const
    magma_int_t status;
 
    status = MFEM_MAGMA_PREFIX(getrf_batched)(
-               n, n, d_A_ptrs, n, d_P_ptrs, info_array.Write(), n_mat,
+               n, n, d_LU_ptrs, n, d_P_ptrs, info_array.Write(), n_mat,
                Magma::Queue());
    MFEM_VERIFY(status == MAGMA_SUCCESS, "");
 

linalg/hypre.cpp

@@ -561,7 +561,8 @@ void CopyMemory(Memory<T> &src, Memory<T> &dst, MemoryClass dst_mc,
     this function. In particular, @a dst should be empty or deleted before
     calling this function. */
 template <typename SrcT, typename DstT>
-void CopyConvertMemory(Memory<SrcT> &src, MemoryClass dst_mc, Memory<DstT> &dst)
+void CopyConvertMemory(const Memory<SrcT> &src, MemoryClass dst_mc,
+                       Memory<DstT> &dst)
 {
    auto capacity = src.Capacity();
    dst.New(capacity, GetMemoryType(dst_mc));
@@ -842,8 +843,8 @@ static int GetPartitioningArraySize(MPI_Comm comm)
 ///
 /// Both @a row and @a col are partitioning arrays, whose length is returned by
 /// GetPartitioningArraySize(), see @ref hypre_partitioning_descr.
-static bool RowAndColStartsAreEqual(MPI_Comm comm, HYPRE_BigInt *rows,
-                                    HYPRE_BigInt *cols)
+static bool RowAndColStartsAreEqual(MPI_Comm comm, const HYPRE_BigInt *rows,
+                                    const HYPRE_BigInt *cols)
 {
    const int part_size = GetPartitioningArraySize(comm);
    bool are_equal = true;
@@ -1131,7 +1132,7 @@ HypreParMatrix::HypreParMatrix(
 HypreParMatrix::HypreParMatrix(MPI_Comm comm,
                                HYPRE_BigInt *row_starts,
                                HYPRE_BigInt *col_starts,
-                               SparseMatrix *sm_a)
+                               const SparseMatrix *sm_a)
 {
    MFEM_ASSERT(sm_a != NULL, "invalid input");
    MFEM_VERIFY(!HYPRE_AssumedPartitionCheck(),
@@ -1145,7 +1146,7 @@ HypreParMatrix::HypreParMatrix(MPI_Comm comm,
 
    hypre_CSRMatrixSetDataOwner(csr_a,0);
    MemoryIJData mem_a;
-   CopyCSR(sm_a, mem_a, csr_a, false);
+   CopyCSR(const_cast<SparseMatrix*>(sm_a), mem_a, csr_a, false);
    hypre_CSRMatrixSetRownnz(csr_a);
 
    // NOTE: this call creates a matrix on host even when device support is
@@ -1307,10 +1308,11 @@ HypreParMatrix::HypreParMatrix(MPI_Comm comm, int id, int np,
 HypreParMatrix::HypreParMatrix(MPI_Comm comm, int nrows,
                                HYPRE_BigInt glob_nrows,
                                HYPRE_BigInt glob_ncols,
-                               int *I, HYPRE_BigInt *J,
-                               real_t *data,
-                               HYPRE_BigInt *rows,
-                               HYPRE_BigInt *cols)
+                               const int *I,
+                               const HYPRE_BigInt *J,
+                               const real_t *data,
+                               const HYPRE_BigInt *rows,
+                               const HYPRE_BigInt *cols)
 {
    Init();
 
@@ -2327,8 +2329,8 @@ void HypreParMatrix::Threshold(real_t threshold)
    /* TODO: GenerateDiagAndOffd() uses an int array of size equal to the number
       of columns in csr_A_wo_z which is the global number of columns in A. This
       does not scale well. */
-   ierr += GenerateDiagAndOffd(csr_A_wo_z,parcsr_A_ptr,
-                               col_start,col_end);
+   ierr += hypre_GenerateDiagAndOffd(csr_A_wo_z,parcsr_A_ptr,
+                                     col_start,col_end);
 
    ierr += hypre_CSRMatrixDestroy(csr_A_wo_z);
 

linalg/hypre.hpp

@@ -25,11 +25,18 @@
 #define HYPRE_TIMING
 
 // hypre header files
+#if MFEM_HYPRE_VERSION < 30000
 #include <seq_mv.h>
 #include <temp_multivector.h>
+#else
+#include <_hypre_seq_mv.h>
+#include <_hypre_lobpcg_temp_multivector.h>
+#endif
 #include <_hypre_parcsr_mv.h>
 #include <_hypre_parcsr_ls.h>
 
+#include <HYPRE_parcsr_ls.h>
+
 #ifdef HYPRE_COMPLEX
 #error "MFEM does not work with HYPRE's complex numbers support"
 #endif
@@ -53,6 +60,10 @@
 #error "MFEM_USE_HIP=YES is required when HYPRE is built with HIP!"
 #endif
 
+#if MFEM_HYPRE_VERSION > 21500
+#define HYPRE_AssumedPartitionCheck() 1
+#endif
+
 namespace mfem
 {
 
@@ -554,7 +565,7 @@ public:
        partitioning arrays @a row_starts and @a col_starts. */
    HypreParMatrix(MPI_Comm comm, HYPRE_BigInt *row_starts,
                   HYPRE_BigInt *col_starts,
-                  SparseMatrix *a); // constructor with 4 arguments, v2
+                  const SparseMatrix *a); // constructor with 4 arguments, v2
 
    /// Creates boolean block-diagonal rectangular parallel matrix.
    /** The new HypreParMatrix does not take ownership of any of the input
@@ -583,9 +594,9 @@ public:
        arrays (so they can be deleted). See @ref hypre_partitioning_descr "here"
        for a description of the partitioning arrays @a rows and @a cols. */
    HypreParMatrix(MPI_Comm comm, int nrows, HYPRE_BigInt glob_nrows,
-                  HYPRE_BigInt glob_ncols, int *I, HYPRE_BigInt *J,
-                  real_t *data, HYPRE_BigInt *rows,
-                  HYPRE_BigInt *cols); // constructor with 9 arguments
+                  HYPRE_BigInt glob_ncols, const int *I, const HYPRE_BigInt *J,
+                  const real_t *data, const HYPRE_BigInt *rows,
+                  const HYPRE_BigInt *cols); // constructor with 9 arguments
 
    /** @brief Copy constructor for a ParCSR matrix which creates a deep copy of
        structure and data from @a P. */

linalg/hypre_parcsr.cpp

@@ -1916,9 +1916,9 @@ hypre_ParCSRMatrixAdd(hypre_ParCSRMatrix *A,
       /* FIXME: GenerateDiagAndOffd() uses an int array of size equal to the
          number of columns in csr_C_temp which is the global number of columns
          in A and B. This does not scale well. */
-      ierr += GenerateDiagAndOffd(csr_C_temp, C,
-                                  hypre_ParCSRMatrixFirstColDiag(A),
-                                  hypre_ParCSRMatrixLastColDiag(A));
+      ierr += hypre_GenerateDiagAndOffd(csr_C_temp, C,
+                                        hypre_ParCSRMatrixFirstColDiag(A),
+                                        hypre_ParCSRMatrixLastColDiag(A));
 
       /* delete CSR version of C */
       ierr += hypre_CSRMatrixDestroy(csr_C_temp);

linalg/hypre_parcsr.hpp

@@ -21,6 +21,10 @@
 // hypre header files
 #include <_hypre_parcsr_mv.h>
 
+#if MFEM_HYPRE_VERSION < 30000
+#define hypre_GenerateDiagAndOffd GenerateDiagAndOffd
+#endif
+
 // Older hypre versions do not define HYPRE_BigInt and HYPRE_MPI_BIG_INT, so we
 // define them here for backward compatibility.
 #if MFEM_HYPRE_VERSION < 21600

linalg/sparsemat.cpp

@@ -46,6 +46,12 @@
 #define MFEM_GPUSPARSE_ALG HIPSPARSE_CSRMV_ALG1
 #endif // defined(MFEM_USE_CUDA)
 
+#if defined(MFEM_USE_SINGLE)
+#define MFEM_CUDA_or_HIP_REAL_T MFEM_CUDA_or_HIP(_R_32F)
+#elif defined(MFEM_USE_DOUBLE)
+#define MFEM_CUDA_or_HIP_REAL_T MFEM_CUDA_or_HIP(_R_64F)
+#endif
+
 namespace mfem
 {
 
@@ -57,8 +63,10 @@ int SparseMatrix::SparseMatrixCount = 0;
 /// @cond Suppress_Doxygen_warnings
 MFEM_cu_or_hip(sparseHandle_t) SparseMatrix::handle = nullptr;
 /// @endcond
+#ifndef MFEM_CUDA_1897_WORKAROUND
 size_t SparseMatrix::bufferSize = 0;
 void * SparseMatrix::dBuffer = nullptr;
+#endif
 #endif // MFEM_USE_CUDA_OR_HIP
 
 void SparseMatrix::InitGPUSparse()
@@ -464,109 +472,67 @@ void SparseMatrix::SortColumnIndices()
    }
 
 #ifdef MFEM_USE_CUDA_OR_HIP
-   if ( Device::Allows( Backend::CUDA_MASK ))
+   if (Device::Allows(Backend::CUDA_MASK) || Device::Allows(Backend::HIP_MASK))
    {
-#if defined(MFEM_USE_CUDA)
-      size_t pBufferSizeInBytes = 0;
-      void *pBuffer = NULL;
-
-      const int n = Height();
-      const int m = Width();
+      const int m = Height();
+      const int n = Width();
       const int nnzA = J.Capacity();
-      real_t * d_a_sorted = ReadWriteData();
-      const int * d_ia = ReadI();
-      int * d_ja_sorted = ReadWriteJ();
-      csru2csrInfo_t sortInfoA;
+      const int *d_ia = ReadI();
+      int *d_ja = ReadWriteJ();
 
-      cusparseMatDescr_t matA_descr;
-      cusparseCreateMatDescr( &matA_descr );
-      cusparseSetMatIndexBase( matA_descr, CUSPARSE_INDEX_BASE_ZERO );
-      cusparseSetMatType( matA_descr, CUSPARSE_MATRIX_TYPE_GENERAL );
+      // Get size of temporary buffer needed to sort the column indices,
+      // allocate the temporary buffer.
+      size_t pBufferSizeInBytes;
+      MFEM_cu_or_hip(sparseXcsrsort_bufferSizeExt)(handle, m, n, nnzA, d_ia,
+                                                   d_ja, &pBufferSizeInBytes);
+      void *pBuffer = MFEM_Cu_or_Hip(MemAlloc)(&pBuffer, pBufferSizeInBytes);
 
-      cusparseCreateCsru2csrInfo( &sortInfoA );
+      // Create matrix descriptor, will have default values
+      // CUSPARSE_INDEX_BASE_ZERO and CUSPARSE_MATRIX_TYPE_GENERAL.
+      MFEM_cu_or_hip(sparseMatDescr_t) matA_descr;
+      MFEM_cu_or_hip(sparseCreateMatDescr)(&matA_descr);
 
-#ifdef MFEM_USE_SINGLE
-      cusparseScsru2csr_bufferSizeExt( handle, n, m, nnzA, d_a_sorted, d_ia,
-                                       d_ja_sorted, sortInfoA,
-                                       &pBufferSizeInBytes);
-#elif defined MFEM_USE_DOUBLE
-      cusparseDcsru2csr_bufferSizeExt( handle, n, m, nnzA, d_a_sorted, d_ia,
-                                       d_ja_sorted, sortInfoA,
-                                       &pBufferSizeInBytes);
-#else
-      MFEM_ABORT("Floating point type undefined");
-#endif
+      // Initialize permutation to identity
+      Array<int> P(nnzA);
+      int *d_P = P.Write();
+      mfem::forall(nnzA, [=] MFEM_HOST_DEVICE (int i) { d_P[i] = i; });
 
-      CuMemAlloc( &pBuffer, pBufferSizeInBytes );
+      // Sort the column indices. The array d_ja will now be sorted. The
+      // permutation required to sort the values will be returned in d_P.
+      MFEM_cu_or_hip(sparseXcsrsort)(handle, m, n, nnzA, matA_descr, d_ia, d_ja,
+                                     d_P, pBuffer);
 
-#ifdef MFEM_USE_SINGLE
-      cusparseScsru2csr( handle, n, m, nnzA, matA_descr, d_a_sorted, d_ia,
-                         d_ja_sorted, sortInfoA, pBuffer);
-#elif defined MFEM_USE_DOUBLE
-      cusparseDcsru2csr( handle, n, m, nnzA, matA_descr, d_a_sorted, d_ia,
-                         d_ja_sorted, sortInfoA, pBuffer);
-#else
-      MFEM_ABORT("Floating point type undefined");
-#endif
+      // Create a copy of the unsorted matrix values.
+      real_t *d_a = ReadWriteData();
+      void *d_a_unsorted = MFEM_Cu_or_Hip(MemAlloc)(&d_a_unsorted,
+                                                    nnzA * sizeof(real_t));
+      MFEM_Cu_or_Hip(MemcpyDtoD)(d_a_unsorted, d_a, nnzA * sizeof(real_t));
 
-      // The above call is (at least in some cases) asynchronous, so we need to
-      // wait for it to finish before we can free device temporaries.
+      // Create the (input) dense vector with the unsorted values.
+      MFEM_cu_or_hip(sparseDnVecDescr_t) d_a_dense;
+      MFEM_cu_or_hip(sparseCreateDnVec)(&d_a_dense, nnzA, d_a_unsorted,
+                                        MFEM_CUDA_or_HIP_REAL_T);
+
+      // Create the (output) sparse vector that will have the sorted values.
+      MFEM_cu_or_hip(sparseSpVecDescr_t) d_a_sparse;
+      MFEM_cu_or_hip(sparseCreateSpVec)(&d_a_sparse, nnzA, nnzA, d_P, d_a,
+                                        MFEM_CU_or_HIP(SPARSE_INDEX_32I),
+                                        MFEM_CU_or_HIP(SPARSE_INDEX_BASE_ZERO),
+                                        MFEM_CUDA_or_HIP_REAL_T);
+
+      // Sort the matrix values using the permutation vector.
+      MFEM_cu_or_hip(sparseGather)(handle, d_a_dense, d_a_sparse);
+
+      // The above calls may be asynchronous, so we need to wait for them to
+      // finish before we can free memory.
       MFEM_STREAM_SYNC;
 
-      cusparseDestroyCsru2csrInfo( sortInfoA );
-      cusparseDestroyMatDescr( matA_descr );
+      MFEM_cu_or_hip(sparseDestroyDnVec)(d_a_dense);
+      MFEM_cu_or_hip(sparseDestroySpVec)(d_a_sparse);
+      MFEM_cu_or_hip(sparseDestroyMatDescr)(matA_descr);
 
-      CuMemFree( pBuffer );
-#endif
-   }
-   else if ( Device::Allows( Backend::HIP_MASK ))
-   {
-#if defined(MFEM_USE_HIP)
-      size_t pBufferSizeInBytes = 0;
-      void *pBuffer = NULL;
-      int *P = NULL;
-
-      const int n = Height();
-      const int m = Width();
-      const int nnzA = J.Capacity();
-      real_t * d_a_sorted = ReadWriteData();
-      const int * d_ia = ReadI();
-      int * d_ja_sorted = ReadWriteJ();
-
-      hipsparseMatDescr_t descrA;
-      hipsparseCreateMatDescr( &descrA );
-      // FIXME: There is not in-place version of csr sort in hipSPARSE currently, so we make
-      //        a temporary copy of the data for gthr, sort that, and then copy the sorted values
-      //        back to the array being returned. Where there is an in-place version available,
-      //        we should use it.
-      Array< real_t > a_tmp( nnzA );
-      real_t *d_a_tmp = a_tmp.Write();
-
-      hipsparseXcsrsort_bufferSizeExt(handle, n, m, nnzA, d_ia, d_ja_sorted,
-                                      &pBufferSizeInBytes);
-
-      HipMemAlloc( &pBuffer, pBufferSizeInBytes );
-      HipMemAlloc( (void**)&P, nnzA * sizeof(int) );
-
-      hipsparseCreateIdentityPermutation(handle, nnzA, P);
-      hipsparseXcsrsort(handle, n, m, nnzA, descrA, d_ia, d_ja_sorted, P, pBuffer);
-
-#if defined(MFEM_USE_SINGLE)
-      hipsparseSgthr(handle, nnzA, d_a_sorted, d_a_tmp, P,
-                     HIPSPARSE_INDEX_BASE_ZERO);
-#elif defined(MFEM_USE_DOUBLE)
-      hipsparseDgthr(handle, nnzA, d_a_sorted, d_a_tmp, P,
-                     HIPSPARSE_INDEX_BASE_ZERO);
-#else
-      MFEM_ABORT("Unsupported floating point type!");
-#endif
-
-      A.CopyFrom( a_tmp.GetMemory(), nnzA );
-      hipsparseDestroyMatDescr( descrA );
-
-      HipMemFree( pBuffer );
-      HipMemFree( P );
-#endif
+      MFEM_Cu_or_Hip(MemFree)(d_a_unsorted);
+      MFEM_Cu_or_Hip(MemFree)(pBuffer);
    }
    else
 #endif // MFEM_USE_CUDA_OR_HIP
@@ -821,27 +787,15 @@ void SparseMatrix::AddMult(const Vector &x, Vector &y, const real_t a) const
             MFEM_CU_or_HIP(SPARSE_INDEX_32I),
             MFEM_CU_or_HIP(SPARSE_INDEX_32I),
             MFEM_CU_or_HIP(SPARSE_INDEX_BASE_ZERO),
-#ifdef MFEM_USE_SINGLE
-            MFEM_CUDA_or_HIP(_R_32F));
-#else
-            MFEM_CUDA_or_HIP(_R_64F));
-#endif
+            MFEM_CUDA_or_HIP_REAL_T);
 
          // Create handles for input/output vectors
          MFEM_cu_or_hip(sparseCreateDnVec)(&vecX_descr,
                                            x.Size(),
                                            const_cast<real_t *>(d_x),
-#ifdef MFEM_USE_SINGLE
-                                           MFEM_CUDA_or_HIP(_R_32F));
-#else
-                                           MFEM_CUDA_or_HIP(_R_64F));
-#endif
+                                           MFEM_CUDA_or_HIP_REAL_T);
          MFEM_cu_or_hip(sparseCreateDnVec)(&vecY_descr, y.Size(), d_y,
-#ifdef MFEM_USE_SINGLE
-                                           MFEM_CUDA_or_HIP(_R_32F));
-#else
-                                           MFEM_CUDA_or_HIP(_R_64F));
-#endif
+                                           MFEM_CUDA_or_HIP_REAL_T);
 #else
          cusparseCreateMatDescr(&matA_descr);
          cusparseSetMatIndexBase(matA_descr, CUSPARSE_INDEX_BASE_ZERO);
@@ -860,11 +814,7 @@ void SparseMatrix::AddMult(const Vector &x, Vector &y, const real_t a) const
          vecX_descr,
          &beta,
          vecY_descr,
-#ifdef MFEM_USE_SINGLE
-         MFEM_CUDA_or_HIP(_R_32F),
-#else
-         MFEM_CUDA_or_HIP(_R_64F),
-#endif
+         MFEM_CUDA_or_HIP_REAL_T,
          MFEM_GPUSPARSE_ALG,
          &newBufferSize);
 
@@ -891,11 +841,7 @@ void SparseMatrix::AddMult(const Vector &x, Vector &y, const real_t a) const
          vecX_descr,
          &beta,
          vecY_descr,
-#ifdef MFEM_USE_SINGLE
-         MFEM_CUDA_or_HIP(_R_32F),
-#else
-         MFEM_CUDA_or_HIP(_R_64F),
-#endif
+         MFEM_CUDA_or_HIP_REAL_T,
          MFEM_GPUSPARSE_ALG,
          dBuffer);
 #else
@@ -4372,6 +4318,14 @@ SparseMatrix::~SparseMatrix()
 #ifdef MFEM_USE_CUDA_OR_HIP
    if (Device::Allows(Backend::CUDA_MASK | Backend::HIP_MASK))
    {
+#ifdef MFEM_CUDA_1897_WORKAROUND
+      if (dBuffer)
+      {
+         MFEM_Cu_or_Hip(MemFree)(dBuffer);
+         dBuffer = nullptr;
+         bufferSize = 0;
+      }
+#endif
       if (SparseMatrixCount==1)
       {
          if (handle)
@@ -4379,12 +4333,14 @@ SparseMatrix::~SparseMatrix()
             MFEM_cu_or_hip(sparseDestroy)(handle);
             handle = nullptr;
          }
+#ifndef MFEM_CUDA_1897_WORKAROUND
          if (dBuffer)
          {
             MFEM_Cu_or_Hip(MemFree)(dBuffer);
             dBuffer = nullptr;
             bufferSize = 0;
          }
+#endif
       }
       SparseMatrixCount--;
    }

linalg/sparsemat.hpp

@@ -98,9 +98,17 @@ protected:
 #ifdef MFEM_USE_CUDA_OR_HIP
    // common for hipSPARSE and cuSPARSE
    static int SparseMatrixCount;
+   mutable bool initBuffers = false;
+
+#if defined(MFEM_USE_CUDA) && CUDA_VERSION >= 12300 && CUDA_VERSION < 12602
+   // Workaround for bug CUSPARSE-1897
+#define MFEM_CUDA_1897_WORKAROUND
+   mutable size_t bufferSize = 0;
+   mutable void *dBuffer = nullptr;
+#else
    static size_t bufferSize;
    static void *dBuffer;
-   mutable bool initBuffers = false;
+#endif
 
 #if defined(MFEM_USE_CUDA)
    cusparseStatus_t status;

linalg/vector.cpp

@@ -92,18 +92,6 @@ struct LpReducer
    }
 };
 
-static Array<real_t>& vector_workspace()
-{
-   static Array<real_t> instance;
-   return instance;
-}
-
-static Array<DevicePair<real_t, real_t>> &Lpvector_workspace()
-{
-   static Array<DevicePair<real_t, real_t>> instance;
-   return instance;
-}
-
 Vector::Vector(const Vector &v)
 {
    const int s = v.Size();
@@ -991,7 +979,7 @@ real_t Vector::Norml2() const
          }
       }
    },
-   L2Reducer{}, UseDevice(), Lpvector_workspace());
+   L2Reducer{}, UseDevice());
    // final answer
    return res.second * sqrt(res.first);
 }
@@ -1006,7 +994,7 @@ real_t Vector::Normlinf() const
    {
       r = fmax(r, fabs(m_data[i]));
    },
-   MaxReducer<real_t> {}, UseDevice(), vector_workspace());
+   MaxReducer<real_t> {}, UseDevice());
    return res;
 }
 
@@ -1020,7 +1008,7 @@ real_t Vector::Norml1() const
    {
       r += fabs(m_data[i]);
    },
-   SumReducer<real_t> {}, UseDevice(), vector_workspace());
+   SumReducer<real_t> {}, UseDevice());
    return res;
 }
 
@@ -1063,7 +1051,7 @@ real_t Vector::Normlp(real_t p) const
             }
          }
       },
-      LpReducer{p}, UseDevice(), Lpvector_workspace());
+      LpReducer{p}, UseDevice());
       // final answer
       return res.second * pow(res.first, 1.0 / p);
    } // end if p < infinity()
@@ -1096,7 +1084,7 @@ real_t Vector::operator*(const Vector &v) const
       {
          r += m_data[i] * v_data[i];
       },
-      SumReducer<real_t> {}, use_dev, vector_workspace());
+      SumReducer<real_t> {}, use_dev);
       return res;
    };
 
@@ -1167,7 +1155,7 @@ real_t Vector::Min() const
       {
          r = fmin(r, m_data[i]);
       },
-      MinReducer<real_t> {}, use_dev, vector_workspace());
+      MinReducer<real_t> {}, use_dev);
       return res;
    };
 
@@ -1213,7 +1201,7 @@ real_t Vector::Max() const
       {
          r = fmax(r, m_data[i]);
       },
-      MaxReducer<real_t> {}, use_dev, vector_workspace());
+      MaxReducer<real_t> {}, use_dev);
       return res;
    };
 
@@ -1248,7 +1236,7 @@ real_t Vector::Sum() const
    {
       r += m_data[i];
    },
-   SumReducer<real_t> {}, UseDevice(), vector_workspace());
+   SumReducer<real_t> {}, UseDevice());
    return res;
 }
 

makefile

@@ -377,7 +377,7 @@ MFEM_CONFIG_VARS = MFEM_CXX MFEM_HOST_CXX MFEM_CPPFLAGS MFEM_CXXFLAGS\
  MFEM_INC_DIR MFEM_TPLFLAGS MFEM_INCFLAGS MFEM_PICFLAG MFEM_FLAGS MFEM_LIB_DIR\
  MFEM_EXT_LIBS MFEM_LIBS MFEM_LIB_FILE MFEM_STATIC MFEM_SHARED MFEM_BUILD_TAG\
  MFEM_PREFIX MFEM_CONFIG_EXTRA MFEM_MPIEXEC MFEM_MPIEXEC_NP MFEM_MPI_NP\
- MFEM_TEST_MK
+ MFEM_TEST_MK MFEM_XLINKER
 
 # Config vars: values of the form @VAL@ are replaced by $(VAL) in config.mk
 MFEM_CPPFLAGS  ?= $(CPPFLAGS)
@@ -394,6 +394,7 @@ MFEM_BUILD_TAG ?= $(shell uname -snm)
 MFEM_PREFIX    ?= $(PREFIX)
 MFEM_INC_DIR   ?= $(if $(CONFIG_FILE_DEF),@MFEM_BUILD_DIR@,@MFEM_DIR@)
 MFEM_LIB_DIR   ?= $(if $(CONFIG_FILE_DEF),@MFEM_BUILD_DIR@,@MFEM_DIR@)
+MFEM_XLINKER   ?= $(XLINKER)
 MFEM_TEST_MK   ?= @MFEM_DIR@/config/test.mk
 # Use "\n" (interpreted by sed) to add a newline.
 MFEM_CONFIG_EXTRA ?= $(if $(CONFIG_FILE_DEF),MFEM_BUILD_DIR ?= @MFEM_DIR@,)

mesh/mesh_readers.cpp

@@ -773,7 +773,7 @@ struct BufferReader : BufferReaderBase
          int header_entry_size = HeaderEntrySize();
          int nblocks = ReadHeaderEntry(header_buf);
          header_buf += header_entry_size;
-         std::vector<int> header(nblocks + 2);
+         std::vector<size_t> header(nblocks + 2);
          for (int i=0; i<nblocks+2; ++i)
          {
             header[i] = ReadHeaderEntry(header_buf);
@@ -792,7 +792,7 @@ struct BufferReader : BufferReaderBase
             dest_ptr += dest_len;
             source_ptr += source_len;
          }
-         MFEM_VERIFY(int(sizeof(F)*n) == (dest_ptr - dest_start),
+         MFEM_VERIFY(size_t(sizeof(F)*n) == (dest_ptr - dest_start),
                      "AppendedData: wrong data size");
          buf = uncompressed_data.data();
 #else

miniapps/autodiff/makefile

@@ -17,10 +17,6 @@ SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/autodiff/,)
 CONFIG_MK = $(or $(wildcard $(MFEM_BUILD_DIR)/config/config.mk),\
    $(wildcard $(MFEM_INSTALL_DIR)/share/mfem/config.mk))
 
-# Include defaults.mk to get XLINKER
-DEFAULTS_MK = $(MFEM_DIR)/config/defaults.mk
-include $(DEFAULTS_MK)
-
 MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 

miniapps/common/makefile

@@ -20,6 +20,7 @@ CONFIG_MK = $(or $(wildcard $(MFEM_BUILD_DIR)/config/config.mk),\
 # Default target
 all: lib-common
 
+# Include defaults.mk to get the definition of BUILD_SOFLAGS
 DEFAULTS_MK = $(MFEM_DIR)/config/defaults.mk
 include $(DEFAULTS_MK)
 
@@ -30,13 +31,7 @@ MFEM_LIB_FILE = mfem_is_not_built
 ifneq (clean,$(MAKECMDGOALS))
    -include $(CONFIG_MK)
 
-   ifeq ($(MFEM_USE_CUDA),YES)
-      XLINKER = $(CUDA_XLINKER)
-   else ifeq ($(MFEM_USE_HIP),YES)
-      XLINKER = $(HIP_XLINKER)
-   else
-      XLINKER = $(CXX_XLINKER)
-   endif
+   XLINKER = $(MFEM_XLINKER)
 
    BUILD_REAL_DIR = $(realpath .)
    BUILD_SOFLAGS := $(subst libmfem.,libmfem-common.,$(BUILD_SOFLAGS))

miniapps/diag-smoothers/makefile

@@ -15,11 +15,9 @@ MFEM_BUILD_DIR ?= ../..
 
 SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/diag-smoothers/,)
 CONFIG_MK = $(MFEM_BUILD_DIR)/config/config.mk
-DEFAULTS_MK = $(MFEM_DIR)/config/defaults.mk
 
 MFEM_LIB_FILE = mfem_is_not_built
 
--include $(DEFAULTS_MK)
 -include $(CONFIG_MK)
 
 DS_COMMON_SRC = ds-common.cpp
@@ -31,8 +29,7 @@ MINIAPPS = $(if $(MFEM_USE_MPI:NO=),$(PAR_MINIAPPS),)
 
 COMMON_LIB = -L$(MFEM_BUILD_DIR)/miniapps/common -lmfem-common
 COMMON_LIB += $(if $(MFEM_SHARED:YES=),,\
-   $(if $(MFEM_USE_CUDA:YES=),$(CXX_XLINKER),$(CUDA_XLINKER))-rpath,$(abspath\
-   $(MFEM_BUILD_DIR)/miniapps/common))
+   $(MFEM_XLINKER)-rpath,$(abspath $(MFEM_BUILD_DIR)/miniapps/common))
 
 APP_DEPS = $(DS_COMMON_OBJ) $(MFEM_LIB_FILE) $(CONFIG_MK)
 APP_LIBS = $(COMMON_LIB) $(MFEM_LIBS)

miniapps/dpg/makefile

@@ -17,10 +17,6 @@ SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/dpg/,)
 CONFIG_MK = $(or $(wildcard $(MFEM_BUILD_DIR)/config/config.mk),\
    $(wildcard $(MFEM_INSTALL_DIR)/share/mfem/config.mk))
 
-# Include defaults.mk to get XLINKER
-DEFAULTS_MK = $(MFEM_DIR)/config/defaults.mk
-include $(DEFAULTS_MK)
-
 MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
@@ -69,8 +65,7 @@ COMMON_LIB = -L$(MFEM_BUILD_DIR)/miniapps/common -lmfem-common
 
 # If MFEM_SHARED is set, add the ../common rpath
 COMMON_LIB += $(if $(MFEM_SHARED:YES=),,\
-   $(if $(MFEM_USE_CUDA:YES=),$(CXX_XLINKER),$(CUDA_XLINKER))-rpath,$(abspath\
-   $(MFEM_BUILD_DIR)/miniapps/common))
+   $(MFEM_XLINKER)-rpath,$(abspath $(MFEM_BUILD_DIR)/miniapps/common))
 
 .SUFFIXES:
 .SUFFIXES: .o .cpp .mk

miniapps/electromagnetics/CMakeLists.txt

@@ -34,6 +34,11 @@ if (MFEM_USE_MPI)
    EXTRA_HEADERS maxwell_solver.hpp ${MFEM_MINIAPPS_COMMON_HEADERS}
    LIBRARIES mfem-common)
 
+ add_mfem_miniapp(lorentz
+   MAIN lorentz.cpp
+   EXTRA_HEADERS ${MFEM_MINIAPPS_COMMON_HEADERS}
+   LIBRARIES mfem-common)
+
   # Add the corresponding tests to the "test" target
   if (MFEM_ENABLE_TESTING)
     add_test(NAME tesla_np=4

miniapps/electromagnetics/lorentz.cpp

@@ -0,0 +1,571 @@
+// 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.
+//
+//           -----------------------------------------------------
+//           Lorentz Miniapp:  Simple Lorentz Force Particle Mover
+//           -----------------------------------------------------
+//
+// This miniapp computes the trajectory of a single charged particle subject to
+// Lorentz forces.
+//
+//                           dp/dt = q (E + v x B)
+//
+// The method used is the explicit Boris algortihm which conserves phase space
+// volume for long term accuracy.
+//
+// The electric and magnetic fields are read from VisItDataCollection objects
+// such as those produced by the Volta and Tesla miniapps. It is notable that
+// these two fields do not need to be defined on the same mesh. Of course, the
+// particle trajectory can only be computed on the intersection of the two
+// domains. The starting point of the path must be chosen within in this
+// intersection and the trajectory will terminate when it leaves the
+// intersection or reaches a specified time duration.
+//
+// Note that the VisItDataCollection objects must have been stored using the
+// parallel format e.g. visit_dc.SetFormat(DataCollection::PARALLEL_FORMAT);.
+// Without this optional format specifier the vector field lookups will fail.
+//
+// Compile with: make lorentz
+//
+// Sample runs:
+//
+//   Free particle moving with constant velocity
+//      mpirun -np 4 lorentz -p0 '1 1 1'
+//
+//   Particle accelerating in a constant electric field
+//      mpirun -np 4 volta -m ../../data/inline-hex.mesh -dbcs '1 6' -dbcv '0 1'
+//      mpirun -np 4 lorentz -er Volta-AMR-Parallel -x0 '0.5 0.5 0.9' -p0 '1 0 0'
+//
+//   Particle accelerating in a constant magnetic field
+//      mpirun -np 4 tesla -m ../../data/inline-hex.mesh -ubbc '0 0 1'
+//      mpirun -np 4 lorentz -br Tesla-AMR-Parallel -x0 '0.1 0.5 0.1' -p0 '0 0.4 0.1' -tf 9
+//
+//   Magnetic mirror effect near a charged sphere and a bar magnet
+//      mpirun -np 4 volta -m ../../data/ball-nurbs.mesh -dbcs 1 -cs '0 0 0 0.1 2e-11' -rs 2 -maxit 4
+//      mpirun -np 4 tesla -m ../../data/fichera.mesh -maxit 4 -rs 3 -bm '-0.1 -0.1 -0.1 0.1 0.1 0.1 0.1 -1e10'
+//      mpirun -np 4 lorentz -er Volta-AMR-Parallel -ec 4 -br Tesla-AMR-Parallel -bc 4 -x0 '0.8 0 0' -p0 '-8 -4 4' -q -10 -tf 0.2 -dt 1e-3 -rf 1e-6
+//
+// This miniapp demonstrates the use of the ParMesh::FindPoints functionality
+// to evaluate field data from stored DataCollection objects.  While this
+// miniapp is far from a full particle-in-cell (PIC) code it does demonstrate
+// some of the building blocks that might be used to construct the particle
+// mover portion of a PIC code.
+
+#include "mfem.hpp"
+#include "../common/fem_extras.hpp"
+#include "../common/pfem_extras.hpp"
+#include "electromagnetics.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+using namespace mfem::common;
+using namespace mfem::electromagnetics;
+
+typedef DataCollection::FieldMapType fields_t;
+
+/// This class implements the Boris algorithm as described in the
+/// article `Why is Boris algorithm so good?` by H. Qin et al in
+/// Physics of Plasmas, Volume 20 Issue 8, August 2013,
+/// https://doi.org/10.1063/1.4818428.
+class BorisAlgorithm
+{
+private:
+   real_t charge_;
+   real_t mass_;
+
+   ParMesh         *E_pmesh_;
+   ParGridFunction *E_field_;
+
+   ParMesh         *B_pmesh_;
+   ParGridFunction *B_field_;
+
+   mutable Array<int>              elem_id_;
+   mutable Array<IntegrationPoint> ip_;
+
+   mutable Vector E_;
+   mutable Vector B_;
+   mutable Vector pxB_;
+   mutable Vector pm_;
+   mutable Vector pp_;
+
+   // Returns true if a usable V has been found. If @a pgf is NULL, V = 0 is
+   // returned as a default value.
+   bool GetValue(ParMesh *pmesh, ParGridFunction *pgf, Vector q, Vector &V)
+   {
+      DenseMatrix point(q.GetData(), 3, 1);
+
+      int pt_found =
+         (pmesh != NULL) ? pmesh->FindPoints(point, elem_id_, ip_, false) : -1;
+
+      // We have a mesh but the point was not found. The path must be outside
+      // the domain of interest.
+      if (pmesh != NULL && pt_found <= 0) { return false; }
+
+      int pt_root = -1;
+
+      if (pt_found > 0 && elem_id_[0] >= 0 && pgf != NULL)
+      {
+         pt_root = pmesh->GetMyRank();
+
+         pgf->GetVectorValue(elem_id_[0], ip_[0], V);
+      }
+      else
+      {
+         pt_root = 0;
+         V = 0.0;
+      }
+
+      // Determine processor which found the field point
+      int glb_pt_root = -1;
+      MPI_Allreduce(&pt_root, &glb_pt_root, 1,
+                    MPI_INT, MPI_MAX, MPI_COMM_WORLD);
+
+      // Send the field value to the root processor
+      if (pmesh != NULL && elem_id_[0] >= 0 && glb_pt_root != 0)
+      {
+         MPI_Send(V.GetData(), 3, MPITypeMap<real_t>::mpi_type,
+                  0, 1030, MPI_COMM_WORLD);
+      }
+
+      // Receive the field value on the root processor
+      if (Mpi::Root() && pmesh != NULL && glb_pt_root != 0)
+      {
+         MPI_Status status;
+         MPI_Recv(V.GetData(), 3, MPITypeMap<real_t>::mpi_type,
+                  glb_pt_root, 1030, MPI_COMM_WORLD, &status);
+      }
+      return true;
+   }
+
+public:
+   BorisAlgorithm(ParGridFunction *E_gf,
+                  ParGridFunction *B_gf,
+                  real_t charge, real_t mass)
+      : charge_(charge), mass_(mass),
+        E_field_(E_gf),
+        B_field_(B_gf),
+        E_(3), B_(3), pxB_(3), pm_(3), pp_(3)
+   {
+      E_pmesh_ = (E_field_) ? E_field_->ParFESpace()->GetParMesh() : NULL;
+      B_pmesh_ = (B_field_) ? B_field_->ParFESpace()->GetParMesh() : NULL;
+   }
+
+   bool Step(Vector &q, Vector &p, real_t &t, real_t &dt)
+   {
+      // Locate current point in each mesh, evaluate the fields, and collect
+      // field values on the root processor.
+      if (!GetValue(E_pmesh_, E_field_, q, E_)) { return false; }
+      if (!GetValue(B_pmesh_, B_field_, q, B_)) { return false; }
+
+      // Compute updated position and momentum using the Boris algorithm
+      if (Mpi::Root())
+      {
+         // Compute half of the contribution from q E
+         add(p, 0.5 * dt * charge_, E_, pm_);
+
+         // Compute the contributiobn from q p x B
+         const real_t B2 = B_ * B_;
+
+         // ... along pm x B
+         const real_t a1 = 4.0 * dt * charge_ * mass_;
+         pm_.cross3D(B_, pxB_);
+         pp_.Set(a1, pxB_);
+
+         // ... along pm
+         const real_t a2 = 4.0 * mass_ * mass_ -
+                           dt * dt * charge_ * charge_ * B2;
+         pp_.Add(a2, pm_);
+
+         // ... along B
+         const real_t a3 = 2.0 * dt * dt * charge_ * charge_ * (B_ * p);
+         pp_.Add(a3, B_);
+
+         // scale by common denominator
+         const real_t a4 = 4.0 * mass_ * mass_ +
+                           dt * dt * charge_ * charge_ * B2;
+         pp_ /= a4;
+
+         // Update the momentum
+         add(pp_, 0.5 * dt * charge_, E_, p);
+
+         // Update the position
+         q.Add(dt / mass_, p);
+      }
+
+      // Update the time
+      t += dt;
+
+      // Broadcast the updated position
+      MPI_Bcast(q.GetData(), 3, MPITypeMap<real_t>::mpi_type,
+                0, MPI_COMM_WORLD);
+
+      // Broadcast the updated momentum
+      MPI_Bcast(p.GetData(), 3, MPITypeMap<real_t>::mpi_type,
+                0, MPI_COMM_WORLD);
+
+      return true;
+   }
+};
+
+// Open the named VisItDataCollection and read the named field.
+// Returns pointers to the two new objects.
+int ReadGridFunction(const char * coll_name, const char * field_name,
+                     int pad_digits_cycle, int pad_digits_rank, int cycle,
+                     VisItDataCollection *&dc, ParGridFunction *& gf);
+
+// By default the initial position will be the center of the intersection
+// of the bounding boxes of the meshes containing the E and B fields.
+void SetInitialPosition(VisItDataCollection *E_dc,
+                        VisItDataCollection *B_dc,
+                        Vector &x_init);
+
+// Build a quadrilateral mesh approximating the trajectory as a
+// ribbon. One edge of the ribbon follows the trajectory of the
+// particle. The opposite edge is offset by the acceleration vector
+// (scaled by a constant called the r_factor).
+Mesh MakeTrajectoryMesh(int step, real_t m, real_t dt, real_t r_factor,
+                        const DenseMatrix &pos_data,
+                        const DenseMatrix &mom_data);
+
+// Prints the program's logo to the given output stream
+void display_banner(ostream & os);
+
+int main(int argc, char *argv[])
+{
+   Mpi::Init(argc, argv);
+   Hypre::Init();
+
+   if ( Mpi::Root() ) { display_banner(cout); }
+
+   const char *E_coll_name = "";
+   const char *E_field_name = "E";
+   int E_cycle = 10;
+   int E_pad_digits_cycle = 6;
+   int E_pad_digits_rank = 6;
+
+   const char *B_coll_name = "";
+   const char *B_field_name = "B";
+   int B_cycle = 10;
+   int B_pad_digits_cycle = 6;
+   int B_pad_digits_rank = 6;
+
+   real_t q = 1.0;
+   real_t m = 1.0;
+   real_t dt = 1e-2;
+   real_t t_init = 0.0;
+   real_t t_final = 1.0;
+   real_t r_factor = -1.0;
+   Vector x_init;
+   Vector p_init;
+   int visport = 19916;
+   bool visualization = true;
+   bool visit = true;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&E_coll_name, "-er", "--e-root-file",
+                  "Set the VisIt data collection E field root file prefix.");
+   args.AddOption(&E_field_name, "-ef", "--e-field-name",
+                  "Set the VisIt data collection E field name");
+   args.AddOption(&E_cycle, "-ec", "--e-cycle",
+                  "Set the E field cycle index to read.");
+   args.AddOption(&E_pad_digits_cycle, "-epdc", "--e-pad-digits-cycle",
+                  "Number of digits in E field cycle.");
+   args.AddOption(&E_pad_digits_rank, "-epdr", "--e-pad-digits-rank",
+                  "Number of digits in E field MPI rank.");
+   args.AddOption(&B_coll_name, "-br", "--b-root-file",
+                  "Set the VisIt data collection B field root file prefix.");
+   args.AddOption(&B_field_name, "-bf", "--b-field-name",
+                  "Set the VisIt data collection B field name");
+   args.AddOption(&B_cycle, "-bc", "--b-cycle",
+                  "Set the B field cycle index to read.");
+   args.AddOption(&B_pad_digits_cycle, "-bpdc", "--b-pad-digits-cycle",
+                  "Number of digits in B field cycle.");
+   args.AddOption(&B_pad_digits_rank, "-bpdr", "--b-pad-digits-rank",
+                  "Number of digits in B field MPI rank.");
+   args.AddOption(&q, "-q", "--charge",
+                  "Particle charge.");
+   args.AddOption(&m, "-m", "--mass",
+                  "Particle mass.");
+   args.AddOption(&dt, "-dt", "--time-step",
+                  "Time Step.");
+   args.AddOption(&t_init, "-ti", "--initial-time",
+                  "Initial Time.");
+   args.AddOption(&t_final, "-tf", "--final-time",
+                  "Final Time.");
+   args.AddOption(&x_init, "-x0", "--initial-position",
+                  "Initial position.");
+   args.AddOption(&p_init, "-p0", "--initial-momentum",
+                  "Initial momentum.");
+   args.AddOption(&r_factor, "-rf", "--ribbon-factor",
+                  "Scale factor for ribbon width (rf * (p1-p0) / (m * dt) "
+                  "where p0 and p1 are computed momenta).");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&visit, "-visit", "--visit", "-no-visit", "--no-visit",
+                  "Enable or disable VisIt visualization.");
+   args.AddOption(&visport, "-p", "--send-port", "Socket for GLVis.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (Mpi::Root())
+      {
+         args.PrintUsage(cout);
+      }
+      return 1;
+   }
+   if (r_factor <= 0.0)
+   {
+      r_factor = dt;
+   }
+   if (Mpi::Root())
+   {
+      args.PrintOptions(cout);
+   }
+
+   VisItDataCollection *E_dc = NULL;
+   ParGridFunction     *E_gf = NULL;
+
+   if (strcmp(E_coll_name, ""))
+   {
+      if (ReadGridFunction(E_coll_name, E_field_name, E_pad_digits_cycle,
+                           E_pad_digits_rank, E_cycle, E_dc, E_gf))
+      {
+         mfem::out << "Error loading E field" << endl;
+         return 1;
+      }
+   }
+
+   VisItDataCollection *B_dc = NULL;
+   ParGridFunction     *B_gf = NULL;
+
+   if (strcmp(B_coll_name, ""))
+   {
+      if (ReadGridFunction(B_coll_name, B_field_name, B_pad_digits_cycle,
+                           B_pad_digits_rank, B_cycle, B_dc, B_gf))
+      {
+         mfem::out << "Error loading B field" << endl;
+         return 1;
+      }
+   }
+
+   if (x_init.Size() < 3)
+   {
+      SetInitialPosition(E_dc, B_dc, x_init);
+   }
+   if (p_init.Size() < 3)
+   {
+      p_init.SetSize(3); p_init = 0.0;
+   }
+   if (Mpi::Root())
+   {
+      mfem::out << "Initial position: "; x_init.Print(mfem::out);
+      mfem::out << "Initial momentum: "; p_init.Print(mfem::out);
+   }
+
+   BorisAlgorithm boris(E_gf, B_gf, q, m);
+   Vector pos(x_init);
+   Vector mom(p_init);
+
+   ofstream ofs("Lorentz.dat");
+   ofs.precision(14);
+
+   int nsteps = 1 + (int)ceil((t_final - t_init) / dt);
+   DenseMatrix pos_data(3, nsteps);
+   DenseMatrix mom_data(3, nsteps + 1);
+   mom_data.SetCol(0, p_init);
+
+   if (Mpi::Root())
+   {
+      mfem::out << "Maximum number of steps: " << nsteps << endl;
+   }
+
+   int step = -1;
+   real_t t = t_init;
+   do
+   {
+      if (Mpi::Root())
+      {
+         ofs << t
+             << '\t' << pos[0] << '\t' << pos[1] << '\t' << pos[2]
+             << '\t' << mom[0] << '\t' << mom[1] << '\t' << mom[2]
+             << '\n';
+      }
+      step++;
+
+      pos_data.SetCol(step, pos);
+      mom_data.SetCol(step + 1, mom);
+   }
+   while (boris.Step(pos, mom, t, dt) && step < nsteps - 1);
+
+   if (Mpi::Root() && (visit || visualization))
+   {
+      Mesh trajectory = MakeTrajectoryMesh(step, m, dt, r_factor,
+                                           pos_data, mom_data);
+
+      L2_FECollection    fec_l2(0, 2);
+      FiniteElementSpace fes_l2(&trajectory, &fec_l2);
+      GridFunction traj_time(&fes_l2);
+      for (int i=0; i<step; i++)
+      {
+         traj_time[i] = dt * i;
+      }
+
+      if (visit)
+      {
+         VisItDataCollection visit_dc("Lorentz", &trajectory);
+         visit_dc.RegisterField("Time", &traj_time);
+         visit_dc.SetCycle(step);
+         visit_dc.SetTime(step * dt);
+         visit_dc.Save();
+      }
+
+      if (visualization)
+      {
+         socketstream traj_sock;
+         traj_sock.precision(8);
+
+         char vishost[] = "localhost";
+
+         int Wx = 0, Wy = 0; // window position
+         int Ww = 350, Wh = 350; // window size
+
+         VisualizeField(traj_sock, vishost, visport,
+                        traj_time, "Trajectory", Wx, Wy, Ww, Wh);
+      }
+   }
+   if (Mpi::Root())
+   {
+      mfem::out << "Number of steps taken: " << step << endl;
+   }
+
+   // Clean up
+   delete E_dc;
+   delete B_dc;
+}
+
+// Print the Lorentz ascii logo to the given ostream
+void display_banner(ostream & os)
+{
+   os << "   ____                                __          "
+      << endl
+      << "  |    |    ___________   ____   _____/  |_________"
+      << endl
+      << "  |    |   /  _ \\_  __ \\_/ __ \\ /    \\   __\\___   /"
+      << endl
+      << "  |    |__(  <_> )  | \\/\\  ___/|   |  \\  |  /    / "
+      << endl
+      << "  |_______ \\____/|__|    \\___  >___|  /__| /_____ \\"
+      << endl
+      << "          \\/                 \\/     \\/           \\/"
+      << endl << flush;
+}
+
+int ReadGridFunction(const char * coll_name, const char * field_name,
+                     int pad_digits_cycle, int pad_digits_rank, int cycle,
+                     VisItDataCollection *&dc, ParGridFunction *& gf)
+{
+   dc = new VisItDataCollection(MPI_COMM_WORLD, coll_name);
+   dc->SetPadDigitsCycle(pad_digits_cycle);
+   dc->SetPadDigitsRank(pad_digits_rank);
+   dc->Load(cycle);
+
+   if (dc->Error() != DataCollection::No_Error)
+   {
+      mfem::out << "Error loading VisIt data collection: "
+                << coll_name << endl;
+      return 1;
+   }
+
+   if (dc->GetMesh()->Dimension() < 3)
+   {
+      mfem::out << "Field must be defined on a three dimensional mesh"
+                << endl;
+      return 1;
+   }
+
+   if (dc->HasField(field_name))
+   {
+      gf = dc->GetParField(field_name);
+   }
+
+   return 0;
+}
+
+void SetInitialPosition(VisItDataCollection *E_dc,
+                        VisItDataCollection *B_dc,
+                        Vector &x_init)
+{
+   x_init.SetSize(3); x_init = 0.0;
+
+   if (E_dc != NULL || B_dc != NULL)
+   {
+      Vector E_p_min(3); E_p_min = -infinity();
+      Vector E_p_max(3); E_p_max = infinity();
+      if (E_dc != NULL)
+      {
+         ParMesh * E_pmesh = dynamic_cast<ParMesh*>(E_dc->GetMesh());
+         E_pmesh->GetBoundingBox(E_p_min, E_p_max);
+      }
+
+      Vector B_p_min(3); B_p_min = -infinity();
+      Vector B_p_max(3); B_p_max = infinity();
+      if (B_dc != NULL)
+      {
+         ParMesh *B_pmesh = dynamic_cast<ParMesh*>(B_dc->GetMesh());
+         B_pmesh->GetBoundingBox(B_p_min, B_p_max);
+      }
+
+      for (int d = 0; d<3; d++)
+      {
+         const real_t p_min = std::max(E_p_min[d], B_p_min[d]);
+         const real_t p_max = std::min(E_p_max[d], B_p_max[d]);
+         x_init[d] = 0.5 * (p_min + p_max);
+      }
+   }
+}
+
+Mesh MakeTrajectoryMesh(int step, real_t m, real_t dt, real_t r_factor,
+                        const DenseMatrix &pos_data,
+                        const DenseMatrix &mom_data)
+{
+   Mesh trajectory(2, 2 * (step + 1), step, 0, 3);
+
+   for (int i=0; i<=step; i++)
+   {
+      trajectory.AddVertex(pos_data(0,i), pos_data(1,i), pos_data(2,i));
+
+      real_t dpx = (mom_data(0, i + 1) - mom_data(0, i)) / (m * dt);
+      real_t dpy = (mom_data(1, i + 1) - mom_data(1, i)) / (m * dt);
+      real_t dpz = (mom_data(2, i + 1) - mom_data(2, i)) / (m * dt);
+
+      trajectory.AddVertex(pos_data(0,i) + r_factor * dpx,
+                           pos_data(1,i) + r_factor * dpy,
+                           pos_data(2,i) + r_factor * dpz);
+   }
+
+   int v[4];
+   for (int i=0; i<step; i++)
+   {
+      v[0] = 2 * i;
+      v[1] = 2 * (i + 1);
+      v[2] = 2 * (i + 1) + 1;
+      v[3] = 2 * i + 1;
+
+      trajectory.AddQuad(v);
+   }
+
+   trajectory.FinalizeQuadMesh(1);
+
+   return trajectory;
+}

miniapps/electromagnetics/makefile

@@ -17,15 +17,11 @@ SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/electromagnetics/,)
 CONFIG_MK = $(or $(wildcard $(MFEM_BUILD_DIR)/config/config.mk),\
    $(wildcard $(MFEM_INSTALL_DIR)/share/mfem/config.mk))
 
-# Include defaults.mk to get XLINKER
-DEFAULTS_MK = $(MFEM_DIR)/config/defaults.mk
-include $(DEFAULTS_MK)
-
 MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
 SEQ_MINIAPPS =
-PAR_MINIAPPS = volta tesla maxwell joule
+PAR_MINIAPPS = volta tesla maxwell joule lorentz
 ifeq ($(MFEM_USE_MPI),NO)
    MINIAPPS = $(SEQ_MINIAPPS)
 else
@@ -41,8 +37,7 @@ COMMON_LIB = -L$(MFEM_BUILD_DIR)/miniapps/common -lmfem-common
 
 # If MFEM_SHARED is set, add the ../common rpath
 COMMON_LIB += $(if $(MFEM_SHARED:YES=),,\
-   $(if $(MFEM_USE_CUDA:YES=),$(CXX_XLINKER),$(CUDA_XLINKER))-rpath,$(abspath\
-   $(MFEM_BUILD_DIR)/miniapps/common))
+   $(MFEM_XLINKER)-rpath,$(abspath $(MFEM_BUILD_DIR)/miniapps/common))
 
 # Remove built-in rules
 %: %.cpp
@@ -56,6 +51,10 @@ all: $(MINIAPPS)
 	$(MFEM_CXX) $(MFEM_LINK_FLAGS) -o $@ $@.o $@_solver.o $(COMMON_LIB) \
 	  $(MFEM_LIBS)
 
+lorentz: %: $(SRC)%.cpp $(MFEM_LIB_FILE) $(CONFIG_MK) | lib-common
+	$(MFEM_CXX) $(MFEM_FLAGS) -c $(<)
+	$(MFEM_CXX) $(MFEM_LINK_FLAGS) -o $@ $@.o $(COMMON_LIB) $(MFEM_LIBS)
+
 # Rules for compiling miniapp dependencies
 $(addsuffix _solver.o,$(MINIAPPS)): \
 %.o: $(SRC)%.cpp $(SRC)%.hpp $(CONFIG_MK)
@@ -86,7 +85,7 @@ include $(MFEM_TEST_MK)
 
 # Testing: Specific execution options
 RUN_MPI = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) $(MFEM_MPI_NP)
-volta-test-par: volta-test-1 volta-test-2
+volta-test-par: volta-test-1 volta-test-2 volta-test-3
 volta-test-1: volta
 	@$(call mfem-test,$<, $(RUN_MPI), Electromagnetic miniapp,\
 	-maxit 2 -dbcs 1 -dbcg -ds '0.0 0.0 0.0 0.2 8.0')
@@ -94,15 +93,29 @@ volta-test-2: volta
 	@$(call mfem-test,$<, $(RUN_MPI), Electromagnetic miniapp,\
 	-maxit 2 -m ../../data/square-disc.mesh \
 	-dbcs '1 2 3 4 5 6 7 8' -dbcv '0 0 0 0 1 1 1 1')
-tesla-test-par: tesla
+volta-test-3: volta
+	@$(call mfem-test,$<, $(RUN_MPI), Electromagnetic miniapp,\
+	-maxit 2 -m ../../data/inline-hex.mesh -dbcs '1 6' -dbcv '0 1')
+tesla-test-par: tesla-test-1 tesla-test-2
+tesla-test-1: tesla
 	@$(call mfem-test,$<, $(RUN_MPI), Electromagnetic miniapp,\
 	-maxit 2 -cr '0 0 -0.2 0 0 0.2 0.2 0.4 1')
+tesla-test-2: tesla
+	@$(call mfem-test,$<, $(RUN_MPI), Electromagnetic miniapp,\
+	-maxit 2 -m ../../data/inline-hex.mesh -ubbc '0 0 1')
 maxwell-test-par: maxwell
 	@$(call mfem-test,$<, $(RUN_MPI), Electromagnetic miniapp,\
 	-abcs '-1' -dp '-0.3 0.0 0.0 0.3 0.0 0.0 0.1 1 .5 .5')
 joule-test-par: joule
 	@$(call mfem-test,$<, $(RUN_MPI), Electromagnetic miniapp,\
 	-m cylinder-hex.mesh -p rod -tf 3)
+lorentz-test-par: lorentz-test-1 lorentz-test-2
+lorentz-test-1: lorentz volta-test-3
+	@$(call mfem-test,$<, $(RUN_MPI), Electromagnetic miniapp,\
+	-er Volta-AMR-Parallel -ec 2 -x0 '0.5 0.5 0.9' -p0 '1 0 0')
+lorentz-test-2: lorentz tesla-test-2
+	@$(call mfem-test,$<, $(RUN_MPI), Electromagnetic miniapp,\
+	-br Tesla-AMR-Parallel -bc 2 -x0 '0.1 0.5 0.1' -p0 '0 0.4 0.1' -tf 9)
 
 # Testing: "test" target and mfem-test* variables are defined in config/test.mk
 
@@ -117,4 +130,4 @@ clean-build:
 	rm -rf *.dSYM *.TVD.*breakpoints
 
 clean-exec:
-	@rm -rf Volta-AMR* Tesla-AMR* Maxwell-Parallel* Joule_*
+	@rm -rf Volta-AMR* Tesla-AMR* Maxwell-Parallel* Joule_* Lorentz*

miniapps/electromagnetics/tesla.cpp

@@ -253,6 +253,7 @@ int main(int argc, char *argv[])
 
    // Initialize VisIt visualization
    VisItDataCollection visit_dc("Tesla-AMR-Parallel", &pmesh);
+   visit_dc.SetFormat(DataCollection::PARALLEL_FORMAT);
 
    if ( visit )
    {

miniapps/electromagnetics/volta.cpp

@@ -266,6 +266,7 @@ int main(int argc, char *argv[])
 
    // Initialize VisIt visualization
    VisItDataCollection visit_dc("Volta-AMR-Parallel", &pmesh);
+   visit_dc.SetFormat(DataCollection::PARALLEL_FORMAT);
 
    if ( visit )
    {

miniapps/gslib/makefile

@@ -17,10 +17,6 @@ SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/gslib/,)
 CONFIG_MK = $(or $(wildcard $(MFEM_BUILD_DIR)/config/config.mk),\
    $(wildcard $(MFEM_INSTALL_DIR)/share/mfem/config.mk))
 
-# Include defaults.mk to get XLINKER
-DEFAULTS_MK = $(MFEM_DIR)/config/defaults.mk
-include $(DEFAULTS_MK)
-
 MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
@@ -42,8 +38,7 @@ COMMON_LIB = -L$(MFEM_BUILD_DIR)/miniapps/common -lmfem-common
 
 # If MFEM_SHARED is set, add the ../common rpath
 COMMON_LIB += $(if $(MFEM_SHARED:YES=),,\
-   $(if $(MFEM_USE_CUDA:YES=),$(CXX_XLINKER),$(CUDA_XLINKER))-rpath,$(abspath\
-   $(MFEM_BUILD_DIR)/miniapps/common))
+   $(MFEM_XLINKER)-rpath,$(abspath $(MFEM_BUILD_DIR)/miniapps/common))
 
 .SUFFIXES:
 .SUFFIXES: .o .cpp .mk

miniapps/meshing/makefile

@@ -17,10 +17,6 @@ SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/meshing/,)
 CONFIG_MK = $(or $(wildcard $(MFEM_BUILD_DIR)/config/config.mk),\
    $(wildcard $(MFEM_INSTALL_DIR)/share/mfem/config.mk))
 
-# Include defaults.mk to get XLINKER
-DEFAULTS_MK = $(MFEM_DIR)/config/defaults.mk
-include $(DEFAULTS_MK)
-
 MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
@@ -39,8 +35,7 @@ COMMON_LIB = -L$(MFEM_BUILD_DIR)/miniapps/common -lmfem-common
 
 # If MFEM_SHARED is set, add the ../common rpath
 COMMON_LIB += $(if $(MFEM_SHARED:YES=),,\
-   $(if $(MFEM_USE_CUDA:YES=),$(CXX_XLINKER),$(CUDA_XLINKER))-rpath,$(abspath\
-   $(MFEM_BUILD_DIR)/miniapps/common))
+   $(MFEM_XLINKER)-rpath,$(abspath $(MFEM_BUILD_DIR)/miniapps/common))
 
 .SUFFIXES:
 .SUFFIXES: .o .cpp .mk

miniapps/mtop/makefile

@@ -17,10 +17,6 @@ SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/mtop/,)
 CONFIG_MK = $(or $(wildcard $(MFEM_BUILD_DIR)/config/config.mk),\
    $(wildcard $(MFEM_INSTALL_DIR)/share/mfem/config.mk))
 
-# Include defaults.mk to get XLINKER
-DEFAULTS_MK = $(MFEM_DIR)/config/defaults.mk
-include $(DEFAULTS_MK)
-
 MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 

miniapps/navier/navier_cht.cpp

@@ -461,7 +461,7 @@ void ConductionOperator::Mult(const Vector &u, Vector &du_dt) const
 
    Kmat.Mult(u, z);
    z.Neg(); // z = -z
-   K->EliminateVDofsInRHS(ess_tdof_list, u, z);
+   K->ParallelEliminateTDofsInRHS(ess_tdof_list, u, z);
 
    M_solver.Mult(z, du_dt);
    du_dt.Print();
@@ -483,7 +483,7 @@ void ConductionOperator::ImplicitSolve(const real_t dt,
    MFEM_VERIFY(dt == current_dt, ""); // SDIRK methods use the same dt
    Kmat.Mult(u, z);
    z.Neg();
-   K->EliminateVDofsInRHS(ess_tdof_list, u, z);
+   K->ParallelEliminateTDofsInRHS(ess_tdof_list, u, z);
 
    T_solver.Mult(z, du_dt);
    du_dt.SetSubVector(ess_tdof_list, 0.0);

miniapps/shifted/distance.cpp

@@ -87,6 +87,8 @@
 //
 //   Problem 4: level set: Union of doughnut and swiss cheese shapes
 //     mpirun -np 4 distance -m ../../data/inline-hex.mesh -rs 3 -o 2 -t 1.0 -p 4
+//   Problem 5: point source in mfem mesh.
+//     mpirun -np 4 distance -m ../../data/mfem.mesh -p 5 -rs 3 -t 300.0
 
 #include <fstream>
 #include <iostream>
@@ -233,7 +235,8 @@ int main(int argc, char *argv[])
                   "1: Circle / sphere level set in 2D / 3D\n\t"
                   "2: 2D sine-looking level set\n\t"
                   "3: Gyroid level set in 2D or 3D\n\t"
-                  "4: Combo of a doughnut and swiss cheese shapes in 3D.");
+                  "4: Combo of a doughnut and swiss cheese shapes in 3D.\n\t"
+                  "5: Point source in MFEM mesh.");
    args.AddOption(&rs_levels, "-rs", "--refine-serial",
                   "Number of times to refine the mesh uniformly in serial.");
    args.AddOption(&order, "-o", "--order",
@@ -299,6 +302,11 @@ int main(int argc, char *argv[])
       ls_coeff = new FunctionCoefficient(doughnut_cheese);
       smooth_steps = 0;
    }
+   else if (problem == 5)
+   {
+      ls_coeff = new DeltaCoefficient(0.0, 0.0, 1000.0);
+      smooth_steps = 0;
+   }
    else { MFEM_ABORT("Unrecognized -problem option."); }
 
    const real_t dx = AvgElementSize(pmesh);
@@ -306,7 +314,7 @@ int main(int argc, char *argv[])
    if (solver_type == 0)
    {
       auto ds = new HeatDistanceSolver(t_param * dx * dx);
-      if (problem == 0)
+      if (problem == 0 || problem == 5)
       {
          ds->transform = false;
       }
@@ -334,7 +342,7 @@ int main(int argc, char *argv[])
    // Smooth-out Gibbs oscillations from the input level set. The smoothing
    // parameter here is specified to be mesh dependent with length scale dx.
    ParGridFunction filt_gf(&pfes_s);
-   if (problem != 0)
+   if (problem != 0 && problem != 5)
    {
       real_t filter_weight = dx;
       // The normalization-based solver needs a more diffused input.

miniapps/shifted/makefile

@@ -17,10 +17,6 @@ SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/shifted/,)
 CONFIG_MK = $(or $(wildcard $(MFEM_BUILD_DIR)/config/config.mk),\
    $(wildcard $(MFEM_INSTALL_DIR)/share/mfem/config.mk))
 
-# Include defaults.mk to get XLINKER
-DEFAULTS_MK = $(MFEM_DIR)/config/defaults.mk
-include $(DEFAULTS_MK)
-
 MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
@@ -46,8 +42,7 @@ COMMON_LIB = -L$(MFEM_BUILD_DIR)/miniapps/common -lmfem-common
 
 # If MFEM_SHARED is set, add the ../common rpath
 COMMON_LIB += $(if $(MFEM_SHARED:YES=),,\
-   $(if $(MFEM_USE_CUDA:YES=),$(CXX_XLINKER),$(CUDA_XLINKER))-rpath,$(abspath\
-   $(MFEM_BUILD_DIR)/miniapps/common))
+   $(MFEM_XLINKER)-rpath,$(abspath $(MFEM_BUILD_DIR)/miniapps/common))
 
 .SUFFIXES:
 .SUFFIXES: .o .cpp .mk

miniapps/spde/makefile

@@ -17,10 +17,6 @@ SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/spde/,)
 CONFIG_MK = $(or $(wildcard $(MFEM_BUILD_DIR)/config/config.mk),\
    $(wildcard $(MFEM_INSTALL_DIR)/share/mfem/config.mk))
 
-# Include defaults.mk to get XLINKER
-DEFAULTS_MK = $(MFEM_DIR)/config/defaults.mk
-include $(DEFAULTS_MK)
-
 MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 

miniapps/tools/makefile

@@ -17,10 +17,6 @@ SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/tools/,)
 CONFIG_MK = $(or $(wildcard $(MFEM_BUILD_DIR)/config/config.mk),\
    $(wildcard $(MFEM_INSTALL_DIR)/share/mfem/config.mk))
 
-# Include defaults.mk to get XLINKER
-DEFAULTS_MK = $(MFEM_DIR)/config/defaults.mk
-include $(DEFAULTS_MK)
-
 MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
@@ -44,8 +40,7 @@ COMMON_LIB = -L$(MFEM_BUILD_DIR)/miniapps/common -lmfem-common
 
 # If MFEM_SHARED is set, add the ../common rpath
 COMMON_LIB += $(if $(MFEM_SHARED:YES=),,\
-   $(if $(MFEM_USE_CUDA:YES=),$(CXX_XLINKER),$(CUDA_XLINKER))-rpath,$(abspath\
-   $(MFEM_BUILD_DIR)/miniapps/common))
+   $(MFEM_XLINKER)-rpath,$(abspath $(MFEM_BUILD_DIR)/miniapps/common))
 
 all: $(MINIAPPS)
 

miniapps/toys/makefile

@@ -17,10 +17,6 @@ SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/toys/,)
 CONFIG_MK = $(or $(wildcard $(MFEM_BUILD_DIR)/config/config.mk),\
    $(wildcard $(MFEM_INSTALL_DIR)/share/mfem/config.mk))
 
-# Include defaults.mk to get XLINKER
-DEFAULTS_MK = $(MFEM_DIR)/config/defaults.mk
-include $(DEFAULTS_MK)
-
 MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
@@ -41,8 +37,7 @@ COMMON_LIB = -L$(MFEM_BUILD_DIR)/miniapps/common -lmfem-common
 
 # If MFEM_SHARED is set, add the ../common rpath
 COMMON_LIB += $(if $(MFEM_SHARED:YES=),,\
-   $(if $(MFEM_USE_CUDA:YES=),$(CXX_XLINKER),$(CUDA_XLINKER))-rpath,$(abspath\
-   $(MFEM_BUILD_DIR)/miniapps/common))
+   $(MFEM_XLINKER)-rpath,$(abspath $(MFEM_BUILD_DIR)/miniapps/common))
 
 all: $(MINIAPPS)
 

tests/unit/general/test_reduction.cpp

@@ -22,7 +22,6 @@ using namespace mfem;
 
 TEST_CASE("Reduce Sum", "[Reduction],[GPU]")
 {
-   Array<int> workspace;
    Array<int> a(1000);
    a.HostReadWrite();
    for (int i = 0; i < a.Size(); ++i)
@@ -36,7 +35,7 @@ TEST_CASE("Reduce Sum", "[Reduction],[GPU]")
       int res = 0;
       mfem::reduce(
       a.Size(), res, [=] MFEM_HOST_DEVICE(int i, int &r) { r += dptr[i]; },
-      SumReducer<int> {}, use_dev, workspace);
+      SumReducer<int> {}, use_dev);
       // correct for even-length summations
       int expected = (AsConst(a)[0] + AsConst(a)[a.Size() - 1]) * a.Size() / 2;
       CAPTURE(use_dev);
@@ -46,7 +45,6 @@ TEST_CASE("Reduce Sum", "[Reduction],[GPU]")
 
 TEST_CASE("Reduce Mult", "[Reduction],[GPU]")
 {
-   Array<long long> workspace;
    Array<long long> a(64);
    a.HostReadWrite();
    for (int i = 0; i < a.Size(); ++i)
@@ -64,7 +62,7 @@ TEST_CASE("Reduce Mult", "[Reduction],[GPU]")
          mfem::reduce(
             a.Size(), res,
          [=] MFEM_HOST_DEVICE(int i, long long &r) { r *= dptr[i]; },
-         MultReducer<long long> {}, use_dev, workspace);
+         MultReducer<long long> {}, use_dev);
          long long expected = 0;
          CAPTURE(use_dev);
          REQUIRE(res == expected);
@@ -76,7 +74,7 @@ TEST_CASE("Reduce Mult", "[Reduction],[GPU]")
          mfem::reduce(
             a.Size(), res,
          [=] MFEM_HOST_DEVICE(int i, long long &r) { r *= dptr[i]; },
-         MultReducer<long long> {}, use_dev, workspace);
+         MultReducer<long long> {}, use_dev);
          long long expected = 21936950640377856;
          CAPTURE(use_dev);
          REQUIRE(res == expected);
@@ -86,7 +84,6 @@ TEST_CASE("Reduce Mult", "[Reduction],[GPU]")
 
 TEST_CASE("Reduce BAnd", "[Reduction],[GPU]")
 {
-   Array<unsigned> workspace;
    Array<unsigned> a(10);
    SECTION("{ Bit unset }")
    {
@@ -108,7 +105,7 @@ TEST_CASE("Reduce BAnd", "[Reduction],[GPU]")
          mfem::reduce(
             a.Size(), res,
          [=] MFEM_HOST_DEVICE(int i, unsigned &r) { r &= dptr[i]; },
-         BAndReducer<unsigned> {}, use_dev, workspace);
+         BAndReducer<unsigned> {}, use_dev);
          CAPTURE(use_dev);
          REQUIRE(res == ((~1u) & ~(1u << unset_bit)));
          REQUIRE((res & (1u << unset_bit)) == 0);
@@ -132,7 +129,7 @@ TEST_CASE("Reduce BAnd", "[Reduction],[GPU]")
          mfem::reduce(
             a.Size(), res,
          [=] MFEM_HOST_DEVICE(int i, unsigned &r) { r &= dptr[i]; },
-         BAndReducer<unsigned> {}, use_dev, workspace);
+         BAndReducer<unsigned> {}, use_dev);
          CAPTURE(use_dev);
          REQUIRE(res == (1u << set_bit));
       }
@@ -141,7 +138,6 @@ TEST_CASE("Reduce BAnd", "[Reduction],[GPU]")
 
 TEST_CASE("Reduce BOr", "[Reduction],[GPU]")
 {
-   Array<unsigned> workspace;
    Array<unsigned> a(0x210);
    a.HostReadWrite();
    for (int i = 0; i < a.Size(); ++i)
@@ -157,7 +153,7 @@ TEST_CASE("Reduce BOr", "[Reduction],[GPU]")
       mfem::reduce(
          a.Size(), res,
       [=] MFEM_HOST_DEVICE(int i, unsigned &r) { r |= dptr[i]; },
-      BOrReducer<unsigned> {}, use_dev, workspace);
+      BOrReducer<unsigned> {}, use_dev);
       CAPTURE(use_dev);
       REQUIRE(res == 0x3ffu);
    }
@@ -165,7 +161,6 @@ TEST_CASE("Reduce BOr", "[Reduction],[GPU]")
 
 TEST_CASE("Reduce Min", "[Reduction],[GPU]")
 {
-   Array<int> workspace;
    Array<int> a(1000);
    auto hptr = a.HostReadWrite();
    for (int i = 0; i < a.Size(); ++i)
@@ -190,7 +185,7 @@ TEST_CASE("Reduce Min", "[Reduction],[GPU]")
             r = dptr[i];
          }
       },
-      MinReducer<int> {}, use_dev, workspace);
+      MinReducer<int> {}, use_dev);
       CAPTURE(use_dev);
       REQUIRE(res == -10);
    }
@@ -198,7 +193,6 @@ TEST_CASE("Reduce Min", "[Reduction],[GPU]")
 
 TEST_CASE("Reduce Max", "[Reduction],[GPU]")
 {
-   Array<int> workspace;
    Array<int> a(1000);
    auto hptr = a.HostReadWrite();
    for (int i = 0; i < a.Size(); ++i)
@@ -223,7 +217,7 @@ TEST_CASE("Reduce Max", "[Reduction],[GPU]")
             r = dptr[i];
          }
       },
-      MaxReducer<int> {}, use_dev, workspace);
+      MaxReducer<int> {}, use_dev);
       CAPTURE(use_dev);
       REQUIRE(res == 999 - 10);
    }
@@ -231,7 +225,6 @@ TEST_CASE("Reduce Max", "[Reduction],[GPU]")
 
 TEST_CASE("Reduce MinMax", "[Reduction],[GPU]")
 {
-   Array<DevicePair<int, int>> workspace;
    Array<int> a(1000);
    auto hptr = a.HostReadWrite();
    for (int i = 0; i < a.Size(); ++i)
@@ -262,7 +255,7 @@ TEST_CASE("Reduce MinMax", "[Reduction],[GPU]")
             r.second = dptr[i];
          }
       },
-      MinMaxReducer<int> {}, use_dev, workspace);
+      MinMaxReducer<int> {}, use_dev);
       CAPTURE(use_dev);
       REQUIRE(res.first == -10);
       REQUIRE(res.second == a.Size() - 11);
@@ -271,7 +264,6 @@ TEST_CASE("Reduce MinMax", "[Reduction],[GPU]")
 
 TEST_CASE("Reduce ArgMin", "[Reduction],[GPU]")
 {
-   Array<DevicePair<double, int>> workspace;
    Array<double> a(1000);
    auto hptr = a.HostReadWrite();
    for (int i = 0; i < a.Size(); ++i)
@@ -297,7 +289,7 @@ TEST_CASE("Reduce ArgMin", "[Reduction],[GPU]")
             r.second = i;
          }
       },
-      ArgMinReducer<double, int> {}, use_dev, workspace);
+      ArgMinReducer<double, int> {}, use_dev);
       CAPTURE(use_dev);
       REQUIRE(res.first == -10);
       REQUIRE(res.second >= 0);
@@ -308,7 +300,6 @@ TEST_CASE("Reduce ArgMin", "[Reduction],[GPU]")
 
 TEST_CASE("Reduce ArgMax", "[Reduction],[GPU]")
 {
-   Array<DevicePair<double, int>> workspace;
    Array<double> a(1000);
 
    auto hptr = a.HostReadWrite();
@@ -337,7 +328,7 @@ TEST_CASE("Reduce ArgMax", "[Reduction],[GPU]")
             r.second = i;
          }
       },
-      ArgMaxReducer<double, int> {}, use_dev, workspace);
+      ArgMaxReducer<double, int> {}, use_dev);
       CAPTURE(use_dev);
       REQUIRE(res.first == a.Size() - 11);
       REQUIRE(res.second >= 0);
@@ -348,7 +339,6 @@ TEST_CASE("Reduce ArgMax", "[Reduction],[GPU]")
 
 TEST_CASE("Reduce ArgMinMax", "[Reduction],[GPU]")
 {
-   Array<MinMaxLocScalar<double, int>> workspace;
    Array<double> a(1000);
    auto hptr = a.HostReadWrite();
    for (int i = 0; i < a.Size(); ++i)
@@ -383,7 +373,7 @@ TEST_CASE("Reduce ArgMinMax", "[Reduction],[GPU]")
             r.max_loc = i;
          }
       },
-      ArgMinMaxReducer<double, int> {}, use_dev, workspace);
+      ArgMinMaxReducer<double, int> {}, use_dev);
       CAPTURE(use_dev);
       REQUIRE(res.min_val == -10);
       REQUIRE(res.min_loc >= 0);

tests/unit/linalg/test_matrix_dense.cpp

@@ -373,6 +373,7 @@ TEST_CASE("Batched Linear Algebra",
    const int n_rhs = 2;
 
    DenseTensor A_batch(n, n, n_mat);
+   DenseTensor A_inv_batch(n, n, n_mat);
    Vector x_batch(n * n_rhs * n_mat), y_batch(n * n_rhs * n_mat);
    std::vector<DenseMatrix> As;
    std::vector<DenseMatrix> xs, ys;
@@ -404,6 +405,7 @@ TEST_CASE("Batched Linear Algebra",
       ys.back() = 0.0;
       AddMult_a(1.5, As.back(), xs.back(), ys.back());
       A_batch(i) = As.back();
+      A_inv_batch(i) = As.back();
    }
 
    // Test batched matrix-vector products
@@ -463,6 +465,33 @@ TEST_CASE("Batched Linear Algebra",
          }
       }
    }
+
+   // Test batched matrix inverse
+   BatchedLinAlg::Get(backend).Invert(A_inv_batch);
+   A_inv_batch.HostReadWrite();
+   Vector output_col(n);
+   Vector col;
+   for (int i = 0; i < n_mat; ++i)
+   {
+      DenseMatrix Ai_inv(A_inv_batch(i));
+      for (int j = 0; j < n; ++j)
+      {
+         output_col = 0.0;
+         As[i].GetColumnReference(j, col);
+         Ai_inv.Mult(col, output_col);
+         for (int k = 0; k < n; ++k)
+         {
+            if (j == k)
+            {
+               REQUIRE(output_col(k) == MFEM_Approx(1.0));
+            }
+            else
+            {
+               REQUIRE(output_col(k) == MFEM_Approx(0.0));
+            }
+         }
+      }
+   }
 }
 
 TEST_CASE("DenseTensor copy", "[DenseMatrix][DenseTensor]")

tests/unit/linalg/test_matrix_sparse.cpp

@@ -240,4 +240,54 @@ TEST_CASE("SparseMatrix printing", "[SparseMatrix]")
    }
 }
 
+TEST_CASE("SparseMatrix cuSPARSE Bug", "[SparseMatrix][GPU]")
+{
+   // This test case ensures that we have a functioning workaround for the bug
+   // CUSPARSE-1897. In versions of cuSPARSE before 12.8, the internal buffer
+   // used for cusparseSpMV must be the same when it is called with the same
+   // matrix.
+   //
+   // By default, MFEM uses one buffer, that is shared by all sparse matrices.
+   // In the code below, a buffer is created for A, then modified for B, then
+   // used again for A. Without the workaround, this fails with cuSPARSE version
+   // earlier than 12.8 (confirmed to fail with 12.4).
+
+   const int n = 100;
+   SparseMatrix A(n, n);
+   Vector d(n);
+   d.Randomize(1);
+   for (int i = 0; i < n; ++i)
+   {
+      A.Set(i, i, d[i]);
+   }
+   A.Finalize();
+
+   Vector x(n);
+   x = 1.0;
+
+   Vector y(n);
+   A.Mult(x, y);
+
+   {
+      SparseMatrix B(20, 20);
+      for (int i = 0; i < 20; ++i)
+      {
+         for (int j = 0; j < 20; ++j)
+         {
+            B.Set(i, j, 1.0);
+         }
+      }
+      B.Finalize();
+      Vector u(20);
+      u = 1.0;
+      Vector v(20);
+      B.Mult(u, v);
+   }
+
+   A.Mult(x, y);
+
+   y -= d;
+   REQUIRE(y.Normlinf() == MFEM_Approx(0.0));
+}
+
 } // namespace mfem