Skip to content

Fix the implicit bug #796

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

Merged
landinjm merged 5 commits into from
Feb 24, 2026
Merged

Fix the implicit bug #796

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
b86a451
fix the implicit bug
jijn Feb 5, 2026
7473374
this does not work
jijn Feb 14, 2026
24b1d37
update self_nonlinear
jijn Feb 23, 2026
3c40140
disable cavity flow and linear solve old
jijn Feb 23, 2026
55df6b6
update regression test
jijn Feb 23, 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
33 changes: 29 additions & 4 deletions src/core/solution_handler.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -415,10 +415,14 @@ SolutionHandler<dim, number>::update(FieldSolveType field_solve_type,
Types::Index variable_index)
{
// Helper function to swap vectors for all dependency types
auto swap_all_dependency_vectors = [this](Types::Index index, auto &new_vector)
auto swap_all_dependency_vectors =
[this](Types::Index index, auto &new_vector, bool is_implicit = false)
{
// Always swap the Normal dependency
new_vector->swap(*(solution_set.at(std::make_pair(index, DependencyType::Normal))));
if (!is_implicit)
{
new_vector->swap(
*(solution_set.at(std::make_pair(index, DependencyType::Normal))));
}

// Swap old dependency types if they exist
const std::array<DependencyType, 4> old_types = {
Expand All @@ -435,6 +439,22 @@ SolutionHandler<dim, number>::update(FieldSolveType field_solve_type,
new_vector->swap(*(solution_set.at(std::make_pair(index, dep_type))));
}
}

if (is_implicit)
{
if (solution_set.contains(std::make_pair(index, DependencyType::OldOne)))
{
*new_vector =
*(solution_set.at(std::make_pair(index, DependencyType::OldOne)));
new_vector->swap(
*(solution_set.at(std::make_pair(index, DependencyType::Normal))));
}
else // auxiliary don't have old
{
new_vector->swap(
*(solution_set.at(std::make_pair(index, DependencyType::Normal))));
}
}
};

// Loop through the solutions and swap them
Expand Down Expand Up @@ -466,13 +486,18 @@ SolutionHandler<dim, number>::update(FieldSolveType field_solve_type,
swap_all_dependency_vectors(index, new_vector);
break;
case FieldSolveType::NonexplicitLinear:
if (variable_index == index)
{
swap_all_dependency_vectors(index, new_vector, true);
}
break;
case FieldSolveType::NonexplicitAuxiliary:
case FieldSolveType::NonexplicitSelfnonlinear:
case FieldSolveType::NonexplicitCononlinear:
// For Nonexplicit types we swap all dependency types if the index matches
if (variable_index == index)
{
swap_all_dependency_vectors(index, new_vector);
swap_all_dependency_vectors(index, new_vector, false);
}
break;
default:
Expand Down
28 changes: 6 additions & 22 deletions src/solvers/sequential_co_nonlinear_solver.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -107,18 +107,6 @@ SequentialCoNonlinearSolver<dim, degree, number>::solve()
return;
}

// Grab the old solution fields in temporary variables
std::map<Types::Index, typename SolutionHandler<dim, number>::VectorType> old_solutions;
for (const auto &[index, variable] : this->get_subset_attributes())
{
if (variable.get_pde_type() != PDEType::Auxiliary)
{
old_solutions[index] =
*(this->get_solution_handler().get_solution_vector(index,
DependencyType::Normal));
}
}

// Set the convergence bool and iteration counter
bool unconverged = true;
unsigned int iteration = 0;
Expand Down Expand Up @@ -204,18 +192,14 @@ SequentialCoNonlinearSolver<dim, degree, number>::solve()
if (variable.get_pde_type() != PDEType::Auxiliary)
{
// The solve will have updated the "old solution" vector with the newton update
// so it's technically the new solution. In order to update the solutions and
// preserve the old states we copy the old solution from above and swap.
// so it's technically the new solution.
*(this->get_solution_handler().get_new_solution_vector(index)) =
old_solutions[index];
this->get_solution_handler()
.get_solution_vector(index, DependencyType::Normal)
->swap(*this->get_solution_handler().get_new_solution_vector(index));
*(this->get_solution_handler().get_solution_vector(index,
DependencyType::Normal));
this->get_solution_handler().update(this->get_field_solve_type(),
this->get_solve_block(),
index);
}

this->get_solution_handler().update(this->get_field_solve_type(),
this->get_solve_block(),
index);
}

// Update the ghosts
Expand Down
29 changes: 13 additions & 16 deletions src/solvers/sequential_self_nonlinear_solver.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -84,10 +84,6 @@ SequentialSelfNonlinearSolver<dim, degree, number>::solve()
// Solve each field
for (const auto &[index, variable] : this->get_subset_attributes())
{
// Grab the old solution in a temporary variable
const auto old_solution = *(
this->get_solution_handler().get_solution_vector(index, DependencyType::Normal));

// Set the convergence bool, iteration counter, and step length
bool unconverged = true;
unsigned int iteration = 0;
Expand Down Expand Up @@ -150,18 +146,19 @@ SequentialSelfNonlinearSolver<dim, degree, number>::solve()
iteration++;
}

// The solve will have updated the "old solution" vector with the newton update so
// it's technically the new solution. In order to update the solutions and preserve
// the old states we copy the old solution from above and swap.
*(this->get_solution_handler().get_new_solution_vector(index)) = old_solution;
this->get_solution_handler()
.get_solution_vector(index, DependencyType::Normal)
->swap(*this->get_solution_handler().get_new_solution_vector(index));

// Update the solutions
this->get_solution_handler().update(this->get_field_solve_type(),
this->get_solve_block(),
index);
if (variable.get_pde_type() != PDEType::Auxiliary)
{
// The solve will have updated the "old solution" vector with the newton update
// so it's technically the new solution. In order to update the solutions and
// preserve the old states we copy the old solution to the new solution and
// update.
*(this->get_solution_handler().get_new_solution_vector(index)) =
*(this->get_solution_handler().get_solution_vector(index,
DependencyType::Normal));
this->get_solution_handler().update(this->get_field_solve_type(),
this->get_solve_block(),
index);
}

// Update the ghosts
Timer::start_section("Update ghosts");
Expand Down
5 changes: 5 additions & 0 deletions src/solvers/sequential_solver.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -393,6 +393,11 @@ SequentialSolver<dim, degree, number>::solve_linear_solver(
// Grab the global field index
const Types::Index global_field_index = variable.get_field_index();

// Zero out the ghosts
Timer::start_section("Zero ghosts");
this->get_solution_handler().zero_out_ghosts();
Timer::end_section("Zero ghosts");
Copy link
Copy Markdown
Contributor Author

@jijn jijn Feb 14, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

should we do it here or in sequential_co_nonlinear_solver.cc|solve()?


// Skip if the field type is ImplicitTimeDependent and the increment is 0.
if (variable.get_pde_type() == PDEType::ImplicitTimeDependent &&
this->get_user_inputs().get_temporal_discretization().get_increment() == 0)
Expand Down
4 changes: 2 additions & 2 deletions tests/regression_tests/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,10 +11,10 @@ set(
allen_cahn_implicit
cahn_hilliard_explicit
cahn_hilliard_implicit
cavity_flow
# cavity_flow
fracture
heat_equation_steady_state
linear_solve_old_solution
# linear_solve_old_solution
mechanics
precipitate_explicit
solution_blocks
Expand Down
76 changes: 38 additions & 38 deletions tests/regression_tests/allen_cahn_implicit/gold_output.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,58 +32,58 @@ Iteration: 0


Iteration: 300
Solution index 0 l2-norm: 109.817 integrated value: 1911.34
Solution index 1 l2-norm: 9.88843 integrated value: 94.0071
Solution index 2 l2-norm: 4.16914 integrated value: 32.4592

Iteration: 600
Solution index 0 l2-norm: 111.095 integrated value: 1942.03
Solution index 1 l2-norm: 9.09367 integrated value: 81.6845
Solution index 2 l2-norm: 3.82926 integrated value: 28.0514

Iteration: 900
Solution index 0 l2-norm: 112.331 integrated value: 1969.05
Solution index 1 l2-norm: 8.72713 integrated value: 72.3072
Solution index 2 l2-norm: 3.62388 integrated value: 24.0953

Iteration: 1200
Iteration: 600
Solution index 0 l2-norm: 112.913 integrated value: 1982.79
Solution index 1 l2-norm: 8.61869 integrated value: 69.1128
Solution index 2 l2-norm: 3.57043 integrated value: 23.0561

Iteration: 900
Solution index 0 l2-norm: 113.642 integrated value: 2003.2
Solution index 1 l2-norm: 8.41 integrated value: 65.3784
Solution index 2 l2-norm: 3.48117 integrated value: 21.7874

Iteration: 1200
Solution index 0 l2-norm: 114.299 integrated value: 2022.76
Solution index 1 l2-norm: 8.22981 integrated value: 62.4847
Solution index 2 l2-norm: 3.40591 integrated value: 20.8132

Iteration: 1500
Solution index 0 l2-norm: 113.295 integrated value: 1993.24
Solution index 1 l2-norm: 8.5108 integrated value: 67.0845
Solution index 2 l2-norm: 3.52369 integrated value: 22.366
Solution index 0 l2-norm: 114.93 integrated value: 2042.11
Solution index 1 l2-norm: 8.07034 integrated value: 60.019
Solution index 2 l2-norm: 3.3396 integrated value: 19.987



+-----------------------------------------------+------------+------------+
| Total wallclock time elapsed since start | 14.7s | |
| Total wallclock time elapsed since start | 23.3s | |
| | | |
| Section | no. calls | wall time | % of total |
+-----------------------------------+-----------+------------+------------+
| Auxiliary solver | 1 | 6.85e-06s | 0% |
| Create FESystem | 1 | 0.00479s | 0% |
| Create constraints | 1 | 0.00476s | 0% |
| Explicit solver | 1500 | 0.00236s | 0% |
| Generate mesh | 1 | 0.0332s | 0.23% |
| Initialization | 1 | 0.355s | 2.4% |
| Nonexplicit auxiliary solver | 1500 | 0.00178s | 0% |
| Nonexplicit co-nonlinear solver | 1501 | 0.00262s | 0% |
| Nonexplicit linear solver | 1501 | 0.0017s | 0% |
| Nonexplicit self-nonlinear solver | 1501 | 13.3s | 90% |
| Output | 6 | 0.402s | 2.7% |
| Postprocess solver | 6 | 0.0296s | 0.2% |
| Solve Increment | 1500 | 13.3s | 90% |
| Solver initialization | 1 | 0.0752s | 0.51% |
| Update ghosts | 4212 | 0.00982s | 0% |
| Update time-dependent constraints | 1500 | 0.00217s | 0% |
| Zero ghosts | 6 | 3.04e-05s | 0% |
| compute element volumes | 1 | 0.00606s | 0% |
| reinitialize DoFHandlers | 1 | 0.00749s | 0% |
| reinitialize element volumes | 1 | 0.000799s | 0% |
| reinitialize invm | 1 | 0.0148s | 0.1% |
| reinitialize matrix-free objects | 1 | 0.054s | 0.37% |
| reinitialize solution set | 1 | 0.000382s | 0% |
| Auxiliary solver | 1 | 3.81e-06s | 0% |
| Create FESystem | 1 | 0.00868s | 0% |
| Create constraints | 1 | 0.011s | 0% |
| Explicit solver | 1500 | 0.00251s | 0% |
| Generate mesh | 1 | 0.0573s | 0.25% |
| Initialization | 1 | 0.394s | 1.7% |
| Nonexplicit auxiliary solver | 1500 | 0.00192s | 0% |
| Nonexplicit co-nonlinear solver | 1501 | 0.00337s | 0% |
| Nonexplicit linear solver | 1501 | 0.00163s | 0% |
| Nonexplicit self-nonlinear solver | 1501 | 21.8s | 94% |
| Output | 6 | 0.364s | 1.6% |
| Postprocess solver | 6 | 0.0257s | 0.11% |
| Solve Increment | 1500 | 21.8s | 94% |
| Solver initialization | 1 | 0.0671s | 0.29% |
| Update ghosts | 4962 | 0.0129s | 0% |
| Update time-dependent constraints | 1500 | 0.00208s | 0% |
| Zero ghosts | 6 | 3.45e-05s | 0% |
| compute element volumes | 1 | 0.00189s | 0% |
| reinitialize DoFHandlers | 1 | 0.0172s | 0% |
| reinitialize element volumes | 1 | 0.000857s | 0% |
| reinitialize invm | 1 | 0.0159s | 0% |
| reinitialize matrix-free objects | 1 | 0.0608s | 0.26% |
| reinitialize solution set | 1 | 0.000637s | 0% |
+-----------------------------------+-----------+------------+------------+
82 changes: 41 additions & 41 deletions tests/regression_tests/cahn_hilliard_implicit/gold_output.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
Using the input parameter file: parameters.prm
Using the input parameter file: ../parameters.prm
Number of constants: 2
Number of variables: 3
number of processes: 1
Expand Down Expand Up @@ -27,63 +27,63 @@ solving 0th timestep for solve block 0 ...
outputting initial condition...
Iteration: 0
Solution index 0 l2-norm: 109.458 integrated value: 1886.53
Solution index 1 l2-norm: 0 integrated value: 0
Solution index 1 l2-norm: 22.013 integrated value: 1.86144e-14
Solution index 2 l2-norm: 4.81343 integrated value: 34.5905


Iteration: 100
Solution index 0 l2-norm: 109.451 integrated value: 1886.53
Solution index 1 l2-norm: 12.6513 integrated value: -9.27182e-15
Solution index 2 l2-norm: 4.66535 integrated value: 33.852

Iteration: 200
Solution index 0 l2-norm: 109.437 integrated value: 1886.52
Solution index 1 l2-norm: 12.1533 integrated value: 3.2889e-14
Solution index 1 l2-norm: 12.1533 integrated value: -8.234e-15
Solution index 2 l2-norm: 4.60564 integrated value: 33.6349

Iteration: 200
Solution index 0 l2-norm: 109.413 integrated value: 1886.52
Solution index 1 l2-norm: 11.6105 integrated value: -2.08338e-14
Solution index 2 l2-norm: 4.53104 integrated value: 33.3667

Iteration: 300
Solution index 0 l2-norm: 109.424 integrated value: 1886.52
Solution index 1 l2-norm: 11.8406 integrated value: 4.23915e-14
Solution index 2 l2-norm: 4.56374 integrated value: 33.4845
Solution index 0 l2-norm: 109.392 integrated value: 1886.52
Solution index 1 l2-norm: 11.2774 integrated value: 7.01883e-14
Solution index 2 l2-norm: 4.48123 integrated value: 33.1852

Iteration: 400
Solution index 0 l2-norm: 109.413 integrated value: 1886.52
Solution index 1 l2-norm: 11.6105 integrated value: 2.74139e-14
Solution index 2 l2-norm: 4.53104 integrated value: 33.3667
Solution index 0 l2-norm: 109.373 integrated value: 1886.52
Solution index 1 l2-norm: 11.0371 integrated value: -5.2234e-15
Solution index 2 l2-norm: 4.44364 integrated value: 33.0455

Iteration: 500
Solution index 0 l2-norm: 109.402 integrated value: 1886.52
Solution index 1 l2-norm: 11.4283 integrated value: -4.50475e-14
Solution index 2 l2-norm: 4.50413 integrated value: 33.2691
Solution index 0 l2-norm: 109.356 integrated value: 1886.52
Solution index 1 l2-norm: 10.8503 integrated value: -7.42764e-14
Solution index 2 l2-norm: 4.41351 integrated value: 32.9311



+-----------------------------------------------+------------+------------+
| Total wallclock time elapsed since start | 29.2s | |
| Total wallclock time elapsed since start | 50.1s | |
| | | |
| Section | no. calls | wall time | % of total |
+-----------------------------------+-----------+------------+------------+
| Auxiliary solver | 1 | 1.61e-06s | 0% |
| Create FESystem | 1 | 0.00522s | 0% |
| Create constraints | 1 | 0.00357s | 0% |
| Explicit solver | 500 | 0.000771s | 0% |
| Generate mesh | 1 | 0.0133s | 0% |
| Initialization | 1 | 0.231s | 0.79% |
| Nonexplicit auxiliary solver | 500 | 0.00049s | 0% |
| Nonexplicit co-nonlinear solver | 501 | 28.5s | 98% |
| Nonexplicit linear solver | 501 | 0.000471s | 0% |
| Nonexplicit self-nonlinear solver | 501 | 0.000427s | 0% |
| Output | 6 | 0.217s | 0.74% |
| Postprocess solver | 6 | 0.00267s | 0% |
| Solve Increment | 500 | 28.5s | 98% |
| Solver initialization | 1 | 0.0326s | 0.11% |
| Update ghosts | 33621 | 0.0456s | 0.16% |
| Update time-dependent constraints | 500 | 0.000501s | 0% |
| Zero ghosts | 16563 | 0.0421s | 0.14% |
| compute element volumes | 1 | 0.000823s | 0% |
| reinitialize DoFHandlers | 1 | 0.00987s | 0% |
| reinitialize element volumes | 1 | 0.000534s | 0% |
| reinitialize invm | 1 | 0.0113s | 0% |
| reinitialize matrix-free objects | 1 | 0.0356s | 0.12% |
| reinitialize solution set | 1 | 0.00293s | 0% |
| Auxiliary solver | 1 | 2.44e-06s | 0% |
| Create FESystem | 1 | 0.0049s | 0% |
| Create constraints | 1 | 0.00378s | 0% |
| Explicit solver | 500 | 0.000602s | 0% |
| Generate mesh | 1 | 0.0131s | 0% |
| Initialization | 1 | 0.224s | 0.45% |
| Nonexplicit auxiliary solver | 500 | 0.000514s | 0% |
| Nonexplicit co-nonlinear solver | 501 | 49.4s | 99% |
| Nonexplicit linear solver | 501 | 0.000498s | 0% |
| Nonexplicit self-nonlinear solver | 501 | 0.000431s | 0% |
| Output | 6 | 0.243s | 0.49% |
| Postprocess solver | 6 | 0.00266s | 0% |
| Solve Increment | 500 | 49.4s | 99% |
| Solver initialization | 1 | 0.0264s | 0% |
| Update ghosts | 62169 | 0.0763s | 0.15% |
| Update time-dependent constraints | 500 | 0.000502s | 0% |
| Zero ghosts | 61668 | 0.158s | 0.32% |
| compute element volumes | 1 | 0.00259s | 0% |
| reinitialize DoFHandlers | 1 | 0.00983s | 0% |
| reinitialize element volumes | 1 | 0.000706s | 0% |
| reinitialize invm | 1 | 0.00965s | 0% |
| reinitialize matrix-free objects | 1 | 0.038s | 0% |
| reinitialize solution set | 1 | 0.00306s | 0% |
+-----------------------------------+-----------+------------+------------+
1 change: 1 addition & 0 deletions tests/regression_tests/fracture/ICs_and_BCs.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,7 @@
#include <prismspf/config.h>

#include <cmath>
#include <numbers>

PRISMS_PF_BEGIN_NAMESPACE

Expand Down