Merge pull request #292 from awslabs/simlapointe/transient-adapt-dt

Adaptive time-stepping for transient solver using SUNDIALS

Author: simlapointe

CHANGELOG.md

@@ -42,6 +42,9 @@ The format of this changelog is based on
   - Nonconformal adaptation is now supported for WavePort boundary conditions. This was
     achieved through a patch applied to MFEM to support `mfem::ParSubMesh` on external
     nonconformal surface subdomains.
+  - Added adaptive time-stepping capability for transient simulations. The new ODE integrators
+    rely on the SUNDIALS library and can be specified by setting the
+    `config["Solver"]["Transient"]["Type"]` option to `"CVODE"` or `"ARKODE"`.
 
 ## [0.13.0] - 2024-05-20
 

CMakeLists.txt

@@ -37,6 +37,8 @@ set(PALACE_WITH_GSLIB ON CACHE BOOL "Build with GSLIB library for high-order fie
 set(PALACE_WITH_STRUMPACK_BUTTERFLYPACK OFF CACHE BOOL "Build with ButterflyPACK support for STRUMPACK solver")
 set(PALACE_WITH_STRUMPACK_ZFP OFF CACHE BOOL "Build with ZFP support for STRUMPACK solver")
 
+set(PALACE_WITH_SUNDIALS ON CACHE BOOL "Build with SUNDIALS differential/algebraic equations solver")
+
 set(ANALYZE_SOURCES_CLANG_TIDY OFF CACHE BOOL "Run static analysis checks using clang-tidy")
 set(ANALYZE_SOURCES_CPPCHECK OFF CACHE BOOL "Run static analysis checks using cppcheck")
 

cmake/ExternalGitTags.cmake

@@ -265,3 +265,17 @@ set(EXTERN_EIGEN_URL
   "https://gitlab.com/libeigen/eigen/-/archive/3.4.0/eigen-3.4.0.tar.gz" CACHE STRING
   "URL for external Eigen build"
 )
+
+# SUNDIALS
+set(EXTERN_SUNDIALS_URL
+  "https://github.com/LLNL/sundials.git" CACHE STRING
+  "URL for external SUNDIALS build"
+)
+set(EXTERN_SUNDIALS_GIT_BRANCH
+  "main" CACHE STRING
+  "Git branch for external SUNDIALS build"
+)
+set(EXTERN_SUNDIALS_GIT_TAG
+  "aaeab8d907c6b7dfca86041401fdc1448f35f826" CACHE STRING
+  "Git tag for external SUNDIALS build"
+)

cmake/ExternalMFEM.cmake

@@ -17,6 +17,9 @@ if(PALACE_BUILD_EXTERNAL_DEPS)
   if(PALACE_WITH_SUPERLU)
     list(APPEND MFEM_DEPENDENCIES superlu_dist)
   endif()
+  if(PALACE_WITH_SUNDIALS)
+    list(APPEND MFEM_DEPENDENCIES sundials)
+  endif()
 else()
   set(MFEM_DEPENDENCIES)
 endif()
@@ -87,6 +90,7 @@ list(APPEND MFEM_OPTIONS
   "-DMFEM_USE_LIBUNWIND=${PALACE_MFEM_WITH_LIBUNWIND}"
   "-DMFEM_USE_METIS_5=YES"
   "-DMFEM_USE_CEED=NO"
+  "-DMFEM_USE_SUNDIALS=${PALACE_WITH_SUNDIALS}"
 )
 if(PALACE_WITH_STRUMPACK OR PALACE_WITH_MUMPS)
   list(APPEND MFEM_OPTIONS
@@ -297,6 +301,30 @@ Intel C++ compiler for MUMPS and STRUMPACK dependencies")
       "-DMUMPS_REQUIRED_LIBRARIES=${SCALAPACK_LIBRARIES}$<SEMICOLON>${STRUMPACK_MUMPS_GFORTRAN_LIBRARY}"
     )
   endif()
+
+  # Configure SUNDIALS
+  if(PALACE_WITH_SUNDIALS)
+    set(SUNDIALS_REQUIRED_PACKAGES "LAPACK" "BLAS" "MPI")
+    if(PALACE_WITH_OPENMP)
+      list(APPEND SUNDIALS_REQUIRED_PACKAGES "OpenMP")
+    endif()
+    if(PALACE_WITH_CUDA)
+      list(APPEND SUNDIALS_REQUIRED_PACKAGES "CUDAToolkit")
+      list(APPEND SUNDIALS_REQUIRED_LIBRARIES ${SUPERLU_STRUMPACK_CUDA_LIBRARIES})
+    endif()
+    string(REPLACE ";" "$<SEMICOLON>" SUNDIALS_REQUIRED_PACKAGES "${SUNDIALS_REQUIRED_PACKAGES}")
+    string(REPLACE ";" "$<SEMICOLON>" SUNDIALS_REQUIRED_LIBRARIES "${SUNDIALS_REQUIRED_LIBRARIES}")
+    list(APPEND MFEM_OPTIONS
+      "-DSUNDIALS_DIR=${CMAKE_INSTALL_PREFIX}"
+      "-DSUNDIALS_REQUIRED_PACKAGES=${SUNDIALS_REQUIRED_PACKAGES}"
+    )
+    if(NOT "${SUNDIALS_REQUIRED_LIBRARIES}" STREQUAL "")
+      list(APPEND MFEM_OPTIONS
+        "-DSUNDIALS_REQUIRED_LIBRARIES=${SUNDIALS_REQUIRED_LIBRARIES}"
+      )
+    endif()
+  endif()
+
 else()
   # Help find dependencies for the internal MFEM build
   # If we trust MFEM's Find<PACKAGE>.cmake module, we can just set <PACKAGE>_DIR and, if
@@ -309,6 +337,7 @@ else()
     "SuperLUDist"
     "STRUMPACK"
     "MUMPS"
+    "SUNDIALS"
   )
   foreach(DEP IN LISTS PALACE_MFEM_DEPS)
     set(${DEP}_DIR "" CACHE STRING "Path to ${DEP} build or installation directory")

cmake/ExternalSUNDIALS.cmake

@@ -0,0 +1,66 @@
+# Copyright Amazon.com, Inc. or its affiliates. All Rights Reserved.
+# SPDX-License-Identifier: Apache-2.0
+
+#
+# Build SUNDIALS
+#
+
+# Force build order
+set(SUNDIALS_DEPENDENCIES)
+
+set(SUNDIALS_OPTIONS ${PALACE_SUPERBUILD_DEFAULT_ARGS})
+list(APPEND SUNDIALS_OPTIONS
+  "-DCMAKE_CXX_COMPILER=${CMAKE_CXX_COMPILER}"
+  "-DCMAKE_CXX_FLAGS=${CMAKE_CXX_FLAGS}"
+  "-DCMAKE_C_COMPILER=${CMAKE_C_COMPILER}"
+  "-DCMAKE_C_FLAGS=${CMAKE_C_FLAGS}"
+  "-DEXAMPLES_ENABLE_C=OFF"
+  "-DEXAMPLES_ENABLE_CXX=OFF"
+  "-DEXAMPLES_ENABLE_CUDA=OFF"
+  "-DEXAMPLES_INSTALL=OFF"
+  "-DENABLE_MPI=ON"
+)
+
+if(PALACE_WITH_OPENMP)
+  list(APPEND SUNDIALS_OPTIONS
+    "-DENABLE_OPENMP=ON"
+  )
+else()
+  list(APPEND SUNDIALS_OPTIONS
+    "-DENABLE_OPENMP=OFF"
+  )
+endif()
+
+# Configure LAPACK dependency
+if(NOT "${BLAS_LAPACK_LIBRARIES}" STREQUAL "")
+  list(APPEND SUNDIALS_OPTIONS
+    "-DENABLE_LAPACK=ON"
+    "-DLAPACK_LIBRARIES=${BLAS_LAPACK_LIBRARIES}"
+    "-DLAPACK_WORKS:BOOL=TRUE"
+  )
+endif()
+
+if(PALACE_WITH_CUDA)
+  list(APPEND SUNDIALS_OPTIONS
+    "-DCMAKE_CUDA_ARCHITECTURES=${CMAKE_CUDA_ARCHITECTURES}"
+    "-DENABLE_CUDA=ON"
+  )
+endif()
+
+string(REPLACE ";" "; " SUNDIALS_OPTIONS_PRINT "${SUNDIALS_OPTIONS}")
+message(STATUS "SUNDIALS_OPTIONS: ${SUNDIALS_OPTIONS_PRINT}")
+
+
+include(ExternalProject)
+ExternalProject_Add(sundials
+  DEPENDS           ${SUNDIALS_DEPENDENCIES}
+  GIT_REPOSITORY    ${EXTERN_SUNDIALS_URL}
+  GIT_TAG           ${EXTERN_SUNDIALS_GIT_TAG}
+  SOURCE_DIR        ${CMAKE_BINARY_DIR}/extern/sundials
+  BINARY_DIR        ${CMAKE_BINARY_DIR}/extern/sundials-build
+  INSTALL_DIR       ${CMAKE_INSTALL_PREFIX}
+  PREFIX            ${CMAKE_BINARY_DIR}/extern/sundials-Cmake
+  UPDATE_COMMAND    ""
+  CONFIGURE_COMMAND ${CMAKE_COMMAND} <SOURCE_DIR> "${SUNDIALS_OPTIONS}"
+  TEST_COMMAND      ""
+)

docs/src/config/solver.md

@@ -208,7 +208,10 @@ error tolerance.
     "ExcitationWidth": <float>,
     "MaxTime": <float>,
     "TimeStep": <float>,
-    "SaveStep": <int>
+    "SaveStep": <int>,
+    "Order": <int>,
+    "RelTol": <float>,
+    "AbsTol": <float>
 }
 ```
 
@@ -218,14 +221,12 @@ with
 the second-order system of differential equations. The available options are:
 
   - `"GeneralizedAlpha"` :  The second-order implicit generalized-``\alpha`` method with
-    ``\rho_\inf = 1.0``. This scheme is unconditionally stable.
-  - `"NewmarkBeta"` :  The second-order implicit Newmark-``\beta`` method with
-    ``\beta = 1/4`` and ``\gamma = 1/2``. This scheme is unconditionally stable.
-  - `"CentralDifference"` :  The second-order explicit central difference method, obtained
-    by setting ``\beta = 0`` and ``\gamma = 1/2`` in the Newmark-``\beta`` method. In this
-    case, the maximum eigenvalue of the system operator is estimated at the start of the
-    simulation and used to restrict the simulation time step to below the maximum stability
-    time step.
+    ``\rho_{\inf} = 1.0``. This scheme is unconditionally stable.
+  - `"ARKODE"` :  SUNDIALS ARKode implicit Runge-Kutta scheme applied to the first-order
+    ODE system for the electric field with adaptive time-stepping. This option is only available when *Palace* has been [built with SUNDIALS support](../install.md#Configuration-options).
+  - `"CVODE"` :  SUNDIALS CVODE implicit multistep method scheme applied to the first-order
+    ODE system for the electric field with adaptive time-stepping. This option is only available when *Palace* has been [built with SUNDIALS support](../install.md#Configuration-options).
+  - `"RungeKutta"` : Two stage, singly diagonal implicit Runge-Kutta (SDIRK) method. Second order and L-stable.
   - `"Default"` :  Use the default `"GeneralizedAlpha"` time integration scheme.
 
 `"Excitation" [None]` :  Controls the time dependence of the source excitation. The
@@ -260,6 +261,16 @@ disk for [visualization with ParaView](../guide/postprocessing.md#Visualization)
 saved in the `paraview/` directory under the directory specified by
 [`config["Problem"]["Output"]`](problem.md#config%5B%22Problem%22%5D).
 
+`"Order" [2]` :  Order of the adaptive Runge-Kutta integrators or maximum order of the
+multistep method, must be within `[2,5]`. Should only be specified if `"Type"` is `"ARKODE"`
+or `"CVODE"`.
+
+`"RelTol" [1e-4]` :  Relative tolerance used in adaptive time-stepping schemes. Should only
+be specified if `"Type"` is `"ARKODE"` or `"CVODE"`.
+
+`"AbsTol" [1e-9]` :  Absolute tolerance used in adaptive time-stepping schemes. Should only
+be specified if `"Type"` is `"ARKODE"` or `"CVODE"`.
+
 ## `solver["Electrostatic"]`
 
 ```json

docs/src/install.md

@@ -133,6 +133,7 @@ Additional build options are (with default values in brackets):
   - `PALACE_WITH_LIBXSMM [ON]` :  Build with LIBXSMM backend for libCEED
   - `PALACE_WITH_MAGMA [ON]` :  Build with MAGMA backend for libCEED
   - `PALACE_WITH_GSLIB [ON]` :  Build with GSLIB library for high-order field interpolation
+  - `PALACE_WITH_SUNDIALS [ON]` : Build with SUNDIALS ODE solver library
 
 The build step is invoked by running (for example with 4 `make` threads)
 
@@ -207,6 +208,7 @@ source code for these dependencies is downloaded during the build process:
     `PALACE_WITH_ARPACK=ON`)
   - [LIBXSMM](https://github.com/libxsmm/libxsmm) (optional, when `PALACE_WITH_LIBXSMM=ON`)
   - [MAGMA](https://icl.utk.edu/magma/) (optional, when `PALACE_WITH_MAGMA=ON`)
+  - [SUNDIALS](https://github.com/LLNL/sundials) (optional, when `PALACE_WITH_SUNDIALS=ON`)
   - [nlohmann/json](https://github.com/nlohmann/json)
   - [fmt](https://fmt.dev/latest)
   - [Eigen](https://eigen.tuxfamily.org)

docs/src/reference.md

@@ -94,9 +94,19 @@ boundary condition which is written for a lumped resistive port boundary, for ex
     = \bm{U}^{inc} \,,\; \bm{x}\in\Gamma_{Z} \,.
 ```
 
-The second-order electric field formulation is chosen to take advantage of unconditionally
-stable implicit time-integration schemes without the expense of a coupled block system
-solve for ``\bm{E}(\bm{x},t)`` and ``\bm{B}(\bm{x},t)``. It offers the additional benefit
+The second-order electric field differential equation is transformed into a first-order
+ODE system which is solved along with the equation for the magnetic flux density
+
+```math
+\left(\begin{matrix} \varepsilon_r & 0 & 0 \\ 0 & I & 0 \\ 0 & 0 & I\end{matrix}\right)
+  \left(\begin{matrix} \ddot{\bm{E}} \\ \dot{\bm{E}} \\ \dot{\bm{B}}\end{matrix} \right)
+  = \left(\begin{matrix} -\sigma & -\nabla\times\mu_r^{-1}\nabla\times & 0\\ I & 0 & 0 \\ 0 & -\nabla\times & 0\end{matrix}\right)
+    \left(\begin{matrix}\dot{\bm{E}}\\ \bm{E} \\ \bm{B} \end{matrix}\right) \,.
+```
+
+The first-order ODE system formulation is chosen to take advantage of implicit adaptive
+time-stepping integration schemes. The ``3 \times 3`` system can be block-eliminated to
+avoid an expensive coupled block system solve. It offers the additional benefit
 of sharing many similarities in the spatial discretization as the frequency domain
 formulation outlined above.
 

extern/CMakeLists.txt

@@ -72,3 +72,9 @@ include(ExternalFmt)
 # Add Eigen
 message(STATUS "===================== Configuring Eigen dependency =====================")
 include(ExternalEigen)
+
+# Add SUNDIALS
+if(PALACE_WITH_SUNDIALS)
+  message(STATUS "===================== Configuring SUNDIALS dependency =====================")
+  include(ExternalSUNDIALS)
+endif()

palace/drivers/transientsolver.cpp

@@ -30,15 +30,8 @@ TransientSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
   std::function<double(double)> dJdt_coef = GetTimeExcitation(true);
   SpaceOperator space_op(iodata, mesh);
   TimeOperator time_op(iodata, space_op, dJdt_coef);
+
   double delta_t = iodata.solver.transient.delta_t;
-  if (time_op.isExplicit())
-  {
-    // Stability limited time step.
-    const double dt_max = time_op.GetMaxTimeStep();
-    const double dts_max = iodata.DimensionalizeValue(IoData::ValueType::TIME, dt_max);
-    Mpi::Print(" Maximum stable time step: {:.6e} ns\n", dts_max);
-    delta_t = std::min(delta_t, 0.95 * dt_max);
-  }
   int n_step = GetNumSteps(0.0, iodata.solver.transient.max_t, delta_t);
   SaveMetadata(space_op.GetNDSpaces());
 
@@ -138,6 +131,7 @@ TransientSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
     step++;
   }
   BlockTimer bt1(Timer::POSTPRO);
+  time_op.PrintStats();
   SaveMetadata(time_op.GetLinearSolver());
   return {indicator, space_op.GlobalTrueVSize()};
 }