Merge pull request #47 from awslabs/sjg/driven-port-debug

Refactor linear algebra functionality for complex-valued problems

.github/workflows/build-and-test-macos.yml

@@ -16,14 +16,14 @@ jobs:
       fail-fast: false
       matrix:
         include:  # Pairwise testing
-          # - compiler: clang
-          #   mpi: openmpi
-          #   math-libs: openblas
-          #   build-shared: shared
-          #   with-64bit-int: int32
-          #   with-openmp: serial
-          #   with-solver: superlu
-          #   with-eigensolver: slepc
+          - compiler: clang
+            mpi: openmpi
+            math-libs: openblas
+            build-shared: shared
+            with-64bit-int: int32
+            with-openmp: serial
+            with-solver: superlu
+            with-eigensolver: slepc
 
           - compiler: gcc
             mpi: openmpi

CHANGELOG.md

@@ -13,16 +13,27 @@ The format of this changelog is based on
 
 ## In progress
 
+  - Changed implementation of complex-valued linear algebra to use new `ComplexVector` and
+    `ComplexOperator` types, which are based on the underlying `mfem::Vector` and
+    `mfem::Operator` classes, instead of PETSc. PETSc is now fully optional and only
+    required when SLEPc eigenvalue solver support is requested. Krylov solvers for real- and
+    complex-valued linear systems are implemented via the built-in `IterativeSolver`
+    classes.
+  - Changed implementation of PROMs for adaptive fast frequency sweep to use the Eigen
+    library for sequential dense linear algebra.
+  - Changed implementation of numeric wave ports to use MFEM's `SubMesh` functionality. As
+    of [#3379](https://github.com/mfem/mfem/pull/3379) in MFEM, this has full ND and RT
+    basis support. For now, support for nonconforming mesh boundaries is limited.
   - Added Apptainer/Singularity container build definition for Palace.
   - Added build dependencies on [libCEED](https://github.com/CEED/libCEED) and
     [LIBXSMM](https://github.com/libxsmm/libxsmm) to support operator partial assembly (CPU-
     based for now).
 
 ## [0.11.2] - 2023-07-14
 
-  - Changed layout and names of `palace/` source directory for better organization.
   - Fixed a regression bug affecting meshes which have domain elements which are not
     assigned material properties in the configuration file.
+  - Changed layout and names of `palace/` source directory for better organization.
   - Added many updates to build system: Removed use of Git submodules to download
     dependencies relying instead directly on CMake's ExternalProject, patch GSLIB dependency
     for shared library builds, add CI tests with ARPACK-NG instead of SLEPc, update all

CMakeLists.txt

@@ -57,12 +57,6 @@ if(NOT DEFINED BUILD_SHARED_LIBS)
   set(BUILD_SHARED_LIBS OFF CACHE BOOL "Global flag to cause add_library() to create shared libraries if ON")
 endif()
 
-# For now, SLEPc is always required
-if(NOT PALACE_WITH_SLEPC)
-  message(STATUS "Building with SLEPc eigenvalue solver as it is required")
-  set(PALACE_WITH_SLEPC ON CACHE BOOL "Build with SLEPc eigenvalue solver" FORCE)
-endif()
-
 # Add extra CMake modules
 list(APPEND CMAKE_MODULE_PATH "${CMAKE_SOURCE_DIR}/cmake")
 

docs/src/config/solver.md

@@ -299,11 +299,16 @@ directory specified by [`config["Problem"]["Output"]`]
     "Tol": <float>,
     "MaxIts": <int>,
     "MaxSize": <int>,
-    "UseGMG": <bool>,
-    "UsePCShifted": <bool>,
+    "UsePCMatShifted": <bool>,
+    "PCSide": <string>,
+    "UseMultigrid": <bool>,
+    "MGAuxiliarySmoother": <bool>,
     "MGCycleIts": <int>,
     "MGSmoothIts": <int>,
-    "MGSmoothOrder": <int>
+    "MGSmoothOrder": <int>,
+    "DivFreeTol": <float>,
+    "DivFreeMaxIts": <float>,
+    "GSOrthogonalization": <string>
 }
 ```
 
@@ -353,46 +358,66 @@ equations arising for each simulation type. The available options are:
     definite (SPD) and the preconditioned conjugate gradient method (`"CG"`) is used as the
     Krylov solver.
 
-`"Tol" [1.0e-6]` :  Relative (preconditioned) residual convergence tolerance for the
-iterative linear solver.
+`"Tol" [1.0e-6]` :  Relative residual convergence tolerance for the iterative linear solver.
 
 `"MaxIts" [100]` :  Maximum number of iterations for the iterative linear solver.
 
 `"MaxSize" [0]` :  Maximum Krylov space size for the GMRES and FGMRES solvers. A value less
 than 1 defaults to the value specified by `"MaxIts"`.
 
-`"UseGMG" [true]` :  Enable or not [geometric multigrid solver]
-(https://en.wikipedia.org/wiki/Multigrid_method) which uses h- and p-multigrid coarsening as
-available to construct the multigrid hierarchy. The solver specified by `"Type"` is used on
-the coarsest level. A Hiptmair smoother is applied to all other levels.
-
-`"UsePCShifted" [false]` :  When set to `true`, constructs the preconditioner for frequency
+`"UsePCMatShifted" [false]` :  When set to `true`, constructs the preconditioner for frequency
 domain problems using a real SPD approximation of the system matrix, which can help
 performance at high frequencies (relative to the lowest nonzero eigenfrequencies of the
 model).
 
+`"PCSide" ["Default"]` :  Side for preconditioning. Not all options are available for all
+iterative solver choices, and the default choice depends on the iterative solver used.
+
+  - `"Left"`
+  - `"Right"`
+  - `"Default"`
+
+`"UseMultigrid" [true]` :  Chose whether to enable [geometric multigrid preconditioning]
+(https://en.wikipedia.org/wiki/Multigrid_method) which uses p- and h-multigrid coarsening as
+available to construct the multigrid hierarchy. The solver specified by `"Type"` is used on
+the coarsest level. Relaxation on the fine levels is performed with Chebyshev smoothing.
+
+`"MGAuxiliarySmoother"` :  Activate hybrid smoothing from Hiptmair for multigrid levels when
+`"UseMultigrid"` is `true`. For non-singular problems involving curl-curl operators, this
+option is `true` by default.
+
 `"MGCycleIts" [1]` :  Number of V-cycle iterations per preconditioner application for
-multigrid preconditioners (when `"UseGMG"` is `true` or `"Type"` is `"AMS"` or
+multigrid preconditioners (when `"UseMultigrid"` is `true` or `"Type"` is `"AMS"` or
 `"BoomerAMG"`).
 
 `"MGSmoothIts" [1]` :  Number of pre- and post-smooth iterations used for multigrid
-preconditioners (when `"UseGMG"` is `true` or `"Type"` is `"AMS"` or `"BoomerAMG"`).
+preconditioners (when `"UseMultigrid"` is `true` or `"Type"` is `"AMS"` or `"BoomerAMG"`).
 
 `"MGSmoothOrder" [3]` :  Order of polynomial smoothing for geometric multigrid
-preconditioning (when `"UseGMG"` is `true`).
+preconditioning (when `"UseMultigrid"` is `true`).
+
+`"DivFreeTol" [1.0e-12]` :  Relative tolerance for divergence-free cleaning used in the
+eigenmode simulation type.
+
+`"DivFreeMaxIts" [100]` :  Maximum number of iterations for divergence-free cleaning use in
+the eigenmode simulation type.
+
+`"GSOrthogonalization" ["MGS"]` :  Gram-Schmidt variant used to explicitly orthogonalize
+vectors in Krylov subspace methods or other parts of the code.
+
+  - `"MGS"` :  Modified Gram-Schmidt
+  - `"CGS"` :  Classical Gram-Schmidt
+  - `"CGS2"` :  Two-step classical Gram-Schmidt with reorthogonalization
 
 ### Advanced linear solver options
 
-  - `"Type"`: `"STRUMPACK-MP"`
-  - `"KSPType"`: `"MINRES"`, `"CGSYM"`, `"FCG"`, `"BCGS"`, `"BCGSL"`, `"FBCGS"`, `"QMRCGS"`,
-    `"TFQMR"`
-  - `"UseMGS" [false]`
-  - `"UseCGS2" [false]`
-  - `"UseKSPPiped" [false]`
-  - `"UseLOR" [false]`
-  - `"PrecondSide" ["Default"]`: `"Left"`, `"Right"`, `"Default"`
-  - `"Reordering" ["Default"]`: `"METIS"`, `"ParMETIS"`, `"Default"`
-  - `"STRUMPACKCompressionType" ["None"]`: `"None"`, `"BLR"`, `"HSS"`, `"HODLR"`
+  - `"UseInitialGuess" [true]`
+  - `"UsePartialAssembly" [false]`
+  - `"UseLowOrderRefined" [false]`
+  - `"Reordering" ["Default"]` :  `"METIS"`, `"ParMETIS"`,`"Scotch"`, `"PTScotch"`,
+    `"Default"`
+  - `"STRUMPACKCompressionType" ["None"]` :  `"None"`, `"BLR"`, `"HSS"`, `"HODLR"`, `"ZFP"`,
+    `"BLR-HODLR"`, `"ZFP-BLR-HODLR"`
   - `"STRUMPACKCompressionTol" [1.0e-3]`
   - `"STRUMPACKLossyPrecision" [16]`
   - `"STRUMPACKButterflyLevels" [1]`

docs/src/guide/model.md

@@ -40,7 +40,7 @@ or region-based refinement, specified using the [`config["Model"]["Refinement"]`
 uniform refinement levels as well as local refinement regions which refines the elements
 inside of a certain box or sphere-shaped region. For simplex meshes, the refinement
 maintains a conforming mesh but meshes containing hexahedra, prism, or pyramid elements
-will be non-conforming after local refinement (this is not supported at this time).
+will be nonconforming after local refinement (this is not supported at this time).
 
 [Adaptive mesh refinement (AMR)](https://en.wikipedia.org/wiki/Adaptive_mesh_refinement)
 according to error estimates in the computed solution is a work in progress for all

palace/drivers/basesolver.cpp

@@ -6,6 +6,7 @@
 #include <complex>
 #include <mfem.hpp>
 #include <nlohmann/json.hpp>
+#include "linalg/ksp.hpp"
 #include "models/domainpostoperator.hpp"
 #include "models/postoperator.hpp"
 #include "models/surfacepostoperator.hpp"
@@ -51,10 +52,9 @@ void WriteMetadata(const std::string &post_dir, const json &meta)
 
 }  // namespace
 
-BaseSolver::BaseSolver(const IoData &iodata_, bool root_, int size, int num_thread,
+BaseSolver::BaseSolver(const IoData &iodata, bool root, int size, int num_thread,
                        const char *git_tag)
-  : iodata(iodata_), post_dir(GetPostDir(iodata_.problem.output)), root(root_),
-    table(8, 9, 6)
+  : iodata(iodata), post_dir(GetPostDir(iodata.problem.output)), root(root), table(8, 9, 6)
 {
   // Create directory for output.
   if (root && !std::filesystem::exists(post_dir))
@@ -100,7 +100,8 @@ void BaseSolver::SaveMetadata(const mfem::ParFiniteElementSpace &fespace) const
   }
 }
 
-void BaseSolver::SaveMetadata(int ksp_mult, int ksp_it) const
+template <typename SolverType>
+void BaseSolver::SaveMetadata(const SolverType &ksp) const
 {
   if (post_dir.length() == 0)
   {
@@ -109,8 +110,8 @@ void BaseSolver::SaveMetadata(int ksp_mult, int ksp_it) const
   if (root)
   {
     json meta = LoadMetadata(post_dir);
-    meta["LinearSolver"]["TotalSolves"] = ksp_mult;
-    meta["LinearSolver"]["TotalIts"] = ksp_it;
+    meta["LinearSolver"]["TotalSolves"] = ksp.NumTotalMult();
+    meta["LinearSolver"]["TotalIts"] = ksp.NumTotalMultIterations();
     WriteMetadata(post_dir, meta);
   }
 }
@@ -556,4 +557,7 @@ void BaseSolver::PostprocessFields(const PostOperator &postop, int step, double
   Mpi::Barrier();
 }
 
+template void BaseSolver::SaveMetadata<KspSolver>(const KspSolver &) const;
+template void BaseSolver::SaveMetadata<ComplexKspSolver>(const ComplexKspSolver &) const;
+
 }  // namespace palace

palace/drivers/basesolver.hpp

@@ -45,9 +45,7 @@ class BaseSolver
     int p;   // Floating point precision for data
     int w1;  // First column width = precision + 7 extra
     int p1;  // Floating point precision for first column
-    Table(int sp_, int p_, int p1_) : w(sp_ + p_ + 7), sp(sp_), p(p_), w1(p1_ + 7), p1(p1_)
-    {
-    }
+    Table(int sp, int p, int p1) : w(sp + p + 7), sp(sp), p(p), w1(p1 + 7), p1(p1) {}
   };
   const Table table;
 
@@ -71,7 +69,7 @@ class BaseSolver
   void PostprocessFields(const PostOperator &postop, int step, double time) const;
 
 public:
-  BaseSolver(const IoData &iodata_, bool root_, int size = 0, int num_thread = 0,
+  BaseSolver(const IoData &iodata, bool root, int size = 0, int num_thread = 0,
              const char *git_tag = nullptr);
   virtual ~BaseSolver() = default;
 
@@ -80,7 +78,8 @@ class BaseSolver
 
   // These methods write different simulation metadata to a JSON file in post_dir.
   void SaveMetadata(const mfem::ParFiniteElementSpace &fespace) const;
-  void SaveMetadata(int ksp_mult, int ksp_it) const;
+  template <typename SolverType>
+  void SaveMetadata(const SolverType &ksp) const;
   void SaveMetadata(const Timer &timer) const;
 };