construct TrilinosWrappers::SparseMatrix from TrilinosWrappers::BlockSparseMatrix to use Amesos direct solver

[Simon]

Dear all,

there are several discussions, e.g.

https://groups.google.com/g/dealii/c/KU7g_YCNVRE/m/h7ab-2gRCAAJ?hl=en

or

https://groups.google.com/g/dealii/c/DdMiTjuPx6Q/m/UeH4fmmBAwAJ?hl=en

about the usage of TrilinosWrappers::BlockSparseMatrix in the context of direct solvers

in a distributed framework.

Still, only TrilinosWrappers::SparseMatrix can be feeded to TrilinosWrappers::SparseDirect solvers. 

I am solving a 2x2 block system with iterative and direct solvers (decision is made at run-time via .prm file). 

Clearly, I do not need a BlockSparseMatrix for a direct solve, however,

to have a similar code base also for the iterative solve (where I need a block preconditioner), 

I decided to use block quantities for each solver type.

It is not a big deal to copy a BlockVector into a Vector, but what is an efficient scheme (my problem is non-linear) to

copy a BlockSparseMatrix into a SparseMatrix as needed for the direct solvers coming with trilinos?

Conceptually, I am looking for a constructor

SparseMatrix (const BlockSparseMatrix &)

or a copy_from function.

For instance, the std::memcpy function is used in
SparseMatrix::reinit (const Epetra_CrsMatrix & , copy_values = true) ,

but I am not sure whether it can be used when translating between Block and standard SparseMatrix.

Can you think of an elegant solution or must the non-zero values be copied via a loop

over the individual blocks?

Best regards,

Simon

[Wolfgang Bangerth]

On 12/18/23 04:47, Simon wrote:
>
> Can you think of an elegant solution or must the non-zero values be copied via
> a loop
> over the individual blocks?

This. I don't see a different way of doing things. Though I might just always
create both objects and then call distribute_local_to_global() with one or the
other matrix object, depending on your run-time flags.

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Simon]

"This. I don't see a different way of doing things. Though I might just always
create both objects and then call distribute_local_to_global() with one or the

other matrix object, depending on your run-time flags. "

Your approach sounds more reasonable than copying data back and forth at every Newton step

for a possibly large linear system.

I use a BlockSolver for the solution and right hand side of the linear system to compute the various norms (non-linear solver) 

of the the individual blocks.

For a direct solver, I would now use a SparseMatrix as matrix type and a BlockVector as vector type.

At each solve, I convert the BlockVector to a standard Vector, which I guess is not too expensive.

This mix (SparseMatrix and BlockVector) should not cause any undue problems, right?

Or would you not recommend this?

Best

Simon

[Wolfgang Bangerth]

On 12/18/23 12:44, Simon wrote:
>
> I use a BlockSolver for the solution and right hand side of the
> linear system to compute the various norms (non-linear solver) of the
> the individual blocks. For a direct solver, I would now use a
> SparseMatrix as matrix type and a BlockVector as vector type. At each
> solve, I convert the BlockVector to a standard Vector, which I guess
> is not too expensive. This mix (SparseMatrix and BlockVector) should
> not cause any undue problems, right? Or would you not recommend
> this?

No. This seems like a reasonable approach to me.
Best
W.

[Jan Philipp Thiele]

I have some experience with the direct solvers of Trilinos and it requires very few modifications to change the assembly to a SparseMatrix instead of a BlockSparseMatrix for vector valued problems.
Especially the assembly with distribute_local_to_global does not actually care apart from the matrix type provided to the method.

So you could try to work with some form of pointer, e.g. shared_ptr (which allows for checking initialization etc.) together with a switch for direct vs. iterative.
Then give both pointers to the assembly class and depending on the switch call distribute_local_to_global on one or the other.
This way only one matrix would actually need memory. Same would go for the vector classes.

The main difference is that you don't extract the views of the local and relevant dof index sets into a vector of sets

but just keep it as is in the sparsity pattern etc.

If you have any further questions on this don't hesitate to ask!

Best

Philipp

[Timo Heister]

What I have done in several codes (and what is used in ASPECT), is to always create block matrices and vectors but to put all DoFs into a single block when a direct solver is used. This decision can be made at rintime. Inside the solve routine you can just access block_matrix.block(0,0), which is a SparseMatrix. Same for the vector. This way, there are no different code paths, different types, or copies involved.

--

Timo Heister

http://www.math.clemson.edu/~heister/

---

From: dea...@googlegroups.com <dea...@googlegroups.com> on behalf of Wolfgang Bangerth <bang...@colostate.edu>
Sent: Monday, December 18, 2023 4:37:34 PM
To: dea...@googlegroups.com <dea...@googlegroups.com>
Subject: Re: [deal.II] construct TrilinosWrappers::SparseMatrix from TrilinosWrappers::BlockSparseMatrix to use Amesos direct solver

 

This Message Is From An External Sender: Use caution when opening links or attachments if you do not recognize the sender.

--
The deal.II project is located at https://nam12.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.dealii.org%2F&data=05%7C02%7Cheister%40clemson.edu%7Cb26c1543ea3f4a770d1a08dc00118fcd%7C0c9bf8f6ccad4b87818d49026938aa97%7C0%7C0%7C638385322953674220%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000%7C%7C%7C&sdata=FCQhyskq54Y%2B9Uc7oyauIeXMsI5CLzO12RE%2FCDG1w3Y%3D&reserved=0
For mailing list/forum options, see https://nam12.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgroups.google.com%2Fd%2Fforum%2Fdealii%3Fhl%3Den&data=05%7C02%7Cheister%40clemson.edu%7Cb26c1543ea3f4a770d1a08dc00118fcd%7C0c9bf8f6ccad4b87818d49026938aa97%7C0%7C0%7C638385322953674220%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000%7C%7C%7C&sdata=QKshKRkLcezx3WBxZW4FajigvudlynOE6GClmS9nNnY%3D&reserved=0
---
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
To view this discussion on the web visit https://nam12.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgroups.google.com%2Fd%2Fmsgid%2Fdealii%2F15994792-761a-4899-b7e6-af2fc6ba2932%2540colostate.edu&data=05%7C02%7Cheister%40clemson.edu%7Cb26c1543ea3f4a770d1a08dc00118fcd%7C0c9bf8f6ccad4b87818d49026938aa97%7C0%7C0%7C638385322953674220%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000%7C%7C%7C&sdata=Lmjub869SyEbNQ5DTAN4RNnzG41%2B2XEWpPs6s97c9OU%3D&reserved=0.

[Simon]

Timo,

" What I have done in several codes (and what is used in ASPECT), is to always create block matrices and vectors but to put all DoFs into a single block when a direct solver is used. This decision can be made at rintime. Inside the solve routine you can just access block_matrix.block(0,0), which is a SparseMatrix. Same for the vector. This way, there are no different code paths, different types, or copies involved. "

In terms of readability, this sounds like the best solution.

My motivation to use block structures even for direct solvers is to compute

norms of the pressure or the velocity, if we talk about the Stokes system.

Or in general, having access to the individual blocks to write them to files,...

If I am not mistaken, you create a nested finite element system with one block

FESystem<3>    feVelocity( FE_Q<3>(degree), 3 );

FE_DGQ<3>       fePressure(degree-1);

FESystem<3>    nestedFE ( feVelocity, 1, fePressure, 1 );

resulting in a 1x1 block system.

But it is not clear to me how to compute the norms of velocity and pressure in that case.

So far, I used a Block indices object as returned by DoFTools::count_dofs_per_fe_block()

since I compute the l2-norm only of the unconstrained DoFs.

That said, how do you access the velocity and pressure block when working with direct solvers?

Thank you,

Simon

[Wolfgang Bangerth]

On 12/20/23 04:35, Simon wrote:
> But it is not clear to me how to compute the norms of velocity and pressure in
> that case.
> So far, I used a Block indices object as returned by
> DoFTools::count_dofs_per_fe_block()
> since I compute the l2-norm only of the unconstrained DoFs.
>
> That said, how do you access the velocity and pressure block when working with
> direct solvers?

You can still renumber things based on what variable a DoF represents, so that
they occupy consecutive vector indices -- you just don't split the vector into
blocks.

But even if you didn't want to do that, you can always use
DoFTools::extract_dofs() to figure out which DoF corresponds to which variable.