Add matrix coloring to beuler solver (and fix imexbdf2) #2454

bendudson · 2021-10-15T20:17:38Z

Enables large reductions in the number of iterations needed to calculate the Jacobian elements. Having the Jacobian then enables good preconditioning.

Simple conduction example, nout=5 and timestep=10:

PVODE:

1.000e+01        644       1.25e-01    74.9    0.0    0.7   13.4   10.9
2.000e+01        161       4.41e-02    50.9    0.0    0.4   29.4   19.2

beuler, SNESType = anderson

1.000e+01        962       2.13e-01    65.8    0.0    0.7    8.2   25.4
2.000e+01        347       8.20e-02    61.3    0.0    0.6   14.0   24.0

beuler, SNESType = newtonls, matrix free (no Jacobian)

1.000e+01       3389       5.68e-01    89.2    0.0    0.8    3.3    6.6
2.000e+01       2325       3.55e-01    89.1    0.0    0.8    3.3    6.8

beuler, SNESType = newtonls, with Jacobian, no coloring

1.000e+01        833       1.53e-01    78.8    0.0    0.7   12.1    8.4
2.000e+01        521       9.81e-02    76.6    0.0    0.8   12.3   10.4

beuler, SNESType = newtonls, with Jacobian coloring

1.000e+01         85       3.78e-02    33.3    0.0    0.3   46.5   19.9
2.000e+01         57       3.16e-02    26.4    0.0    0.2   48.9   24.6

Also includes some reformatting and drive-by fixes for things flagged by Apple Clang. These are in separate commits (mostly) for reviewing.

Solvers can pass a flag to say whether the function can be linearised (e.g. inside inner linear solve, or Jacobian calculation). Physics models don't have to use the argument: rhs functions with a single (time) argument continue to work.

Should use const reference to avoid copying index

Deprecated, flagged as warning by Apple Clang

Enables large reductions in the number of iterations needed to calculate the Jacobian elements. Having the Jacobian then enables good preconditioning. Simple conduction example, nout=5 and timestep=10: PVODE: ``` 1.000e+01 644 1.25e-01 74.9 0.0 0.7 13.4 10.9 2.000e+01 161 4.41e-02 50.9 0.0 0.4 29.4 19.2 ``` beuler, SNESType = anderson ``` 1.000e+01 962 2.13e-01 65.8 0.0 0.7 8.2 25.4 2.000e+01 347 8.20e-02 61.3 0.0 0.6 14.0 24.0 ``` beuler, SNESType = newtonls, matrix free (no Jacobian) ``` 1.000e+01 3389 5.68e-01 89.2 0.0 0.8 3.3 6.6 2.000e+01 2325 3.55e-01 89.1 0.0 0.8 3.3 6.8 ``` beuler, SNESType = newtonls, with Jacobian, no coloring ``` 1.000e+01 833 1.53e-01 78.8 0.0 0.7 12.1 8.4 2.000e+01 521 9.81e-02 76.6 0.0 0.8 12.3 10.4 ``` beuler, SNESType = newtonls, with Jacobian coloring ``` 1.000e+01 85 3.78e-02 33.3 0.0 0.3 46.5 19.9 2.000e+01 57 3.16e-02 26.4 0.0 0.2 48.9 24.6 ```

Was destroying the iscoloring object before setup, and then reverting to brute-force calculation of the Jacobian.

github-actions

clang-tidy made some suggestions

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

github-actions · 2021-10-15T20:25:48Z

src/solver/impls/imex-bdf2/imex-bdf2.cxx


 // Redundent definition because < C++17
 constexpr int IMEXBDF2::MAX_SUPPORTED_ORDER;

-IMEXBDF2::IMEXBDF2(Options *opt)
+IMEXBDF2::IMEXBDF2(Options* opt)


warning: constructor does not initialize these fields: maxOrder, out_timestep, nsteps, timestep, ninternal, mxstep, adaptive, nadapt, mxstepAdapt, scaleCushUp, scaleCushDown, adaptRtol, dtMin, dtMax, dtMinFatal, dtImp, nlocal, neq, implicit_gamma, implicit_curtime, predictor, diagnose, verbose, linear_fails, nonlinear_fails, have_constraints, fdcoloring [cppcoreguidelines-pro-type-member-init]

IMEXBDF2::IMEXBDF2(Options* opt) ^

src/solver/impls/imex-bdf2/imex-bdf2.cxx

github-actions · 2021-10-15T20:25:50Z

src/solver/impls/imex-bdf2/imex-bdf2.cxx

-                o_nnz[localIndex+i] += (n3d + n2d);
+              for (int i = 0; i < n3d; i++) {
+                d_nnz[localIndex + i] -= (n3d + n2d);
+                o_nnz[localIndex + i] += (n3d + n2d);


warning: do not use pointer arithmetic [cppcoreguidelines-pro-bounds-pointer-arithmetic]

o_nnz[localIndex + i] += (n3d + n2d); ^

github-actions · 2021-10-15T20:25:50Z

src/solver/impls/imex-bdf2/imex-bdf2.cxx

-          o_nnz[localIndex+i] += (n3d + n2d);
+        for (int i = 0; i < n2d + n3d; i++) {
+          // d_nnz[localIndex+i] -= (n3d + n2d);
+          o_nnz[localIndex + i] += (n3d + n2d);


warning: do not use pointer arithmetic [cppcoreguidelines-pro-bounds-pointer-arithmetic]

o_nnz[localIndex + i] += (n3d + n2d); ^

github-actions · 2021-10-15T20:25:50Z

src/solver/impls/imex-bdf2/imex-bdf2.cxx

-            o_nnz[localIndex+i] += (n3d + n2d);
+          for (int i = 0; i < n3d; i++) {
+            // d_nnz[localIndex+i] -= (n3d + n2d);
+            o_nnz[localIndex + i] += (n3d + n2d);


warning: do not use pointer arithmetic [cppcoreguidelines-pro-bounds-pointer-arithmetic]

o_nnz[localIndex + i] += (n3d + n2d); ^

github-actions · 2021-10-15T20:25:50Z

src/solver/impls/imex-bdf2/imex-bdf2.cxx

-          o_nnz[localIndex+i] += (n3d + n2d);
+        for (int i = 0; i < n2d + n3d; i++) {
+          // d_nnz[localIndex+i] -= (n3d + n2d);
+          o_nnz[localIndex + i] += (n3d + n2d);


warning: do not use pointer arithmetic [cppcoreguidelines-pro-bounds-pointer-arithmetic]

o_nnz[localIndex + i] += (n3d + n2d); ^

github-actions · 2021-10-15T20:25:51Z

src/solver/impls/imex-bdf2/imex-bdf2.cxx

-            o_nnz[localIndex+i] += (n3d + n2d);
+          for (int i = 0; i < n3d; i++) {
+            // d_nnz[localIndex+i] -= (n3d + n2d);
+            o_nnz[localIndex + i] += (n3d + n2d);


warning: do not use pointer arithmetic [cppcoreguidelines-pro-bounds-pointer-arithmetic]

o_nnz[localIndex + i] += (n3d + n2d); ^

Describe all the options, and how to use newtonls with Jacobian coloring.

[skip ci]

- Defaults to using coloring to calculate a Jacobian - Iterations and lag jacobian set to values that work for SD1D tests so far. - Diagnostic outputs now include number of linear iterations and number of failures.

Reflect changes in defaults, and some notes on when these are likely to fail.

github-actions

clang-tidy made some suggestions

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

github-actions · 2021-10-19T18:41:52Z

src/solver/impls/imex-bdf2/imex-bdf2.cxx

-        for(int i=0;i<n2d+n3d;i++) {
-          o_nnz[localIndex+i] -= (n3d + n2d);
+        for (int i = 0; i < n2d + n3d; i++) {
+          o_nnz[localIndex + i] -= (n3d + n2d);


warning: do not use pointer arithmetic [cppcoreguidelines-pro-bounds-pointer-arithmetic]

o_nnz[localIndex + i] -= (n3d + n2d); ^

github-actions · 2021-10-19T18:41:52Z

src/solver/impls/imex-bdf2/imex-bdf2.cxx

-          for(int i=0;i<n3d;i++) {
-            o_nnz[localIndex+i] -= (n3d + n2d);
+          for (int i = 0; i < n3d; i++) {
+            o_nnz[localIndex + i] -= (n3d + n2d);


warning: do not use pointer arithmetic [cppcoreguidelines-pro-bounds-pointer-arithmetic]

o_nnz[localIndex + i] -= (n3d + n2d); ^

github-actions · 2021-10-19T18:41:52Z

src/solver/impls/imex-bdf2/imex-bdf2.cxx

-        for(int i=0;i<n2d+n3d;i++) {
-          o_nnz[localIndex+i] -= (n3d + n2d);
+        for (int i = 0; i < n2d + n3d; i++) {
+          o_nnz[localIndex + i] -= (n3d + n2d);


warning: do not use pointer arithmetic [cppcoreguidelines-pro-bounds-pointer-arithmetic]

o_nnz[localIndex + i] -= (n3d + n2d); ^

github-actions · 2021-10-19T18:41:53Z

src/solver/impls/imex-bdf2/imex-bdf2.cxx

-          for(int i=0;i<n3d;i++) {
-            o_nnz[localIndex+i] -= (n3d + n2d);
+          for (int i = 0; i < n3d; i++) {
+            o_nnz[localIndex + i] -= (n3d + n2d);


warning: do not use pointer arithmetic [cppcoreguidelines-pro-bounds-pointer-arithmetic]

o_nnz[localIndex + i] -= (n3d + n2d); ^

github-actions · 2021-10-19T18:41:53Z

src/solver/impls/imex-bdf2/imex-bdf2.cxx

-      const int xoffset[5] = {0,-1, 1, 0, 0};
-      const int yoffset[5] = {0, 0, 0,-1, 1};
-
+      const int xoffset[5] = {0, -1, 1, 0, 0};


warning: do not declare C-style arrays, use std::array<> instead [cppcoreguidelines-avoid-c-arrays]

const int xoffset[5] = {0, -1, 1, 0, 0}; ^

Not touched in this PR, but this is a good suggestion we should do at some point

github-actions · 2021-10-19T18:41:54Z

src/solver/impls/imex-bdf2/imex-bdf2.cxx

          // An error occurred. If adaptive, reduce timestep
-          if(!adaptive)
+          if (!adaptive)


warning: statement should be inside braces [readability-braces-around-statements]

if (!adaptive) ^ {

github-actions · 2021-10-19T18:41:54Z

src/solver/impls/imex-bdf2/imex-bdf2.cxx

+        // Now we can calculate the error and decide what we want to do
+        if (checkingErr) {
+          // Now we want to find the actual (abs) error
+          BoutReal errTot[3] = {0, 0, 0};


warning: do not declare C-style arrays, use std::array<> instead [cppcoreguidelines-avoid-c-arrays]

BoutReal errTot[3] = {0, 0, 0}; ^

github-actions · 2021-10-19T18:41:55Z

src/solver/impls/imex-bdf2/imex-bdf2.cxx

+        if (checkingErr) {
+          // Now we want to find the actual (abs) error
+          BoutReal errTot[3] = {0, 0, 0};
+          BoutReal errGlobTot[3] = {0, 0, 0};


warning: do not declare C-style arrays, use std::array<> instead [cppcoreguidelines-avoid-c-arrays]

BoutReal errGlobTot[3] = {0, 0, 0}; ^

github-actions · 2021-10-19T18:41:55Z

src/solver/impls/imex-bdf2/imex-bdf2.cxx

-      //Increment order if we're not at the maximum requested
-      if(order<maxOrder) order++;
+      // Increment order if we're not at the maximum requested
+      if (order < maxOrder)


warning: statement should be inside braces [readability-braces-around-statements]

if (order < maxOrder) ^ {

github-actions · 2021-10-19T18:41:55Z

src/solver/impls/imex-bdf2/imex-bdf2.cxx


    iteration++; // Advance iteration number

    /// Call the monitor function

-    if(call_monitors(simtime, s, nsteps)) {
+    if (call_monitors(simtime, s, nsteps)) {


warning: implicit conversion int -> bool [readability-implicit-bool-conversion]

if (call_monitors(simtime, s, nsteps)) { ^ != 0

Clang tidy suggestions

If timestep falls below this threshold, reset the Jacobian and preconditioner, and try to take a large timestep. The idea is that sometimes the solver gets stuck in a local minimum, and the observation that sometimes stopping the simulation and restarting can help.

github-actions

clang-tidy made some suggestions

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

github-actions · 2021-10-19T23:16:18Z

src/solver/impls/imex-bdf2/imex-bdf2.cxx

@@ -1152,97 +1176,102 @@ PetscErrorCode IMEXBDF2::solve_implicit(BoutReal curtime, BoutReal gamma) {
  implicit_gamma = gamma;

  // Set initial guess at the solution
-  BoutReal *xdata;
+  BoutReal* xdata;


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

BoutReal* xdata; ^ = nullptr

github-actions · 2021-10-19T23:16:18Z

src/solver/impls/imex-bdf2/imex-bdf2.cxx

-    for(int i=0;i<nlocal;i++) {
-      xdata[i] = uV[0][i];     // Use previous solution
+    for (int i = 0; i < nlocal; i++) {
+      xdata[i] = uV[0][i]; // Use previous solution


warning: do not use pointer arithmetic [cppcoreguidelines-pro-bounds-pointer-arithmetic]

xdata[i] = uV[0][i]; // Use previous solution ^

github-actions · 2021-10-19T23:16:19Z

src/solver/impls/imex-bdf2/imex-bdf2.cxx

-    for(int i=0;i<nlocal;i++) {
-      xdata[i] = 2.*uV[0][i] - uV[1][i];
+    for (int i = 0; i < nlocal; i++) {
+      xdata[i] = 2. * uV[0][i] - uV[1][i];


warning: do not use pointer arithmetic [cppcoreguidelines-pro-bounds-pointer-arithmetic]

xdata[i] = 2. * uV[0][i] - uV[1][i]; ^

github-actions · 2021-10-19T23:16:19Z

src/solver/impls/imex-bdf2/imex-bdf2.cxx

-    for(int i=0;i<nlocal;i++) {
-      xdata[i] = 3.*uV[0][i] - 3.*uV[1][i] + uV[2][i];
+    for (int i = 0; i < nlocal; i++) {
+      xdata[i] = 3. * uV[0][i] - 3. * uV[1][i] + uV[2][i];


warning: do not use pointer arithmetic [cppcoreguidelines-pro-bounds-pointer-arithmetic]

xdata[i] = 3. * uV[0][i] - 3. * uV[1][i] + uV[2][i]; ^

github-actions · 2021-10-19T23:16:19Z

src/solver/impls/imex-bdf2/imex-bdf2.cxx

-    for(int i=0;i<nlocal;i++) {
-      xdata[i] = rhs[i];   // If G = 0
+    for (int i = 0; i < nlocal; i++) {
+      xdata[i] = rhs[i]; // If G = 0


warning: do not use pointer arithmetic [cppcoreguidelines-pro-bounds-pointer-arithmetic]

xdata[i] = rhs[i]; // If G = 0 ^

github-actions · 2021-10-19T23:16:20Z

src/solver/impls/imex-bdf2/imex-bdf2.cxx

-              op.run(jx, jy, jz, u); ++u;
+      if (mesh->firstX() && !mesh->periodicX) {
+        for (int jx = 0; jx < mesh->xstart; ++jx)
+          for (int jy = mesh->ystart; jy <= mesh->yend; ++jy)


warning: statement should be inside braces [readability-braces-around-statements]

for (int jy = mesh->ystart; jy <= mesh->yend; ++jy) ^ {

src/solver/impls/imex-bdf2/imex-bdf2.cxx

github-actions · 2021-10-19T23:16:21Z

src/solver/impls/imex-bdf2/imex-bdf2.cxx

-            for(int jz=0; jz < mesh->LocalNz; ++jz) {
-              op.run(jx, jy, jz, u); ++u;
+      if (mesh->lastX() && !mesh->periodicX) {
+        for (int jx = mesh->xend + 1; jx < mesh->LocalNx; ++jx)


warning: statement should be inside braces [readability-braces-around-statements]

for (int jx = mesh->xend + 1; jx < mesh->LocalNx; ++jx) ^ {

github-actions · 2021-10-19T23:16:21Z

src/solver/impls/imex-bdf2/imex-bdf2.cxx

-              op.run(jx, jy, jz, u); ++u;
+      if (mesh->lastX() && !mesh->periodicX) {
+        for (int jx = mesh->xend + 1; jx < mesh->LocalNx; ++jx)
+          for (int jy = mesh->ystart; jy <= mesh->yend; ++jy)


warning: statement should be inside braces [readability-braces-around-statements]

for (int jy = mesh->ystart; jy <= mesh->yend; ++jy) ^ {

src/solver/impls/imex-bdf2/imex-bdf2.cxx

- If resetting doesn't work once, quit rather than get stuck in an infinite loop - Turn off the predictor when resetting. The predictor seems to make convergence slower when nearly in steady state.

github-actions

clang-tidy made some suggestions

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

github-actions · 2021-10-19T23:39:33Z

src/solver/impls/imex-bdf2/imex-bdf2.cxx

-          for(int jz=0; jz < mesh->LocalNz; ++jz) {
-            op.run(*xi, jy, jz, u); ++u;
+      for (RangeIterator xi = mesh->iterateBndryLowerY(); !xi.isDone(); ++xi) {
+        for (int jy = 0; jy < mesh->ystart; ++jy)


warning: statement should be inside braces [readability-braces-around-statements]

for (int jy = 0; jy < mesh->ystart; ++jy) ^ {

src/solver/impls/imex-bdf2/imex-bdf2.cxx

github-actions · 2021-10-19T23:39:34Z

src/solver/impls/imex-bdf2/imex-bdf2.cxx

-          for(int jz=0; jz < mesh->LocalNz; ++jz) {
-            op.run(*xi, jy, jz, u); ++u;
+      for (RangeIterator xi = mesh->iterateBndryUpperY(); !xi.isDone(); ++xi) {
+        for (int jy = mesh->yend + 1; jy < mesh->LocalNy; ++jy)


warning: statement should be inside braces [readability-braces-around-statements]

for (int jy = mesh->yend + 1; jy < mesh->LocalNy; ++jy) ^ {

github-actions · 2021-10-19T23:39:34Z

src/solver/impls/imex-bdf2/imex-bdf2.cxx

+        for (int jy = mesh->yend + 1; jy < mesh->LocalNy; ++jy)
+          for (int jz = 0; jz < mesh->LocalNz; ++jz) {
+            op.run(*xi, jy, jz, u);
+            ++u;


warning: do not use pointer arithmetic [cppcoreguidelines-pro-bounds-pointer-arithmetic]

++u; ^

github-actions · 2021-10-19T23:39:34Z

src/solver/impls/imex-bdf2/imex-bdf2.cxx

-      for(int jy=mesh->ystart; jy <= mesh->yend; ++jy)
-        for(int jz=0; jz < mesh->LocalNz; ++jz) {
-          op.run(jx, jy, jz, u); ++u;
+    for (int jx = mesh->xstart; jx <= mesh->xend; ++jx)


warning: statement should be inside braces [readability-braces-around-statements]

for (int jx = mesh->xstart; jx <= mesh->xend; ++jx) ^ {

src/solver/impls/snes/snes.cxx

github-actions · 2021-10-19T23:39:36Z

src/solver/impls/snes/snes.cxx

+      MatSetFromOptions(Jmf);
+
+      PetscInt *d_nnz, *o_nnz;
+      PetscMalloc((localN) * sizeof(PetscInt), &d_nnz);


warning: do not use C-style cast to convert between unrelated types [cppcoreguidelines-pro-type-cstyle-cast]

PetscMalloc((localN) * sizeof(PetscInt), &d_nnz); ^ /usr/lib/petsc/include/petscsys.h:453:99: note: expanded from macro 'PetscMalloc' #define PetscMalloc(a,b) ((*PetscTrMalloc)((a),PETSC_FALSE,__LINE__,PETSC_FUNCTION_NAME,__FILE__,(void**)(b))) ^

The docs suggest using PetscNew or PetscMalloc1 instead.

Could we even use std::vector<PetscInt> d_nnz(localN) here instead? I can't quite get my head around all the alignment business and if it's required, but I think this should be fine? This would eliminate the calls to PetscFree, and should hopefully avoid all the warnings about pointer arithmetic

ZedThree · 2021-10-20T07:59:43Z

I think we probably need to turn off cppcoreguidelines-pro-bounds-pointer-arithmetic, as it appears we use pointer arithmetic too much

Missing on for loops in imex-bdf2 code

github-actions

clang-tidy made some suggestions

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

src/solver/impls/snes/snes.cxx

dschwoerer · 2021-10-22T07:59:09Z

Shouldn't `output.print("\r")` be the proper replacement?

Thanks to @dschworer suggestion. [skip ci]

dschwoerer · 2021-10-25T08:33:17Z

For me the old defaults give better results 😞
Would it make sense to have a meta-option, to switch between different defaults? That would allow users to try different things, without having to worry about the various options initially ...

github-actions

clang-tidy made some suggestions

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

src/solver/impls/snes/snes.cxx

github-actions · 2021-11-05T16:27:36Z

src/solver/impls/snes/snes.cxx

+      // Mark non-zero entries
+
+      // Offsets for a 5-point pattern
+      const int xoffset[5] = {0, -1, 1, 0, 0};


warning: do not declare C-style arrays, use std::array<> instead [cppcoreguidelines-avoid-c-arrays]

const int xoffset[5] = {0, -1, 1, 0, 0}; ^

As this is new code, I think it's worth implementing this suggestion:

Suggested change

const int xoffset[5] = {0, -1, 1, 0, 0};

constexpr std::array<int, 5> xoffset = {0, -1, 1, 0, 0};

github-actions · 2021-11-05T16:27:36Z

src/solver/impls/snes/snes.cxx

+
+      // Offsets for a 5-point pattern
+      const int xoffset[5] = {0, -1, 1, 0, 0};
+      const int yoffset[5] = {0, 0, 0, -1, 1};


warning: do not declare C-style arrays, use std::array<> instead [cppcoreguidelines-avoid-c-arrays]

const int yoffset[5] = {0, 0, 0, -1, 1}; ^

Suggested change

const int yoffset[5] = {0, 0, 0, -1, 1};

constexpr std::array<int, 5> yoffset = {0, 0, 0, -1, 1};

src/solver/impls/snes/snes.cxx

If a step which would have ended at an output time fails, the looping variable should be reset to true so that the retry occurs. Thanks to Mike Kryjak for the report

Sets the linear solver and preconditioner to use

github-actions

clang-tidy made some suggestions

src/solver/impls/snes/snes.cxx

github-actions · 2021-11-09T23:05:33Z

src/solver/solver.cxx

@@ -1246,13 +1246,13 @@ int Solver::run_rhs(BoutReal t) {

    save_vars(tmp.begin()); // Copy variables into tmp
    pre_rhs(t);
-    status = model->runConvective(t);
+    status = model->runConvective(t, linear);


warning: Value stored to status is never read [clang-analyzer-deadcode.DeadStores]

status = model->runConvective(t, linear); ^ src/solver/solver.cxx:1249:5: note: Value stored to 'status' is never read

beuler sometimes appears to reach a steady state, with timesteps continuously increasing, and zero nonlinear iterations e.g. ``` Time: 2189781.780775979, timestep: 9781.780775979214, nl iter: 0, lin iter: 1 Time: 2190000.0, timestep: 10759.958853577136, nl iter: 0, lin iter: 1 2.190e+06 31 1.94e+00 94.9 0.0 0.0 3.0 2.0 Time: 2200000.0, timestep: 10759.958853577136, nl iter: 0, lin iter: 1 2.200e+06 16 1.05e+00 91.0 0.0 0.0 6.3 2.7 Time: 2210000.0, timestep: 10759.958853577136, nl iter: 0, lin iter: 1 2.210e+06 16 1.08e+00 91.9 0.0 0.0 5.5 2.6 Time: 2220000.0, timestep: 10759.958853577136, nl iter: 0, lin iter: 1 ``` Restarting from this state however shows that it is not in steady state and even doesn't converge. To try to prevent this, force SNES to take at least one iteration, add more checks and reporting.

Will sometimes hit a state where snes stops after zero iterations (due to stol tolerance), and continues "converging" until the end of the simulation. Taking a small Euler step seems to help.

Add ability to change line search type. When SNES fails to converge, print diagnostic information on the fields and their time derivatives.

ZedThree

LGTM. There's a few bits that could be polished perhaps. Lots of noise from clang-tidy because PETSc.

There's a fair bit of repeated code between IMEX-BDF2 and SNES for the colouring -- is it possible to pull this out into a shared helper function/class?

src/solver/impls/imex-bdf2/imex-bdf2.cxx

ZedThree · 2021-11-22T09:20:13Z

src/solver/impls/imex-bdf2/imex-bdf2.cxx

-      const int xoffset[5] = {0,-1, 1, 0, 0};
-      const int yoffset[5] = {0, 0, 0,-1, 1};
-
+      const int xoffset[5] = {0, -1, 1, 0, 0};


Not touched in this PR, but this is a good suggestion we should do at some point

src/solver/impls/snes/snes.cxx

ZedThree · 2021-11-22T09:33:45Z

src/solver/impls/snes/snes.cxx

+      MatSetFromOptions(Jmf);
+
+      PetscInt *d_nnz, *o_nnz;
+      PetscMalloc((localN) * sizeof(PetscInt), &d_nnz);


The docs suggest using PetscNew or PetscMalloc1 instead.

Could we even use std::vector<PetscInt> d_nnz(localN) here instead? I can't quite get my head around all the alignment business and if it's required, but I think this should be fine? This would eliminate the calls to PetscFree, and should hopefully avoid all the warnings about pointer arithmetic

ZedThree · 2021-11-22T09:37:11Z

src/solver/impls/snes/snes.cxx

+            } else {
+              // Only 3D fields
+              for (int i = 0; i < n3d; i++) {
+                d_nnz[localIndex + i] -= (n3d + n2d);


Just flagging this to be double-checked: is the rhs supposed to match the loop limit? i.e. should this be

Suggested change

d_nnz[localIndex + i] -= (n3d + n2d);

d_nnz[localIndex + i] -= n3d;

ZedThree · 2021-11-22T09:40:48Z

src/solver/impls/snes/snes.cxx

+            if (z == 0) {
+              // All 2D and 3D fields
+              for (int i = 0; i < n2d + n3d; i++) {
+                d_nnz[localIndex + i] -= (n3d + n2d);
+              }
+            } else {
+              // Only 3D fields
+              for (int i = 0; i < n3d; i++) {
+                d_nnz[localIndex + i] -= (n3d + n2d);
+              }


Assuming the loop body is supposed to be identical between the branches, here's a more concise way of writing this that avoids the repetition in the loop body:

Suggested change

if (z == 0) {

// All 2D and 3D fields

for (int i = 0; i < n2d + n3d; i++) {

d_nnz[localIndex + i] -= (n3d + n2d);

}

} else {

// Only 3D fields

for (int i = 0; i < n3d; i++) {

d_nnz[localIndex + i] -= (n3d + n2d);

}

const auto num_fields = (z == 0) ? n2d + n3d : n3d;

for (int i = 0; i < num_fields; i++) {

d_nnz[localIndex + i] -= (n3d + n2d);

}

I also wonder if these loop bodies could be wrapped up into functions and reused? Might cut down on this section a fair bit

ZedThree · 2021-11-22T09:46:12Z

src/solver/impls/snes/snes.cxx

+      // Mark non-zero entries
+
+      // Offsets for a 5-point pattern
+      const int xoffset[5] = {0, -1, 1, 0, 0};


As this is new code, I think it's worth implementing this suggestion:

Suggested change

const int xoffset[5] = {0, -1, 1, 0, 0};

constexpr std::array<int, 5> xoffset = {0, -1, 1, 0, 0};

ZedThree · 2021-11-22T09:46:39Z

src/solver/impls/snes/snes.cxx

+
+      // Offsets for a 5-point pattern
+      const int xoffset[5] = {0, -1, 1, 0, 0};
+      const int yoffset[5] = {0, 0, 0, -1, 1};


Suggested change

const int yoffset[5] = {0, 0, 0, -1, 1};

constexpr std::array<int, 5> yoffset = {0, 0, 0, -1, 1};

ZedThree · 2021-11-22T09:47:59Z

src/solver/impls/snes/snes.cxx

+            PetscInt row = ind0 + i;
+
+            // Loop through each point in the 5-point stencil
+            for (int c = 0; c < 5; c++) {


One day C++ will have zip!

ZedThree · 2021-11-22T09:55:09Z

src/solver/impls/snes/snes.hxx

+  /// @param[out] f  The vector for the result f(x)
+  /// @param[in] linear  Specifies that the SNES solver is in a linear (KSP) inner loop,
+  ///                    so the operator should be linearised if possible
+  PetscErrorCode snes_function(Vec x, Vec f, bool linear); ///< Nonlinear function


Suggested change

PetscErrorCode snes_function(Vec x, Vec f, bool linear); ///< Nonlinear function

PetscErrorCode snes_function(Vec x, Vec f, bool linear);

johnomotani · 2021-11-22T10:25:48Z

Sorry, didn't have time to look through the code yet (will try to get to it).

One thought though - this might be something for a separate PR, but I don't think these solvers are using the Petsclib features for setting options that were added in #1795 (like for example LaplaceXY and LaplacePetsc3dAmg do). The PETSc options feature was backported to 4.4 too. Petsclib might need a small update to support setting options for SNES (maybe an extra method - I haven't checked or thought about it...) but using it has the advantage that all PETSc options are available without having to add a BOUT++ option and call the setter function - including new options in future PETSc versions!

If this is a good thing to do, it might be good to add sooner rather than later, since it would change the input file structure - i.e. use a beuler:petsc subsection instead of options in beuler.

github-actions

clang-tidy made some suggestions

src/solver/impls/snes/snes.cxx

Fails to compile on github due to using petsc 3.7.7

equation_form switches between: - Pseudo-transient (like UEDGE) - A rearranged backward Euler - The original backward Euler

ZedThree

LGTM, but a couple of things to double check, and a small number of things to fix

src/solver/impls/snes/snes.cxx

ZedThree · 2021-12-08T16:12:47Z

src/solver/impls/snes/snes.cxx

+
+          // Only 3D fields
+          for (int i = 0; i < n3d; i++) {
+            // d_nnz[localIndex+i] -= (n3d + n2d);


Suggested change

// d_nnz[localIndex+i] -= (n3d + n2d);

ZedThree · 2021-12-08T16:13:18Z

src/solver/impls/snes/snes.cxx

+        localIndex = ROUND(index(x, mesh->yend, 0));
+        // All 2D and 3D fields
+        for (int i = 0; i < n2d + n3d; i++) {
+          // d_nnz[localIndex+i] -= (n3d + n2d);


Suggested change

// d_nnz[localIndex+i] -= (n3d + n2d);

ZedThree · 2021-12-08T16:13:35Z

src/solver/impls/snes/snes.cxx

+
+          // Only 3D fields
+          for (int i = 0; i < n3d; i++) {
+            // d_nnz[localIndex+i] -= (n3d + n2d);


Suggested change

// d_nnz[localIndex+i] -= (n3d + n2d);

ZedThree · 2021-12-08T16:17:18Z

src/solver/impls/snes/snes.cxx

+            if (z == 0) {
+              // All 2D and 3D fields
+              for (int i = 0; i < n2d + n3d; i++) {
+                d_nnz[localIndex + i] -= (n3d + n2d);
+              }
+            } else {
+              // Only 3D fields
+              for (int i = 0; i < n3d; i++) {
+                d_nnz[localIndex + i] -= (n3d + n2d);
+              }


I also wonder if these loop bodies could be wrapped up into functions and reused? Might cut down on this section a fair bit

ZedThree · 2021-12-08T16:18:28Z

src/solver/impls/snes/snes.cxx

+              for (int j = 0; j < n2d; j++) {
+                PetscInt col = ind2 + j;
+
+                MatSetValues(Jmf, 1, &row, 1, &col, &val, INSERT_VALUES);


Is it worth looking at using our Petsc Matrix wrapper for this new code?

ZedThree · 2021-12-08T16:55:23Z

src/solver/impls/snes/snes.cxx


+  equation_form = (*options)["equation_form"]


Could equation_form be a BOUT_ENUM_CLASS? Then users could use the names directly, rather than integers

In y boundaries the X index can be out of the domain, leading to negative indices. This caused out of bounds memory access, and a lot of slowdown in the Jacobian coloring setup. Also added some progress output

github-actions

clang-tidy made some suggestions

src/solver/impls/snes/snes.cxx

Also add some notes on preconditioner options

…ev into beuler-jacobian-color

Co-authored-by: Peter Hill <zed.three@gmail.com>

Limits how large the timestep can be made, to try and prevent repeated increases and failures

beuler solver additions

* next: (37 commits) Merge in Solver and PhysicsModel changes from next SNES solver merges Fix contributor's name Update release date Update changelog Fix ambiguous visit call Revert SONAME/SOVERSION to 4.4.0 Fix test-bout-override-default-option for cmake Update DOI and release date for 4.4.1 Update changelog for 4.4.1 Update translation files for 4.4.1 Bump version to 4.4.1 Enable setting SNES solver PETSc options from input file Try to consolidate some loops/branches in SNESSolver Use `BOUT_ENUM_CLASS` for `SNESSolver::equation_form` Use `std::vector/array` instead of C style arrays Remove `__FUNCT__` macros for PETSc callbacks Update backport of beuler/snes solver Fix typo in docs Remove IDL section and recommend xBOUT for analysis ...

ZedThree · 2022-02-08T15:08:10Z

The black check failed because we're running it on both pull requests and pushes. There's a few CI things to clean up, I'll do so in another PR

bendudson added 7 commits October 15, 2021 10:32

Fix petsc interface loops

de3e6bc

Should use const reference to avoid copying index

Replace finite with std::isfinite

06fa323

Deprecated, flagged as warning by Apple Clang

Remove unused private variables in PETSc solvers

422a5fd

Fix imexbdf2 coloring

ba17afb

Was destroying the iscoloring object before setup, and then reverting to brute-force calculation of the Jacobian.

Clang format imexbdf2

3200e4d

github-actions bot reviewed Oct 15, 2021

View reviewed changes

bendudson added 4 commits October 15, 2021 14:35

Update beuler manual section

39e500b

Describe all the options, and how to use newtonls with Jacobian coloring.

Clang tidy, add couple of braces

89e3a67

[skip ci]

SNES / beuler solver change defaults, diagnostics

34f1a8b

- Defaults to using coloring to calculate a Jacobian - Iterations and lag jacobian set to values that work for SD1D tests so far. - Diagnostic outputs now include number of linear iterations and number of failures.

Update documentation on snes solver

436ee6d

Reflect changes in defaults, and some notes on when these are likely to fail.

github-actions bot reviewed Oct 19, 2021

View reviewed changes

bendudson added 2 commits October 19, 2021 15:36

imex-bdf2 add missing braces

4f42716

Clang tidy suggestions

github-actions bot reviewed Oct 19, 2021

View reviewed changes

beuler/snes recovery improvements

e5902cf

- If resetting doesn't work once, quit rather than get stuck in an infinite loop - Turn off the predictor when resetting. The predictor seems to make convergence slower when nearly in steady state.

github-actions bot reviewed Oct 19, 2021

View reviewed changes

Clang tidy some braces

4516729

Missing on for loops in imex-bdf2 code

github-actions bot reviewed Oct 20, 2021

View reviewed changes

dschwoerer reviewed Oct 21, 2021

View reviewed changes

src/solver/impls/snes/snes.cxx Outdated Show resolved Hide resolved

Replace std::cout with output.print

6746a74

Thanks to @dschworer suggestion. [skip ci]

github-actions bot reviewed Nov 5, 2021

View reviewed changes

bendudson added 2 commits November 5, 2021 16:34

beuler fix: Keep looping on failure

c2b60cc

If a step which would have ended at an output time fails, the looping variable should be reset to true so that the retry occurs. Thanks to Mike Kryjak for the report

Add beuler solver:pc_type and solver:ksp_type settings

e4dd1bf

Sets the linear solver and preconditioner to use

github-actions bot reviewed Nov 9, 2021

View reviewed changes

bendudson mentioned this pull request Nov 11, 2021

Back-port beuler/snes improvements to master #2459

Merged

bendudson added 3 commits November 17, 2021 17:17

If SNES thinks it has converged, try an Euler step

ced8530

Will sometimes hit a state where snes stops after zero iterations (due to stol tolerance), and continues "converging" until the end of the simulation. Taking a small Euler step seems to help.

Add beuler line_search_type and more diagnostics

b6e9918

Add ability to change line search type. When SNES fails to converge, print diagnostic information on the fields and their time derivatives.

ZedThree previously approved these changes Nov 22, 2021

View reviewed changes

bendudson dismissed ZedThree’s stale review via b6e9918 November 24, 2021 00:07

github-actions bot reviewed Nov 24, 2021

View reviewed changes

src/solver/impls/snes/snes.cxx Outdated Show resolved Hide resolved

src/solver/impls/snes/snes.cxx Outdated Show resolved Hide resolved

src/solver/impls/snes/snes.cxx Outdated Show resolved Hide resolved

bendudson added 2 commits November 30, 2021 22:19

SNESSetForceIteration was added in PETSc 3.8

44cfd58

Fails to compile on github due to using petsc 3.7.7

Add different forms of the nonlinear equation

41a07b1

equation_form switches between: - Pseudo-transient (like UEDGE) - A rearranged backward Euler - The original backward Euler

ZedThree requested changes Dec 8, 2021

View reviewed changes

bendudson and others added 2 commits December 14, 2021 15:19

SNES/beuler solver: Fix out of bounds indices

88652cd

In y boundaries the X index can be out of the domain, leading to negative indices. This caused out of bounds memory access, and a lot of slowdown in the Jacobian coloring setup. Also added some progress output

[skip ci] Apply black changes

22c3e0b

github-actions bot reviewed Dec 14, 2021

View reviewed changes

src/solver/impls/snes/snes.cxx Outdated Show resolved Hide resolved

src/solver/impls/snes/snes.cxx Outdated Show resolved Hide resolved

src/solver/impls/snes/snes.cxx Outdated Show resolved Hide resolved

bendudson and others added 8 commits December 14, 2021 17:48

Fix upper Y boundary too

160970b

Also add some notes on preconditioner options

Merge branch 'beuler-jacobian-color' of github.com:boutproject/BOUT-d…

58ae3be

…ev into beuler-jacobian-color

Update src/solver/impls/snes/snes.cxx

7803a35

Co-authored-by: Peter Hill <zed.three@gmail.com>

Merge branch 'next' into beuler-jacobian-color

374ecdd

Add max_timestep option to beuler solver

3be4293

Limits how large the timestep can be made, to try and prevent repeated increases and failures

Merge pull request #2489 from boutproject/beuler-additions

7938e9c

beuler solver additions

[skip ci] Apply black changes

3a42624

ZedThree approved these changes Feb 8, 2022

View reviewed changes

ZedThree merged commit d746cda into next Feb 9, 2022

ZedThree deleted the beuler-jacobian-color branch February 9, 2022 11:09

	const int xoffset[5] = {0, -1, 1, 0, 0};
	constexpr std::array<int, 5> xoffset = {0, -1, 1, 0, 0};

	const int yoffset[5] = {0, 0, 0, -1, 1};
	constexpr std::array<int, 5> yoffset = {0, 0, 0, -1, 1};

	d_nnz[localIndex + i] -= (n3d + n2d);
	d_nnz[localIndex + i] -= n3d;

	PetscErrorCode snes_function(Vec x, Vec f, bool linear); ///< Nonlinear function
	PetscErrorCode snes_function(Vec x, Vec f, bool linear);

Add matrix coloring to beuler solver (and fix imexbdf2) #2454

Add matrix coloring to beuler solver (and fix imexbdf2) #2454

Conversation

bendudson commented Oct 15, 2021

github-actions bot left a comment

Choose a reason for hiding this comment

github-actions bot Oct 15, 2021

Choose a reason for hiding this comment

github-actions bot Oct 15, 2021

Choose a reason for hiding this comment

github-actions bot Oct 15, 2021

Choose a reason for hiding this comment

github-actions bot Oct 15, 2021

Choose a reason for hiding this comment

github-actions bot Oct 15, 2021

Choose a reason for hiding this comment

github-actions bot Oct 15, 2021

Choose a reason for hiding this comment

github-actions bot left a comment

Choose a reason for hiding this comment

github-actions bot Oct 19, 2021

Choose a reason for hiding this comment

github-actions bot Oct 19, 2021

Choose a reason for hiding this comment

github-actions bot Oct 19, 2021

Choose a reason for hiding this comment

github-actions bot Oct 19, 2021

Choose a reason for hiding this comment

github-actions bot Oct 19, 2021

Choose a reason for hiding this comment

Choose a reason for hiding this comment

github-actions bot Oct 19, 2021

Choose a reason for hiding this comment

github-actions bot Oct 19, 2021

Choose a reason for hiding this comment

github-actions bot Oct 19, 2021

Choose a reason for hiding this comment

github-actions bot Oct 19, 2021

Choose a reason for hiding this comment

github-actions bot Oct 19, 2021

Choose a reason for hiding this comment

github-actions bot left a comment

Choose a reason for hiding this comment

github-actions bot Oct 19, 2021

Choose a reason for hiding this comment

github-actions bot Oct 19, 2021

Choose a reason for hiding this comment

github-actions bot Oct 19, 2021

Choose a reason for hiding this comment

github-actions bot Oct 19, 2021

Choose a reason for hiding this comment

github-actions bot Oct 19, 2021

Choose a reason for hiding this comment

github-actions bot Oct 19, 2021

Choose a reason for hiding this comment

github-actions bot Oct 19, 2021

Choose a reason for hiding this comment

github-actions bot Oct 19, 2021

Choose a reason for hiding this comment

github-actions bot left a comment

Choose a reason for hiding this comment

github-actions bot Oct 19, 2021

Choose a reason for hiding this comment

github-actions bot Oct 19, 2021

Choose a reason for hiding this comment

github-actions bot Oct 19, 2021

Choose a reason for hiding this comment

github-actions bot Oct 19, 2021

Choose a reason for hiding this comment

github-actions bot Oct 19, 2021

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ZedThree commented Oct 20, 2021

github-actions bot left a comment

Choose a reason for hiding this comment

dschwoerer commented Oct 22, 2021 via email

dschwoerer commented Oct 25, 2021

github-actions bot left a comment

Choose a reason for hiding this comment

github-actions bot Nov 5, 2021