[WIP] Harmonize solvers

docs/source/USAGE.rst

@@ -278,6 +278,14 @@ package.
    print a variety of diagnostic information.  All other solvers should use
    `Logging`_ and :envvar:`FIPY_LOG_CONFIG`.
 
+.. envvar:: FIPY_DEFAULT_CRITERION
+
+   Changes the default solver :ref:`CONVERGENCE` criterion to the specified
+   value.  Valid choices are "``legacy``", "``unscaled``", "``RHS``",
+   "``matrix``", "``initial``", "``solution``", "``preconditioned``",
+   "``natural``", "``default``".  A value of "``default``" is admittedly
+   circular, but it works.
+
 .. envvar:: FIPY_VIEWER
 
    Forces the use of the specified viewer. Valid values are any
@@ -303,7 +311,7 @@ package.
    :ref:`PETSC` solvers in order to see what options are possible.  Ignored
    if solver is not :ref:`PETSC`.
 
-.. _PETSc configuration options: https://docs.petsc.org/en/latest/manual/other/#sec-options
+.. _PETSc configuration options: https://petsc.org/main/manual/other/#runtime-options
 
 .. _PARALLEL:
 
@@ -1063,13 +1071,13 @@ Thanks to the future_ package and to the contributions of pya_ and
 woodscn_, :term:`FiPy` runs under both :term:`Python 3` and :term:`Python`
 2.7, without conversion or modification.
 
-Because :term:`Python` itself will `drop support for Python 2.7 on January
-1, 2020`_ and many of the prerequisites for :term:`FiPy` have `pledged to
-drop support for Python 2.7 no later than 2020`_, we have prioritized adding
+Because :term:`Python` itself `dropped support for Python 2.7 on January
+1, 2020`_ and many of the prerequisites for :term:`FiPy` `pledged to
+drop support for Python 2.7 no later than 2020`_, we prioritized adding
 support for better :term:`Python 3` solvers, starting with
 :term:`petsc4py`.
 
-Because the faster :term:`PySparse` and :term:`Trilinos` solvers are not
+Because the faster :term:`PySparse` solvers are not
 available under :term:`Python 3`, we will maintain :term:`Python` 2.x
 support as long as practical. Be aware that the conda-forge_ packages that
 :term:`FiPy` depends upon are not well-maintained on :term:`Python` 2.x
@@ -1081,7 +1089,7 @@ become available on conda-forge_.
 .. _future: http://python-future.org
 .. _pya: https://github.com/pya
 .. _woodscn: https://github.com/pya
-.. _drop support for Python 2.7 on January 1, 2020: https://www.python.org/dev/peps/pep-0373/#update
+.. _dropped support for Python 2.7 on January 1, 2020: https://www.python.org/dev/peps/pep-0373/#update
 .. _pledged to drop support for Python 2.7 no later than 2020: https://python3statement.org
 
 ------

docs/source/_static/residuals.csv

@@ -0,0 +1,10 @@
+criterion,residual,scale,PETSc_ [#KSP_Convergence_Tests]_,pyamgx_ [#AMGX_convergence]_,PySparse_,SciPy_ [#SciPy_Convergence_Test]_,Trilinos_ [#AztecOO_convergence]_
+``unscaled``,:math:`\|\mathsf{L}\vec{x} - \vec{b}\|_2`,:math:`1`,[#FiPy_Convergence_Test]_,``ABSOLUTE``,[#FiPy_Convergence_Test]_,[#FiPy_Convergence_Test]_,``AZ_noscaled``
+``RHS``,:math:`\|\mathsf{L}\vec{x} - \vec{b}\|_2`,:math:`\|\vec{b}\|_2`,|KSP_NORM_UNPRECONDITIONED|_,[#FiPy_Convergence_Test]_,"`cgs <https://sourceforge.net/p/pysparse/git/ci/master/tree/pysparse/itsolvers/src/cgs.c#l90>`_, `pcg <https://sourceforge.net/p/pysparse/git/ci/master/tree/pysparse/itsolvers/src/pcg.c#l154>`_, `qmres <https://sourceforge.net/p/pysparse/git/ci/master/tree/pysparse/itsolvers/src/qmrs.c#l139>`_, or [#FiPy_Convergence_Test]_",``default``,``AZ_rhs``
+``matrix``,:math:`\|\mathsf{L}\vec{x} - \vec{b}\|_2`,:math:`\|\mathsf{L}\|_\infty`,[#FiPy_Convergence_Test]_,[#FiPy_Convergence_Test]_,[#FiPy_Convergence_Test]_,[#FiPy_Convergence_Test]_,``AZ_Anorm``
+``initial``,:math:`\|\mathsf{L}\vec{x} - \vec{b}\|_2`,:math:`\|\mathsf{L}\vec{x} - \vec{b}\|_2^{(0)}`,[#FiPy_Convergence_Test]_,``RELATIVE_INI_CORE``,"`bicgstab <https://sourceforge.net/p/pysparse/git/ci/master/tree/pysparse/itsolvers/src/bicgstab.c#l287>`_, `gmres <https://sourceforge.net/p/pysparse/git/ci/master/tree/pysparse/itsolvers/src/gmres.c#l159>`_, `minres <https://sourceforge.net/p/pysparse/git/ci/master/tree/pysparse/itsolvers/src/minres.c#l195>`_, or [#FiPy_Convergence_Test]_",[#FiPy_Convergence_Test]_,``AZ_r0``
+``solution``,:math:`\|\mathsf{L}\vec{x} - \vec{b}\|_\infty`,:math:`\|\mathsf{L}\|_\infty  * \|\vec{x}\|_1 + \|\vec{b}\|_\infty`,,,,,``AZ_sol``
+``preconditioned``,:math:`\left\|\mathsf{P}^{-1}(\mathsf{L}\vec{x} - \vec{b})\right\|_2`,:math:`\left\|\vec{b}\right\|_2`,|KSP_NORM_PRECONDITIONED|_,,,,
+``natural``,:math:`\sqrt{(\mathsf{L}\vec{x} - \vec{b})\mathsf{P}^{-1}(\mathsf{L}\vec{x} - \vec{b})}`,:math:`\left\|\vec{b}\right\|_2`,|KSP_NORM_NATURAL|_,,,,
+``legacy``,,,|KSP_NORM_DEFAULT|_ (``RHS`` or ``preconditioned``),``initial``,``RHS`` or ``initial``,``RHS``,``initial``
+``default``,,,``RHS``,``RHS``,``RHS``,``RHS``,``RHS``

docs/source/_static/solver_convergence.csv

@@ -0,0 +1,9 @@
+,,PETSc_,pyamgx_,PySparse_,SciPy_,Trilinos_
+:class:`~fipy.solvers.convergence.Convergence`,Convergence criteria met.,,|AMGX_SOLVE_SUCCESS|_,,|SciPy_0|_,|AZ_normal|_
+:class:`~fipy.solvers.convergence.IterationConvergence`,Requested iterations complete (and no residual calculated).,|KSP_CONVERGED_ITS|_,,,,
+:class:`~fipy.solvers.convergence.AbsoluteToleranceConvergence`,"Converged, residual is as small as seems reasonable on this machine.",|KSP_CONVERGED_ATOL|_,,|PySparse_2|_,,
+:class:`~fipy.solvers.convergence.RHSZeroConvergence`,"Converged, :math:`\mathbf{b} = 0`, so the exact solution is :math:`\mathbf{x} = 0`.",,,|PySparse_1|_,,
+:class:`~fipy.solvers.convergence.RelativeToleranceConvergence`,"Converged, relative error appears to be less than tolerance.",|KSP_CONVERGED_RTOL|_,,|PySparse_0|_,,
+:class:`~fipy.solvers.convergence.HappyBreakdownConvergence`,"""Exact"" solution found and more iterations will just make things worse.",|KSP_CONVERGED_HAPPY_BREAKDOWN|_,,,,
+:class:`~fipy.solvers.convergence.LossOfAccuracyConvergence`,The iterative solver has terminated due to a lack of accuracy in the recursive residual (caused by rounding errors).,,,,,|AZ_loss|_
+:class:`~fipy.solvers.convergence.IteratingConvergence`,Solve still in progress.,|KSP_CONVERGED_ITERATING|_,,,,

docs/source/_static/solver_divergence.csv

@@ -0,0 +1,10 @@
+,,PETSc_,pyamgx_,PySparse_,SciPy_,Trilinos_
+:class:`~fipy.solvers.convergence.BreakdownDivergence`,Illegal input or the iterative solver has broken down.,|KSP_DIVERGED_BREAKDOWN|_,|AMGX_SOLVE_FAILED|_,,|SciPy_lt0|_,|AZ_breakdown|_
+:class:`~fipy.solvers.convergence.IterationDivergence`,Maximum number of iterations was reached.,|KSP_DIVERGED_ITS|_,|AMGX_SOLVE_DIVERGED|_,|PySparse_neg1|_,|SciPy_gt0|_,|AZ_maxits|_
+:class:`~fipy.solvers.convergence.PreconditioningDivergence`,The system involving the preconditioner was ill-conditioned.,|KSP_DIVERGED_PC_FAILED|_,,|PySparse_neg2|_,,
+:class:`~fipy.solvers.convergence.IllConditionedPreconditionerDivergence`,"An inner product of the form :math:`\mathbf{x}^T \mathsf{P}^{-1} \mathbf{x}` was not positive, so the preconditioning matrix :math:`\mathsf{P}` does not appear to be positive definite.",|KSP_DIVERGED_INDEFINITE_PC|_,,|PySparse_neg3|_,,
+:class:`~fipy.solvers.convergence.IllConditionedDivergence`,The matrix :math:`\mathsf{L}` appears to be ill-conditioned.,|KSP_DIVERGED_INDEFINITE_MAT|_,,|PySparse_neg4|_,,|AZ_ill_cond|_
+:class:`~fipy.solvers.convergence.StagnatedDivergence`,The method stagnated.,,,|PySparse_neg5|_,,
+:class:`~fipy.solvers.convergence.OutOfRangeDivergence`,A scalar quantity became too small or too large to continue computing.,|KSP_DIVERGED_NANORINF|_,,|PySparse_neg6|_,,
+:class:`~fipy.solvers.convergence.NullDivergence`,Breakdown when solving the Hessenberg system within GMRES.,|KSP_DIVERGED_NULL|_,,,,
+:class:`~fipy.solvers.convergence.ToleranceDivergence`,The residual norm increased by a factor of ``divtol``.,|KSP_DIVERGED_DTOL|_,,,,

docs/source/_static/widetable.css

@@ -0,0 +1,38 @@
+/* 
+ * Wide tables are unclickable, apparently because they're under the sidebar.
+ * 
+ * Adapted from:
+ * https://groups.google.com/g/sphinx-users/c/ZV8gGObVrQU
+ * https://blogs.perficient.com/2021/01/18/freezing-row-and-column-in-html-table-using-css/
+ */
+
+/* 
+ * Make wide table horizontally scrollable.
+ * This has significant UI issues.
+ */
+table.widescroll {
+    display: block;
+    overflow: auto;
+}
+
+table.widescroll th.stub {
+    left: 0;
+    position: sticky;
+    z-index: 50;
+}
+
+table.widescroll caption span.caption-text {
+    left: 0;
+    position: sticky;
+}
+
+/* Move wide table above whatever is blocking it. */
+table.wideshow {
+	/* 
+	 * "The `position: static` property prevents z-index from having an effect."
+	 * but where does `position: static` come from?
+	 */
+    position: relative;
+    z-index: 25;
+}
+

docs/source/conf.py

@@ -130,8 +130,8 @@
 # [@MadPhysicist](https://stackoverflow.com/users/2988730/mad-physicist)
 # [CC BY-SA 4.0](https://creativecommons.org/licenses/by-sa/4.0/)
 # https://stackoverflow.com/a/66182779
-napoleon_use_param = True
-# napoleon_preprocess_types = True
+# Corrected for https://github.com/sphinx-doc/sphinx/issues/10963
+napoleon_preprocess_types = True
 napoleon_type_aliases = {
     'array-like': ':term:`array-like <array_like>`',
     'array_like': ':term:`array_like`'
@@ -244,6 +244,8 @@
         \definecolor{bluish}{rgb}{0.216,0.188,0.533}
 
         \newcommand{\fipylogo}{\scalebox{10}{\rotatebox{4}{\textcolor{redish}{\( \varphi \)}}\kern-.70em\raisebox{-.15em}{\textcolor{bluish}{\( \pi\)}}}}
+
+        \usepackage{pdflscape}
     """,
     'printindex': r'\footnotesize\raggedright\printindex',
 }
@@ -332,6 +334,7 @@ def autodoc_skip_member(app, what, name, obj, skip, options):
 
 def setup(app):
     app.connect('autodoc-skip-member', autodoc_skip_member)
+    app.add_css_file('widetable.css')
 
 # lifted from astropy/astropy@e68ca1a1
 

examples/cahnHilliard/mesh2D.py

@@ -84,10 +84,7 @@
 
 >>> dexp = -5
 >>> elapsed = 0.
->>> if __name__ == "__main__":
-...     duration = 1000.
-... else:
-...     duration = 1000.
+>>> duration = 1000.
 
 >>> while elapsed < duration:
 ...     dt = min(100, numerix.exp(dexp))
@@ -118,4 +115,3 @@
     exec(fipy.tests.doctestPlus._getScript())
 
     input('finished')
-

examples/cahnHilliard/tanh1D.py

@@ -103,11 +103,11 @@
    single: LinearLUSolver
    single: DefaultSolver
 
->>> import fipy.solvers.solver
->>> if fipy.solvers.solver_suite in ['pysparse', 'pyamgx']:
+>>> from fipy.solvers import solver_suite
+>>> if solver_suite in ['pysparse']:
 ...     solver = LinearLUSolver(tolerance=1e-15, iterations=100)
 ... else:
-...     solver = DefaultSolver()
+...     solver = DefaultSolver(tolerance=1e-15, iterations=100)
 
 The solution to this 1D problem over an infinite domain is given by,
 

examples/convection/exponential1D/tri2D.py

@@ -34,8 +34,10 @@
 >>> eq = (DiffusionTerm(coeff=diffCoeff)
 ...       + ExponentialConvectionTerm(coeff=convCoeff))
 
->>> eq.solve(var = var,
-...          solver=DefaultAsymmetricSolver(iterations=10000))
+It should be possible to drive this steady-state solution to high accuracy.
+
+>>> solver = DefaultAsymmetricSolver(tolerance=1e-10, iterations=10000)
+>>> eq.solve(var=var, solver=solver)
 
 The analytical solution test for this problem is given by:
 
@@ -45,8 +47,8 @@
 >>> DD = 1. - numerix.exp(-convCoeff[axis] * L / diffCoeff)
 >>> analyticalArray = CC / DD
 
->>> print(var.allclose(analyticalArray, rtol = 1e-6, atol = 1e-6))
-1
+>>> print(var.allclose(analyticalArray, rtol=1e-6, atol=1e-6))
+True
 
 >>> if __name__ == '__main__':
 ...     viewer = Viewer(vars = var)

examples/convection/exponential1DSource/tri2D.py

@@ -40,7 +40,7 @@
 ...       - ExponentialConvectionTerm(coeff = convCoeff))
 
 >>> eq.solve(var=var,
-...          solver=DefaultAsymmetricSolver(tolerance=1.e-15, iterations=10000))
+...          solver=DefaultAsymmetricSolver(tolerance=1.e-8, iterations=10000))
 
 The analytical solution test for this problem is given by:
 

examples/convection/exponential2D/tri2D.py

@@ -27,7 +27,10 @@
 >>> eq = (DiffusionTerm(coeff=diffCoeff)
 ...       + ExponentialConvectionTerm(coeff=convCoeff))
 
->>> eq.solve(var = var)
+It should be possible to drive this steady-state solution to high accuracy.
+
+>>> solver = eq.getDefaultSolver(tolerance=1e-10)
+>>> eq.solve(var=var, solver=solver)
 
 The analytical solution test for this problem is given by:
 
@@ -36,8 +39,8 @@
 >>> CC = 1. - numerix.exp(-convCoeff[axis] * x / diffCoeff)
 >>> DD = 1. - numerix.exp(-convCoeff[axis] * L / diffCoeff)
 >>> analyticalArray = CC / DD
->>> print(var.allclose(analyticalArray, rtol = 1e-10, atol = 1e-10))
-1
+>>> print(var.allclose(analyticalArray, rtol=1e-10, atol=1e-10))
+True
 
 >>> if __name__ == '__main__':
 ...     viewer = Viewer(vars = var)

examples/diffusion/mesh1D.py

@@ -708,7 +708,7 @@
 and see that :math:`\phi` dissipates to the expected average value of 0.2 with
 reasonable accuracy.
 
->>> print(numerix.allclose(phi, 0.2, atol=1e-5))
+>>> print(numerix.allclose(phi, 0.2, atol=4e-5))
 True
 
 If we reset the initial condition
@@ -788,7 +788,13 @@
 >>> if __name__ == '__main__':
 ...     viewer.plot()
 
->>> (TransientTerm() == DiffusionTerm(D)).solve(var=phi, dt=1e6*dt)
+>>> eq = (TransientTerm() == DiffusionTerm(D))
+
+The initial residual is much larger than the norm of the right-hand-side
+vector, so we use `"initial"` tolerance scaling.
+
+>>> solver = eq.getDefaultSolver(criterion="initial")
+>>> eq.solve(var=phi, dt=1e6*dt, solver=solver)
 >>> if __name__ == '__main__':
 ...     viewer.plot()
 >>> from fipy import input

examples/diffusion/nthOrder/input4thOrder1D.py

@@ -15,7 +15,7 @@
 .. index::
    single: Grid1D
 
->>> from fipy import CellVariable, Grid1D, NthOrderBoundaryCondition, DiffusionTerm, Viewer, GeneralSolver
+>>> from fipy import CellVariable, Grid1D, NthOrderBoundaryCondition, DiffusionTerm, Viewer, LinearLUSolver
 
 >>> nx = 500
 >>> dx = L / nx
@@ -56,11 +56,7 @@
 
 >>> eq = DiffusionTerm(coeff=(1, 1)) == 0
 
->>> import fipy.solvers.solver
->>> if fipy.solvers.solver_suite  == 'petsc':
-...     solver = GeneralSolver(precon='lu')
-... else:
-...     solver = GeneralSolver()
+>>> solver = LinearLUSolver()
 
 We perform one implicit timestep to achieve steady state
 

examples/diffusion/steadyState/mesh1D/inputPeriodic.py

@@ -33,15 +33,21 @@
 solver has no fixed value and can become unstable.
 
 >>> eq = TransientTerm(coeff=1e-8) - DiffusionTerm()
->>> eq.solve(var=var, dt=1.)
+
+The initial residual is much larger than the norm of the right-hand-side
+vector, so we use `"initial"` tolerance scaling with a tolerance that will
+drive to an accurate solution.
+
+>>> solver = eq.getDefaultSolver(criterion="initial", tolerance=1e-8)
+>>> eq.solve(var=var, dt=1., solver=solver)
 
 >>> if __name__ == '__main__':
 ...     viewer.plot()
 
 The result of the calculation will be the average value over the domain.
 
 >>> print(var.allclose((valueLeft + valueRight) / 2., rtol = 1e-5))
-1
+True
 
 """
 from __future__ import unicode_literals

examples/diffusion/steadyState/mesh1D/tri2Dinput.py

@@ -13,9 +13,12 @@
 :func:`~fipy.meshes.factoryMeshes.Grid1D` object.
 
 Here, one time step is executed to implicitly find the steady state
-solution.
+solution.  We increase the solver tolerance from the default
+:math:`10^{-5}` in order to achieve a good solution.
 
-    >>> DiffusionTerm().solve(var)
+    >>> eq = DiffusionTerm()
+    >>> solver = eq.getDefaultSolver(tolerance=1e-8)
+    >>> eq.solve(var, solver=solver)
 
 To test the solution, the analytical result is required. The `x`
 coordinates from the mesh are gathered and the length of the domain,
@@ -30,7 +33,7 @@
 tolerance of `1e-10`.
 
     >>> print(var.allclose(analyticalArray))
-    1
+    True
 
 """
 from __future__ import print_function

examples/diffusion/steadyState/mesh20x20/tri2Dinput.py

@@ -1,16 +1,20 @@
 """
 
 This input file again solves a 2D diffusion problem on a triangular mesh.
+We increase the solver tolerance from the default :math:`10^{-5}` in order
+to achieve a good solution.
 
-    >>> DiffusionTerm().solve(var)
+    >>> eq = DiffusionTerm()
+    >>> solver = eq.getDefaultSolver(tolerance=1e-10)
+    >>> eq.solve(var, solver=solver)
 
 The result is again tested in the same way:
 
     >>> Lx = nx * dx
     >>> x = mesh.cellCenters[0]
     >>> analyticalArray = valueLeft + (valueRight - valueLeft) * x / Lx
     >>> print(var.allclose(analyticalArray, rtol = 1e-8))
-    1
+    True
 
 """
 from __future__ import unicode_literals

examples/diffusion/steadyState/mesh50x50/input.py

@@ -6,12 +6,14 @@
 
 The result is again tested in the same way:
 
-    >>> DiffusionTerm().solve(var)
+    >>> eq = DiffusionTerm()
+    >>> solver = eq.getDefaultSolver(tolerance=1e-10)
+    >>> eq.solve(var, solver=solver)
     >>> Lx = nx * dx
     >>> x = mesh.cellCenters[0]
     >>> analyticalArray = valueLeft + (valueRight - valueLeft) * x / Lx
     >>> print(var.allclose(analyticalArray, rtol = 1e-9))
-    1
+    True
 
 """
 from __future__ import unicode_literals

examples/diffusion/steadyState/mesh50x50/tri2Dinput.py

@@ -12,12 +12,14 @@
 
 The result is again tested in the same way:
 
-    >>> DiffusionTerm().solve(var)
+    >>> eq = DiffusionTerm()
+    >>> solver = eq.getDefaultSolver(tolerance=1e-10)
+    >>> eq.solve(var, solver=solver)
     >>> Lx = nx * dx
     >>> x = mesh.cellCenters[0]
     >>> analyticalArray = valueLeft + (valueRight - valueLeft) * x / Lx
     >>> print(var.allclose(analyticalArray, atol = 1e-7))
-    1
+    True
 
 """
 from __future__ import unicode_literals