%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (c) 2003-2018 by The University of Queensland
% http://www.uq.edu.au
%
% Primary Business: Queensland, Australia
% Licensed under the Apache License, version 2.0
% http://www.apache.org/licenses/LICENSE-2.0
%
% Development until 2012 by Earth Systems Science Computational Center (ESSCC)
% Development 2012-2013 by School of Earth Sciences
% Development from 2014 by Centre for Geoscience Computing (GeoComp)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\chapter{Using Trilinos}
\label{TRILINOS}
Trilinos has a number of packages and provides a large collection of both direct and indirect solvers.  We refer the reader to \cite{TrilinosWeb} for details.  Escript needs to be installed with Trilinos to be able to use the Trilinos solvers.  See the install guide for details.  We show a few examples for the Trilinos options with a simple example.

Consider Laplacian in domain ($\Omega\in \mathbb{R}^3$) with a simple right hand side,
\begin{align}
 -\nabla^t\; \nabla u &= 1,  &&\text{ in } \Omega, \label{CONST1a}\\
 u &= 0, &&\text{ on } \Gamma_D,\label{BC1}\\
 \mathbf{n}^t \; \nabla u &= 0, &&\text{ on }\Gamma_N,\,\label{BC2}
\end{align} 
with $\mathbf{n}$ the outward normal, with $\Gamma_D$ the left boundary and $\Gamma=\partial\Omega\backslash\Gamma_D$.  Gravity forward weak form of PDE (\ref{CONST1a})-(\ref{BC2}), where $(~,~)$ is the standard $L^2$ inner product on $\Omega$, is
\begin{equation}\label{weak}
(\nabla u ~,~\nabla v ) = -(f ~,~v ),   
\end{equation}
for all admissible potential functions $v$. 


For this example, we just consider a simple, unstructured, 3D domain.  The mesh is created with GMSH using simplemesh.geo (file \ref{simplemesh}) see /escript/doc/examples/usersguide/simplemesh.geo.  To test AMG, a finer mesh is created using    /escript/doc/examples/usersguide/simplemeshfine.geo.

\begin{python}[caption=simplemesh.geo,label=simplemesh]
// dimensions and mesh size
xdim = 100.;
ydim = 200.;
zdim = 50.;
mtop = 2.;
mbase = 5.;

//Points
Point(1) = {0., 0., 0., mbase};
Point(2) = {xdim, 0., 0., mbase};
Point(3) = {0., ydim, 0., mbase};
Point(4) = {xdim, ydim, 0., mbase};
Point(5) = {0., 0., zdim, mtop};
Point(6) = {xdim, 0., zdim, mtop};
Point(7) = {0., ydim, zdim, mtop};
Point(8) = {xdim, ydim, zdim, mtop};

//Lines and surfaces
Line(1) = {1, 2};
Line(2) = {3, 4};
Line(3) = {1, 3};
Line(4) = {2, 4};
Line(5) = {5, 6};
Line(6) = {7, 8};
Line(7) = {5, 7};
Line(8) = {6, 8};
Line(9) = {1, 5};
Line(10) = {3, 7};
Line(11) = {2, 6};
Line(12) = {4, 8};
Line Loop(1) = {-1, 3, 2, -4};
Plane Surface(1) = {1};
Line Loop(2) = {5, 8, -6, -7};
Plane Surface(2) = {2};
Line Loop(3) = {1, 11, -5, -9};
Plane Surface(3) = {3};
Line Loop(4) = {-2, 10, 6, -12};
Plane Surface(4) = {4};
Line Loop(5) = {-3, 9, 7, -10};
Plane Surface(5) = {5};
Line Loop(6) = {4, 12, -8, -11};
Plane Surface(6) = {6};

// domain
Surface Loop(1) = {1:6};
Volume(1) = {1};
\end{python}

The most basic escript script to solve this pde using escript defaults is in program \ref{basic}.  The computed solution is saved in a silo file "asimple.silo" and can be visualized using \VisIt.
\begin{python}[caption=basic solve using defaults only, label=basic ]
from esys.escript import *
from esys.weipa import saveVTK, saveSilo
from esys.escript.linearPDEs import LinearSinglePDE
from esys.finley import ReadGmsh
  
domain=ReadGmsh("simplemesh.msh", 3,  optimize=True )
       
pde = LinearSinglePDE(domain, isComplex=False)
pde.setSymmetryOn()
x = domain.getX()
pde.setValue(A=kronecker(3), Y = 1., q = whereZero(x[0]-inf(x[0])))

u=pde.getSolution()    
saveSilo("asimple", u=u)    
\end{python}

This script uses default Trilinos solver PCG with Jacobi preconditioner. (It takes 271 iterations to reach the default tolerance of $10^{-8}$).  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()        
\end{python}



\section{Direct Solvers}
There are a number of options that can be set for solving the pde.  The simplest way to choose a direct solver is to use the default Trilinos direct solver KLU2.  Escript can only use this feature if it is available in Trilinos.

To choose the default Trilinos direct solver, set the tolerance for the PDE solver to $10^{-8}$,  use the Trilinos suite of solvers and output information about the solvers, we add code to listing \ref{basic}.
\begin{python}[caption=default Direct , label=defaultdirect ]
options = pde.getSolverOptions()  
options.setSolverMethod(SolverOptions.DIRECT)         
options.setPackage(SolverOptions.TRILINOS)
options.setVerbosityOn()        
\end{python}
The default Trilinos Direct solver is the Amesos2 LU factorisation KLU. It is a serial, unsymmetric sparse, partial-pivoting, direct matrix solver.  The last line above ensures that the output includes details of the methods used.


To use SUPERLU the code that needs to be added to listing \ref{basic} is 
\begin{python}[caption=SuperLU, label=superLU ]
options = pde.getSolverOptions()  
options.setTolerance(1e-8)        
options.setPackage(SolverOptions.TRILINOS)
options.setSolverMethod(SolverOptions.DIRECT_SUPERLU)
options.setVerbosityOn()        
\end{python}
Trilinos is the default solver package so the setPackage line is not really necessary.

\section{Preconditioned Conjugate Gradient}
If we want to use preconditioned conjugate gradient, then we remover the DIRECT solver line and add the line
\begin{python}[caption=Preconditioned conjugate gradient defaults, label=PCG]
options.setSolverMethod(SolverOptions.PCG)        
\end{python}
This takes 271 iterations to reach tolerance using a Jacobi preconditioner, BELOS Pseudo Block CG with Ifpack2.  To change the preconditioner to Gauss-Siedel, we add the line
\begin{python}
options.setPreconditioner(SolverOptions.GAUSS_SEIDEL)
\end{python}
This takes 120 iterations to reach tolerance.



\section{Multigrid}
Geometric multigrid methods were introduced for structured grids to maximise the advantages of iterative methods for solving matrix equations derived from discretised partial differential equations.  Iterative methods for these problems are effective in reducing oscillatory error but stall for smooth error and smooth error appears oscillatory when restricted to a coarser grid.  The idea is to have a succession of grids from fine to coarse and to remove error by iterating on each of the grids in turn reducing the oscillatory error on that grid.  The coarsest grid is chosen small enough so that it can be quickly solved directly or iteratively. If $n$ is the size of the problem then a multigrid algorithm is order $n$.

The matrix equation,  derived from the PDE, is
\begin{equation}\label{matrixEQ}
    \mathbf{A}\mathbf{u}=\mathbf{f},
\end{equation}
where $\mathbf{u}$ is the unknown $n\times 1$ vector, $\mathbf{A}$ is the $n\times n$ matrix and $\mathbf{f}$ is the nown right hand side $n\times 1$ vector.  The residual $\mathbf{r}$ is defined
\begin{equation}\label{res}
    \mathbf{r}=\mathbf{f}-\mathbf{A}\mathbf{u}%=\mathbf{A}\mathbf{u}^*-\mathbf{A}\mathbf{u}=\mathbf{A}\mathbf{u}^*-\mathbf{u}),
\end{equation}
and the residual equation is defined 
\begin{equation}\label{resEQ}
    \mathbf{A}\mathbf{e}=\mathbf{r},
\end{equation}
where $\mathbf{e}=\mathbf{u}^*-\mathbf{u}$ is the error and $\mathbf{u}^*$ is the exact solution.  Relaxing on \ref{resEQ}) with the residual as the right hand side is the same as relaxing on (\ref{matrixEQ}) with the original right and side.  We use superscripts on these terms to indicate the discretization representing element length, $h$, for a fine grid and $H$ for one level coarser discretization, $h<H$.  For a structured mesh, $H=\frac{h}{2}$ and $\mathbf{A}^h$ is an $n\times n$ matrix and $\mathbf{A}^H$ is an $N\times N$ matrix obtained in the coarsening process with $N<n$.  Prolongation operator, $\mathbf{P}^h_H$, an $n\times N$ matrix, is used to interpolate the error from the coarse grid to the fine grid and restriction operator $\mathbf{P}_h^H$, an $N\times n$ matrix, is used to restrict the residual from the fine grid to the coarse grid.  It is not necessary for the restriction operator to be the transpose of the interpolation operator.  The coarse grid matrix operator $\mathbf{A}^H$ is computed by multiplying the fine grid matrix operator $\mathbf{A}^h$ on the left by the restriction operator and the right by the interpolation operator resulting in an $N\times N$ matrix.  

The MG algorithm is best described using a recursion algorithm.  After iteration on a fine grid, the fine grid error is smooth and the residual can be restricted to a coarse grid. The error on the coarse grid appears more oscillatory and is (smoothed) reduced by iteration.  If this is the coarsest grid then the matrix equation is solved using a direct method or sufficient iterations of the solver to get within discretization error, otherwise the algorithm is called again with the coarse grid replacing the fine grid and the next coarser grid as the coarse grid.  The coarse grid error, is interpolated to the fine grid where it corrects the fine grid approximation and post-smoothing iteration smooths the error.  For structured grids, the choices for interpolating from a coarse grid to a fine grid or restricting from a fine grid to a coarse grid are reasonably obvious, see \cite{Briggs2000}. 


\begin{table}\center
\begin{tabular}{|l|l|ll}
\hline
AMG($\mathbf{u}^h;\mathbf{A}^h,\mathbf{f}^h)$ &\\ \hline
    \quad $\mathbf{A}^h\mathbf{u}^h=\mathbf{f}^h$ & \textbf{p steps pre-smoothing iteration}\\
	\quad $\mathbf{r}^h=\mathbf{f}^h-\mathbf{A}^h \mathbf{u}^h$ & \textbf{fine grid residual}\\
	\quad $\mathbf{r}^H=\mathbf{R}_h^H \mathbf{r}^h$ & \textbf{restrict residual} \\
	\quad $\mathbf{A}^H=\mathbf{R}_h^{H} \mathbf{A}^h \mathbf{P}_H^h$& \textbf{restrict operator}\\
	\quad if not coarsest  &\\
	  \quad\quad AMG($\mathbf{e}^H;\mathbf{A}^H,\mathbf{r}^H$) &\textbf{recursion}\\
	  \quad else &\\
	  \quad\quad $\mathbf{A}^H\mathbf{e}^H=\mathbf{r}^H$ & \textbf{coarsest solve}\\
	\quad$\mathbf{e}^h=\mathbf{P}_H^h \mathbf{e}^H$ & \textbf{interpolate error} \\
	\quad$\mathbf{u}^h=\mathbf{u}^h+\mathbf{e}^h$ & \textbf{fine grid correction}\\
	\quad$\mathbf{A}^h\mathbf{u}^h=\mathbf{f}^h$ & \textbf{q steps post-smoothing iteration}\\
	\quad return $\mathbf{u}^h$ &\\
\hline
\end{tabular}\caption{An AMG V(p.q) cycle algorithm to solve $\mathbf{A}\mathbf{u}=\mathbf{f}$, where $h$ and $H$ represent grid sizes with $h<H$.}
\end{table}


Algebraic multigrid methods were developed for unstructured grids and do not reference the grid but instead use interpolation and restriction operators derived from the matrix (see \cite{Briggs2000, Stuben2001281, Vanek1996, Tuminaro2000}).  The terms "coarse grid" and "fine grid" are still used but do not refer to actual grids.  "Smooth error" is defined to be the error not reduced by iteration and "oscillatory error" is the error reduced by iteration.  Coarse levels are chosen from the relative sizes of the off diagonal terms in the fine matrix.
Once the coarse grid is chosen the restriction and interpolation operators are computed. Restriction operators need to be chosen so that "smooth error" on a fine grid will appear "oscillatory" on a coarse grid ensuring that it can be reduced by iterating on this grid. There are a number of algorithm options available in Trilinos to compute the coarse grid and the choice will depend on the original PDE and smoothing options. For any multigrid method, there are basic choices: 
\begin{itemize}%[topsep=0pt,itemsep=-1ex,partopsep=1ex,parsep=1ex]
    \item pre-smoothing iterative solver and number of iterations 
    \item post-smoothing iterative solver and number of iterations
    \item choosing coarse grids    
    \item number of coarse grids or size of coarsest grid
    \item interpolation operator
    \item restriction operator
    \item coarsest grid solver
    \item cycle type 
\end{itemize}
We use Trilinos solvers and it is possible to access Trilinos options either within the escript script or, if more complicated control is needed, in an XML file. There are many Trilinos packages that can be used by escript including 
MueLu   - setup of AMG, Belos   - linear solvers - Pseudo Block CG, Ifpack2 - iterative solvers (Jacobi, Gauss-Siedel), Amesos2 - direct solvers for coarse level and Voltan or Voltan2 - repartitioning for caorse grids


\subsection{Default MueLu Trilinos options for algebraic multigrid preconditioned conjugate gradient (AMG-PCG)}
More detail on the various options can be found in the MUELU user guide and other Trilinos user guides \cite{TrilinosWeb}. MUELU uses other Trilinos packages and to access these parameters the XML file must be used.  
Recall $\mathbf{R}_h^H$ is the restriction operator and $\mathbf{P}^h_H$ is the interpolation (prolongation) operator.  The default values used are shown in Table \ref{defaultAMG}.  To use AMG-PCG with default paramerters we use script \ref{basicAMG}. 

\begin{table}\center
\begin{tabular}{|r|l|}
\hline
    number of equations & 1\\
    problem: symmetric & True, $\mathbf{R}_h^H = (\mathbf{P}^h_H)^t$ \\
    pre-smoothing iterative solver & Symmetric Gauss-Seidel\\
    post-smoothing iterative solver & Symmetric Gauss-Seidel\\
    pre-smoothing iterations & 1 \\
    post-smoothing iterations & 1 \\
    minimum aggregate size & 2 \\
    maximum aggregate size & unlimited \\
    aggregation & uncoupled \\
    maximum number of levels & 10\\
    maximum size of coarsest grid & 2000\\
    choosing coarse grids & classical smoothed aggregation\\
    coarsest grid solver &  SuperLU\\
    cycle type & V(1,1) \\
    \hline
\end{tabular}\caption{Default parameters for AMG-PCG}\label{defaultAMG}
\end{table}


\begin{python}[caption=basic Trilinos PCG-AMG defaults only script, label=basicAMG ]
from esys.escript import *
from esys.weipa import saveVTK, saveSilo
from esys.escript.linearPDEs import LinearSinglePDE
from esys.finley import ReadGmsh

domain=ReadGmsh("simplemesh.msh", 3,  optimize=True )

pde = LinearSinglePDE(domain, isComplex=False)
pde.setSymmetryOn()
x = domain.getX()
pde.setValue(A=kronecker(3), Y=1, q=whereZero(x[0]-inf(x[0])))

options = pde.getSolverOptions()
options.setPackage(SolverOptions.TRILINOS)
options.setSolverMethod(SolverOptions.PCG)
options.setPreconditioner(SolverOptions.AMG)

u=pde.getSolution()    
saveSilo("asimple",u=u)    
\end{python}

Only two grids were used for the simple mesh.  The first coarse grid, replaces, if possible, 27 fine grid nodes represented by  one coarse grid node.  In 2D this would be 9 fine grid nodes represented with one coarse grid node. The ratio of the fine to coarse grid in this example is 19.21.  For a finer mesh, with mtop = 2 and mbase=1 in the geo file, 3 levels of grids and the coarsest grid is solved with a direct method.  AMG-PCG took 11 iterations to reach tolerance.

\subsection{Altering MUELU parameters}
To access MUELU parameters in the python script we use the general form  
\begin{python}
options.setTrilinosParameter( "A", "B")        
\end{python}
where "A" is the Trilinos parameter and "B" is its string value. It is extremely important to have correct spaces in the strings.  

\subsubsection{Debug output}
Options are "none", "low", "medium", "high" and "extreme".
\begin{python}
options.setTrilinosParameter("verbosity", "low")  
\end{python}
Options are\\
\var{"low"} - setup time,\\
\var{"medium"} - basic AMG data, mesh sizes, smoothers + "low" \\
\var{"high"} - input data, relaxation solvers and data, aggregate data + "medium"\\
\var{"extreme"} - may include solver details that MueLu calls + "high "\\

\subsubsection{Problem type}
Changes default multigrid algorithm, block size and smoother.  
Options are 
\begin{itemize}
    \item \var{"unknown"}: default
    \item \var{"Poisson-2D" or "Poisson-3D"} : using smoothed aggregation, Chebyshev smoother and block size of 1,
    \item \var{"Elasticity-2D" and "Elasticity-3D"}: using smoothed aggregation, Chebyshev smoother and block size of 2 or 3 respectively,
    \item \var{"Poisson-2D-complex" and "Poisson-3D-complex"}: using smoothed aggregation, symmetric Gauss-Seidel and block size 1,
    \item \var{"Elasticity-2D-complex" and "Elasticity-3D-complex"}: using smoothed aggregation, symmetric Gauss-Seidel and block size 2 and 3 respectively,
    \item \var{"ConvectionDiffusion"}: using Petrov-Galerkin AMG, Gauss-Seidel and 1 block,
    \item \var{"MHD"}: using unsmoothed aggregation and Additive Schwarts method with one level of overlap and ILU(0) as a subdomain solver.
\end{itemize}
\begin{python}
options.setTrilinosParameter("problem:type", "Poisson-3D")    
\end{python}

\subsubsection{number of equations}
Number of PDE equations at each grid node.
\begin{python}
options.setTrilinosParameter("number of equations", 1)        
\end{python}

\subsubsection{AMG algorithm}
The multigrid algorithm for computing the coarse levels and interpolation and restriction operators is controlled with \var{"multigrid algorithm"}. The default value is smoothed aggregation and is selected with \var{"sa"} and a damping factor can be imposed.  The other options are \var{"unsmoothed"}, no Jacobi prolongation improvement step; \var{"pg"}, $\mathbf{A}$ prolongation smoothing and $\mathbf{A}^T$ restriction smoothing; \var{"emin"} basis functions for grid transfer using energy constrained minimisation; \var{"interp"}, piecewise constant ("interpolation order" set to 0) or linear interpolation (\var{"interpolation order"} set to 1) from coarse to fine and is only possible with structured aggregation; and \var{"semicoarsen"}, coarsen fully in z direction (will need to set rate in this direction).   It is also possible to use an implicit transpose for the restriction operator.
\begin{python}
# smoothed aggregation
options.setTrilinosParameter("multigrid algorithm", "sa")
options.setTrilinosParameter("sa: damping factor", 1.3)
options.setTrilinosParameter("sa: use filtered matrix", True)
options.setTrilinosParameter("filtered matrix: use lumping", True)
options.setTrilinosParameter("filtered matrix: reuse eigenvalue", True)
# unsmoothed
options.setTrilinosParamter("multigrid algorithm", "unsmoothed")
# pg
options.setTrilinosParameter("multigrid algorithm", "pg")
# interpolation
options.setTrilinosParameter("multigrid algorithm", "interp")
options.setTrilinosParameter("interp: interpolation order", 1)    
                                                          # 0, 1
options.setTrilinosParameter("interp: build coarse coordinates", True)
# emin
options.setTrilinosParameter("multigrid algorithm", "emin")
options.setTrilinosParameter("emin: iterative method", "cg") 
                                                     # "cg", "gmres", "sd"
options.setTrilinosParameter("emin: num iterations", 2)
options.setTrilinosParameter("emin: num reuse iterations", 1)
options.setTrilinosParameter("emin: pattern", "AkPtent")
options.setTrilinosParameter("emin: pattern order", 1)
# semicoarsen
options.setTrilinosParameter("multigrid algorithm", "semicoarsen")
options.setTrilinosParameter("semicoarsen: coarsen rate", 3)
#
options.setTrilinosParameter("transpose: use implicit", False) 
\end{python}

\subsubsection{Maximum levels, coarse mesh, coarse solver}
It is possible to limit the size of the coarsest level as well as limit the number of levels.  The default for the size of the coarsest level is 2000.  So once the size of the coarse level is less than 2000 then no more coarse levels are created.  Additionally, it is possible to limit the number of levels by setting "max levels" in the hierarchy.  This includes the fine grid.  The default value is 10 but depending on fine grid size, changing this could improve performance of the algorithm.  The coarsest level can be solved using a direct solver.  Possibilities are KLU, KLU2, SuperLU, SuperLU\_dist, Umfpack and Mumps

\begin{python}
options.setTrilinosParameter("max levels", 10)         
options.setTrilinosParameter("coarse: max size", 2000)
options.setTrilinosParameter("coarse: type", "SuperLU")
\end{python}


\subsubsection{Aggregation}
It is possible to influence aggregation options.  If the fine mesh is a structured grid then aggregates can be created in a  "structured" way and the aggregation attempts to form hexahedral coarse levels.  This uses a default coarsening rate of 3 in each direction. The option "hybrid" allows user determined "structured" or "unstructured" aggregation for each level, To get optimal size coarse mesh ($3^d$ in $d$ dimensions) "uncoupled" or "coupled" is used with "coupled" allowing aggregates to span processors.  It is suggested that "coupled" should be used with care.  "brick" attempts to make rectangular aggregates. Some of the options are below with more detail and more options in the MUELU user guide.   
\begin{python}
options.setTrilinosParameter("aggregation: type", "structured")
options.setTrilinosParameter("aggregation: ordering", "natural")
                                                 # "natural", "graph", "random"
options.setTrilinosParameter("aggregation: drop scheme", "classical")
                                                 # "classical", "distance laplacian" 
options.setTrilinosParameter("aggregation: drop tol", 0.0)
options.setTrilinosParameter("aggregation: min agg size", 2)
options.setTrilinosParameter("aggregation: max agg size", -1) 
                                                # -1 means unlimited    
options.setTrilinosParameter("aggregation: Dirichlet threshold", 1e-5)
\end{python}
    
\subsubsection{Relationship between $\mathbf{R}_h^H$  and $\mathbf{P}^h_H$}
For $\mathbf{R}_h^H=(\mathbf{P}^h_H)^t$
\begin{python}
options.setTrilinosParameter("problem: symmetric", True)
\end{python}
this is the default.

\subsubsection{Smoothers}
In the escript script it is possible to choose smoother type, "RELAXATION", "CHEBYSHEV" and "ILUT" or "RILUT" but for more specific control the XML file needs to be used.  It is possible to use different pre and post smoothers.  "RELAXATION" could use Jacobi, Gauss-Seidel, symmetric Gauss-Seidel, multithreaded Gauss-Seidel.  To specify which one the XML file must be used.  Some examples for this are 
\begin{python}
options.setTrilinosParameter("smoother: pre or post", "both")
options.setTrilinosParameter("smoother: type", "RELAXATION")
options.setTrilinosParameter("smoother: pre type", "CHEBYSHEV")
options.setTrilinosParameter("smoother: post type", "RELAXATION")
\end{python}

\subsubsection{Cycle type}
Allowable cycle types are "V" and "W".  The default is a "V" cycle.
\begin{python}
options.setTrilinosParameter("cycle type", "V")
\end{python}

\subsubsection{reuse}
If multiple PDEs are being solved the reuse strategy can use elements of previous computations.  The level of reuse varies from none to full.  Options are \var{"none"}; \var{"S"}, symbolic coarse levels information; \var{"tP"}, reuse tentative prolongation operator; \var{"emin"}, reuse old prolongator for initial guess; \var{"RP"}, reuse smoothed restrictor and prolongator; \var{"RAP"}, compute only fine level smoothers and reuse all other operators, and \var{"full"}, reuse everything.
\begin{python}
options.setTrilinosParameter("reuse: type", "full")
\end{python}

\subsubsection{repartitioning}
If there are multiple processors it might be benificial to repartition as the mesh are coarsened including perhaps using only one processor for the coarsest grid.  This is to reduce communication costs for the caorser grids.
\begin{python}
options.setTrilinosParameter("repartition: enable", False)
options.setTrilinosParameter("repartition: start level", 2)
options.setTrilinosParameter("repartition: min rows per proc", 800)
options.setTrilinosParameter("repartition: max imbalance", 1.2)
options.setTrilinosParameter("repartition: remap parts", True)
options.setTrilinosParameter("repartition: rebalance P and R", False)
\end{python}


\subsection{commands in XML file}
All the previous commands can be placed into an XML file.  The XML file option allows the user to choose parameters for the programs that MUELU calls, so more control is possible on the iterative solvers.
\begin{python}
from esys.escript import *
from esys.weipa import saveVTK, saveSilo
from esys.escript.linearPDEs import LinearSinglePDE
from esys.finley import ReadGmsh

domain=ReadGmsh("simplemesh.msh", 3,  optimize=True )

pde = LinearSinglePDE(domain, isComplex=False)
pde.setSymmetryOn()
x = domain.getX()
pde.setValue(A=kronecker(3), Y=1, q=whereZero(x[0]-inf(x[0])))

options = pde.getSolverOptions()
options.setPackage(SolverOptions.TRILINOS)
options.setSolverMethod(SolverOptions.PCG)
options.setPreconditioner(SolverOptions.AMG)
options.setTrilinosParameter("xml parameter file", "simplebob.xml")

u=pde.getSolution()    
saveSilo("asimple",u=u) 
\end{python}

It is possible to specify how many sweeps of the iterative solvers for pre and post smoothing and we could choose cycles with different numbers of pre and post sweeps.  Amesos2 provides direct solvers including superLU and Mumps. MueLu passes parameters directly to solver library.  To specify CHEBYSHEV parameters, for example, an XML file must be used.

Ifpack2 or Ifpack provides iterative matrix solvers Jacobi, Gauss Seidel, polynomial, distribution relaxation, domain decomposition solvers and incomplete factorizations.


A very simple XML file is in file \ref{simplebob}
\begin{python}[caption=simplebob.xml,label=simplebob]
<ParameterList name="MueLu"> 
  <Parameter name="verbosity"               type="string"    value="high"/> 
  <Parameter name="max levels"              type="int"       value="4"/>
  <Parameter name="coarse: max size"        type="int"       value="200"/>
  <Parameter name="multigrid algorithm"     type="string"    value="sa"/>
  <Parameter name="reuse: type"             type="string"    value="full"/>
  <Parameter name="transpose: use implicit" type="bool"      value="true"/>
  <Parameter name="sa: damping factor"      type="double"    value="0.1"/> 
  <Parameter name="sa: use filtered matrix" type="bool"      value="true"/>
</ParameterList>
\end{python}

A more complicated example that controls the number of pre and post sweeps is in the listing \ref{complicatedbob}
\begin{python}[caption=complicatedbob.xml,label=complicatedbob]
<ParameterList name="MueLu">
  <!--    General    -->
  <Parameter name="verbosity"               type="string"    value="high"/> 
  <Parameter name="max levels"              type="int"       value="4"/>
  <Parameter name="coarse: max size"        type="int"       value="200"/>
  <Parameter name="multigrid algorithm"     type="string"    value="sa"/>
  <Parameter name="reuse: type"             type="string"    value="full"/>
  <Parameter name="transpose: use implicit" type="bool"      value="true"/>
  <Parameter name="sa: damping factor"      type="double"    value="0.1"/> 
  <Parameter name="sa: use filtered matrix" type="bool"      value="true"/>

  <!-- Smoothing -->
  <Parameter name="smoother: pre or post"        type="string"  value="both"/>

  <Parameter name="smoother: pre type"           type="string"  value="CHEBYSHEV"/>
  <ParameterList name="smoother: pre params">
    <Parameter name="relaxation: type"           type="string"  value="Symmetric Gauss-Seidel"/>
    <Parameter name="relaxation: sweeps"         type="int"     value="5"/>
    <Parameter name="relaxation: damping factor" type="double"  value="0.9"/>
  </ParameterList>
  
  <ParameterList name="smoother: params">
    <Parameter name="chebyshev: degree"           type="int"     value="3"/>
    <Parameter name="chebyshev: ratio eigenvalue" type="double"  value="15"/>
  </ParameterList>

  <Parameter name="smoother: post type"           type="string"  value="RELAXATION"/>
  <ParameterList name="smoother: post params">
    <Parameter name="relaxation: type"           type="string"  value="Symmetric Gauss-Seidel"/>
    <Parameter name="relaxation: sweeps"         type="int"     value="5"/>
    <Parameter name="relaxation: damping factor" type="double"  value="0.9"/>
  </ParameterList>

  <!-- Aggregation -->
  <Parameter name="aggregation: type"           type="string"  value="uncoupled"/>
  <Parameter name="aggregation: min agg size"   type="int"     value="3"/>
  <Parameter name="aggregation: max agg size"   type="int"     value="27"/>

  <!--  for different level parameter list -->
  <ParameterList name="level 2">
    <Parameter name="smoother: type" type="string" value="CHEBYSHEV"/>
  </ParameterList>

</ParameterList>
\end{python}