Skip to content
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.

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

Remove the local copy of Jacobian matrix when doing LU decomposition #646

Merged
merged 1 commit into from
Sep 4, 2024
Merged
Changes from all commits
Commits
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
9 changes: 6 additions & 3 deletions include/micm/solver/rosenbrock.inl
Original file line number Diff line number Diff line change
Expand Up @@ -245,10 +245,9 @@ namespace micm
singular = false;
while (true)
{
auto jacobian = state.jacobian_;
double alpha = 1 / (H * gamma);
static_cast<Derived*>(this)->AlphaMinusJacobian(jacobian, alpha);
linear_solver_.Factor(jacobian, state.lower_matrix_, state.upper_matrix_, singular);
static_cast<Derived*>(this)->AlphaMinusJacobian(state.jacobian_, alpha);
linear_solver_.Factor(state.jacobian_, state.lower_matrix_, state.upper_matrix_, singular);
stats.decompositions_ += 1;

// if we are checking for singularity and the matrix is not singular, we can break the loop
Expand All @@ -260,6 +259,10 @@ namespace micm
if (++n_consecutive > 5)
break;
H /= 2;
// Reconstruct the Jacobian matrix if a substepping is performed here
state.jacobian_.Fill(0.0);
rates_.SubtractJacobianTerms(state.rate_constants_, number_densities, state.jacobian_);
stats.jacobian_updates_ += 1;
}
}

Expand Down
Loading