MueLu: pass blocked nullspace vectors through "user data" list

[mayrmt]

Enhancement

@trilinos/muelu

Background

When applying MueLu to a blocked matrix, the nullspace can be provided for every block. This is done via "Nullspace1", "Nullspace2", ..., "Nullspace9".

In contrast to the regular "Nullspace" for non-blocked matrices, the individual block nullspaces cannot be passed through the parameter list "user data", but have to be set. on the fine level of the hierarchy manually via

RCP<Hierarchy> hierarchy = ...;
RCP<MultiVector> nspVector1 = ...;

hierarchy->GetLevel(0)->Set("Nullspace1", nspVector1);

It is desired to pass them through the "user data" parameter list in order to enable the use of MueLu::CreateXpetraPreconditioner().

Implementation Details

Passing data through the "user data" parameter list requires the handling of non-serialisable data in MueLu::HierarchyUtils::AddNonSerializableDataToHierarchy() and MueLu::Utilities::ExtractNonSerializableData(). Right now, adding new user data requires to manually keep two long if-clauses in sync (spread across two classes). This is tedious and error prone.

Idea/question: Can we define a list/struct/... of non-serializable data a a central place and then just loop over that list/struct/... to extract/add this data where necessary. This would remove the task of coherent maintenance of one list a two places. Yet, I'm no sure about the technical details. Any ideas are welcome!

Definition of Done

• [ ] allow to pass the nullspace for each block through the "user data" parameter list.
• [ ] maybe rework the design of this mechanism as briefly sketched out above.

[mayrmt]

Another bit of code in MueLu_NullspaceFactory_def.hpp to think about:

    // TODO not very elegant.
    // 1/20/2016: we could add a sublist (e.g. "Nullspaces" which is excluded from parameter validation)
    validParamList->set< RCP<const FactoryBase> >("Nullspace1", Teuchos::null, "Generating factory of the fine level null space associated with the (1,1) block in your n x n block matrix.");
    validParamList->set< RCP<const FactoryBase> >("Nullspace2", Teuchos::null, "Generating factory of the fine level null space associated with the (2,2) block in your n x n block matrix.");
    validParamList->set< RCP<const FactoryBase> >("Nullspace3", Teuchos::null, "Generating factory of the fine level null space associated with the (3,3) block in your n x n block matrix.");
    validParamList->set< RCP<const FactoryBase> >("Nullspace4", Teuchos::null, "Generating factory of the fine level null space associated with the (4,4) block in your n x n block matrix.");
    validParamList->set< RCP<const FactoryBase> >("Nullspace5", Teuchos::null, "Generating factory of the fine level null space associated with the (5,5) block in your n x n block matrix.");
    validParamList->set< RCP<const FactoryBase> >("Nullspace6", Teuchos::null, "Generating factory of the fine level null space associated with the (6,6) block in your n x n block matrix.");
    validParamList->set< RCP<const FactoryBase> >("Nullspace7", Teuchos::null, "Generating factory of the fine level null space associated with the (7,7) block in your n x n block matrix.");
    validParamList->set< RCP<const FactoryBase> >("Nullspace8", Teuchos::null, "Generating factory of the fine level null space associated with the (8,8) block in your n x n block matrix.");
    validParamList->set< RCP<const FactoryBase> >("Nullspace9", Teuchos::null, "Generating factory of the fine level null space associated with the (9,9) block in your n x n block matrix.");

[csiefer2]

I never liked the Nullspace1...Nullspace9 goop...

[mayrmt]

@trilinos/muelu @maxfirmbach

For blocked operators, we just pass the operator via

hierarchy->GetLevel(0)->Set("A", blockedOperator)

to the fine level and then use the SubBlockAFactory to extract individual blocks. Can't we follow a similar concept here as well?

In detail:

• create the null space vectors for each block
• create a blocked-type multi vector that just collects all multi vectors from the individual blocks
• use a SubBlockNullspaceFactory (or SubBlockVectorFactory if we want to be more generic) to extract the null space of any block

We would then pass the null space via

hierarchy->GetLevel(0)->Set("Nullspace", blockedNullspace)

and define the SubBlockNullspaceFactory objects in the xml-file/parameter list.

To account for different numbers of null space vectors in each block, maybe let's create the global null space MultiVector with max(nullspaceVectors) columns and fill unused columns in a block with less null space vectors with zeros, e.g. for 2D incompressible flow:

        | 1  0 |
        | 0  1 |
        | 1  0 |
        | 0  1 | 
        | .  . |
        | .  . |
        | -  - |
nsp =   | 1  0 |
        | 1  0 |
        | .  . |
        | .  . |

The number of columns to be extracted from nsp for a given block i can be specified on the parameter list of the SubBlockNullspaceFactory.

What do you think?

[jhux2]

@mayrmt Presumably a subblock's nullspace multivector (MV) would have the length of that subblock's domain map. Would an array mvArray of MVs with length equal to the number of subblocks be suitable for your needs? For subblock i, mvArray[i] is the corresponding nullspace. Maybe this is too simplistic, and I don't fully understand what's required.

[mayrmt]

@jhux2 That sounds like a much simpler, but totally sufficient approach! I just want to pass the nullspace "all at once" on the fine level and need some extraction mechanism. During coarsening, each block does take care of its nullspace independent form the other blocks (which is already part of MueLu).

I'll look into your idea a bit more, but it sounds promising and should help us to remove the manual and tedious specification of "nullspace1", "nullspace2", ..., "nullspace9". Thanks!

/cc @maxfirmbach

[github-actions[bot]]

This issue has had no activity for 365 days and is marked for closure. It will be closed after an additional 30 days of inactivity. If you would like to keep this issue open please add a comment and/or remove the MARKED_FOR_CLOSURE label. If this issue should be kept open even with no activity beyond the time limits you can add the label DO_NOT_AUTOCLOSE. If it is ok for this issue to be closed, feel free to go ahead and close it. Please do not add any comments or change any labels or otherwise touch this issue unless your intention is to reset the inactivity counter for an additional year.

[mayrmt]

Let's keep this issue alive for a little longer...

[nasseralkmim]

This feature would also be useful for users of Stratimikos/Thyra interface. It is straightforward to modify the parameter list and pass that to Stratimikos::DefaultLinearSolverBuilder.

However, it is not clear how to access the MueLu hierarchy from the Thyra::LinearOpWithSolveFactoryBase. Does anyone has an idea on how to access the MueLu hierarchy from the Thyra factory?

[cgcgcg]

@nasseralkmim Have a look at how the reuse code path of MueLu's Stratimikos adapter works: https://github.com/trilinos/Trilinos/blob/4efc8d8ba3fa49a497ba92c5b0e8611f9ee3f702/packages/muelu/adapters/stratimikos/Thyra_MueLuPreconditionerFactory_def.hpp#L612-L613 https://github.com/trilinos/Trilinos/blob/4efc8d8ba3fa49a497ba92c5b0e8611f9ee3f702/packages/muelu/adapters/stratimikos/Thyra_MueLuPreconditionerFactory_def.hpp#L643-L645

[nasseralkmim]

Hi @cgcgcg, I made an attempt but still rough and not fully functional. Here is a description of what I did (using Trilinos 13.4.0):

With the LOWS-factory, we can get the preconditioner factory. With this factory, I create a new preconditioner object.

RCP<Thyra::PreconditionerFactoryBase<Scalar>> precFact = lowsFactory->getPreconditionerFactory();
RCP<Thyra::PreconditionerBase<Scalar>> thyra_prec = precFact->createPrec();

Then, I need to initialize this preconditioner (otherwise I can not extract the Xpetra operator[^1]).

Thyra::initializePrec<Scalar>(*precFact, fwdOp, thyra_prec.ptr());

But, I have not passed any information to the hierarchy yet. This results in the MueLu has been successfully set up. Now, I perform the multiple cast like in the file you showed, to finally arrive at the hierarchy:

RCP<Thyra::LinearOpBase<Scalar>> thyra_precOp = thyra_prec->getNonconstUnspecifiedPrecOp();

using XpLinOp = Thyra::XpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, DeviceNode>;
RCP<XpLinOp> thyra_xp_precOp = Teuchos::rcp_dynamic_cast<XpLinOp>(thyra_precOp, true);

using MueLuXpOp = MueLu::XpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, DeviceNode>;
RCP<MueLuXpOp> xp_precOp =
Teuchos::rcp_dynamic_cast<MueLuXpOp>(thyra_xp_precOp->getXpetraOperator(), true);

using Hierarchy = MueLu::Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, DeviceNode>;
RCP<Hierarchy> hierarchy = Teuchos::rcp_dynamic_cast<Hierarchy>(xp_precOp->GetHierarchy());

With this hierarchy object, I set the "Nullspace1", information. I assume that I have to initialize the preconditioner again, so it can be set up with the new information. Is that correct?

In the end, I initialize the preconditioned operator and solve the system,

th_invA = lowsFactory->createOp();
Thyra::initializePreconditionedOp<Scalar>(*lowsFactory, fwdOp, thyra_prec, th_invA.ptr());

For some reason, it is not registering the null space information that I just passed. It just set up the hierarchy twice. Do you have any ideas to improve on that?

[^1] Here is the debugger backtrace,

#0  0x0000555555936492 in Teuchos::RCP<Xpetra::Operator<double, int, int, Kokkos::Compat::KokkosDeviceWrapperNode<Kokkos::Serial, Kokkos::HostSpace> > const>::get (this=0x38) at /Trilinos/build/install/include/Teuchos_RCP.hpp:391
#1  0x0000555555939dfe in Teuchos::ConstNonconstObjectContainer<Xpetra::Operator<double, int, int, Kokkos::Compat::KokkosDeviceWrapperNode<Kokkos::Serial, Kokkos::HostSpace> > >::getNonconstObj (this=0x38) at /Trilinos/build/install/include/Teuchos_ConstNonconstObjectContainer.hpp:327
#2  0x000055555592977c in Thyra::XpetraLinearOp<double, int, int, Kokkos::Compat::KokkosDeviceWrapperNode<Kokkos::Serial, Kokkos::HostSpace> >::getXpetraOperator (this=0x0) at /Trilinos/build/install/include/Thyra_XpetraLinearOp_def.hpp:94

[cgcgcg]

Have another look further down in the reuse code path: https://github.com/trilinos/Trilinos/blob/4efc8d8ba3fa49a497ba92c5b0e8611f9ee3f702/packages/muelu/adapters/stratimikos/Thyra_MueLuPreconditionerFactory_def.hpp#L642-L676 In particular, you might need to call SetupRe on the hierarchy.

But let me back up a bit: what are you actually trying to achieve? Modify an already set up Stratimikos-MueLu preconditioner? Or is this purely a question of how to get a blocked nullspace through the Stratimikos interface to MueLu?

[nasseralkmim]

In particular, you might need to call SetupRe on the hierarchy.

Thanks, I will probably need to setup the hierarchy from scratch and not rely on Thyra::initializePrec.

But let me back up a bit: what are you actually trying to achieve? Modify an already set up Stratimikos-MueLu preconditioner? Or is this purely a question of how to get a blocked nullspace through the Stratimikos interface to MueLu?

Yes, it is a little of both. Ideally, I want to pass the null space before the Stratimikos-MueLu preconditioner is initialized.

Since it is not possible to pass it via the "user data" in the parameter list, I'm attempting to reuse the information already in the factory (Thyra::PreconditionerFactoryBase) to recreate the preconditioner manually with the null space information passed to the hierarchy.

[cgcgcg]

Ok. I'd think that fixing this and allow to pass the block nullspaces via parameterlist is actually quite simple.