diff --git a/.envrc b/.envrc new file mode 100644 index 0000000000..3550a30f2d --- /dev/null +++ b/.envrc @@ -0,0 +1 @@ +use flake diff --git a/GridKit/Model/PhasorDynamics/Bus/BusSignalVoltageOut/BusSignalVoltageOut.hpp b/GridKit/Model/PhasorDynamics/Bus/BusSignalVoltageOut/BusSignalVoltageOut.hpp index e951eb9b24..0e32ff0ee8 100644 --- a/GridKit/Model/PhasorDynamics/Bus/BusSignalVoltageOut/BusSignalVoltageOut.hpp +++ b/GridKit/Model/PhasorDynamics/Bus/BusSignalVoltageOut/BusSignalVoltageOut.hpp @@ -23,15 +23,15 @@ namespace GridKit * - publishes _Vr_ and _Vi_ on signal outlets `vr` and `vi`, and * - sets its residuals f[0] and f[1] to the current injections read from * signal inlets `ir` and `ii`, respectively. Both inlets are - * mandatory: verify() reports an error if either is not connected to a linked - * signal, and no default value is ever used. + * mandatory: verify() reports an error if either is not connected to a + * linked signal, and no default value is ever used. * * Components attached to the bus directly (without signals) keep adding * their currents to the residuals after the bus residual is evaluated, * exactly as they do for @ref Bus. * - * @note Signal ports have to be connected before allocate() is called, since - * the output signals are linked to the bus variables there. + * @note Signal ports have to be connected before allocate() is called, + * since the output signals are linked to the bus variables there. */ template class BusSignalVoltageOut : public BusBase @@ -67,6 +67,7 @@ namespace GridKit virtual ~BusSignalVoltageOut(); virtual int setBusID(IdxT) override final; + virtual int setFault(bool status, RealT R, RealT X) override final; virtual int allocate() override final; virtual int tagDifferentiable() override final; virtual int setAbsoluteTolerance(RealT rel_tol) override final; @@ -192,6 +193,8 @@ namespace GridKit private: ScalarT Vr0_{0.0}; ScalarT Vi0_{0.0}; + RealT fault_g_{0.0}; + RealT fault_b_{0.0}; /// Signal ports SignalPortsT ports_; diff --git a/GridKit/Model/PhasorDynamics/Bus/BusSignalVoltageOut/BusSignalVoltageOutEnzyme.cpp b/GridKit/Model/PhasorDynamics/Bus/BusSignalVoltageOut/BusSignalVoltageOutEnzyme.cpp index de7559b94e..f70d4fbdb1 100644 --- a/GridKit/Model/PhasorDynamics/Bus/BusSignalVoltageOut/BusSignalVoltageOutEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/Bus/BusSignalVoltageOut/BusSignalVoltageOutEnzyme.cpp @@ -12,9 +12,8 @@ namespace GridKit /** * @brief Jacobian evaluation for the Enzyme build. * - * The first four entries are zero placeholders for the bus voltage - * columns, as in Bus. They provide the indices for entries that other - * components contribute to and that are later deduplicated. + * This sets the fault admittance terms of the bus which are zero when + * cleared. * * One additional entry per connected signal inlet holds the derivative of * the residual with respect to the signal variable, which is one. When @@ -58,10 +57,10 @@ namespace GridKit J_cols_buffer_[1] = variable_indices_.at(1); J_cols_buffer_[2] = variable_indices_.at(0); J_cols_buffer_[3] = variable_indices_.at(1); - J_vals_buffer_[0] = 0.0; - J_vals_buffer_[1] = 0.0; - J_vals_buffer_[2] = 0.0; - J_vals_buffer_[3] = 0.0; + J_vals_buffer_[0] = -fault_g_; + J_vals_buffer_[1] = fault_b_; + J_vals_buffer_[2] = -fault_b_; + J_vals_buffer_[3] = -fault_g_; IdxT k = num_bus_entries; @@ -71,8 +70,9 @@ namespace GridKit { return; } - const IdxT signal_index = port.linked() ? port.signalVariableIndex() : INVALID_INDEX; - J_rows_buffer_[k] = residual_indices_.at(residual); + const IdxT signal_index = + port.linked() ? port.signalVariableIndex() : INVALID_INDEX; + J_rows_buffer_[k] = residual_indices_.at(residual); if (signal_index != INVALID_INDEX) { J_cols_buffer_[k] = signal_index; diff --git a/GridKit/Model/PhasorDynamics/Bus/BusSignalVoltageOut/BusSignalVoltageOutImpl.hpp b/GridKit/Model/PhasorDynamics/Bus/BusSignalVoltageOut/BusSignalVoltageOutImpl.hpp index 053dbbbc8c..d535c654ca 100644 --- a/GridKit/Model/PhasorDynamics/Bus/BusSignalVoltageOut/BusSignalVoltageOutImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Bus/BusSignalVoltageOut/BusSignalVoltageOutImpl.hpp @@ -46,8 +46,7 @@ namespace GridKit */ template BusSignalVoltageOut::BusSignalVoltageOut(const ModelDataT& data) - : Vr0_(data.Vr0), - Vi0_(data.Vi0) + : Vr0_(data.Vr0), Vi0_(data.Vi0) { bus_id_ = data.bus_id; size_ = 2; @@ -57,8 +56,9 @@ namespace GridKit { return Vr(); }); monitor_->set(Variable::Vi, [this] { return Vi(); }); - monitor_->set(Variable::Vm, [this] - { return std::sqrt(Vr() * Vr() + Vi() * Vi()); }); + monitor_->set( + Variable::Vm, [this] + { return std::sqrt(Vr() * Vr() + Vi() * Vi()); }); monitor_->set(Variable::Va, [this] { return std::atan2(Vi(), Vr()); }); } @@ -144,11 +144,15 @@ namespace GridKit const auto& port = ports_.in.template port(); if (!port.connected()) { - checks.fail(std::format("{} signal inlet is not connected; a default current is not allowed", name)); + checks.fail(std::format( + "{} signal inlet is not connected; a default current is not " + "allowed", + name)); } else if (!port.linked()) { - checks.fail(std::format("{} signal attached with no linked source", name)); + checks.fail( + std::format("{} signal attached with no linked source", name)); } }; @@ -157,7 +161,10 @@ namespace GridKit const auto& port = ports_.out.template port(); if (port.connected() && !port.linked()) { - checks.fail(std::format("{} signal attached but not linked; connect signal ports before allocate()", name)); + checks.fail(std::format( + "{} signal attached but not linked; connect signal ports before " + "allocate()", + name)); } }; @@ -179,6 +186,28 @@ namespace GridKit return 0; } + template + int BusSignalVoltageOut::setFault(bool status, RealT R, RealT X) + { + if (status && !(ZERO < R * R + X * X)) + { + Log::error() + << "BusSignalVoltageOut: fault impedance R + jX must be nonzero\n"; + return 1; + } + + fault_g_ = ZERO; + fault_b_ = ZERO; + + if (status) + { + fault_g_ = R / (R * R + X * X); + fault_b_ = -X / (R * R + X * X); + } + + return 0; + } + /*! * @brief Bus variables are algebraic. */ @@ -235,9 +264,9 @@ namespace GridKit * * Residuals f[0] and f[1] are set to the values read from the `ir` and * `ii` signal inlets, respectively. Both inlets are mandatory; verify() - * reports an error if either is not connected to a linked signal, and no default - * value is used here. Components attached to the bus add their currents - * afterwards. + * reports an error if either is not connected to a linked signal, and no + * default value is used here. Components attached to the bus add their + * currents afterwards. * * @warning This implementation assumes bus residuals are always evaluated * _before_ component model residuals. @@ -247,8 +276,10 @@ namespace GridKit { auto* f = f_.getData(); - f[0] = ports_.in.template port().readSignal(); - f[1] = ports_.in.template port().readSignal(); + f[0] = ports_.in.template port().readSignal() + - (fault_g_ * Vr() - fault_b_ * Vi()); + f[1] = ports_.in.template port().readSignal() + - (fault_b_ * Vr() + fault_g_ * Vi()); f_.setDataUpdated(); return 0; diff --git a/GridKit/Model/PhasorDynamics/Bus/BusSignalVoltageOut/README.md b/GridKit/Model/PhasorDynamics/Bus/BusSignalVoltageOut/README.md index 83642e8d0d..773909c587 100644 --- a/GridKit/Model/PhasorDynamics/Bus/BusSignalVoltageOut/README.md +++ b/GridKit/Model/PhasorDynamics/Bus/BusSignalVoltageOut/README.md @@ -1,16 +1,16 @@ # BusSignalVoltageOut -`BusSignalVoltageOut` is a bus with signal ports. It owns the same voltage variables -and current-balance residuals as `Bus`, publishes the voltage components on -output signals, and adds current injections received on input signals to its -residuals. Devices attached directly to the bus still add their currents after -the bus residual is evaluated. +`BusSignalVoltageOut` is a bus with signal ports. It owns the same voltage +variables and current-balance residuals as `Bus`, publishes the voltage +components on output signals, and adds current injections received on input +signals to its residuals. Devices attached directly to the bus still add their +currents after the bus residual is evaluated. ## Notes -- Signal ports must be connected before `allocate()` is called. The signal outlets - are linked to the bus voltage variables and their system indices in - `allocate()`. +- Signal ports must be connected before `allocate()` is called. The signal + outlets are linked to the bus voltage variables and their system indices + in `allocate()`. - Both current inlets are mandatory. `verify()` reports an error if a current inlet is not connected or not linked, or if a connected outlet is not linked. No default current is ever used. @@ -22,12 +22,14 @@ Same as `Bus`. ## Model Ports -Signal | Direction | Units | Description | Note ------|-----------|--------|----------------------------------------------------|----- -`vr` | out | [p.u.] | Bus voltage, real component $V_r$ | -`vi` | out | [p.u.] | Bus voltage, imaginary component $V_i$ | -`ir` | in | [p.u.] | Current injection, real component $I_r^{s}$ | Required, sets $f_0$ -`ii` | in | [p.u.] | Current injection, imaginary component $I_i^{s}$ | Required, sets $f_1$ +Signal | Direction | Units | Description +-------|-----------|-------|----------------------------------------- +`vr` | out | [p.u.] | Bus voltage, real component $V_r$ +`vi` | out | [p.u.] | Bus voltage, imaginary component $V_i$ +`ir` | in | [p.u.] | Current injection, real component $I_r^{s}$ +`ii` | in | [p.u.] | Current injection, imaginary component $I_i^{s}$ + +Both current inputs are required. `ir` sets $f_0$, and `ii` sets $f_1$. ## Model Variables @@ -55,17 +57,18 @@ $I_i^{s}$ | [p.u.] | Current injection on signal inlet `ii` | #### Algebraic -Let $\mathcal{D}$ denote the set of devices attached directly to the bus. +Let $\mathcal{D}$ denote the set of devices attached directly to the bus +and $G_f$, $B_f$ be the conductance and susceptance of a fault applied +to the bus, respectively. During a fault, $G_f + jB_f = 1/(R_f + jX_f)$. +Otherwise, both are zero. ```math \begin{aligned} -0 &= I_r^{s} + \sum_{d \in \mathcal{D}} I_{r,d} \\ -0 &= I_i^{s} + \sum_{d \in \mathcal{D}} I_{i,d} +0 &= I_r^{s} - (G_f V_r - B_f V_i) + \sum_{d \in \mathcal{D}} I_{r,d} \\ +0 &= I_i^{s} - (B_f V_r + G_f V_i) + \sum_{d \in \mathcal{D}} I_{i,d} \end{aligned} ``` -Both signal inlets must be connected; see `verify()`. - ## Initialization Same as `Bus`. @@ -76,7 +79,7 @@ Same as `Bus`. ## Testing -Unit tests in `tests/UnitTests/PhasorDynamics/BusSignalVoltageOutTests.hpp` cover -construction, output signal linking, residual evaluation with input signals, -verification of unlinked inputs, dependency-tracking derivatives, and (when -Enzyme is enabled) the sparse Jacobian entries. +Unit tests in `tests/UnitTests/PhasorDynamics/BusSignalVoltageOutTests.hpp` +cover construction, output signal linking, residual evaluation with input +signals, verification of unlinked inputs, dependency-tracking derivatives, +and (when Enzyme is enabled) the sparse Jacobian entries. diff --git a/GridKit/Solver/CMakeLists.txt b/GridKit/Solver/CMakeLists.txt index febd6808b3..e9b4b145e9 100644 --- a/GridKit/Solver/CMakeLists.txt +++ b/GridKit/Solver/CMakeLists.txt @@ -4,6 +4,8 @@ # - Slaven Peles #]] +install(FILES SundialsException.hpp DESTINATION include/GridKit/Solver) + if(TARGET SUNDIALS::kinsol) add_subdirectory(SteadyState) endif() diff --git a/GridKit/Solver/Dynamic/Ida.cpp b/GridKit/Solver/Dynamic/Ida.cpp index 0d3d4118a8..0a9eaa8873 100644 --- a/GridKit/Solver/Dynamic/Ida.cpp +++ b/GridKit/Solver/Dynamic/Ida.cpp @@ -33,9 +33,10 @@ namespace AnalysisManager } /** - * @brief Destroy the Ida< Scalar T, Idx T>:: Ida object + * @brief Destroy the @ref Ida object * - * @note if sysmodel is freed before this will fail. May want something agnostic to this + * @note if sysmodel is freed before this will fail + * @todo consider making this agnostic to the above * * @tparam ScalarT * @tparam IdxT @@ -96,7 +97,14 @@ namespace AnalysisManager retval = IDASetId(solver_, tag_); checkOutput(retval, "IDASetId"); - setIDAOptions(solver_, time_step_, rel_tol_, abs_tol_override_, max_steps_, max_order_, suppress_alg_); + setIDAOptions( + solver_, + time_step_, + rel_tol_, + abs_tol_override_, + max_steps_, + max_order_, + suppress_alg_); // Set up linear solver return this->configureLinearSolver(); @@ -116,7 +124,7 @@ namespace AnalysisManager { int retval = 0; -/// Todo - Implement a cleaner way to handle the Sparse versus Dense solvers +// TODO: Implement a cleaner way to handle the Sparse versus Dense solvers #ifdef GRIDKIT_ENABLE_SUNDIALS_SPARSE if (model_->hasJacobian()) { @@ -127,10 +135,11 @@ namespace AnalysisManager this->configureLinearSolverDense(); } #else - /// Todo - Improve error handling capabilities and hasJacobian_ ownership + // TODO: Improve error handling capabilities and hasJacobian_ ownership if (model_->hasJacobian()) { - Log::warning() << "SUNDIALS is not configured with KLU, but the model has a (sparse) Jacobian. " + Log::warning() << "SUNDIALS is not configured with KLU, but the model " + "has a (sparse) Jacobian. " << "Falling back to dense Jacobian.\n"; } this->configureLinearSolverDense(); @@ -153,14 +162,11 @@ namespace AnalysisManager { int retval = 0; - sunindextype n = static_cast(model_->size()); - sunindextype nnz = static_cast(model_->getCsrJacobian()->getNnz()); + sunindextype n = static_cast(model_->size()); + sunindextype nnz = + static_cast(model_->getCsrJacobian()->getNnz()); - JacobianMat_ = SUNSparseMatrix(n, - n, - nnz, - CSR_MAT, - context_); + JacobianMat_ = SUNSparseMatrix(n, n, nnz, CSR_MAT, context_); checkAllocation((void*) JacobianMat_, "SUNSparseMatrix"); linearSolver_ = SUNLinSol_KLU(yy_, JacobianMat_, context_); @@ -191,9 +197,10 @@ namespace AnalysisManager { int retval = 0; - JacobianMat_ = SUNDenseMatrix(static_cast(model_->size()), - static_cast(model_->size()), - context_); + JacobianMat_ = SUNDenseMatrix( + static_cast(model_->size()), + static_cast(model_->size()), + context_); checkAllocation((void*) JacobianMat_, "SUNDenseMatrix"); linearSolver_ = SUNLinSol_Dense(yy_, JacobianMat_, context_); @@ -229,7 +236,8 @@ namespace AnalysisManager * @tparam IdxT */ template - int Ida::initializeSimulation(RealT t0, bool findConsistent) + int Ida::initializeSimulation( + RealT t0, bool findConsistent, bool printMonitor) { int retval = 0; @@ -243,7 +251,8 @@ namespace AnalysisManager if (findConsistent) { const int consistentICType = getIDAConsistentICType(); - retval = IDACalcIC(solver_, consistentICType, t0 + 0.1); + + retval = IDACalcIC(solver_, consistentICType, t0 + 0.1); checkOutput(retval, "IDACalcIC"); retval = IDAGetConsistentIC(solver_, yy_, yp_); @@ -253,7 +262,7 @@ namespace AnalysisManager copyVec(yp_, model_->yp()); } - if (model_->monitoring()) + if (printMonitor && model_->monitoring()) { updateModelState(t0); model_->printMonitoredVariables(); @@ -277,10 +286,10 @@ namespace AnalysisManager return 1; } - const RealT n_est = (tf - t_init_) / dt; - const RealT epsilon = std::numeric_limits::epsilon() - * std::max({std::abs(t_init_), std::abs(tf), RealT(1.0)}) - / dt; + const RealT n_est = (tf - t_init_) / dt; + const RealT epsilon = + std::numeric_limits::epsilon() + * std::max({std::abs(t_init_), std::abs(tf), RealT(1.0)}) / dt; return static_cast(std::ceil(n_est - epsilon)); } @@ -319,20 +328,23 @@ namespace AnalysisManager default: GridKit::Utilities::Logger::error() << "Invalid IDA consistent initial condition type " - << static_cast(consistent_ic_type_) - << ".\n"; + << static_cast(consistent_ic_type_) << ".\n"; throw SundialsException(); } } /** - * @brief Run the IDA solver and optionally produce monitor output every `dt_monitor`. + * @brief Run the IDA solver and optionally produce monitor output every + * `dt_monitor`. * * When `dt_monitor` is zero, the simulation runs directly to the final * time. The final time is always solved and monitored. */ template - int Ida::runSimulation(RealT tf, RealT dt_monitor, const std::optional step_callback) + int Ida::runSimulation( + RealT tf, + RealT dt_monitor, + const std::optional step_callback) { int retval = 0; int nsteps = getStepCount(tf, dt_monitor); @@ -349,21 +361,23 @@ namespace AnalysisManager template int Ida::runSimulationStep( - RealT tf, RealT dt_monitor, int step, int nsteps, const std::optional callback) + RealT tf, + RealT dt_monitor, + int step, + int nsteps, + const std::optional callback, + bool printMonitor) { const RealT tout = getMonitorTime(tf, dt_monitor, step, nsteps); RealT tret; auto retval = IDASolve(solver_, tout, &tret, yy_, yp_, IDA_NORMAL); checkOutput(retval, "IDASolve"); + // always allow observation of the updated values even without a callback + updateModelState(tret); if (callback.has_value() || model_->monitoring()) { - // The callback may try to observe updated values in the model, so we - // should update them here (At this point, the model's values are one - // internal integrator step out of date) - updateModelState(tret); - - if (model_->monitoring()) + if (printMonitor && model_->monitoring()) { model_->printMonitoredVariables(); } @@ -408,7 +422,8 @@ namespace AnalysisManager int retval = 0; // Create and initialize quadratures - q_ = N_VNew_Serial(static_cast(model_->sizeQuadrature()), context_); + q_ = N_VNew_Serial( + static_cast(model_->sizeQuadrature()), context_); checkAllocation((void*) q_, "N_VNew_Serial"); // Set integrand function and allocate quadrature workspace @@ -416,7 +431,8 @@ namespace AnalysisManager checkOutput(retval, "IDAQuadInit"); // Set tolerances and error control for quadratures - setQuadratureTolerance(solver_, quadrature_rel_tol_, quadrature_abs_tol_override_); + setQuadratureTolerance( + solver_, quadrature_rel_tol_, quadrature_abs_tol_override_); // Include quadrature in error checking retval = IDASetQuadErrCon(solver_, SUNTRUE); @@ -511,7 +527,8 @@ namespace AnalysisManager ypB_ = N_VClone(yyB_); checkAllocation((void*) ypB_, "N_VClone"); - qB_ = N_VNew_Serial(static_cast(model_->sizeParams()), context_); + qB_ = N_VNew_Serial( + static_cast(model_->sizeParams()), context_); checkAllocation((void*) qB_, "N_VNew_Serial"); return 0; @@ -556,16 +573,18 @@ namespace AnalysisManager checkOutput(retval, "IDACreateB"); // IDAInitB must be called after forward simulation run. - retval = IDAInitB(solver_, backwardID_, this->adjointResidual, tf, yyB_, ypB_); + retval = + IDAInitB(solver_, backwardID_, this->adjointResidual, tf, yyB_, ypB_); checkOutput(retval, "IDAInitB"); - setIDAOptions(IDAGetAdjIDABmem(solver_, backwardID_), - backward_time_step_, - backward_rel_tol_, - backward_abs_tol_override_, - backward_max_steps_, - backward_max_order_, - backward_suppress_alg_); + setIDAOptions( + IDAGetAdjIDABmem(solver_, backwardID_), + backward_time_step_, + backward_rel_tol_, + backward_abs_tol_override_, + backward_max_steps_, + backward_max_order_, + backward_suppress_alg_); retval = IDASetUserDataB(solver_, backwardID_, model_); checkOutput(retval, "IDASetUserDataB"); @@ -573,9 +592,10 @@ namespace AnalysisManager // Allocate Jacobian matrix, if not already if (JacobianMatB_ == nullptr) { - JacobianMatB_ = SUNDenseMatrix(static_cast(model_->size()), - static_cast(model_->size()), - context_); + JacobianMatB_ = SUNDenseMatrix( + static_cast(model_->size()), + static_cast(model_->size()), + context_); checkAllocation((void*) JacobianMatB_, "SUNDenseMatrix"); } @@ -587,14 +607,18 @@ namespace AnalysisManager } // Setup linear solver (only dense supported at this time) - retval = IDASetLinearSolverB(solver_, backwardID_, linearSolverB_, JacobianMatB_); + retval = IDASetLinearSolverB( + solver_, backwardID_, linearSolverB_, JacobianMatB_); checkOutput(retval, "IDASetLinearSolverB"); // Also reinitialize quadratures. retval = IDAQuadInitB(solver_, backwardID_, this->adjointIntegrand, qB_); checkOutput(retval, "IDAQuadInitB"); - setQuadratureTolerance(IDAGetAdjIDABmem(solver_, backwardID_), backward_quadrature_rel_tol_, backward_quadrature_abs_tol_override_); + setQuadratureTolerance( + IDAGetAdjIDABmem(solver_, backwardID_), + backward_quadrature_rel_tol_, + backward_quadrature_abs_tol_override_); // Include quadratures in error control retval = IDASetQuadErrConB(solver_, backwardID_, SUNTRUE); @@ -617,9 +641,10 @@ namespace AnalysisManager int retval = 0; // Create Jacobian matrix - JacobianMatB_ = SUNDenseMatrix(static_cast(model_->size()), - static_cast(model_->size()), - context_); + JacobianMatB_ = SUNDenseMatrix( + static_cast(model_->size()), + static_cast(model_->size()), + context_); checkAllocation((void*) JacobianMatB_, "SUNDenseMatrix"); // Create linear solver @@ -627,14 +652,16 @@ namespace AnalysisManager checkAllocation((void*) linearSolverB_, "SUNLinSol_Dense"); // Attach linear solver to IDA - retval = IDASetLinearSolverB(solver_, backwardID_, linearSolverB_, JacobianMatB_); + retval = IDASetLinearSolverB( + solver_, backwardID_, linearSolverB_, JacobianMatB_); checkOutput(retval, "IDASetLinearSolverB"); return retval; } /** - * @brief Run forward simulation for adjoint analysis with optional output every `dt_monitor`. + * @brief Run forward simulation for adjoint analysis with optional output + * every `dt_monitor`. * * When `dt_monitor` is zero, the simulation runs directly to the final * time. The final time is always solved. @@ -685,7 +712,8 @@ namespace AnalysisManager checkOutput(retval, "IDASolveB"); IDAGetNumSteps(IDAGetAdjIDABmem(solver_, backwardID_), &nstB); - Log::summary() << "Backward integration for adjoint analysis done (nst = " << nstB << ")\n"; + Log::summary() << "Backward integration for adjoint analysis done (nst = " + << nstB << ")\n"; retval = IDAGetB(solver_, backwardID_, &time, yyB_, ypB_); checkOutput(retval, "IDAGetB"); @@ -738,9 +766,11 @@ namespace AnalysisManager * @tparam IdxT */ template - int Ida::Residual(RealT tres, N_Vector yy, N_Vector yp, N_Vector rr, void* user_data) + int Ida::Residual( + RealT tres, N_Vector yy, N_Vector yp, N_Vector rr, void* user_data) { - GridKit::Model::Evaluator* model = static_cast*>(user_data); + GridKit::Model::Evaluator* model = + static_cast*>(user_data); copyVec(yy, model->y()); copyVec(yp, model->yp()); @@ -761,9 +791,20 @@ namespace AnalysisManager * @tparam IdxT */ template - int Ida::Jac(RealT t, RealT cj, N_Vector yy, N_Vector yp, N_Vector, SUNMatrix J, void* user_data, N_Vector, N_Vector, N_Vector) + int Ida::Jac( + RealT t, + RealT cj, + N_Vector yy, + N_Vector yp, + N_Vector, + SUNMatrix J, + void* user_data, + N_Vector, + N_Vector, + N_Vector) { - GridKit::Model::Evaluator* model = static_cast*>(user_data); + GridKit::Model::Evaluator* model = + static_cast*>(user_data); copyVec(yy, model->y()); copyVec(yp, model->yp()); @@ -803,9 +844,11 @@ namespace AnalysisManager * @tparam IdxT */ template - int Ida::Integrand(RealT tt, N_Vector yy, N_Vector yp, N_Vector rhsQ, void* user_data) + int Ida::Integrand( + RealT tt, N_Vector yy, N_Vector yp, N_Vector rhsQ, void* user_data) { - GridKit::Model::Evaluator* model = static_cast*>(user_data); + GridKit::Model::Evaluator* model = + static_cast*>(user_data); copyVec(yy, model->y()); copyVec(yp, model->yp()); @@ -824,9 +867,17 @@ namespace AnalysisManager * @tparam IdxT */ template - int Ida::adjointResidual(RealT tt, N_Vector yy, N_Vector yp, N_Vector yyB, N_Vector ypB, N_Vector rrB, void* user_data) + int Ida::adjointResidual( + RealT tt, + N_Vector yy, + N_Vector yp, + N_Vector yyB, + N_Vector ypB, + N_Vector rrB, + void* user_data) { - GridKit::Model::Evaluator* model = static_cast*>(user_data); + GridKit::Model::Evaluator* model = + static_cast*>(user_data); copyVec(yy, model->y()); copyVec(yp, model->yp()); @@ -847,9 +898,17 @@ namespace AnalysisManager * @tparam IdxT */ template - int Ida::adjointIntegrand(RealT tt, N_Vector yy, N_Vector yp, N_Vector yyB, N_Vector ypB, N_Vector rhsQB, void* user_data) + int Ida::adjointIntegrand( + RealT tt, + N_Vector yy, + N_Vector yp, + N_Vector yyB, + N_Vector ypB, + N_Vector rhsQB, + void* user_data) { - GridKit::Model::Evaluator* model = static_cast*>(user_data); + GridKit::Model::Evaluator* model = + static_cast*>(user_data); copyVec(yy, model->y()); copyVec(yp, model->yp()); @@ -876,8 +935,9 @@ namespace AnalysisManager const auto ysize = static_cast(y.getSize()); if (xsize != ysize) { - GridKit::Utilities::Logger::error() << "N_Vector size (" << xsize << ") does not match vector size (" - << y.getSize() << ").\n"; + GridKit::Utilities::Logger::error() + << "N_Vector size (" << xsize << ") does not match vector size (" + << y.getSize() << ").\n"; throw SundialsException(); } @@ -899,8 +959,9 @@ namespace AnalysisManager const auto xsize = static_cast(x.getSize()); if (xsize != ysize) { - GridKit::Utilities::Logger::error() << "vector size (" << x.getSize() << ") does not match N_Vector size (" - << ysize << ").\n"; + GridKit::Utilities::Logger::error() + << "vector size (" << x.getSize() + << ") does not match N_Vector size (" << ysize << ").\n"; throw SundialsException(); } @@ -921,8 +982,9 @@ namespace AnalysisManager const auto ysize = static_cast(N_VGetLength(y)); if (x.size() != ysize) { - GridKit::Utilities::Logger::error() << "std::vector size (" << x.size() << ") does not match N_Vector size (" - << ysize << ").\n"; + GridKit::Utilities::Logger::error() + << "std::vector size (" << x.size() + << ") does not match N_Vector size (" << ysize << ").\n"; throw SundialsException(); } @@ -988,24 +1050,25 @@ namespace AnalysisManager } /** - * @brief Accumulate another stats object into this one, allowing for stats to be kept - * across multiple simulations with IDA + * @brief Accumulate another stats object into this one, allowing for stats + * to be kept across multiple simulations with IDA */ IdaStats& IdaStats::operator+=(const IdaStats& other) { - num_steps_ += other.num_steps_; - num_residual_evals_ += other.num_residual_evals_; - num_linear_decompositions_ += other.num_linear_decompositions_; - num_error_test_fails_ += other.num_error_test_fails_; - num_nonlinear_iters_ += other.num_nonlinear_iters_; - num_nonlinear_convergence_fails_ += other.num_nonlinear_convergence_fails_; + num_steps_ += other.num_steps_; + num_residual_evals_ += other.num_residual_evals_; + num_linear_decompositions_ += other.num_linear_decompositions_; + num_error_test_fails_ += other.num_error_test_fails_; + num_nonlinear_iters_ += other.num_nonlinear_iters_; + num_nonlinear_convergence_fails_ += + other.num_nonlinear_convergence_fails_; return *this; } /** - * @brief Generate a string containing all of the stats in a formatted report. - * All columns are aligned. To change the width of a column, + * @brief Generate a string containing all of the stats in a formatted + * report. All columns are aligned. To change the width of a column, * modify `label_width` or `stat_width`. */ std::string IdaStats::report() const @@ -1014,20 +1077,28 @@ namespace AnalysisManager int stat_width = 12; std::stringstream out; - out << std::setw(label_width) << "Steps" << " : " << std::setw(stat_width) << num_steps_ << '\n' - << std::setw(label_width) << "Residual evals" << " : " << std::setw(stat_width) << num_residual_evals_ << '\n' - << std::setw(label_width) << "Linear decompositions" << " : " << std::setw(stat_width) << num_linear_decompositions_ << '\n' - << std::setw(label_width) << "Error test failures" << " : " << std::setw(stat_width) << num_error_test_fails_ << '\n' - << std::setw(label_width) << "Nonlinear iterations" << " : " << std::setw(stat_width) << num_nonlinear_iters_ << '\n' - << std::setw(label_width) << "Nonlinear convergence failures" << " : " << std::setw(stat_width) << num_nonlinear_convergence_fails_; + out << std::setw(label_width) << "Steps" << " : " << std::setw(stat_width) + << num_steps_ << '\n' + << std::setw(label_width) << "Residual evals" << " : " + << std::setw(stat_width) << num_residual_evals_ << '\n' + << std::setw(label_width) << "Linear decompositions" << " : " + << std::setw(stat_width) << num_linear_decompositions_ << '\n' + << std::setw(label_width) << "Error test failures" << " : " + << std::setw(stat_width) << num_error_test_fails_ << '\n' + << std::setw(label_width) << "Nonlinear iterations" << " : " + << std::setw(stat_width) << num_nonlinear_iters_ << '\n' + << std::setw(label_width) << "Nonlinear convergence failures" << " : " + << std::setw(stat_width) << num_nonlinear_convergence_fails_; return out.str(); } /** - * @brief Construct and return an `IdaStats` object containing the statistics of the current IDA workspace. - * Several statistics returned by IDA are ignored because they are about the current state of IDA, - * rather than about the simulation at large. + * @brief Construct and return an `IdaStats` object containing the + * statistics of the current IDA workspace. + * + * Several statistics returned by IDA are ignored because they are about the + * current state of IDA, rather than about the simulation at large. */ template IdaStats Ida::getStats() const @@ -1038,20 +1109,24 @@ namespace AnalysisManager int dummy; sunrealtype dummy2; - int retval = IDAGetIntegratorStats(solver_, - &stats.num_steps_, - &stats.num_residual_evals_, - &stats.num_linear_decompositions_, - &stats.num_error_test_fails_, - &dummy, - &dummy, - &dummy2, - &dummy2, - &dummy2, - &dummy2); + int retval = IDAGetIntegratorStats( + solver_, + &stats.num_steps_, + &stats.num_residual_evals_, + &stats.num_linear_decompositions_, + &stats.num_error_test_fails_, + &dummy, + &dummy, + &dummy2, + &dummy2, + &dummy2, + &dummy2); checkOutput(retval, "IDAGetIntegratorStats"); - retval = IDAGetNonlinSolvStats(solver_, &stats.num_nonlinear_iters_, &stats.num_nonlinear_convergence_fails_); + retval = IDAGetNonlinSolvStats( + solver_, + &stats.num_nonlinear_iters_, + &stats.num_nonlinear_convergence_fails_); checkOutput(retval, "IDAGetNonlinSolvStats"); return stats; @@ -1068,7 +1143,9 @@ namespace AnalysisManager { if (v == NULL) { - GridKit::Utilities::Logger::error() << "Function " << functionName << " failed -- returned NULL pointer!\n"; + GridKit::Utilities::Logger::error() + << "Function " << functionName + << " failed -- returned NULL pointer!\n"; throw SundialsException(); } } @@ -1085,7 +1162,9 @@ namespace AnalysisManager if (retval < 0) { char* msg = IDAGetReturnFlagName(retval); - GridKit::Utilities::Logger::error() << "Function " << functionName << " failed with flag " << msg << "!\n"; + GridKit::Utilities::Logger::error() + << "Function " << functionName << " failed with flag " << msg + << "!\n"; std::free(msg); throw SundialsException(); } @@ -1130,8 +1209,8 @@ namespace AnalysisManager * @tparam IdxT Index data type */ template - void Ida::setTolerance(ScalarT rel_tol, - ScalarT abs_tol_override) + void Ida::setTolerance( + ScalarT rel_tol, ScalarT abs_tol_override) { rel_tol_ = rel_tol; abs_tol_override_ = abs_tol_override; @@ -1149,8 +1228,8 @@ namespace AnalysisManager * @tparam IdxT Index data type */ template - void Ida::setBackwardTolerance(ScalarT rel_tol, - ScalarT abs_tol_override) + void Ida::setBackwardTolerance( + ScalarT rel_tol, ScalarT abs_tol_override) { backward_rel_tol_ = rel_tol; backward_abs_tol_override_ = abs_tol_override; @@ -1168,8 +1247,7 @@ namespace AnalysisManager * @tparam IdxT Index data type */ template - void Ida::setQuadratureTolerance(ScalarT rel_tol, - ScalarT abs_tol_override) + void Ida::setQuadratureTolerance(ScalarT rel_tol, ScalarT abs_tol_override) { quadrature_rel_tol_ = rel_tol; quadrature_abs_tol_override_ = abs_tol_override; @@ -1187,15 +1265,15 @@ namespace AnalysisManager * @tparam IdxT Index data type */ template - void Ida::setBackwardQuadratureTolerance(ScalarT rel_tol, - ScalarT abs_tol_override) + void Ida::setBackwardQuadratureTolerance(ScalarT rel_tol, ScalarT abs_tol_override) { backward_quadrature_rel_tol_ = rel_tol; backward_quadrature_abs_tol_override_ = abs_tol_override; } /** - * @brief Set whether IDA suppresses local error tests on algebraic variables + * @brief Set whether IDA suppresses local error tests on algebraic + * variables * * @param suppress If true, algebraic variables are excluded from IDA's * local error test @@ -1276,7 +1354,8 @@ namespace AnalysisManager } /** - * @brief Set the maximum integration method order for the backward simulation + * @brief Set the maximum integration method order for the backward + * simulation * * @param max_order The maximum integration method order * @tparam ScalarT Scalar data type @@ -1305,13 +1384,14 @@ namespace AnalysisManager * @tparam IdxT Index data type */ template - void Ida::setIDAOptions(void* mem, - ScalarT time_step, - ScalarT rel_tol, - ScalarT abs_tol_override, - IdxT max_steps, - int max_order, - bool suppress_alg) + void Ida::setIDAOptions( + void* mem, + ScalarT time_step, + ScalarT rel_tol, + ScalarT abs_tol_override, + IdxT max_steps, + int max_order, + bool suppress_alg) { int retval = 0; retval = IDASetMinStep(mem, time_step); @@ -1320,7 +1400,8 @@ namespace AnalysisManager checkOutput(retval, "IDASetMaxStep"); retval = IDASetMaxNumSteps(mem, static_cast(max_steps)); checkOutput(retval, "IDASetMaxNumSteps"); - retval = IDASetMaxOrd(mem, time_step == 0 ? max_order : std::min(max_order, 2)); + retval = + IDASetMaxOrd(mem, time_step == 0 ? max_order : std::min(max_order, 2)); checkOutput(retval, "IDASetMaxOrd"); retval = IDASetSuppressAlg(mem, suppress_alg ? SUNTRUE : SUNFALSE); checkOutput(retval, "IDASetSuppressAlg"); @@ -1334,7 +1415,8 @@ namespace AnalysisManager /* Enable more nonlinear iterations because a failed nonlinear solve * causes a failed integration with fixed steps */ static constexpr int FIXED_STEP_NONLIN_ITRS = 16; - retval = IDASetMaxNonlinIters(mem, FIXED_STEP_NONLIN_ITRS); + + retval = IDASetMaxNonlinIters(mem, FIXED_STEP_NONLIN_ITRS); checkOutput(retval, "IDASetMaxNonlinIters"); // Set a large tolerance so the error test will never fail @@ -1349,12 +1431,14 @@ namespace AnalysisManager static constexpr RealT DEFAULT_NONLIN_CONV_COEF_IC = 0.0033; - retval = IDASetNonlinConvCoefIC(mem, DEFAULT_NONLIN_CONV_COEF_IC / FIXED_STEP_TOL_FAC); + retval = IDASetNonlinConvCoefIC( + mem, DEFAULT_NONLIN_CONV_COEF_IC / FIXED_STEP_TOL_FAC); checkOutput(retval, "IDASetNonlinConvCoefIC"); // Scale the IC line-search threshold with the same error weights. - const RealT step_tol_ic = std::pow(std::numeric_limits::epsilon(), 2.0 / 3.0); - retval = IDASetStepToleranceIC(mem, step_tol_ic / FIXED_STEP_TOL_FAC); + const RealT step_tol_ic = + std::pow(std::numeric_limits::epsilon(), 2.0 / 3.0); + retval = IDASetStepToleranceIC(mem, step_tol_ic / FIXED_STEP_TOL_FAC); checkOutput(retval, "IDASetStepToleranceIC"); } } @@ -1373,16 +1457,15 @@ namespace AnalysisManager * @tparam IdxT Index data type */ template - void Ida::setTolerance(void* mem, - ScalarT rel_tol, - ScalarT abs_tol_override, - ScalarT tol_fac) + void Ida::setTolerance( + void* mem, ScalarT rel_tol, ScalarT abs_tol_override, ScalarT tol_fac) { int retval = 0; if (abs_tol_override > 0) { - retval = IDASStolerances(mem, tol_fac * rel_tol, tol_fac * abs_tol_override); + retval = + IDASStolerances(mem, tol_fac * rel_tol, tol_fac * abs_tol_override); checkOutput(retval, "IDASStolerances"); return; } @@ -1414,9 +1497,8 @@ namespace AnalysisManager * @tparam IdxT Index data type */ template - void Ida::setQuadratureTolerance(void* mem, - ScalarT rel_tol, - ScalarT abs_tol_override) + void Ida::setQuadratureTolerance( + void* mem, ScalarT rel_tol, ScalarT abs_tol_override) { int retval = 0; @@ -1438,10 +1520,10 @@ namespace AnalysisManager N_VDestroy(abs_tol_vec); } - // Compiler will prevent building modules with data type incompatible with sunrealtype + // Compiler will prevent building modules with data type incompatible with + // sunrealtype template class Ida; template class Ida; template class Ida; - } // namespace Sundials } // namespace AnalysisManager diff --git a/GridKit/Solver/Dynamic/Ida.hpp b/GridKit/Solver/Dynamic/Ida.hpp index 9a782f5d93..435b642bc0 100644 --- a/GridKit/Solver/Dynamic/Ida.hpp +++ b/GridKit/Solver/Dynamic/Ida.hpp @@ -8,8 +8,8 @@ #include #include -#include /* access to dense linear solver */ -#include /* access to sparse SUNMatrix */ +#include /* dense linear solver */ +#include /* sparse SUNMatrix */ #include @@ -19,6 +19,7 @@ #include #include +#include #include namespace AnalysisManager @@ -67,16 +68,22 @@ namespace AnalysisManager #endif int configureLinearSolverDense(); int getDefaultInitialCondition(); - int initializeSimulation(RealT t0, bool findConsistent = true); + int initializeSimulation( + RealT t0, bool findConsistent = true, bool printMonitor = true); int getStepCount(RealT tf, RealT dt) const; - int runSimulation(RealT tf, RealT dt_monitor = 0, std::optional step_callback = {}); - int runSimulationStep(RealT tf, - RealT dt_monitor, - int step, - int nsteps, - std::optional callback = {}); + int runSimulation( + RealT tf, + RealT dt_monitor = 0, + std::optional step_callback = {}); + int runSimulationStep( + RealT tf, + RealT dt_monitor, + int step, + int nsteps, + std::optional callback = {}, + bool printMonitor = true); int deleteSimulation(); int configureQuadrature(); @@ -141,10 +148,8 @@ namespace AnalysisManager using DynamicSolver::setTolerance; void setTolerance(ScalarT rel_tol, ScalarT abs_tol_override) override; void setBackwardTolerance(ScalarT rel_tol, ScalarT abs_tol_override = 0); - void setQuadratureTolerance(ScalarT rel_tol, - ScalarT abs_tol_override = 0); - void setBackwardQuadratureTolerance(ScalarT rel_tol, - ScalarT abs_tol_override = 0); + void setQuadratureTolerance(ScalarT rel_tol, ScalarT abs_tol_override = 0); + void setBackwardQuadratureTolerance(ScalarT rel_tol, ScalarT abs_tol_override = 0); void setSuppressAlgebraicErrors(bool suppress); void setBackwardSuppressAlgebraicErrors(bool suppress); void setConsistentICType(IdaConsistentICType consistent_ic_type); @@ -156,44 +161,39 @@ namespace AnalysisManager IdaStats getStats() const; private: - static int Residual(RealT t, - N_Vector yy, - N_Vector yp, - N_Vector rr, - void* user_data); - - static int Jac(RealT t, - RealT cj, - N_Vector yy, - N_Vector yp, - N_Vector resvec, - SUNMatrix J, - void* user_data, - N_Vector tmp1, - N_Vector tmp2, - N_Vector tmp3); - - static int Integrand(RealT t, - N_Vector yy, - N_Vector yp, - N_Vector rhsQ, - void* user_data); - - static int adjointResidual(RealT t, - N_Vector yy, - N_Vector yp, - N_Vector yyB, - N_Vector ypB, - N_Vector rrB, - void* user_data); - - static int adjointIntegrand(RealT t, - N_Vector yy, - N_Vector yp, - N_Vector yyB, - N_Vector ypB, - N_Vector rhsQB, - void* user_data); + static int Residual(RealT t, N_Vector yy, N_Vector yp, N_Vector rr, void* user_data); + + static int Jac( + RealT t, + RealT cj, + N_Vector yy, + N_Vector yp, + N_Vector resvec, + SUNMatrix J, + void* user_data, + N_Vector tmp1, + N_Vector tmp2, + N_Vector tmp3); + + static int Integrand(RealT t, N_Vector yy, N_Vector yp, N_Vector rhsQ, void* user_data); + + static int adjointResidual( + RealT t, + N_Vector yy, + N_Vector yp, + N_Vector yyB, + N_Vector ypB, + N_Vector rrB, + void* user_data); + + static int adjointIntegrand( + RealT t, + N_Vector yy, + N_Vector yp, + N_Vector yyB, + N_Vector ypB, + N_Vector rhsQB, + void* user_data); RealT getMonitorTime(RealT tf, RealT dt_monitor, int step, int nsteps) const; int getIDAConsistentICType() const; @@ -256,29 +256,20 @@ namespace AnalysisManager static void checkAllocation(void* v, const char* functionName); static void checkOutput(int retval, const char* functionName); - void setIDAOptions(void* mem, - ScalarT time_step, - ScalarT rel_tol, - ScalarT abs_tol_override, - IdxT max_steps, - int max_order, - bool suppress_alg); - void setTolerance(void* mem, - ScalarT rel_tol, - ScalarT abs_tol_override, - ScalarT tol_fac = 1); - void setQuadratureTolerance(void* mem, - ScalarT rel_tol, - ScalarT abs_tol_override); - }; - - /// Simple exception to use within Ida class. - class SundialsException : public std::exception - { - virtual const char* what() const throw() - { - return "Method in Ida class failed!\n"; - } + void setIDAOptions( + void* mem, + ScalarT time_step, + ScalarT rel_tol, + ScalarT abs_tol_override, + IdxT max_steps, + int max_order, + bool suppress_alg); + void setTolerance( + void* mem, + ScalarT rel_tol, + ScalarT abs_tol_override, + ScalarT tol_fac = 1); + void setQuadratureTolerance(void* mem, ScalarT rel_tol, ScalarT abs_tol_override); }; } // namespace Sundials diff --git a/GridKit/Solver/SteadyState/CMakeLists.txt b/GridKit/Solver/SteadyState/CMakeLists.txt index 576d796e4c..39e6f68da2 100644 --- a/GridKit/Solver/SteadyState/CMakeLists.txt +++ b/GridKit/Solver/SteadyState/CMakeLists.txt @@ -9,4 +9,10 @@ gridkit_add_library( solvers_steady SOURCES Kinsol.cpp HEADERS Kinsol.hpp SteadyStateSolver.hpp - LINK_LIBRARIES PUBLIC SUNDIALS::kinsol GridKit::dense_vector) + LINK_LIBRARIES + PUBLIC + SUNDIALS::kinsol + PRIVATE + GridKit::dense_vector + PRIVATE + GridKit::utilities_logger) diff --git a/GridKit/Solver/SteadyState/Kinsol.cpp b/GridKit/Solver/SteadyState/Kinsol.cpp index 9cb3907a66..ac018c0577 100644 --- a/GridKit/Solver/SteadyState/Kinsol.cpp +++ b/GridKit/Solver/SteadyState/Kinsol.cpp @@ -1,17 +1,19 @@ - /** * @file Kinsol.cpp * @author Slaven Peles * * Contains definition of interface to KINSOL nonlinear solver from * SUNDIALS library. - * */ #include "Kinsol.hpp" +#include +#include #include #include +#include +#include #include // access to KINSOL func., consts. #include // access to serial N_Vector @@ -19,17 +21,31 @@ #include // access to dense SUNMatrix #include -#include namespace AnalysisManager { - namespace Sundials { - template Kinsol::Kinsol(GridKit::Model::Evaluator* model) - : SteadyStateSolver(model) + : Kinsol(model->size()) + { + model_ = model; + residual_ = + [model](std::span values, std::span residual) + { + std::copy(values.begin(), values.end(), model->y().getData()); + model->y().setDataUpdated(); + const int status = model->evaluateResidual(); + std::copy_n( + model->getResidual().getData(), residual.size(), residual.begin()); + return status; + }; + } + + template + Kinsol::Kinsol(IdxT size) + : SteadyStateSolver(nullptr), size_(size) { int retval = 0; @@ -51,26 +67,35 @@ namespace AnalysisManager template int Kinsol::configureSimulation() { + if (model_) + { + size_ = model_->size(); + } + if (size_ <= 0) + { + throw std::invalid_argument("KINSOL requires a positive problem size"); + } + int retval = 0; // Allocate solution vectors - yy_ = N_VNew_Serial(static_cast(model_->size()), context_); - checkAllocation((void*) yy_, "N_VNew_Serial"); + yy_ = N_VNew_Serial(static_cast(size_), context_); + checkAllocation(context_, (void*) yy_, "N_VNew_Serial"); // Allocate scaling vector scale_ = N_VClone(yy_); - checkAllocation((void*) scale_, "N_VClone"); + checkAllocation(context_, (void*) scale_, "N_VClone"); // Create vectors to store restart initial condition yy0_ = N_VClone(yy_); - checkAllocation((void*) yy0_, "N_VClone"); + checkAllocation(context_, (void*) yy0_, "N_VClone"); // Allocate and initialize KIN workspace retval = KINInit(solver_, this->Residual, yy_); checkOutput(retval, "KINInit"); - // Set pointer to model data - retval = KINSetUserData(solver_, model_); + // Set the callback context + retval = KINSetUserData(solver_, this); checkOutput(retval, "KINSetUserData"); // Set up linear solver @@ -83,13 +108,14 @@ namespace AnalysisManager int retval = 0; // Set up linear solver - JacobianMat_ = SUNDenseMatrix(static_cast(model_->size()), - static_cast(model_->size()), - context_); - checkAllocation((void*) JacobianMat_, "SUNDenseMatrix"); + JacobianMat_ = SUNDenseMatrix( + static_cast(size_), + static_cast(size_), + context_); + checkAllocation(context_, (void*) JacobianMat_, "SUNDenseMatrix"); linearSolver_ = SUNLinSol_Dense(yy_, JacobianMat_, context_); - checkAllocation((void*) linearSolver_, "SUNLinSol_Dense"); + checkAllocation(context_, (void*) linearSolver_, "SUNLinSol_Dense"); retval = KINSetLinearSolver(solver_, linearSolver_, JacobianMat_); checkOutput(retval, "KINSetLinearSolver"); @@ -100,6 +126,12 @@ namespace AnalysisManager template int Kinsol::getDefaultInitialCondition() { + if (!model_) + { + throw std::logic_error( + "KINSOL requires an explicit initial guess when used with a callback"); + } + model_->initialize(); copyVec(model_->y(), yy_); @@ -110,9 +142,56 @@ namespace AnalysisManager template int Kinsol::runSimulation() { - int retval = 0; N_VConst(1.0, scale_); - retval = KINSol(solver_, yy_, KIN_LINESEARCH, scale_, scale_); + return solve(); + } + + template + int Kinsol::runSimulation(std::span values, + ResidualFunction residual, + bool line_search) + { + if (values.size() != static_cast(size_) || !residual) + { + throw std::invalid_argument( + "KINSOL requires a matching initial guess and a residual callback"); + } + + // ensure the model callback is restored on scope exit + struct Restore + { + ResidualFunction& active; + ResidualFunction previous; + + ~Restore() + { + active.swap(previous); + } + } restore{residual_, std::move(residual_)}; + + residual_ = std::move(residual); + std::copy(values.begin(), values.end(), N_VGetArrayPointer(yy_)); + N_VConst(1.0, scale_); + const int status = solve(line_search); + std::copy_n(N_VGetArrayPointer(yy_), values.size(), values.begin()); + return status; + } + + template + int Kinsol::solve(bool line_search) + { + if (!residual_) + { + throw std::logic_error("KINSOL requires a residual callback"); + } + + callback_error_ = nullptr; + const int retval = KINSol( + solver_, yy_, line_search ? KIN_LINESEARCH : KIN_NONE, scale_, scale_); + if (callback_error_) + { + std::rethrow_exception(callback_error_); + } checkOutput(retval, "KINSol"); return retval; } @@ -132,15 +211,18 @@ namespace AnalysisManager template int Kinsol::Residual(N_Vector yy, N_Vector rr, void* user_data) { - GridKit::Model::Evaluator* model = - static_cast*>(user_data); - - copyVec(yy, model->y()); - - model->evaluateResidual(); - copyVec(model->getResidual(), rr); - - return 0; + auto& self = *static_cast(user_data); + try + { + return self.residual_( + {N_VGetArrayPointer(yy), static_cast(self.size_)}, + {N_VGetArrayPointer(rr), static_cast(self.size_)}); + } + catch (...) + { + self.callback_error_ = std::current_exception(); + return -1; + } } template @@ -165,7 +247,7 @@ namespace AnalysisManager sunrealtype* yval = N_VGetArrayPointer(yy_); std::cout << std::setprecision(5) << std::setw(7); - for (IdxT i = 0; i < model_->size(); ++i) + for (IdxT i = 0; i < size_; ++i) { std::cout << yval[i] << " "; } @@ -194,11 +276,14 @@ namespace AnalysisManager } template - void Kinsol::checkAllocation(void* v, const char* functionName) + void Kinsol::checkAllocation(SUNContext c, void* v, const char* functionName) { if (v == NULL) { - GridKit::Utilities::Logger::error() << "Function " << functionName << " failed -- returned NULL pointer!\n"; + // NOTE: this will only work if called immediately after the function + auto errorCode = SUNContext_GetLastError(c); + Logger::error() << std::format( + "Function {} failed: {}\n", functionName, SUNGetErrMsg(errorCode)); throw SundialsException(); } } @@ -208,7 +293,8 @@ namespace AnalysisManager { if (retval < 0) { - GridKit::Utilities::Logger::error() << "Function " << functionName << " failed with flag " << retval << "!\n"; + Logger::error() << std::format( + "Function {} failed: {}\n", functionName, SUNGetErrMsg(retval)); throw SundialsException(); } } @@ -216,17 +302,36 @@ namespace AnalysisManager template void Kinsol::setTolerance(ScalarT tol) { - int retval = KINSetFuncNormTol(solver_, tol); + setTolerance(tol, tol); + } + + template + void Kinsol::setTolerance(ScalarT residual_tol, ScalarT step_tol) + { + int retval = KINSetFuncNormTol(solver_, residual_tol); checkOutput(retval, "KINSetFuncNormTol"); - retval = KINSetScaledStepTol(solver_, tol); + retval = KINSetScaledStepTol(solver_, step_tol); checkOutput(retval, "KINSetScaledStepTol"); } - // Compiler will prevent building modules with data type incompatible with sunrealtype + template + void Kinsol::setFunctionError(RealT error) + { + checkOutput(KINSetRelErrFunc(solver_, error), "KINSetRelErrFunc"); + } + + template + void Kinsol::setJacobianUpdateInterval(long int iterations) + { + checkOutput( + KINSetMaxSetupCalls(solver_, iterations), "KINSetMaxSetupCalls"); + } + + // Compiler will prevent building modules with data type incompatible with + // sunrealtype template class Kinsol; template class Kinsol; template class Kinsol; - } // namespace Sundials } // namespace AnalysisManager diff --git a/GridKit/Solver/SteadyState/Kinsol.hpp b/GridKit/Solver/SteadyState/Kinsol.hpp index ef992d44c4..adfc834cbd 100644 --- a/GridKit/Solver/SteadyState/Kinsol.hpp +++ b/GridKit/Solver/SteadyState/Kinsol.hpp @@ -5,20 +5,23 @@ * * Contains declaration of interface to KINSOL nonlinear solver from * SUNDIALS library. - * */ + #pragma once #include +#include #include +#include #include -#include /* access to sparse SUNMatrix */ -// #include /* access to KLU linear solver */ -#include /* access to dense linear solver */ +#include +#include #include #include +#include +#include namespace AnalysisManager { @@ -28,14 +31,22 @@ namespace AnalysisManager class Kinsol : public SteadyStateSolver { using SteadyStateSolver::model_; + using Logger = ::GridKit::Utilities::Logger; + using RealT = typename GridKit::ScalarTraits::RealT; + using VectorT = + typename GridKit::Model::Evaluator::VectorT; - using RealT = typename GridKit::ScalarTraits::RealT; - using VectorT = typename GridKit::Model::Evaluator::VectorT; - - static_assert(std::is_same_v, "RealT must be the same type as sunrealtype"); + static_assert( + std::is_same_v, + "RealT must be the same type as sunrealtype"); public: + using ResidualFunction = + std::function, std::span)>; + Kinsol(GridKit::Model::Evaluator* model); + explicit Kinsol(IdxT size); + ~Kinsol(); int configureSimulation(); @@ -44,8 +55,33 @@ namespace AnalysisManager // int setIntegrationTime(RealT t_init, RealT t_final, int nout); // int initializeSimulation(); int runSimulation(); + + /// Use the given values as the initial guess and return the accepted + /// solution in the same vector. + /// + /// @param[in,out] values initial guess, overwritten with the solution + /// @param residual function that fills the residual for trial values + /// and returns zero on success + /// @param line_search whether or not to use KINSOL's backtracking line + /// search + int runSimulation( + std::span values, + ResidualFunction residual, + bool line_search = true); + int deleteSimulation(); + void setTolerance(ScalarT tol) override; + void setTolerance(ScalarT residual_tol, ScalarT step_tol); + + /// Estimated relative error of residual evaluations. + /// + /// This is used to choose finite-difference perturbations. + void setFunctionError(RealT error); + + /// Maximum number of Newton iterations between Jacobian evaluations. + void setJacobianUpdateInterval(long int iterations); + // int configureQuadrature(); // int initializeQuadrature(); // int runSimulationQuadrature(RealT tf, int nout=1); @@ -112,6 +148,7 @@ namespace AnalysisManager private: static int Residual(N_Vector yy, N_Vector rr, void* user_data); + int solve(bool line_search = true); // static int Integrand(sunrealtype t, // N_Vector yy, N_Vector yp, @@ -133,6 +170,10 @@ namespace AnalysisManager SUNMatrix JacobianMat_{}; SUNLinearSolver linearSolver_{}; + IdxT size_; + ResidualFunction residual_; + std::exception_ptr callback_error_; + N_Vector yy_{}; ///< Solution vector N_Vector scale_{}; ///< Scaling vector N_Vector q_{}; ///< Integrand vector @@ -146,19 +187,8 @@ namespace AnalysisManager static void copyVec(const std::vector& x, N_Vector y); // int check_flag(void *flagvalue, const char *funcname, int opt); - static inline void checkAllocation(void* v, const char* functionName); + static inline void checkAllocation(SUNContext ctx, void* v, const char* functionName); static inline void checkOutput(int retval, const char* functionName); - - void setTolerance(ScalarT tol) override; - }; - - /// Simple exception to use within Kinsol class. - class SundialsException : public std::exception - { - virtual const char* what() const throw() - { - return "Method in Kinsol class failed!\n"; - } }; } // namespace Sundials diff --git a/GridKit/Solver/SundialsException.hpp b/GridKit/Solver/SundialsException.hpp new file mode 100644 index 0000000000..ca1bb2c007 --- /dev/null +++ b/GridKit/Solver/SundialsException.hpp @@ -0,0 +1,28 @@ +#pragma once + +#include +#include + +namespace AnalysisManager +{ + namespace Sundials + { + /// Exception wrapping an error that occurs within SUNDIALS. + class SundialsException : public std::exception + { + std::string message_; + + public: + explicit SundialsException( + const char* message = "A SUNDIALS class wrapper failed!\n") + : message_(message) + { + } + + const char* what() const noexcept override + { + return message_.c_str(); + } + }; + } // namespace Sundials +} // namespace AnalysisManager diff --git a/GridKit/Utilities/CliArgs/CliArgs.cpp b/GridKit/Utilities/CliArgs/CliArgs.cpp index 77a109e8fe..239728a5e8 100644 --- a/GridKit/Utilities/CliArgs/CliArgs.cpp +++ b/GridKit/Utilities/CliArgs/CliArgs.cpp @@ -4,6 +4,7 @@ #include #include #include +#include #include #include diff --git a/GridKit/Utilities/CliArgs/CliArgs.hpp b/GridKit/Utilities/CliArgs/CliArgs.hpp index 58a40b01b4..75610902a1 100644 --- a/GridKit/Utilities/CliArgs/CliArgs.hpp +++ b/GridKit/Utilities/CliArgs/CliArgs.hpp @@ -19,7 +19,6 @@ namespace GridKit { namespace Utilities { - /** * @brief Table with which to specify and parse command-line arguments * @@ -28,12 +27,11 @@ namespace GridKit * * ``` * CliArgs args{{.name = {"--opt1"}, - * .help = "an option with an argument"}, + * .help = "An option with an argument"}, * * {.name = {"--opt2", "-o"}, - * .help = "An option with an argument (and a short-hand)", - * .type = ArgType::Integer, - * .defaults = 0}, + * .help = "An option with an argument (and shorthand)", + * .type = ArgType::Integer, .defaults = 0}, * * {.name = {"--flag1", "-f"}, * .help = "A flag option (no argument)", @@ -41,12 +39,14 @@ namespace GridKit * ``` * * @note The `{"--help", "-h"}` argument is added automatically. When it is - * parsed from the command line. The help description will be printed and - * the application will exit. + * parsed from the command line. The help description will be printed + * and the application will exit. * * @note The usage syntax is intended to be similar to * [ArgParse.jl](https://carlobaldassi.github.io/ArgParse.jl/stable/). * + * @note `repeatable = true` allows a named option to occur multiple times. + * * @todo Allow option grouping for a more expressive help description. */ class CliArgs @@ -124,14 +124,13 @@ namespace GridKit * @brief Stream insertion operator for printing the argument table * * @note This is useful mainly for debugging purposes. It prints both the - * list of options and the option name map, so the same Option will be - * shown multiple times. Each Option is printed with its associated parsed - * argument values. + * list of options and the option name map, so the same Option will + * be shown multiple times. Each Option is printed with its + * associated parsed argument values. */ friend std::ostream& operator<<(std::ostream& os, const CliArgs& args); std::unique_ptr pImpl_; }; - } // namespace Utilities } // namespace GridKit diff --git a/GridKit/Utilities/CliArgs/CliArgsImpl.hpp b/GridKit/Utilities/CliArgs/CliArgsImpl.hpp index 44e890a3d2..2de5adb602 100644 --- a/GridKit/Utilities/CliArgs/CliArgsImpl.hpp +++ b/GridKit/Utilities/CliArgs/CliArgsImpl.hpp @@ -1,3 +1,5 @@ +#include + #include #include @@ -7,34 +9,28 @@ namespace GridKit { using Log = ::GridKit::Utilities::Logger; - /** - * @brief Associate an Option with its parsed arguments - */ + /// Associate an @ref Option with its parsed arguments. struct OptionArgumentPair { - /** - * @brief Construct with an Option spec and an empty set of values - */ + /// Construct an @ref OptionArgumentPair with a specification and an empty + /// set of values. explicit OptionArgumentPair(Option opt) : option(std::move(opt)) { } - // OptionArgumentPair(OptionArgumentPair + /// Reference to an existing option. + Option option; - /// Reference to existing option - Option option; - /// Parsed values + /// Parsed values. ArgVector values{}; }; - /** - * @brief Implementation of CliArgs - */ + /// Implementation of @ref CliArgs. struct CliArgsImpl { ///@{ - /// @brief Implementation of CliArgs + /// Implementation of CliArgs CliArgsImpl(std::initializer_list