Add domain-level default solver package (issue #130)

Addresses the 'odd construct' where the solver package in SolverOptions
affected the system matrix type, meaning a change in solver method could
inadvertently change the matrix backend and trigger an unnecessary rebuild.

New API:
  domain.setDefaultSolverPackage(SolverOptions.TRILINOS)  # or PASO/DEFAULT
  domain.getDefaultSolverPackage()

The domain default is applied when SolverOptions.package == SO_DEFAULT,
after the hard constraints (complex/direct → Trilinos) but before the
build-time fallback. An explicit package in SolverOptions still overrides
the domain default for individual solves.

Changed files:
- escriptcore/src/AbstractContinuousDomain.{h,cpp}: new member + getter/setter
- escriptcore/src/escriptcpp.cpp: Python binding
- finley/src/FinleyDomain.cpp: honour domain default in getSystemMatrixTypeId
- ripley/src/RipleyDomain.cpp: same
- oxley/src/OxleyDomain.cpp: same
- doc/user/linearPDE.tex: document new domain methods
- doc/user/trilinos.tex: document domain-level package selection

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

Author: Lutz Gross

doc/user/linearPDE.tex

@@ -594,6 +594,31 @@ \subsubsection*{class SolverOptions}
 \paragraph{SolverOptions.getPackage()}
 returns the solver package key.
 
+\paragraph{domain.setDefaultSolverPackage(package)}
+\label{DOMAIN DEFAULT PACKAGE}
+sets the default solver package for the domain.
+The value of \var{package} must be one of \member{SolverOptions.DEFAULT},
+\member{SolverOptions.PASO} or \member{SolverOptions.TRILINOS}.
+Setting the package at the domain level keeps the matrix type stable regardless
+of which solver method is selected per solve.  In particular, changing the
+iterative solver method (e.g.\ from \PCG to \GMRES) will never trigger an
+unnecessary rebuild of the system matrix.
+
+If \member{SolverOptions.DEFAULT} is used (the default), the package is chosen
+automatically: \TRILINOS is used for direct and complex-valued problems; \PASO
+is used otherwise (when available).
+
+Example usage:
+\begin{python}
+from esys.escript import SolverOptions
+from esys.finley import Rectangle
+domain = Rectangle(20, 20, 1)
+domain.setDefaultSolverPackage(SolverOptions.TRILINOS)
+\end{python}
+
+\paragraph{domain.getDefaultSolverPackage()}
+returns the default solver package set for the domain, or
+\member{SolverOptions.DEFAULT} if none has been explicitly set.
 
 \paragraph{SolverOptions.setTrilinosParameter(\optional{ name , value)}
 }

doc/user/trilinos.tex

@@ -138,12 +138,27 @@ \section{Example Problem}
 \end{itemize}
 Tolerance and solver output can be controlled by adding to the basic listing.
  \begin{python}[caption=tolerance , label=tolerance ]
-options = pde.getSolverOptions()  
-options.setTolerance(1e-8)        
-options.setVerbosityOn()        
+options = pde.getSolverOptions()
+options.setTolerance(1e-8)
+options.setVerbosityOn()
 \end{python}
 
+\subsection{Setting the Solver Package at the Domain Level}
+When both \PASO and \TRILINOS are available, the default iterative solver uses
+\PASO.  To use Trilinos for all solves on a given domain---regardless of the
+iterative method selected per solve---set the default solver package on the
+domain itself:
+\begin{python}
+from esys.escript import SolverOptions
+domain.setDefaultSolverPackage(SolverOptions.TRILINOS)
+\end{python}
+This fixes the solver backend for the domain so that changing the iterative
+method (e.g.\ from \PCG to \GMRES) never triggers an unnecessary rebuild of
+the system matrix.  The package set in \verb|SolverOptions.setPackage()| still
+overrides the domain default for individual solves.
 
+Use \verb|domain.getDefaultSolverPackage()| to query the current setting.
+See also Section~\ref{DOMAIN DEFAULT PACKAGE} in Chapter~\ref{SEC LinearPDE}.
 
 \section{Solver Options Summary}
 Table~\ref{tab:solver_summary} provides a quick reference for the most commonly used Trilinos solver options in escript.

escriptcore/src/AbstractContinuousDomain.cpp

@@ -141,11 +141,21 @@ int AbstractContinuousDomain::getSystemMatrixTypeId(const boost::python::object&
    return 0;
 }
 
-int AbstractContinuousDomain::getTransportTypeId(int solver, int precondioner, int package, bool symmetry) const 
+int AbstractContinuousDomain::getTransportTypeId(int solver, int precondioner, int package, bool symmetry) const
 {
    return 0;
 }
 
+void AbstractContinuousDomain::setDefaultSolverPackage(int package)
+{
+    m_defaultSolverPackage = package;
+}
+
+int AbstractContinuousDomain::getDefaultSolverPackage() const
+{
+    return m_defaultSolverPackage;
+}
+
 void AbstractContinuousDomain::addPDEToSystem(
                      AbstractSystemMatrix& mat, escript::Data& rhs,
                      const escript::Data& A, const escript::Data& B, const escript::Data& C, 

escriptcore/src/AbstractContinuousDomain.h

@@ -302,6 +302,22 @@ class ESCRIPT_DLL_API AbstractContinuousDomain : public AbstractDomain
                       const escript::FunctionSpace& functionspace,
                       const int type) const;
 
+  /**
+     \brief
+     Set the default solver package for this domain.
+     The package is used when the solver options leave the package at SO_DEFAULT.
+     Valid values are SO_DEFAULT (auto), SO_PACKAGE_PASO, SO_PACKAGE_TRILINOS.
+     This allows the solver backend to be fixed at the domain level rather than
+     per-solve, avoiding unexpected operator rebuilds when the solver method changes.
+  */
+  void setDefaultSolverPackage(int package);
+
+  /**
+     \brief
+     Returns the default solver package for this domain, or SO_DEFAULT if not set.
+  */
+  int getDefaultSolverPackage() const;
+
   /**
      \brief
      Return the number of data points summed across all MPI processes
@@ -326,6 +342,10 @@ class ESCRIPT_DLL_API AbstractContinuousDomain : public AbstractDomain
      \param full
   */
   virtual void Print_Mesh_Info(const bool full=false) const;
+
+ private:
+  /// Default solver package; SO_DEFAULT means use the build-time fallback.
+  int m_defaultSolverPackage = 0; // 0 == SO_DEFAULT
 };
 
 } // end of namespace

escriptcore/src/escriptcpp.cpp

@@ -270,6 +270,19 @@ args("options"),
         ":rtype: int")
        .def("getTransportTypeId",&escript::AbstractContinuousDomain::getTransportTypeId,
 args("solver", "preconditioner", "package", "symmetry"))
+       .def("setDefaultSolverPackage",&escript::AbstractContinuousDomain::setDefaultSolverPackage,
+args("package"),
+        "Set the default solver package for this domain.\n\n"
+        "The package is used when SolverOptions leaves the package at SO_DEFAULT. "
+        "Valid values are ``SolverOptions.DEFAULT``, ``SolverOptions.PASO``, "
+        "``SolverOptions.TRILINOS``. Setting the package at the domain level avoids "
+        "unexpected operator rebuilds when the solver method changes between solves.\n\n"
+        ":param package: solver package identifier\n"
+        ":type package: `int`")
+       .def("getDefaultSolverPackage",&escript::AbstractContinuousDomain::getDefaultSolverPackage,
+        "Returns the default solver package set for this domain, or ``SolverOptions.DEFAULT`` "
+        "if none has been set.\n\n"
+        ":rtype: `int`")
 
       .def("addPDEToSystem",&escript::AbstractContinuousDomain::addPDEToSystem,
 args("mat", "rhs","A", "B", "C", "D", "X", "Y", "d", "y", "d_contact", "y_contact", "d_dirac", "y_dirac"),

finley/src/FinleyDomain.cpp

@@ -1745,6 +1745,12 @@ int FinleyDomain::getSystemMatrixTypeId(const bp::object& options) const
         }
     }
 #endif
+    // honour the domain-level default before falling back to build-time defaults
+    if (package == escript::SO_DEFAULT) {
+        int domDefault = getDefaultSolverPackage();
+        if (domDefault != escript::SO_DEFAULT)
+            package = domDefault;
+    }
 #ifdef ESYS_HAVE_PASO
     if (package == escript::SO_DEFAULT)
         package = escript::SO_PACKAGE_PASO;

oxley/src/OxleyDomain.cpp

@@ -1127,6 +1127,12 @@ namespace oxley {
             }
         }
 #endif
+        // honour the domain-level default before falling back to build-time defaults
+        if (package == escript::SO_DEFAULT) {
+            int domDefault = getDefaultSolverPackage();
+            if (domDefault != escript::SO_DEFAULT)
+                package = domDefault;
+        }
 #ifdef ESYS_HAVE_PASO
         if (package == escript::SO_DEFAULT)
             package = escript::SO_PACKAGE_PASO;

ripley/src/RipleyDomain.cpp

@@ -947,6 +947,12 @@ int RipleyDomain::getSystemMatrixTypeId(const bp::object& options) const
         }
     }
 #endif
+    // honour the domain-level default before falling back to build-time defaults
+    if (package == escript::SO_DEFAULT) {
+        int domDefault = getDefaultSolverPackage();
+        if (domDefault != escript::SO_DEFAULT)
+            package = domDefault;
+    }
 #ifdef ESYS_HAVE_PASO
     if (package == escript::SO_DEFAULT)
         package = escript::SO_PACKAGE_PASO;