Skip to content

Fixed Step IDA #332

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
Steven-Roberts wants to merge 3 commits into
base: develop
Choose a base branch
from
Draft

Fixed Step IDA #332

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
de77286
Add fixed time step for IDA forward simulation
Steven-Roberts Jan 9, 2026
d6e668c
Add fixed steps to adjoint
Steven-Roberts Jan 9, 2026
391c6d0
Apply pre-commmit fixes
Steven-Roberts Jan 9, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
87 changes: 64 additions & 23 deletions GridKit/Solver/Dynamic/Ida.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,8 +16,8 @@ namespace AnalysisManager
{

template <class ScalarT, typename IdxT>
Ida<ScalarT, IdxT>::Ida(GridKit::Model::Evaluator<ScalarT, IdxT>* model)
: DynamicSolver<ScalarT, IdxT>(model)
Ida<ScalarT, IdxT>::Ida(GridKit::Model::Evaluator<ScalarT, IdxT>* model, RealT dt)
: DynamicSolver<ScalarT, IdxT>(model), dt_(dt)
{
int retval = 0;

Expand Down Expand Up @@ -82,20 +82,7 @@ namespace AnalysisManager
retval = IDASetUserData(solver_, model_);
checkOutput(retval, "IDASetUserData");

// Set tolerances
RealT rel_tol;
RealT abs_tol;

model_->setTolerances(rel_tol, abs_tol); ///< \todo Function name should be "getTolerances"!
retval = IDASStolerances(solver_, rel_tol, abs_tol);
checkOutput(retval, "IDASStolerances");

IdxT msa;
model_->setMaxSteps(msa);

/// \todo Need to set max number of steps based on user input!
retval = IDASetMaxNumSteps(solver_, static_cast<long>(msa));
checkOutput(retval, "IDASetMaxNumSteps");
configureIDA(solver_);

// Tag differential variables
std::vector<bool>& tag = model_->tag();
Expand Down Expand Up @@ -155,6 +142,65 @@ namespace AnalysisManager
return retval;
}

/**
* @brief Configure IDA for forward or backward simulation
*
* @tparam ScalarT
* @tparam IdxT
*/
template <class ScalarT, typename IdxT>
void Ida<ScalarT, IdxT>::configureIDA(void* ida)
{
int retval = 0;

// Set tolerances
RealT rel_tol;
RealT abs_tol;

model_->setTolerances(rel_tol, abs_tol); ///< \todo Function name should be "getTolerances"!

if (dt_ > 0)
{
retval = IDASetMinStep(ida, dt_);
checkOutput(retval, "IDASetMinStep");
retval = IDASetMaxStep(ida, dt_);
checkOutput(retval, "IDASetMaxStep");

/* Since the starting procedure is first order, the maximum global order
* of convergence is two */
IDASetMaxOrd(ida, 2);
checkOutput(retval, "IDASetMaxOrd");

/* Enable more nonlinear iterations because a failed nonlinear solve
* causes a failed integration with fixed steps */
retval = IDASetMaxNonlinIters(ida, 100);
checkOutput(retval, "IDASetMaxNonlinIters");

// Set a large tolerance so the error test will never fail
static constexpr RealT FIXED_STEP_TOL_FAC = 1e10;
retval = IDASStolerances(ida, FIXED_STEP_TOL_FAC, FIXED_STEP_TOL_FAC * abs_tol / rel_tol);
checkOutput(retval, "IDASStolerances");

/* We want the nonlinear solver tolerance to be ~rel_tol, but the with
* the large tolerances set above, we need to choose this tolerance to
* "undo" the fac scaling. */
retval = IDASetNonlinConvCoef(ida, rel_tol / FIXED_STEP_TOL_FAC);
checkOutput(retval, "IDASetNonlinConvCoef");
}
else
{
retval = IDASStolerances(ida, rel_tol, abs_tol);
checkOutput(retval, "IDASStolerances");
}

IdxT msa;
model_->setMaxSteps(msa);

/// \todo Need to set max number of steps based on user input!
retval = IDASetMaxNumSteps(ida, static_cast<long>(msa));
checkOutput(retval, "IDASetMaxNumSteps");
}

#ifdef GRIDKIT_ENABLE_SUNDIALS_SPARSE
/**
* @brief Configure a sparse linear solver
Expand Down Expand Up @@ -559,16 +605,10 @@ namespace AnalysisManager
retval = IDAInitB(solver_, backwardID_, this->adjointResidual, tf, yyB_, ypB_);
checkOutput(retval, "IDAInitB");

model_->setTolerances(rel_tol, abs_tol);
retval = IDASStolerancesB(solver_, backwardID_, rel_tol, abs_tol);
checkOutput(retval, "IDASStolerancesB");

retval = IDASetUserDataB(solver_, backwardID_, model_);
checkOutput(retval, "IDASetUserDataB");

/// \todo Need to set max number of steps based on user input!
retval = IDASetMaxNumStepsB(solver_, backwardID_, 2000);
checkOutput(retval, "IDASetMaxNumSteps");
configureIDA(IDAGetAdjIDABmem(solver_, backwardID_));

// Allocate Jacobian matrix, if not already
if (JacobianMatB_ == nullptr)
Expand All @@ -595,6 +635,7 @@ namespace AnalysisManager
checkOutput(retval, "IDAQuadInitB");

// retval = IDAQuadSStolerancesB(solver_, backwardID_, rel_tol*1.1, abs_tol*1.1);
model_->setTolerances(rel_tol, abs_tol);
retval = IDAQuadSStolerancesB(solver_, backwardID_, rel_tol * 0.1, abs_tol * 0.1);
checkOutput(retval, "IDAQuadSStolerancesB");

Expand Down
6 changes: 5 additions & 1 deletion GridKit/Solver/Dynamic/Ida.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,7 @@ namespace AnalysisManager
using RealT = typename GridKit::ScalarTraits<ScalarT>::RealT;

public:
Ida(GridKit::Model::Evaluator<ScalarT, IdxT>* model);
Ida(GridKit::Model::Evaluator<ScalarT, IdxT>* model, RealT dt = 0);
~Ida();

int configureSimulation();
Expand Down Expand Up @@ -156,6 +156,8 @@ namespace AnalysisManager
void* user_data);

private:
RealT dt_;

void* solver_{};
SUNContext context_{};
SUNMatrix JacobianMat_{};
Expand All @@ -182,6 +184,8 @@ namespace AnalysisManager
int backwardID_{};

private:
void configureIDA(void* ida);

// static void copyMat(Model::Evaluator::Mat& J, SlsMat Jida);
static void copyVec(const N_Vector x, std::vector<ScalarT>& y);
static void copyVec(const std::vector<ScalarT>& x, N_Vector y);
Expand Down