Best practices for solving linear systems associated with adjoint equations

[josearnal]

Hello everyone,

I am interested in solving large sparse linear systems that arise from the adjoint equations of a CFD solver, and I am wondering what would be the best Trilinos stack to do this.

Physics
Consider a fully discrete nonlinear system of equations arising from the discretization of fluid-dynamics equations (e.g Euler, Navier-Stokes, MHD, etc.). This system can be represented by $\mathbf{R}(\mathbf{u}) = \mathbf{0}$, where $\mathbf{R}$ is the residual vector, and $\mathbf{u}$ is the solution vector. If we wish to compute the sensitivity of a scalar function $f(\mathbf{u})$ with respect to the solution vector, we must solve the following so-called adjoint equation

$\frac{\partial \mathbf{R}}{ \partial \mathbf{u}}^T \mathbf{\psi} = \frac{\partial f}{\partial \mathbf{u}}$

Numerics
We can therefore write this as $Ax = b$. In practice, I don't store $\frac{\partial \mathbf{R}}{ \partial \mathbf{u}}^T$, but instead I have a computer program that can compute the matrix-vector-product of $\frac{\partial \mathbf{R}}{ \partial \mathbf{u}}^T$ with an arbitrary vector. However, I do explicitly store a preconditioner matrix, M, which closely resembles $\frac{\partial \mathbf{R}}{ \partial \mathbf{u}}^T$ but is easier to invert. I want to use this matrix to do right preconditioning:

$A M^{-1} M x = b$.

Questions
My question is therefore, what would be the best Trilinos packages to do the following:

1. Build the matrix-vector operation $Ax$, without explicitly storing $A$
2. Build and store the preconditioner matrix $M$
3. Use 1. and 2. to iteratively solve $Ax = b$ with GMRES
4. Make use of additive Schwartz for further global preconditioning together with block incomplete lower-upper local preconditioning

Furthermore, are there any examples / tutorials / open source code that show these steps or at least subsets of them?

[jhux2]

@josearnal I converted your issue into a discussion.

[josearnal]

Thank you!

[cgcgcg]

Some questions about what you envision:

• For 1: Do I understand correctly that you already have methods for applying $A$ and $M^{-1}$, but you want to now wrap them to be used with Trilinos?
• For 4: what kind of interface do you have to build the subdomain operators?

[josearnal]

I’ll try to explain more clearly what I want to achieve.

Since I don't store $A$ explicitly (I only have a matrix-free method for computing $A x$), it’s challenging to build a preconditioner directly for $A$ as you point out. However, I do store a matrix $M$ that is an approximate representation of $A$ (i.e., $M \approx A$). My plan is to apply additive Schwarz preconditioning to $M$, rather than to $A$ directly.

More specifically, I want to construct a global preconditioner in the following way:

Ingredient 1: Construct a Lower-Fidelity Approximation $M \approx A$
Use a lower-order discretization to build the matrix $M$, which approximates $A$. This matrix is stored explicitly.

Ingredient 2: Apply Additive Schwarz Domain Decomposition
Decompose the global domain into overlapping subdomains. On each subdomain, restrict $M$ to form a local matrix $M_i$. The global preconditioner will then be expressed as:
$M^{-1} \approx \sum R_i^T \bar M_i^{-1} R_i$

where $R_i$ is the restriction operator for subdomain $i$, and the sum is taken over all subdomains (e.g., blocks in our CFD grid).

Ingredient 3: Block Incomplete LU (ILU) Factorization on Each Subdomain
Instead of fully inverting each local matrix $M_i$, I want to compute an approximate inversion using an incomplete LU factorization. For each subdomain $i$, we have:

$M_i \approx \bar M_i = L_i U_i,$

where $L_i$ and $U_i$ come from the block ILU factorization of $M_i$. The action of $(M_i)^{-1}$ is then approximated by performing forward and backward substitutions with $L_i$ and $U_i$, respectively.

Ingredient 4: Assemble the Global Preconditioner and Use it in GMRES
Combine the contributions from all subdomains to form the global preconditioner, and use this within an iterative solver like GMRES. Specifically, I would use right preconditioning, solving:
$A M^{-1} y = b$ , with $x=M^{-1} y,$

where the preconditioner $M^{-1}$ is applied using the additive Schwarz method with block ILU on each subdomain.

[cgcgcg]

Oh, you're essentially running a high-order matrix-free problem and you want to precondition using a low-order discretization on a refined mesh? Since you do have this discretization available in matrix form, Trilinos can build an additive Schwarz preconditioner for you. At the moment, we have additive Schwarz in two packages: Ifpack2 (one-level Schwarz only) and FROSch (multilevel Schwarz). I would recommend the following approach:

• Wrap your $A$ into a Tpetra::Operator.
• Construct a Tpetra::CrsMatrix for your $M$.
• Wrap both as Thyra operators and pass them into Stratimikos. This gives you access to all solvers/preconditioners.

It probably makes sense to first only work with $A$ and get a simple Krylov solver running. And then add the preconditioner.

[josearnal]

Awesome, thanks for your patience and help!

Would Ifpack2 or FROSch be responsible for doing the block ILU factorization? Or would that take place with something else?

[cgcgcg]

Yes, both of these will take care of the factorization for you. It's a simple parameterlist option, so you will have the choice between direct subdomain solves, incomplete factorizations and iterative methods with preconditioning.

[cgcgcg]

This might be a good starting point:
https://github.com/trilinos/Trilinos/blob/master/packages/stratimikos/example/tpetra/simple_tpetra_stratimikos_example.cpp

[josearnal]

Please let me know if there is a better avenue for me to post this followup question.

The way my code is structured right now, I know how to build M one column at a time. However, my understanding from this tutorial is that I should fill the CrsMatrix one row at a time. What would be the best practice in this case? Should I build the transpose, $M^T$ one row at a time and then manually transpose it (with this or this) afterwards to obtain $M$?

[cgcgcg]

In order to get started you could assemble $M^T$ and then transpose it. In the long run I would definitely recommend to directly assemble $M$ to avoid the cost of the transpose. You could also skip the explicit transpose and set up a preconditioner that uses transpose applies $x \mapsto M^T x$. That would involve writing a small wrapper.