Backward Euler solver improvements #3009

bendudson · 2024-10-24T06:17:55Z

Replaces #2631

Running 1D-recycling, the first 200 output timesteps to t=1e6. BOUT++ compiled with PETSc 3.22.0, using one MPI process. Using PETSc ILU preconditioner.

Baseline: next branch: 104,905 calls. 8 minutes 23 seconds

Improvements (from most to least successful):

Use matrix free Jacobian-vector operator, keep using finite difference Jacobian to calculate the preconditoner: 23,206 calls. 2 minutes 10 seconds. This option matrix_free_operator works so well that it is set to true by default.
Scaling the time derivatives (Jacobian row scaling). By itself it improves performance over baseline (40,876 calls, 3 min 42 secs) but with (1) it reduces performance slightly (26,786 calls, 2 min 35 secs).
Prune the Jacobian by removing elements with small values. Intended to reduce work needed to form the Jacobian and calculate the preconditioner. Tried various strategies, not successful yet: With (1) takes 162,249 calls, 11 mins 35 secs. Depending on strategy, either the preconditioner is significantly worse or the Jacobian needs to be recalculated more often.
Scale variables by their values. Seems to work poorly, perhaps due to bad choice of normalisation or implementation bug. Reduces performance so the solver either failed or I lost patience before finishing.

Scales the RHS (time derivatives) in the Jacobian calculation. Speeds up the 1D-recycling example by about x2.

Uses the magnitude of each variable to update its normalisation. Doesn't improve convergence for 1D-recycling test case (makes worse)

When using coloring, identify near-zero values and remove them from the Jacobian. This could improve efficiency by reducing the effort needed to recompute the Jacobian and preconditioner. Currently over-prunes based on initial conditions, resulting in a worse preconditioner. Need to occasionally reset and add back in the removed elements, maybe following a convergence failure.

Works well, true by default. Configures SNES so that a matrix-free approximation is used for the Jacobian-vector operator. The finite difference Jacobian is used to construct the preconditioner.

github-actions

clang-tidy made some suggestions

There were too many comments to post at once. Showing the first 25 out of 49. Check the log or trigger a new build to see more.

src/solver/impls/snes/snes.cxx

github-actions · 2024-10-24T06:21:10Z

src/solver/impls/snes/snes.cxx

+  // Get the the SNESSolver pointer from the function call context
+  SNESSolver* fctx;
+  err = MatFDColoringGetFunction(static_cast<MatFDColoring>(ctx), nullptr,
+                                 reinterpret_cast<void**>(&fctx));


warning: do not use reinterpret_cast [cppcoreguidelines-pro-type-reinterpret-cast]

reinterpret_cast<void**>(&fctx)); ^

I believe this might be required because PETSc uses the wrong type (or that might be a similar but unrelated situation).

Two options:

on line above: // NOLINTNEXTLINE(cppcoreguidelines-pro-type-reinterpret-cast)

if we need to use this function more than once, make a wrapper with the correct type and the NOLINT comment inside?

github-actions · 2024-10-24T06:21:12Z

src/solver/impls/snes/snes.cxx


 int SNESSolver::run() {
  TRACE("SNESSolver::run()");
+  int ierr;


warning: variable 'ierr' is not initialized [cppcoreguidelines-init-variables]

Suggested change

int ierr;

int ierr = 0;

github-actions · 2024-10-24T06:21:12Z

src/solver/impls/snes/snes.cxx

+        if (loop_count % 100 == 0) {
+          // Rescale state (snes_x) so that all quantities are around 1
+          // If quantities are near zero then RTOL is used
+          int istart, iend;


warning: variable 'istart' is not initialized [cppcoreguidelines-init-variables]

int istart, iend; ^

this fix will not be applied because it overlaps with another fix

github-actions · 2024-10-24T06:21:12Z

src/solver/impls/snes/snes.cxx

+        if (loop_count % 100 == 0) {
+          // Rescale state (snes_x) so that all quantities are around 1
+          // If quantities are near zero then RTOL is used
+          int istart, iend;


warning: variable 'iend' is not initialized [cppcoreguidelines-init-variables]

int istart, iend; ^

this fix will not be applied because it overlaps with another fix

src/solver/impls/snes/snes.cxx

github-actions · 2024-10-24T06:21:12Z

src/solver/impls/snes/snes.cxx

+          ierr = VecGetArray(snes_x, &snes_x_data);
+          CHKERRQ(ierr);
+          PetscScalar* x1_data;
+          ierr = VecGetArray(x1, &x1_data);


warning: variable 'var_scaling_factors_data' is not initialized [cppcoreguidelines-init-variables]

Suggested change

ierr = VecGetArray(x1, &x1_data);

PetscScalar* var_scaling_factors_data = nullptr;

ZedThree

LGTM, would be good to address some of the clang-tidy comments

src/solver/impls/snes/snes.cxx

ZedThree · 2024-10-24T08:15:59Z

src/solver/impls/snes/snes.cxx

+  // Get the the SNESSolver pointer from the function call context
+  SNESSolver* fctx;
+  err = MatFDColoringGetFunction(static_cast<MatFDColoring>(ctx), nullptr,
+                                 reinterpret_cast<void**>(&fctx));


I believe this might be required because PETSc uses the wrong type (or that might be a similar but unrelated situation).

Two options:

on line above: // NOLINTNEXTLINE(cppcoreguidelines-pro-type-reinterpret-cast)

if we need to use this function more than once, make a wrapper with the correct type and the NOLINT comment inside?

src/solver/impls/snes/snes.cxx

src/solver/impls/snes/snes.hxx

src/solver/impls/snes/snes.cxx

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> Co-authored-by: Peter Hill <peter.hill@york.ac.uk>

github-actions

clang-tidy made some suggestions

There were too many comments to post at once. Showing the first 25 out of 32. Check the log or trigger a new build to see more.

github-actions · 2024-10-24T15:47:55Z

src/solver/impls/snes/snes.cxx

+          PetscScalar* snes_x_data = nullptr;
+          ierr = VecGetArray(snes_x, &snes_x_data);
+          CHKERRQ(ierr);
+          PetscScalar* x1_data;


warning: variable 'x1_data' is not initialized [cppcoreguidelines-init-variables]

Suggested change

PetscScalar* x1_data;

PetscScalar* x1_data = nullptr;

github-actions · 2024-10-24T15:47:56Z

src/solver/impls/snes/snes.cxx

+          PetscScalar* x1_data;
+          ierr = VecGetArray(x1, &x1_data);
+          CHKERRQ(ierr);
+          PetscScalar* var_scaling_factors_data;


warning: variable 'var_scaling_factors_data' is not initialized [cppcoreguidelines-init-variables]

Suggested change

PetscScalar* var_scaling_factors_data;

PetscScalar* var_scaling_factors_data = nullptr;

github-actions · 2024-10-24T15:47:56Z

src/solver/impls/snes/snes.cxx

+
+          if (diagnose) {
+            // Print maximum and minimum scaling factors
+            PetscReal max_scale, min_scale;


warning: multiple declarations in a single statement reduces readability [readability-isolate-declaration]

Suggested change

PetscReal max_scale, min_scale;

PetscReal max_scale;

PetscReal min_scale;

github-actions · 2024-10-24T15:47:56Z

src/solver/impls/snes/snes.cxx

+
+          if (diagnose) {
+            // Print maximum and minimum scaling factors
+            PetscReal max_scale, min_scale;


warning: variable 'max_scale' is not initialized [cppcoreguidelines-init-variables]

PetscReal max_scale, min_scale; ^

this fix will not be applied because it overlaps with another fix

github-actions · 2024-10-24T15:47:56Z

src/solver/impls/snes/snes.cxx

+
+          if (diagnose) {
+            // Print maximum and minimum scaling factors
+            PetscReal max_scale, min_scale;


warning: variable 'min_scale' is not initialized [cppcoreguidelines-init-variables]

PetscReal max_scale, min_scale; ^

this fix will not be applied because it overlaps with another fix

github-actions · 2024-10-24T15:47:56Z

src/solver/impls/snes/snes.cxx

+      const BoutReal* xdata = nullptr;
+      ierr = VecGetArrayRead(scaled_x, &xdata);
+      CHKERRQ(ierr);
+      load_vars(const_cast<BoutReal*>(xdata));


warning: do not use const_cast [cppcoreguidelines-pro-type-const-cast]

; ^

github-actions · 2024-10-24T15:47:56Z

src/solver/impls/snes/snes.cxx

+    const BoutReal* xdata = nullptr;
+    ierr = VecGetArrayRead(scaled_x, &xdata);
+    CHKERRQ(ierr);
+    load_vars(const_cast<BoutReal*>(xdata));


warning: do not use const_cast [cppcoreguidelines-pro-type-const-cast]

; ^

github-actions · 2024-10-24T15:47:56Z

src/solver/impls/snes/snes.cxx

+    const BoutReal* xdata = nullptr;
+    int ierr = VecGetArrayRead(x, &xdata);
+    CHKERRQ(ierr);
+    load_vars(const_cast<BoutReal*>(xdata));


warning: do not use const_cast [cppcoreguidelines-pro-type-const-cast]

; ^

github-actions · 2024-10-24T15:47:57Z

src/solver/impls/snes/snes.cxx

+  int rstart, rend;
+  MatGetOwnershipRange(B, &rstart, &rend);


warning: multiple declarations in a single statement reduces readability [readability-isolate-declaration]

Suggested change

int rstart, rend;

MatGetOwnershipRange(B, &rstart, &rend);

essor

reint rstart;

int rend;

github-actions · 2024-10-24T15:47:57Z

src/solver/impls/snes/snes.cxx

+
+  // Get index of rows owned by this processor
+  int rstart, rend;
+  MatGetOwnershipRange(B, &rstart, &rend);


warning: variable 'rend' is not initialized [cppcoreguidelines-init-variables]

rend; ^

this fix will not be applied because it overlaps with another fix

- Making more indices const - Renaming Jacobian scaling function, removing SNES prefix to avoid potential collision with PETSc, and making function static.

github-actions

clang-tidy made some suggestions

There were too many comments to post at once. Showing the first 25 out of 27. Check the log or trigger a new build to see more.

github-actions · 2024-10-25T21:09:16Z

src/solver/impls/snes/snes.cxx

+    const BoutReal* xdata = nullptr;
+    ierr = VecGetArrayRead(scaled_x, &xdata);
+    CHKERRQ(ierr);
+    load_vars(const_cast<BoutReal*>(


warning: do not use const_cast [cppcoreguidelines-pro-type-const-cast]

; ^

github-actions · 2024-10-25T21:09:16Z

src/solver/impls/snes/snes.cxx

+  } else {
+    const BoutReal* xdata = nullptr;
+    int ierr = VecGetArrayRead(x, &xdata);
+    CHKERRQ(ierr);


warning: do not use const_cast [cppcoreguidelines-pro-type-const-cast]

; ^

ZedThree · 2024-10-30T10:54:46Z

I've deleted a lot of the clang-tidy comments because they were not on the right lines and were very confusing

PR #3009 included changes to allow matrix-free Jacobian-vector products, while using the finite difference Jacobian for the preconditioner. That gave significant speedups that will hopefully also help now that the coloring stencil is fixed.

bendudson and others added 5 commits October 23, 2024 11:46

beuler solver: scale_rhs option

1126b28

Scales the RHS (time derivatives) in the Jacobian calculation. Speeds up the 1D-recycling example by about x2.

beuler: Add variable scaling (Jacobian column)

4220175

Uses the magnitude of each variable to update its normalisation. Doesn't improve convergence for 1D-recycling test case (makes worse)

beuler: Add matrix_free_operator flag

1df2f34

Works well, true by default. Configures SNES so that a matrix-free approximation is used for the Jacobian-vector operator. The finite difference Jacobian is used to construct the preconditioner.

Apply clang-format changes

393bbe5

github-actions bot reviewed Oct 24, 2024

View reviewed changes

bendudson mentioned this pull request Oct 24, 2024

Improving backward Euler solver #2631

Closed

ZedThree previously approved these changes Oct 24, 2024

View reviewed changes

Apply suggestions from code review

07a1b60

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> Co-authored-by: Peter Hill <peter.hill@york.ac.uk>

bendudson dismissed ZedThree’s stale review via 07a1b60 October 24, 2024 15:44

Apply clang-format changes

adf2311

github-actions bot reviewed Oct 24, 2024

View reviewed changes

bendudson and others added 2 commits October 25, 2024 14:04

beuler solver: Implement tidy suggestions

dec8ae7

- Making more indices const - Renaming Jacobian scaling function, removing SNES prefix to avoid potential collision with PETSc, and making function static.

Apply clang-format changes

89c32b7

github-actions bot reviewed Oct 25, 2024

View reviewed changes

boutproject deleted a comment from github-actions bot Oct 30, 2024

bendudson merged commit 9d8c036 into next Nov 5, 2024

bendudson mentioned this pull request Apr 24, 2025

Feature/snes stencil #3104

Merged

	ierr = VecGetArray(x1, &x1_data);
	PetscScalar* var_scaling_factors_data = nullptr;

	PetscScalar* var_scaling_factors_data;
	PetscScalar* var_scaling_factors_data = nullptr;

	PetscReal max_scale, min_scale;
	PetscReal max_scale;
	PetscReal min_scale;

Backward Euler solver improvements #3009

Backward Euler solver improvements #3009

Uh oh!

Conversation

bendudson commented Oct 24, 2024

Uh oh!

github-actions bot left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

github-actions bot Oct 24, 2024

Choose a reason for hiding this comment

Uh oh!

ZedThree Oct 24, 2024

Choose a reason for hiding this comment

Uh oh!

github-actions bot Oct 24, 2024

Choose a reason for hiding this comment

Uh oh!

github-actions bot Oct 24, 2024

Choose a reason for hiding this comment

Uh oh!

github-actions bot Oct 24, 2024

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

github-actions bot Oct 24, 2024

Choose a reason for hiding this comment

Uh oh!

ZedThree left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

ZedThree Oct 24, 2024

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

github-actions bot left a comment

Choose a reason for hiding this comment

Uh oh!

github-actions bot Oct 24, 2024

Choose a reason for hiding this comment

Uh oh!

github-actions bot Oct 24, 2024

Choose a reason for hiding this comment

Uh oh!

github-actions bot Oct 24, 2024

Choose a reason for hiding this comment

Uh oh!

github-actions bot Oct 24, 2024

Choose a reason for hiding this comment

Uh oh!

github-actions bot Oct 24, 2024

Choose a reason for hiding this comment

Uh oh!

github-actions bot Oct 24, 2024

Choose a reason for hiding this comment

Uh oh!

github-actions bot Oct 24, 2024

Choose a reason for hiding this comment

Uh oh!

github-actions bot Oct 24, 2024

Choose a reason for hiding this comment

Uh oh!

github-actions bot Oct 24, 2024

Choose a reason for hiding this comment

Uh oh!

github-actions bot Oct 24, 2024

Choose a reason for hiding this comment

Uh oh!

github-actions bot left a comment

Choose a reason for hiding this comment

Uh oh!