Update docs, JSON Schema, changelog

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

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>,
+    "Orthogonalization": <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.
+
+`"Orthogonalization" ["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/eigensolver.cpp

@@ -117,9 +117,7 @@ void EigenSolver::Solve(std::vector<std::unique_ptr<mfem::ParMesh>> &mesh,
     }
     slepc->SetProblemType(slepc::SlepcEigenSolver::ProblemType::GEN_NON_HERMITIAN);
     slepc->SetOrthogonalization(
-        iodata.solver.linear.orthog_type == config::LinearSolverData::OrthogType::MGS ||
-            iodata.solver.linear.orthog_type ==
-                config::LinearSolverData::OrthogType::DEFAULT,
+        iodata.solver.linear.orthog_type == config::LinearSolverData::OrthogType::MGS,
         iodata.solver.linear.orthog_type == config::LinearSolverData::OrthogType::CGS2);
     eigen = std::move(slepc);
 #endif

palace/linalg/ksp.cpp

@@ -95,30 +95,27 @@ std::unique_ptr<IterativeSolver<OperType>> ConfigureKrylovSolver(MPI_Comm comm,
   }
 
   // Configure orthogonalization method for GMRES/FMGRES.
-  if (iodata.solver.linear.orthog_type != config::LinearSolverData::OrthogType::DEFAULT)
+  if (type != config::LinearSolverData::KspType::GMRES &&
+      type != config::LinearSolverData::KspType::FGMRES)
   {
-    if (type != config::LinearSolverData::KspType::GMRES &&
-        type != config::LinearSolverData::KspType::FGMRES)
-    {
-      Mpi::Warning(comm, "Orthogonalization method will be ignored for non-GMRES/FGMRES "
-                         "iterative solvers!\n");
-    }
-    else
+    Mpi::Warning(comm, "Orthogonalization method will be ignored for non-GMRES/FGMRES "
+                       "iterative solvers!\n");
+  }
+  else
+  {
+    // Because FGMRES inherits from GMRES, this is OK.
+    auto *gmres = static_cast<GmresSolver<OperType> *>(ksp.get());
+    switch (iodata.solver.linear.orthog_type)
     {
-      // Because FGMRES inherits from GMRES, this is OK.
-      auto *gmres = static_cast<GmresSolver<OperType> *>(ksp.get());
-      switch (iodata.solver.linear.orthog_type)
-      {
-        case config::LinearSolverData::OrthogType::MGS:
-          gmres->SetOrthogonalization(GmresSolver<OperType>::OrthogType::MGS);
-          break;
-        case config::LinearSolverData::OrthogType::CGS:
-          gmres->SetOrthogonalization(GmresSolver<OperType>::OrthogType::CGS);
-          break;
-        case config::LinearSolverData::OrthogType::CGS2:
-          gmres->SetOrthogonalization(GmresSolver<OperType>::OrthogType::CGS2);
-          break;
-      }
+      case config::LinearSolverData::OrthogType::MGS:
+        gmres->SetOrthogonalization(GmresSolver<OperType>::OrthogType::MGS);
+        break;
+      case config::LinearSolverData::OrthogType::CGS:
+        gmres->SetOrthogonalization(GmresSolver<OperType>::OrthogType::CGS);
+        break;
+      case config::LinearSolverData::OrthogType::CGS2:
+        gmres->SetOrthogonalization(GmresSolver<OperType>::OrthogType::CGS2);
+        break;
     }
   }
 

palace/models/romoperator.cpp

@@ -123,8 +123,7 @@ RomOperator::RomOperator(const IoData &iodata, SpaceOperator &spaceop) : spaceop
   // The initial PROM basis is empty. Orthogonalization uses MGS by default, else CGS2.
   dim_V = 0;
   orthog_mgs =
-      (iodata.solver.linear.orthog_type == config::LinearSolverData::OrthogType::DEFAULT ||
-       iodata.solver.linear.orthog_type == config::LinearSolverData::OrthogType::MGS);
+      (iodata.solver.linear.orthog_type == config::LinearSolverData::OrthogType::MGS);
 
   // Seed the random number generator for parameter space sampling.
   engine.seed(std::chrono::system_clock::now().time_since_epoch().count());

palace/utils/configfile.cpp

@@ -1536,8 +1536,7 @@ NLOHMANN_JSON_SERIALIZE_ENUM(LinearSolverData::OrthogType,
                              {{LinearSolverData::OrthogType::INVALID, nullptr},
                               {LinearSolverData::OrthogType::MGS, "MGS"},
                               {LinearSolverData::OrthogType::CGS, "CGS"},
-                              {LinearSolverData::OrthogType::CGS2, "CGS2"},
-                              {LinearSolverData::OrthogType::DEFAULT, "Default"}})
+                              {LinearSolverData::OrthogType::CGS2, "CGS2"}})
 
 void LinearSolverData::SetUp(json &solver)
 {

palace/utils/configfile.hpp

@@ -825,10 +825,9 @@ struct LinearSolverData
     MGS,
     CGS,
     CGS2,
-    DEFAULT,
     INVALID = -1
   };
-  OrthogType orthog_type = OrthogType::DEFAULT;
+  OrthogType orthog_type = OrthogType::MGS;
 
   void SetUp(json &solver);
 };

palace/utils/geodata.cpp

@@ -190,7 +190,7 @@ void RefineMesh(const IoData &iodata, std::vector<std::unique_ptr<mfem::ParMesh>
        mesh[0]->MeshGenerator() & 8))
   {
     // XX TODO: Region-based refinement won't work if the ParMesh has been constructed from
-    //          a conforming mesh, but non-conforming refinement is needed. Unclear if the
+    //          a conforming mesh, but nonconforming refinement is needed. Unclear if the
     //          current mesh distribution scheme will work even for a conforming serial mesh
     //          which is a NCMesh after Mesh::EnsureNCMesh is called.
     MFEM_ABORT("Region-based refinement is currently only supported for simplex meshes!");

scripts/schema/config/solver.json

@@ -100,17 +100,16 @@
         "Tol": { "type": "number", "minimum": 0.0 },
         "MaxIts": { "type": "integer", "exclusiveMinimum": 0 },
         "MaxSize": { "type": "integer", "exclusiveMinimum": 0 },
-        "UseMGS": { "type": "boolean" },
-        "UseCGS2": { "type": "boolean" },
         "UseInitialGuess": { "type": "boolean" },
-        "UseKSPPiped": { "type": "boolean" },
-        "UseGMG": { "type": "boolean" },
-        "UseLOR": { "type": "boolean" },
-        "UsePCShifted": { "type": "boolean" },
+        "UsePartialAssembly": { "type": "boolean" },
+        "UseLowOrderRefined": { "type": "boolean" },
+        "UsePCMatShifted": { "type": "boolean" },
+        "PCSide": { "type": "string" },
+        "UseMultigrid": { "type": "boolean" },
+        "MGAuxiliarySmoother": { "type": "boolean" },
         "MGCycleIts": { "type": "integer", "exclusiveMinimum": 0 },
         "MGSmoothIts": { "type": "integer", "exclusiveMinimum": 0 },
         "MGSmoothOrder": { "type": "integer", "exclusiveMinimum": 0 },
-        "PrecondSide": { "type": "string" },
         "Reordering": { "type": "string" },
         "STRUMPACKCompressionType": { "type": "string" },
         "STRUMPACKCompressionTol": { "type": "number", "minimum": 0.0 },
@@ -119,7 +118,8 @@
         "SuperLU3D": { "type": "boolean" },
         "AMSVector": { "type": "boolean" },
         "DivFreeTol": { "type": "number", "minimum": 0.0 },
-        "DivFreeMaxIts": { "type": "integer", "minimum": 0 }
+        "DivFreeMaxIts": { "type": "integer", "minimum": 0 },
+        "Orthogonalization": { "type": "string" }
       }
     }
   }