Declare location of include headers
In [1]:
.I /opt/install/do-conf-ep-serial/include/
Out[1]:
Load Trilinos libraries
In [2]:
.L libepetra
Out[2]:
In [3]:
.L libxpetra
Out[3]:
In [4]:
.L libmuelu
Out[4]:
In [5]:
.L libteuchoscore
Out[5]:
Include some standard header files to create a communicator and an Xpetra::Map
In [6]:
#include "Kokkos_DefaultNode.hpp"
#include "Epetra_SerialComm.h"
#include "Epetra_Vector.h"
#include "Epetra_CrsMatrix.h"
#include "Teuchos_RCP.hpp"
#include "Teuchos_DefaultSerialComm.hpp"
#include "Xpetra_MapFactory.hpp"
#include "Xpetra_Map.hpp"
#include "Xpetra_EpetraCrsMatrix.hpp"
Out[6]:
Again, we use standard template parameters:
In [7]:
typedef double Scalar;
typedef int LocalOrdinal;
typedef int GlobalOrdinal;
typedef Kokkos::Compat::KokkosSerialWrapperNode EpetraNode;
typedef EpetraNode Node;
typedef Scalar SC;
typedef LocalOrdinal LO;
typedef GlobalOrdinal GO;
typedef Node NO;
Out[7]:
In [8]:
#include "Xpetra_CrsMatrixWrap.hpp"
#include "Xpetra_CrsMatrix.hpp"
#include "Xpetra_MultiVector.hpp"
#include "Xpetra_BlockedMultiVector.hpp"
#include "Xpetra_BlockedVector.hpp"
#include "Xpetra_MultiVectorFactory.hpp"
Out[8]:
Let's create a matrix and wrap it as an Xpetra::CrsMatrixWrap object
In [9]:
int NumMyElements = 100;
Out[9]:
In [10]:
Epetra_SerialComm Comm;
Epetra_Map Map(-1, NumMyElements, 0, Comm);
int NumGlobalElements = Map.NumGlobalElements();
Epetra_CrsMatrix A(Copy, Map, 3);
double negOne = -1.0;
double posTwo = 2.0;
for (int i=0; i<NumMyElements; i++) {
int GlobalRow = A.GRID(i); int RowLess1 = GlobalRow - 1; int RowPlus1 = GlobalRow + 1;
if (RowLess1!=-1) A.InsertGlobalValues(GlobalRow, 1, &negOne, &RowLess1);
if (RowPlus1!=NumGlobalElements) A.InsertGlobalValues(GlobalRow, 1, &negOne, &RowPlus1);
A.InsertGlobalValues(GlobalRow, 1, &posTwo, &GlobalRow);
};
A.FillComplete();
Out[10]:
In [11]:
Teuchos::RCP<Epetra_CrsMatrix> rcpA = Teuchos::rcpFromRef(A);
Xpetra::EpetraCrsMatrixT<GO,NO> xA(rcpA);
Teuchos::RCP<Xpetra::CrsMatrix<SC,LO,GO,NO> > xCrsMat = Teuchos::rcpFromRef(xA);
Xpetra::CrsMatrixWrap<SC,LO,GO,NO> xCrsWrap = Xpetra::CrsMatrixWrap<SC,LO,GO,NO>(xCrsMat);
Out[11]:
The xCrsWrap object contains the matrix object that can be used as input for MueLu.
First, include the necessary header files:
In [13]:
#include "MueLu_Level.hpp"
#include "MueLu_MutuallyExclusiveTime.hpp"
#include "MueLu_TimeMonitor.hpp"
#include "MueLu_Hierarchy.hpp"
Out[13]:
Then, we create a new Hierarchy object. We use the input matrix xCrsWrap as parameter. The constructor automatically generates the first (finest) level and stores the provided matrix object as variable A on the finest level.
Note: there is also an empty constructor that you can use for more special cases. You can add new levels using the ```AddLevel``` routine and fill them with data by hand.
In [14]:
MueLu::Hierarchy<SC,LO,GO,NO> H(Teuchos::rcpFromRef(xCrsWrap));
Out[14]:
We set the maximum size for the coarsest level to $10$ and increase the verbosity output level to get more screen output.
In [15]:
H.setDefaultVerbLevel(Teuchos::VERB_HIGH);
H.SetMaxCoarseSize(10);
Out[15]:
You can access the finest level that has been created in the constructor using the Hierarchy::GetLevel(int) function. If you do not provide the level number it assumes the finest level. The MueLu::Level class also has a verbosity level for the output details related to the Level data.
In [16]:
Teuchos::RCP<MueLu::Level> Finest = H.GetLevel();
Finest->setDefaultVerbLevel(Teuchos::VERB_HIGH);
Finest->print(std::cout, MueLu::Extreme);
Out[16]:
As you can see, the finest level already contains the fine level operator A provided by the user. This is the xCrsWrap operator you provided in the constructor for the Hierarchy.
Note: Instead of providing the ```Xpetra::Matrix``` (or ```Xpetra::CrsMatrixWrap``` object) in the constructor, you can also use the *empty* standard constructor and add the fine level operator by hand using a code line like ```Finest->Set("A",xCrsWrap);```.
The AmalgamationFactory takes the matrix A as input and amalgamates it. That is, if the matrix uses more than one DOFs per node, it builds the corresponding amalgamated node map and stores information on how to unamalgamate node-based information to dof-based information in the UnamalgamationInfo container.
In [17]:
#include "MueLu_AmalgamationFactory.hpp"
Out[17]:
In [18]:
MueLu::AmalgamationFactory<SC,LO,GO,NO> amalgFact;
Out[18]:
The CoalesceDropFactory actually generates the node-based graph given a dof-based matrix A and stores it as Graph in the Level container. It also accepts a threshold parameter for dropping small entries in A. It may use information from the UnAmalgamationInfo container.
In [19]:
#include "MueLu_CoalesceDropFactory.hpp"
Out[19]:
In [20]:
MueLu::CoalesceDropFactory<SC,LO,GO,NO> coalFact;
Out[20]:
In [21]:
Teuchos::RCP<MueLu::CoalesceDropFactory<SC,LO,GO,NO> > rcpCoalFact = Teuchos::rcpFromRef(coalFact);
Out[21]:
Note: In ```cling``` we create all factories on the stack and use ```Teuchos::rcpFromRef``` to define ```RCP``` pointers. We explicitly have to create an ```RCP``` pointer to a factory if we need it in several places.
The UncoupledAggregationFactory builds the aggregates given a Graph and the DofsPerNode variable. DofsPerNode is provided by the CoalesceDropFactory and contains how many DOFs per node the matrix has.
In [22]:
#include "MueLu_UncoupledAggregationFactory.hpp"
Out[22]:
In [23]:
MueLu::UncoupledAggregationFactory<LO,GO,NO> aggFact;
Out[23]:
There are several ways how to define data dependencies between factories. You can (and sometimes have to) directly connect factories: for example, to use Graph and DofsPerNode generated by the previously defined CoalesceDropFactory rcpCoalFact in the AggregationFactory aggFact you use:
In [24]:
aggFact.SetFactory("Graph", rcpCoalFact);
aggFact.SetFactory("DofsPerNode", rcpCoalFact);
Out[24]:
Note: In practice, it might be too complicated to define all necessary factory dependencies by hand. For example, most factories need ```A``` as input, which is generated on the previous level (by a *RAPFactory* object). *MueLu* provides a way to globally define which factory provides which data. That is, you can declare default factories which produce certain data. This functionality is provided by the *FactoryManager* which in detail will be introduced later.
The NullspaceFactory provides constant vectors as default near-nullspace modes for the finest level. On the coarser levels it (usually) just pipes through the Nullspaace which is produced by the TentativePFactory (see below)
In [25]:
#include "MueLu_NullspaceFactory.hpp"
Out[25]:
In [26]:
MueLu::NullspaceFactory<SC,LO,GO,NO> nspFact;
Out[26]:
In [27]:
Teuchos::RCP<MueLu::NullspaceFactory<SC,LO,GO,NO>> rcpNspFact = Teuchos::rcpFromRef(nspFact);
Out[27]:
The TentativePFactory build the tentative prolongation operator (with piecewise constant basis functions) using Aggregates and the fine level Nullspace. It not only calculates and stores P on the coarsest level, but also generates the coarse null space (variable Nullspace on the coarse level).
In [28]:
#include "MueLu_TentativePFactory.hpp"
Out[28]:
In [29]:
MueLu::TentativePFactory<SC,LO,GO,NO> PFact;
Out[29]:
For demonstration purposes we explicitly set aggFact and nspFact to be the generating factories for the Aggregates and the fine level Nullspace.
In [30]:
PFact.SetFactory("Aggregates", Teuchos::rcpFromRef(aggFact));
PFact.SetFactory("Nullspace", rcpNspFact);
Out[30]:
So, it uses the Nullspace provided by the NullspaceFactory (see above), which on the finest levels creates a default set of near-nullspace vectors or uses user-provided null-space vectors. Then, the TentativePFactory creates a new coarse Nullspace on the coarse level (which is supposed to be the fine Nullspace for the next level).
On the next coarser level, NullspaceFactory again is supposed to provide the fine level Nullspace (which actually has been generated by TentativePFactory). For doing so, NullspaceFactory just takes the coarse Nullspace from TentativePFactory and stores it as fine Nullspace generated by NullspaceFactory.
For doing so, we have to tell the NullspaceFactory that it is supposed to use the Nullspace from the TentativePFactory as input. The NullspaceFactory contains some code for the finest level (to generate the default nullspace) to break through the circle dependency.
In [31]:
Teuchos::RCP<MueLu::TentativePFactory<SC,LO,GO,NO>> rcpPFact = Teuchos::rcpFromRef(PFact);
Out[31]:
In [32]:
rcpNspFact->SetFactory("Nullspace", rcpPFact);
Out[32]:
Note: Here, for ```cling``` using the C++ interface, we can easily declare the circle dependency between *TentativePFactory* and *NullspaceFactory*. For a more general XML based input deck, we need a *FactoryManager*. Details will be discussed later.
The TransPFactory takes a prolongation operator P as input and stores its transposed as R
In [33]:
#include "MueLu_TransPFactory.hpp"
Out[33]:
In [34]:
MueLu::TransPFactory<SC,LO,GO,NO> RFact;
Out[34]:
In [35]:
RFact.SetFactory("P", Teuchos::rcpFromRef(PFact));
Out[35]:
The RAPFactory builds the triple Galerkin product to calculate the coarse level operator $A_c = RAP$. It takes the prolongation and restriction operators P and R as well as the fine level operator A as input and stores the coarse A on the coarse level.
In [36]:
#include "MueLu_RAPFactory.hpp"
Out[36]:
In [37]:
MueLu::RAPFactory<SC,LO,GO,NO> RAPFact;
Out[37]:
One can skip the following lines, if you have only one type of P and R whose associated factories are defined in a global FactoryManager (details later). Nevertheless, directly defining the dependencies does not hurt:
In [38]:
RAPFact.SetFactory("P", Teuchos::rcpFromRef(PFact));
Out[38]:
In [39]:
RAPFact.SetFactory("R", Teuchos::rcpFromRef(RFact));
Out[39]:
The CoarseMapFactory creates the coarse map which is used as domain map for the (tentative) prolongation operator. It is usually only used within the TentativePFactory and an average MueLu developer should never get into touch with it. Nevertheless, it is listed here for the sake of completeness of factories that you need somewhere when doing aggregation-based algebraic multigrid methods.
It uses the Aggregates and Nullspace as input.
Usually, you would not see the CoarseMapFactory at all, since it is usually used internally in the TentativePFactory. Nevertheless, it is there and plays an important role. In the usual case, the FactoryManager (details later) will automatically generate a default CoarseMapFactory for you that is used under the hood.
Note: The *FactoryManager* manages a list of (default) factories for all variables using ```RCP``` pointers. However, with ```cling``` it is not safe to have ```RCP``` objects within other objects, that may not live on the stack. The C++ interpreter may crash since ```RCP``` objects suddenly may have ran out of scope.
In [40]:
#include "MueLu_CoarseMapFactory.hpp"
Out[40]:
In [41]:
MueLu::CoarseMapFactory<SC,LO,GO,NO> coarseMapFact;
Out[41]:
Again, just for demonstration purposes, we directly connect the nspFact and coarseMapFact using a call to SetFactory
In [42]:
coarseMapFact.SetFactory("Nullspace", rcpNspFact);
Out[42]:
Last, but not least, the SmootherFactory handles level smoothers. More details about smoothers in MueLu can be found in a separate notebook.
In [43]:
#include "MueLu_SmootherFactory.hpp"
Out[43]:
In [44]:
MueLu::SmootherFactory<SC,LO,GO,NO> smootherFact;
Out[44]:
In the previous section, we briefly introduced the minimal set of necessary factories to define an aggregation-based algebraic multigrid setup. We have seen, that one can define direct data dependencies using the Factory::SetFactory(VarName, rcpFactory) routine. Calling this routine defines VarName provided by rcpFactory to be used by the current factory.
It is hard to define all data dependencies this way. In addition to SetFactory we can also globally define certain factories to provide certain data. The FactoryManager manages a list of classes associated to certain data. Unless directly defined, a Factory uses then the information from the FactoryManager to find the factory supposed to provide requested data.
First, we create an empty FactoryManager object using
In [45]:
#include "MueLu_FactoryManager.hpp"
Out[45]:
In [46]:
MueLu::FactoryManager<SC,LO,GO,NO> M;
Out[46]:
Then, for each variable name, we declare the (default) factory, which is providing the variable.
In [47]:
M.SetFactory("Smoother", Teuchos::rcpFromRef(smootherFact));
M.SetFactory("CoarseSolver", Teuchos::rcpFromRef(smootherFact));
M.SetFactory("A", Teuchos::rcpFromRef(RAPFact));
M.SetFactory("P", rcpPFact);
M.SetFactory("R", Teuchos::rcpFromRef(RFact));
M.SetFactory("DofsPerNode", rcpCoalFact);
M.SetFactory("Graph", rcpCoalFact);
M.SetFactory("UnAmalgamationInfo", Teuchos::rcpFromRef(amalgFact));
M.SetFactory("Nullspace", rcpPFact);
M.SetFactory("Aggregates", Teuchos::rcpFromRef(aggFact));
M.SetFactory("CoarseMap", Teuchos::rcpFromRef(coarseMapFact));
Out[47]:
For example: if not otherwise provided, the variable Nullspace is provided by the TentativePFactory PFact. Or, even more important, variable A is calculated by the RAPFactory. Here, A primarily describes the coarse matrices A, since A is provided by the user on the finest level.
You may notice, that some factories (like the CoalesceDropFactory) produce several different variables (such as Graph and DofsPerNode).
For simple standard cases, one can get rid of most of the Factory::SetFactory calls. It should be sufficient to just declare all factories in the FactoryManager in one central place. Vice versa, making sure, that all direct dependencies are set using Factory::SetFactory, one can reduce the number of entries in the FactoryManager.
One can ask the FactoryManager to return the associated (default) factory for a given variable using, e.g.
In [48]:
M.GetFactory("A");
Out[48]:
If the variable is not associated with a factory (through a FactoryManager::SetFactory call), the FactoryManager class will try to return hard-coded default factories (for most variable names).
Note: With ```cling``` this is technically not working. That is, in this notebook we have to explicitly creating instances of all the factories we need on the stack. We cannot use auto-generated default factories from the *FactoryManager*
There are two more rather special factories, which a MueLu user and a standard MueLu developer should not see. But they are there and one should be aware of their existence:
It is clear, that most multigrid methods consist of smoothers and transfer operators. Therefore, the variables Smoother and CoarseSolver as well as P and R are special in the sense, that you not only create them during the Setup phase, but you also make heavily use of them during the iteration phase. That is, it makes sense, to keep these data persistent and simplify access to it. This is done by the TopRAPFactory and TopSmootherFactory. These factories make sure, that above variables can be easily accessed on all levels without providing the generating factories. I.e., on all levels you just can use level.Get<Xpetra::Matrix<SC,LO,GO,NO>>("A"); and you will just get the matrix (same for Smoother...)
In [49]:
Teuchos::RCP<MueLu::FactoryManager<SC,LO,GO,NO>> rcpM = Teuchos::rcpFromRef(M);
#include "MueLu_TopSmootherFactory.hpp"
MueLu::TopSmootherFactory<SC,LO,GO,NO> ts(rcpM,"Smoother");
Out[49]:
In [50]:
MueLu::TopSmootherFactory<SC,LO,GO,NO> tsc(rcpM,"CoarseSolver");
Out[50]:
Let's create a TopRAPFactory object:
In [52]:
#include "MueLu_TopRAPFactory.hpp"
Out[52]:
In [53]:
MueLu::TopRAPFactory<SC,LO,GO,NO> tr(rcpM);
Out[53]:
Again, you should never really need to use these factories directly (unless you are modifying the Hierarchy in special ways).
With the FactoryManager set up using all the previously defined Factories, we finally can call the Setup routine to generate the multigrid hierarchy.
For example, with Hierarchy::Setup(M,0,1); we create a 1-level method using the FactoryManager M.
In [54]:
H.Setup(M,0,1);
Out[54]:
Obviously, we built a one-level method with no smoother. We can build a two-level method (starting level = 0, ending level = 1) using:
In [55]:
H.Setup(M,0,2);
Out[55]:
We got a two-level method with no smoothers. Of course, without level smoothers this would be a very inefficient multigrid method.
Why does it not create level smoothers even though we extra defined a SmootherFactory in the FactoryManager?
Let's add level smoothers:
In [58]:
#include "Teuchos_ParameterList.hpp"
Out[58]:
In [59]:
Teuchos::ParameterList smooParams;
Out[59]:
In [60]:
smooParams.set("type", Teuchos::ParameterEntry("RELAXATION"));
Teuchos::ParameterList& subSmooParams = smooParams.sublist("ParamterList");
subSmooParams.set("relaxation: type", "Symmetric Gauss-Seidel");
subSmooParams.set("relaxation: sweeps",1);
subSmooParams.set("relaxation: damping factor",1.0);
Out[60]:
In [61]:
std::cout << smooParams << std::endl;
Out[61]:
In [62]:
#include "MueLu_TrilinosSmoother.hpp"
Out[62]:
In [63]:
MueLu::TrilinosSmoother<SC,LO,GO,NO> trSmoother("RELAXATION", subSmooParams);
Out[63]:
In [64]:
PFact.DisableMultipleCheckGlobally();
Out[64]:
In [65]:
smootherFact.SetSmootherPrototypes(Teuchos::rcpFromRef(trSmoother));
Out[65]:
In [66]:
H.Setup(M,0,2);
Out[66]:
In [67]:
H.Setup(M,0,3);
Out[67]:
In [ ]: