regression in build_pc_ilu in ex09

[gdmcbain]

ex09 is running differently to how I remember. I remember the second solve with solver_iter_pcg, using M=build_pc_ilu(Aint), as being faster than the first, withM=None`. (That was presumably the point of the example.) However now, it seems that the second solve is failing with

UserWarning: Convergence not achieved!

This doesn't fail tests.test_examples because there there is a third preconditioner, based on pyamg.smoothed_aggregation_solver, which overwrites the solution and it's only the last value of ex09.x that is assessed.

I've tested this on Ubuntu 18.04.5 LTS, 20.04.1 LTS, and Microsoft Windows 10.

Perhaps tests.test_examples.TextEx09 should be strengthened to test all preconditioners present?

I note that build_pc_ilu is also exercised in ex30 (Uzawa) and it is still working fine.

[gdmcbain]

pip install scikit-fem==1.0.0 is O. K. (all four preconditioners).

[gdmcbain]

As is 1.2.0, but 2.0.0 introduces the failure. That's in terms of releases; I haven't tried git bisect yet.

[gdmcbain]

This is running with conda install scipy, which gives scipy 1.5.2.

[kinnala]

Diff between 1.2.0 and 2.0.0 in utils.py:

diff --git a/skfem/utils.py b/skfem/utils.py
index d348ce8..aab3f68 100644
--- a/skfem/utils.py
+++ b/skfem/utils.py
@@ -2,8 +2,8 @@
 SciPy linear solvers."""

 import warnings
-from typing import Optional, Union, Tuple,\
-    Callable, Dict
+from typing import Optional, Union, Tuple, Callable, Dict
+from inspect import signature

 import numpy as np
 import scipy.sparse as sp
@@ -12,7 +12,7 @@ import scipy.sparse.linalg as spl
 from numpy import ndarray
 from scipy.sparse import spmatrix

-from skfem.assembly import asm, BilinearForm, LinearForm, Dofs
+from skfem.assembly import asm, BilinearForm, LinearForm, DofsView
 from skfem.assembly.basis import Basis
 from skfem.element import ElementVectorH1

@@ -28,7 +28,7 @@ CondensedSystem = Union[spmatrix,
                         Tuple[spmatrix, ndarray, ndarray],
                         Tuple[spmatrix, ndarray, ndarray, ndarray],
                         Tuple[spmatrix, spmatrix, ndarray, ndarray]]
-DofsCollection = Union[ndarray, Dofs, Dict[str, Dofs]]
+DofsCollection = Union[ndarray, DofsView, Dict[str, DofsView]]

 # preconditioners, e.g. for :func:`skfem.utils.solver_iter_krylov`
@@ -76,6 +76,7 @@ def solver_eigen_scipy(**kwargs) -> EigenSolver:

 def solver_direct_scipy(**kwargs) -> LinearSolver:
+    """The default linear solver of SciPy."""
     def solver(A, b, **solve_time_kwargs):
         kwargs.update(solve_time_kwargs)
         return spl.spsolve(A, b, **kwargs)
@@ -185,10 +186,10 @@ def _flatten_dofs(S: DofsCollection) -> ndarray:
     else:
         if isinstance(S, ndarray):
             return S
-        elif isinstance(S, Dofs):
-            return S.all()
+        elif isinstance(S, DofsView):
+            return S.flatten()
         elif isinstance(S, dict):
-            return np.unique(np.concatenate([S[key].all() for key in S]))
+            return np.unique(np.concatenate([S[key].flatten() for key in S]))
         raise NotImplementedError("Unable to flatten the given set of DOF's.")

@@ -200,7 +201,15 @@ def condense(A: spmatrix,
              expand: bool = True) -> CondensedSystem:
     """Eliminate degrees-of-freedom from a linear system.

-    Supports also generalized eigenvalue problems.
+    The user should provide the linear system ``A`` and ``b``
+    and either the set of DOF's to eliminate (``D``) or the set
+    of DOF's to keep (``I``).  Optionally, nonzero values for
+    the eliminated DOF's can be supplied via ``x``.
+
+    .. note::
+
+        Supports also generalized eigenvalue problems
+        where ``b`` is a matrix.

     Parameters
     ----------
@@ -219,13 +228,14 @@ def condense(A: spmatrix,
     expand
         If `True` (default), returns also `x` and `I`. As a consequence,
         :func:`skfem.utils.solve` will expand the solution vector
-        automatically. Otherwise, r
+        automatically.

     Returns
     -------
     CondensedSystem
-        The condensed linear system (and optionally information about
-        the boundary values).
+        The condensed linear system and (optionally) information about
+        the boundary values.
+
     """
     D = _flatten_dofs(D)
     I = _flatten_dofs(I)
@@ -275,16 +285,17 @@ def rcm(A: spmatrix,
 def adaptive_theta(est, theta=0.5, max=None):
     """For choosing which elements to refine in an adaptive strategy."""
     if max is None:
-        return np.nonzero(theta*np.max(est) < est)[0]
+        return np.nonzero(theta * np.max(est) < est)[0]
     else:
-        return np.nonzero(theta*max < est)[0]
+        return np.nonzero(theta * max < est)[0]

 def project(fun,
 import numpy as np
 import scipy.sparse as sp
@@ -12,7 +12,7 @@ import scipy.sparse.linalg as spl
 from numpy import ndarray
 from scipy.sparse import spmatrix

-from skfem.assembly import asm, BilinearForm, LinearForm, Dofs
+from skfem.assembly import asm, BilinearForm, LinearForm, DofsView
 from skfem.assembly.basis import Basis
 from skfem.element import ElementVectorH1

@@ -28,7 +28,7 @@ CondensedSystem = Union[spmatrix,
                         Tuple[spmatrix, ndarray, ndarray],
                         Tuple[spmatrix, ndarray, ndarray, ndarray],
                         Tuple[spmatrix, spmatrix, ndarray, ndarray]]
-DofsCollection = Union[ndarray, Dofs, Dict[str, Dofs]]
+DofsCollection = Union[ndarray, DofsView, Dict[str, DofsView]]

 # preconditioners, e.g. for :func:`skfem.utils.solver_iter_krylov`
@@ -76,6 +76,7 @@ def solver_eigen_scipy(**kwargs) -> EigenSolver:

 def solver_direct_scipy(**kwargs) -> LinearSolver:
+    """The default linear solver of SciPy."""
     def solver(A, b, **solve_time_kwargs):
         kwargs.update(solve_time_kwargs)
         return spl.spsolve(A, b, **kwargs)
@@ -185,10 +186,10 @@ def _flatten_dofs(S: DofsCollection) -> ndarray:
     else:
         if isinstance(S, ndarray):
             return S
-        elif isinstance(S, Dofs):
-            return S.all()
+        elif isinstance(S, DofsView):
+            return S.flatten()
         elif isinstance(S, dict):
-            return np.unique(np.concatenate([S[key].all() for key in S]))
+            return np.unique(np.concatenate([S[key].flatten() for key in S]))
         raise NotImplementedError("Unable to flatten the given set of DOF's.")

@@ -200,7 +201,15 @@ def condense(A: spmatrix,
              expand: bool = True) -> CondensedSystem:
     """Eliminate degrees-of-freedom from a linear system.

-    Supports also generalized eigenvalue problems.
+    The user should provide the linear system ``A`` and ``b``
+    and either the set of DOF's to eliminate (``D``) or the set
+    of DOF's to keep (``I``).  Optionally, nonzero values for
+    the eliminated DOF's can be supplied via ``x``.
+

Doesn't seem to touch the preconditioner part? Maybe there is another bug somewhere?

[kinnala]

Perhaps tests.test_examples.TextEx09 should be strengthened to test all preconditioners present?

Sounds like a good idea.

[gdmcbain]

According to git bisect

8d6446de63f1bcd6548483f195bc53d9c9575c23 is the first bad commit

[gdmcbain]

Maybe the local reimplementation of MeshTet.init_tensor, previously delegated to skfem.zoo.tet_tensor.build?

https://github.com/kinnala/scikit-fem/commit/8d6446de63f1bcd6548483f195bc53d9c9575c23#diff-7d5249076e12208eb9bff97b642c77ff

[gdmcbain]

Yep; ex09 runs fine in master if the MeshTet is build with meshzoo.cube.

from meshzoo import cube

# p = np.linspace(0, 1, 16)
# m = MeshTet.init_tensor(*(p,) * 3)
points, cells = cube(nx=16, ny=16, nz=16)
m = MeshTet(points.T, cells.T)

[gdmcbain]

I say fine, though it does fail the test (pytest -k ex09 tests/test_examples.py):

FAILED tests/test_examples.py::TestEx09::runTest - AssertionError: 0.05559677587711787 != 0.05528520791811886 within 7 places (0.00031156795899901085 difference)

but then 8d6446d did move the goalposts too.

        -self.assertAlmostEqual(np.max(ex09.x), 0.055596791644282988)
        +self.assertAlmostEqual(np.max(ex09.x), 0.05528520791811886)

[gdmcbain]

The example also runs fine with a refined MeshTet().

m = MeshTet()
m.refine(4)

This isn't exactly the same size but quite close: 20480 tets, cf. 20250 in the current example.

[kinnala]

Default parameters of build_pc_ilu may not be the best ones for all scenarios. But now it seems to converge.