From 6692a0062fac475e6c2a3dcefdfb764aa5810a47 Mon Sep 17 00:00:00 2001 From: Philipp Dumitrescu Date: Thu, 5 Dec 2024 18:43:26 -0800 Subject: [PATCH] Fix max_size meaning in PROM construction --- palace/drivers/drivensolver.cpp | 19 +++++++++---------- palace/models/romoperator.cpp | 13 ++++++++----- palace/models/romoperator.hpp | 2 +- 3 files changed, 18 insertions(+), 16 deletions(-) diff --git a/palace/drivers/drivensolver.cpp b/palace/drivers/drivensolver.cpp index 9f06389e2..9f19dab83 100644 --- a/palace/drivers/drivensolver.cpp +++ b/palace/drivers/drivensolver.cpp @@ -206,16 +206,15 @@ ErrorIndicator DrivenSolver::SweepAdaptive(SpaceOperator &space_op, // Configure PROM parameters if not specified. double offline_tol = iodata.solver.driven.adaptive_tol; - int max_size = iodata.solver.driven.adaptive_max_size; - MFEM_VERIFY(max_size <= 0 || max_size > 2, + int max_size_per_excitation = iodata.solver.driven.adaptive_max_size; + MFEM_VERIFY(max_size_per_excitation <= 0 || max_size_per_excitation > 2, "Adaptive frequency sweep must sample at least two frequency points!"); - if (max_size <= 0) + if (max_size_per_excitation <= 0) { - max_size = 20 * excitation_helper.Size(); // Default value + max_size_per_excitation = 20; // Default value } - max_size = std::min(max_size, - (n_step - step0) * - int(excitation_helper.Size())); // Maximum size dictated by sweep + // Maximum size — no more than nr steps needed + max_size_per_excitation = std::min(max_size_per_excitation, (n_step - step0)); int convergence_memory = iodata.solver.driven.adaptive_memory; // Allocate negative curl matrix for postprocessing the B-field and vectors for the @@ -247,7 +246,7 @@ ErrorIndicator DrivenSolver::SweepAdaptive(SpaceOperator &space_op, " {:d} points for frequency sweep over [{:.3e}, {:.3e}] GHz\n", n_step - step0, omega0 * unit_GHz, (omega0 + (n_step - step0 - 1) * delta_omega) * unit_GHz); - RomOperator prom_op(iodata, space_op, max_size); + RomOperator prom_op(iodata, space_op, max_size_per_excitation); space_op.GetWavePortOp().SetSuppressOutput( true); // Suppress wave port output for offline @@ -317,7 +316,7 @@ ErrorIndicator DrivenSolver::SweepAdaptive(SpaceOperator &space_op, { memory = 0; } - if (it == max_size) + if (it == max_size_per_excitation) { break; } @@ -335,7 +334,7 @@ ErrorIndicator DrivenSolver::SweepAdaptive(SpaceOperator &space_op, } Mpi::Print("\nAdaptive sampling{} {:d} frequency samples:\n" " n = {:d}, error = {:.3e}, tol = {:.3e}, memory = {:d}/{:d}\n", - (it == max_size) ? " reached maximum" : " converged with", it, + (it == max_size_per_excitation) ? " reached maximum" : " converged with", it, prom_op.GetReducedDimension(), max_errors.back(), offline_tol, memory, convergence_memory); utils::PrettyPrint(prom_op.GetSamplePoints(excitation_idx), unit_GHz, diff --git a/palace/models/romoperator.cpp b/palace/models/romoperator.cpp index 5f5446405..83e28c840 100644 --- a/palace/models/romoperator.cpp +++ b/palace/models/romoperator.cpp @@ -313,7 +313,8 @@ std::vector MinimalRationInterpolation::FindMaxError(int N) const return vals; } -RomOperator::RomOperator(const IoData &iodata, SpaceOperator &space_op, int max_size) +RomOperator::RomOperator(const IoData &iodata, SpaceOperator &space_op, + int max_size_per_excitation) : space_op(space_op) { // Construct the system matrices defining the linear operator. PEC boundaries are handled @@ -336,11 +337,14 @@ RomOperator::RomOperator(const IoData &iodata, SpaceOperator &space_op, int max_ ksp = std::make_unique(iodata, space_op.GetNDSpaces(), &space_op.GetH1Spaces()); + auto excitation_helper = space_op.BuildPortExcitationHelper(); + // The initial PROM basis is empty. The provided maximum dimension is the number of sample // points (2 basis vectors per point). Basis orthogonalization method is configured using // GMRES/FGMRES settings. - MFEM_VERIFY(max_size > 0, "Reduced order basis storage must have > 0 columns!"); - V.resize(2 * max_size, Vector()); + MFEM_VERIFY(max_size_per_excitation * excitation_helper.Size() > 0, + "Reduced order basis storage must have > 0 columns!"); + V.resize(2 * max_size_per_excitation * excitation_helper.Size(), Vector()); switch (iodata.solver.linear.gs_orthog_type) { case config::LinearSolverData::OrthogType::MGS: @@ -354,10 +358,9 @@ RomOperator::RomOperator(const IoData &iodata, SpaceOperator &space_op, int max_ break; } // Set up MRI - auto excitation_helper = space_op.BuildPortExcitationHelper(); for (const auto &[excitation_idx, data] : excitation_helper) { - mri.emplace(excitation_idx, MinimalRationInterpolation(max_size)); + mri.emplace(excitation_idx, MinimalRationInterpolation(max_size_per_excitation)); } } diff --git a/palace/models/romoperator.hpp b/palace/models/romoperator.hpp index 31710d192..f9b87061c 100644 --- a/palace/models/romoperator.hpp +++ b/palace/models/romoperator.hpp @@ -82,7 +82,7 @@ class RomOperator std::map mri; public: - RomOperator(const IoData &iodata, SpaceOperator &space_op, int max_size); + RomOperator(const IoData &iodata, SpaceOperator &space_op, int max_size_per_excitation); // Return the HDM linear solver. const ComplexKspSolver &GetLinearSolver() const { return *ksp; }