Merge branch 'master' into main-dev

Remove hypre+umpire check in CMake

.codecov.yml

@@ -12,7 +12,7 @@ coverage:
         threshold: 0%
         base: auto
         branches:
-          - master
+          - main
         if_ci_failed: error
         informational: true
         only_pulls: true
@@ -22,7 +22,7 @@ coverage:
         threshold: 1%  # allows variations around the target
         base: auto
         branches:
-          - master
+          - main
         if_ci_failed: error
         only_pulls: true
 

.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

@@ -40,6 +40,8 @@ Discretization improvements
   provided that neighboring hexahedra are not refined in conflicting directions.
   A new ParMesh method is added to check for such conflicts, before refinement.
 
+- Renamed the default GitHub branch from "master" to "main".
+
 Meshing improvements
 --------------------
 
@@ -91,6 +93,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

@@ -652,7 +652,7 @@ if (MFEM_USE_CUDA)
 endif()
 
 add_subdirectory(config)
-set(MASTER_HEADERS
+set(MAIN_HEADERS
   ${PROJECT_SOURCE_DIR}/mfem.hpp
   ${PROJECT_SOURCE_DIR}/mfem-performance.hpp)
 
@@ -689,7 +689,7 @@ set(MFEM_SOURCE_DIR ${CMAKE_CURRENT_SOURCE_DIR})
 set(MFEM_INSTALL_DIR ${CMAKE_INSTALL_PREFIX})
 
 # Declaring the library
-mfem_add_library(mfem ${SOURCES} ${HEADERS} ${MASTER_HEADERS})
+mfem_add_library(mfem ${SOURCES} ${HEADERS} ${MAIN_HEADERS})
 # message(STATUS "TPL_LIBRARIES = ${TPL_LIBRARIES}")
 target_link_libraries(mfem PUBLIC ${TPL_LIBRARIES} ${TPL_TARGETS})
 if (TPL_TARGETS)
@@ -888,12 +888,12 @@ install(TARGETS ${PROJECT_NAME}
   LIBRARY DESTINATION ${INSTALL_LIB_DIR}
   ARCHIVE DESTINATION ${INSTALL_LIB_DIR})
 
-# Install the master headers
+# Install the main headers
 foreach(Header mfem.hpp mfem-performance.hpp)
   install(FILES ${PROJECT_BINARY_DIR}/InstallHeaders/${Header}
     DESTINATION ${INSTALL_INCLUDE_DIR})
 endforeach()
-install(FILES ${MASTER_HEADERS} DESTINATION ${INSTALL_INCLUDE_DIR}/mfem)
+install(FILES ${MAIN_HEADERS} DESTINATION ${INSTALL_INCLUDE_DIR}/mfem)
 
 # Install the headers (except common miniapp which is installed from its subdir)
 install(DIRECTORY ${MFEM_SOURCE_DIRS}

CONTRIBUTING.md

@@ -3,10 +3,10 @@
 </p>
 
 <p align="center">
-<a href="https://github.com/mfem/mfem/blob/master/LICENSE"><img alt="License" src="https://img.shields.io/badge/License-BSD-brightgreen.svg"></a>
-<a href="https://github.com/mfem/mfem/actions?query=workflow%3Arepo-check+branch%3Amaster"><img alt="Repo check" src="https://github.com/mfem/mfem/actions/workflows/repo-check.yml/badge.svg?branch=master"></a>
-<a href="https://github.com/mfem/mfem/actions?query=workflow%3Abuild-analysis+branch%3Amaster"><img alt="Build Analysis" src="https://github.com/mfem/mfem/actions/workflows/mfem-analysis.yml/badge.svg?branch=master"></a>
-<a href="https://github.com/mfem/mfem/actions?query=workflow%3Abuilds-and-tests+branch%3Amaster"><img alt="Builds and Tests" src="https://github.com/mfem/mfem/actions/workflows/builds-and-tests.yml/badge.svg?branch=master"></a>
+<a href="https://github.com/mfem/mfem/blob/main/LICENSE"><img alt="License" src="https://img.shields.io/badge/License-BSD-brightgreen.svg"></a>
+<a href="https://github.com/mfem/mfem/actions?query=workflow%3Arepo-check+branch%3Amain"><img alt="Repo check" src="https://github.com/mfem/mfem/actions/workflows/repo-check.yml/badge.svg?branch=main"></a>
+<a href="https://github.com/mfem/mfem/actions?query=workflow%3Abuild-analysis+branch%3Amain"><img alt="Build Analysis" src="https://github.com/mfem/mfem/actions/workflows/mfem-analysis.yml/badge.svg?branch=main"></a>
+<a href="https://github.com/mfem/mfem/actions?query=workflow%3Abuilds-and-tests+branch%3Amain"><img alt="Builds and Tests" src="https://github.com/mfem/mfem/actions/workflows/builds-and-tests.yml/badge.svg?branch=main"></a>
 <a href="https://ci.appveyor.com/project/mfem/mfem"><img alt="Build Status" src="https://ci.appveyor.com/api/projects/status/19non9sqm6msi2wy?svg=true"></a>
 <a href="https://docs.mfem.org/html/index.html"><img alt="Doxygen" src="https://img.shields.io/badge/code-documented-brightgreen.svg"></a>
 </p>
@@ -27,7 +27,7 @@ in the MFEM community, you agree to abide by its rules.
 If you plan on contributing to MFEM, consider reviewing the
 [issue tracker](https://github.com/mfem/mfem/issues) first to check if a thread
 already exists for your desired feature or the bug you ran into. Use a pull
-request (PR) toward the `mfem:master` branch to propose your contribution. If
+request (PR) toward the `mfem:main` branch to propose your contribution. If
 you are planning significant code changes or have questions, you may want to
 open an [issue](https://github.com/mfem/mfem/issues) before issuing a PR. In
 addition to technical contributions, we are also interested in your results and
@@ -47,7 +47,7 @@ back to them before issuing pull requests:
   - [Pull Requests](#pull-requests)
   - [MFEM PR Rules](#mfem-pr-rules)
   - [Pull Request Checklist](#pull-request-checklist)
-  - [Master/Next Workflow](#masternext-workflow)
+  - [Main/Next Workflow](#mainnext-workflow)
   - [Releases](#releases)
   - [Release Checklist](#release-checklist)
 - [LLNL Workflow](#llnl-workflow)
@@ -66,12 +66,12 @@ Origin](#developers-certificate-of-origin-11) at the end of this file.*
 ## Quick Summary
 
 - We encourage you to [join the MFEM organization](#mfem-organization) and create
-  development branches off `mfem:master`.
+  development branches off `mfem:main`.
 - Please follow the [developer guidelines](#developer-guidelines), in particular
   with regards to documentation and code styling.
 - Please do not commit large/binary files to the central repository (use a fork
   instead).
-- Pull requests should be issued toward `mfem:master`. Make sure
+- Pull requests should be issued toward `mfem:main`. Make sure
   to check the items off the [Pull Request Checklist](#pull-request-checklist) and
   follow the [MFEM PR Rules](#mfem-pr-rules).
 - When your contribution is fully working and ready to be reviewed, add
@@ -81,8 +81,8 @@ Origin](#developers-certificate-of-origin-11) at the end of this file.*
 - The reviewers have 3 weeks to evaluate the PR and work with the author to
   fix issues and implement improvements.
 - During review there should be no force pushes/rewriting history in the branch.
-- After approval, MFEM developers merge the PR manually in the [mfem:next branch](#masternext-workflow).
-- After a week of testing in `mfem:next`, the original PR is merged in `mfem:master`.
+- After approval, MFEM developers merge the PR manually in the [mfem:next branch](#mainnext-workflow).
+- After a week of testing in `mfem:next`, the original PR is merged in `mfem:main`.
 - We use [milestones](https://github.com/mfem/mfem/milestones) to coordinate the
   work on different PRs toward a release.
 - Don't hesitate to [contact us](#contact-information) if you have any questions.
@@ -289,7 +289,7 @@ Before you can start, you need a GitHub account, here are a few suggestions:
   author, is willing to work on it and be its champion.
 
 - The author creates a branch for the new feature (with suffix `-dev`), off
-  the `master` branch, or another existing feature branch, for example:
+  the `main` branch, or another existing feature branch, for example:
 
   ```
   # Clone assuming you have setup your ssh keys on GitHub:
@@ -298,8 +298,8 @@ Before you can start, you need a GitHub account, here are a few suggestions:
   # Alternatively, clone using the "https" protocol:
   git clone https://github.com/mfem/mfem.git
 
-  # Create a new feature branch starting from "master":
-  git checkout master
+  # Create a new feature branch starting from "main":
+  git checkout main
   git pull
   git checkout -b feature-dev
 
@@ -375,7 +375,7 @@ Before you can start, you need a GitHub account, here are a few suggestions:
 ### Pull Requests
 
 - When your branch is ready for other developers to review / comment on
-  the code, create a pull request towards `mfem:master`.
+  the code, create a pull request towards `mfem:main`.
 
 - Pull request typically have titles like:
 
@@ -411,7 +411,7 @@ Before you can start, you need a GitHub account, here are a few suggestions:
   assigned, the PR is considered under review. To help with the review process
   there should be no force pushes/rewriting history in the branch.
 
-- After approval, the PR is [tested](#masternext-workflow) for a week with
+- After approval, the PR is [tested](#mainnext-workflow) for a week with
   other approved PRs in the `mfem:next` branch.
 
 - Consider manually running the tests in `tests/scripts` before merging in
@@ -447,7 +447,7 @@ The Pull Request (PR) approval process in MFEM is similar to the approval of pap
 
 3. A PR can be (manually) merged in the *next* branch only if 2 of the assigned reviewers have approved it and it has passed internal testing. This merge can be performed by any of the assigned reviewers or by any of the editors.
 
-4. A PR can be merged in the *master* branch only if it has been tested successfully for a week in *next* and an editor has (optionally) taken a final look. This merge can be performed only by one of the editors.
+4. A PR can be merged in the *main* branch only if it has been tested successfully for a week in *next* and an editor has (optionally) taken a final look. This merge can be performed only by one of the editors.
 
 #### Responsibilities of Editors
 
@@ -470,7 +470,7 @@ The current list of MFEM editors is:
 
 5. To remind the reviewers about timely completion of their review.
 
-6. To take a final look and complete the PR merge in *master*. The final look step is optional and shouldn't take more than 3 days.
+6. To take a final look and complete the PR merge in *main*. The final look step is optional and shouldn't take more than 3 days.
 
 7. The assignment of bugfixes should be expedited proportional to their importance, e.g. in some cases the editor can assign much shorter review window.
 
@@ -492,7 +492,7 @@ Everyone on the MFEM team can be asked to serve as a reviewer on a PR in their a
 
 5. To test the PR thoroughly before merging in *next*. The PR author is also encouraged to perform testing and inform the reviewers about the results.
 
-6. To monitor the PR impact on the testing in the *next* branch and alert the editors that the PR is ready for merging in *master*.
+6. To monitor the PR impact on the testing in the *next* branch and alert the editors that the PR is ready for merging in *main*.
 
 7. The review of bugfixes should be expedited proportional to their importance. The review window can be much less than three weeks in such cases.
 
@@ -582,11 +582,11 @@ Before a PR can be merged, it should satisfy the following:
    - [ ] Update internal tests to include the new features.
 
 
-### Master/Next Workflow
+### Main/Next Workflow
 
-MFEM uses a `master`/`next`-branch workflow as described below:
+MFEM uses a `main`/`next`-branch workflow as described below:
 
-- The `master` branch should always be of release quality and changes should not
+- The `main` branch should always be of release quality and changes should not
   be merged until they have been fully tested. This branch is protected, and
   changes can only be made through pull requests.
 
@@ -613,20 +613,20 @@ MFEM uses a `master`/`next`-branch workflow as described below:
 
 - After a week of testing in `next` (excluding bugfixes), both on GitHub, as
   well as [internally](#tests-at-llnl) at LLNL, the original PR is merged into
-  `master` (provided there are no issues).
+  `main` (provided there are no issues).
 
 - After the merge, the feature branch is deleted (unless it is a long-term
   project with periodic PRs).
 
 - The `next` branch is used just for integrated testing of all PRs approved for
-  merging into `master` to verify that each works individually and that all of
+  merging into `main` to verify that each works individually and that all of
   them work as a group. This branch can be discarded at any time, though we
   typically do that only at the end of a [release cycle](#releases).
 
 
 ### Releases
 
-- Releases are just tags in the `master` branch, e.g. https://github.com/mfem/mfem/releases/tag/v3.3.2,
+- Releases are just tags in the `main` branch, e.g. https://github.com/mfem/mfem/releases/tag/v3.3.2,
   and have a version that ends in an even "patch" number, e.g. `v3.2.2` or
   `v3.4` (by convention `v3.4` is the same as `v3.4.0`.)  Between releases, the
   version ends in an odd "patch" number, e.g. `v3.3.3`.
@@ -679,19 +679,19 @@ MFEM uses a `master`/`next`-branch workflow as described below:
 
 ### Mirroring on Bitbucket
 
-- The GitHub `master` and `next` branches are mirrored to the LLNL institutional
-  Bitbucket repository as `gh-master` and `gh-next`.
+- The GitHub `main` and `next` branches are mirrored to the LLNL institutional
+  Bitbucket repository as `gh-main` and `gh-next`.
 
-- `gh-master` is merged into LLNL's internal `master` through pull requests; write
-  permissions to `master` are restricted to ensure this is the only way in which it
+- `gh-main` is merged into LLNL's internal `main` through pull requests; write
+  permissions to `main` are restricted to ensure this is the only way in which it
   gets updated.
 
 - We never push directly from LLNL to GitHub.
 
 - Versions of the code on LLNL's internal server, from most to least stable:
   - MFEM official release on mfem.org -- Most stable, tested in many apps.
-  - `mfem:master` -- Recent development version, guaranteed to work.
-  - `mfem:gh-master` -- Stable development version, passed testing, you can use
+  - `mfem:main` -- Recent development version, guaranteed to work.
+  - `mfem:gh-main` -- Stable development version, passed testing, you can use
      it to build your code between releases.
   - `mfem:gh-next` -- Bleeding-edge development version, may be broken, use at
      your own risk.
@@ -721,14 +721,14 @@ directory.
 
 
 ### Linux and Mac smoke tests
-We use GitHub Actions to drive the default tests on the `master` and `next`
+We use GitHub Actions to drive the default tests on the `main` and `next`
 branches. See the `.github/workflows` files and the logs at
 [https://github.com/mfem/mfem/actions](https://github.com/mfem/mfem/actions).
 
 Testing using GitHub Actions should be kept lightweight, as there is a time
 constraint on jobs. Two virtual machines are configured - Mac (OS X) and Linux.
 
-- Tests on the `master` branch are triggered whenever a PR is issued on this branch.
+- Tests on the `main` branch are triggered whenever a PR is issued on this branch.
 - Tests on the `next` branch are currently scheduled to run each night.
 
 
@@ -744,7 +744,7 @@ and debug build is performed with a simple run of `ex1` to verify the executable
 
 ### Tests at LLNL
 
-- We mirror the `master` and `next` branches internally (to `gh-master` and
+- We mirror the `main` and `next` branches internally (to `gh-main` and
   `gh-next`) and run longer nightly tests via cron. On the weekends, a more
   extensive test is run which extracts and executes all the different sample
   runs from each example and most miniapps.

config/cmake/modules/FindHYPRE.cmake

@@ -197,16 +197,6 @@ if (HYPRE_FOUND AND HYPRE_USING_HIP)
   message(STATUS "Updated HYPRE_LIBRARIES: ${HYPRE_LIBRARIES}")
 endif()
 
-# Hypre+Umpire check
-if (HYPRE_FOUND AND (HYPRE_USING_CUDA OR HYPRE_USING_HIP) AND NOT MFEM_USE_UMPIRE)
-  message(WARNING
-"===============================================================
- Detected GPU-enabled HYPRE build without Umpire support.
- This is not recommended for performance reasons!
- Consider rebuilding HYPRE with Umpire support.
-===============================================================")
-endif()
-
 find_package_handle_standard_args(HYPRE
   REQUIRED_VARS HYPRE_LIBRARIES HYPRE_INCLUDE_DIRS HYPRE_VERSION
 )

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=

config/docker/Dockerfile.base

@@ -33,8 +33,8 @@ RUN mkdir -p /opt/mfem-env \
 RUN cd /opt/mfem-env && \
     . /opt/spack/share/spack/setup-env.sh && \
     spack env activate . && \
-    spack develop --path /code mfem@master+examples+miniapps && \
-    spack add mfem@master+examples+miniapps && \
+    spack develop --path /code mfem@main+examples+miniapps && \
+    spack add mfem@main+examples+miniapps # && \
     spack install
 
 # ensure mfem always on various paths

config/docker/spack.yaml

@@ -1,11 +1,11 @@
 spack:
-  specs: [mfem@master+examples+miniapps]
-  view: 
+  specs: [mfem@main+examples+miniapps]
+  view:
     mfem:
-      root: /opt/mfem-view 
+      root: /opt/mfem-view
       link_type: copy
   concretization: together
   develop:
     mfem:
       path: /code
-      spec: mfem@master+examples+miniapps
+      spec: mfem@main+examples+miniapps

config/githooks/pre-push

@@ -115,7 +115,7 @@ fi
 
 # branch-history
 if [[ "${option}" == "--history" || "${option}" == "" ]]; then
-  git fetch origin master:master
+  git fetch origin main:main
   cd tests/scripts
   if ! ./runtest branch-history; then code=1; fi
   cd -

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

examples/jupyter/README.md

@@ -1,6 +1,6 @@
 # Jupyter Notebooks using xeus-cling
 
-[![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/mfem/mfem/master?filepath=examples%2Fjupyter%2Fex.ipynb)
+[![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/mfem/mfem/main?filepath=examples%2Fjupyter%2Fex.ipynb)
 
 [xeus-cling](https://github.com/jupyter-xeus/xeus-cling) is a C++ Jupyter Kernel based on [cling](https://github.com/root-project/cling),
 which can be used to create interactive C++ MFEM and GLVis notebooks.

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;

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

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

@@ -21,7 +21,7 @@ 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
@@ -51,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)
@@ -81,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')
@@ -89,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
 
@@ -112,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/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.

tests/scripts/README

@@ -89,7 +89,7 @@ The currently available test scripts are:
 - gitignore: Checks for compliance with MFEM's .gitignore rules
 
 - branch-history: Checks if the current branch history (the commits that will be
-  merged in master) contains unusually large files, unusually large number of
+  merged in main) contains unusually large files, unusually large number of
   changes in a commit, unusually large number of commits, etc.
 
 Adding new test scripts is easy and we welcome contributions from users and

tests/scripts/branch-history

@@ -27,7 +27,7 @@ sub usage {
           -h|--help   prints this usage information and exits
 
    This script checks if the current branch history, defined as the commits that
-   will be merged in master (those shown in a GitHub PR), contains unusually
+   will be merged in main (those shown in a GitHub PR), contains unusually
    large files, unusually large number of changes in a commit, unusually large
    number of commits, etc.
 
@@ -72,7 +72,7 @@ my $max_branch_kb = 1000;
 my $status = 0; # Return code
 
 # Get SHA hash of all commits in this branch
-my @commits = split /\n/, `git log --pretty=format:%H master..$branch`;
+my @commits = split /\n/, `git log --pretty=format:%H main..$branch`;
 
 # Check if total number of commits in this branch exceeds the maximum allowable
 my $ncommits = scalar @commits;

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