Form a BlockSparseMatrix based on block matrices from a different BlockSparseMatrix

[Tao Jin]

Dear all,

This is a similar question related to a post two years ago:

https://groups.google.com/g/dealii/c/1L6gLFtvTCo/m/frIyTDXbBwAJ

I am developing a fractional time stepping algorithm based on an operator split. Basically, I can obtain the following monolithc system after the finite element discretization:

                                   A   B   0                        u1                    b1

tangent matrix =     B   C   D  , solution =  u2,       RHS = b2

                                   E   0    F                        u3                    b3

where tangent matrix is a 3 by 3 BlockSparseMatrix. 

During a time interval, I want to use the first split operator

A  B * u1 = b1  to solve u1 and u2.

B  C    u2    b2

Then, plug the solved u1 and u2 into the second split operator

F * u3 = b3

where F contains the information of the newly solved u1 and u2.

My question is, is there an easy way now to extract the sparse matrices A, B, B, C to form a new BlockSparseMatrix, something like the following:

new BlockSparseMatrix =   tangent matrix.block(0, 0), tangent matrix.block(0, 1) 

                                                 tangent matrix.block(1, 0), tangent matrix.block(1, 1)

Thank you

Tao

[Wolfgang Bangerth]

On 4/15/23 22:31, Tao Jin wrote:
>
> During a time interval, I want to use the first split operator

> *A  B * u1 = b1 *to solve*u1 *and*u2.*
> *B  C    u2    b2*

> Then, plug the solved u1 and u2 into the second split operator

> *F * u3 = b3*
> where *F* contains the information of the newly solved *u1* and *u2*.

You will want to look at how step-32 solves the coupled problem. There, it is
a time-dependent one where we first want to solve for the first two variables,
and then the third.

If F depends on u1 and u2, you will have to assemble that matrix separate from
the A,B,C matrices. In that case, you might also want to consider the approach
of step-31.

step-31 and -32 solve the same problem. It might be useful to compare and
contrast.

Best
W.

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

[Tao Jin]

Dear Wolfgang,

Thank you so much for your inputs. Indeed, Step-31 is a good example for what I am looking for. However, I have a follow-up question regarding using the WorkStream class based on TBB during the assembly stage when two different dof_handlers are involved.

To be specific, in Step-31, two different dof_handlers are defined, one for the Stoke sub-problem and the other for the thermal sub-problem. Since the thermal sub-problem depends on the velocity field solved from the Stokes sub-problem, the following is done to ensure that the cell from the two dof_handlers are in sync:

auto cell = stokes_dof_handler.begin_active();

const auto endc = stokes_dof_handler.end();

auto temperature_cell = temperature_dof_handler.begin_active();

for (; cell != endc; ++cell, ++temperature_cell)

{

stokes_fe_values.reinit(cell);

temperature_fe_values.reinit(temperature_cell);

...

}

My question is, is there a way to involve two different dof_handlers when using the WorkStream class, something like the following:

    WorkStream::run(
      m_dof_handler_mechanical.active_cell_iterators(),
      [this](const typename DoFHandler<dim>::active_cell_iterator &cell,
             ScratchData_ASM_mechanical &                                     scratch,
             PerTaskData_ASM_mechanical &                                     data) {
        this->assemble_system_one_cell_mechanical(cell, scratch, data);
      },
      [this](const PerTaskData_ASM_mechanical &data) {
        this->m_constraints_mechanical.distribute_local_to_global(data.m_cell_matrix,
                                                                  data.m_cell_rhs,
                                                                  data.m_local_dof_indices,
                                                                  m_tangent_matrix_mechanical,
                                                                  m_system_rhs_mechanical);
      },
      scratch_data,
      per_task_data);

Thank you so much!

Best,

Tao

[Wolfgang Bangerth]

On 4/20/23 17:27, Tao Jin wrote:
>
> My question is, is there a way to involve two different dof_handlers when
> using the WorkStream class, something like the following:

What you are looking for is an object that looks like a *pair* of iterators
and when you call ++it on it, then it increments both elements of the pair.
The class SynchronousIterator does exactly this.

Alternatively, it is possible to convert an iterator into one DoFHandler to an
iterator into another DoFHandler, like in
typename DoFHandler::active_cell_iterator cell_stokes
= stokes_dof_handler.begin_active();

for (..., ++cell_stokes)
{
typename DoFHandler::active_cell_iterator
cell_T (cell_stokes->get_triangulation(),
cell_stokes->level(),
cell_stokes->index(),
&temperature_dof_handler);
...do something with both iterators...

[Tao Jin]

Dear Wolfgang,

Thank you so much for your guidance. Step-35 is the example where  "SynchronousIterators" is used.

""

    using IteratorTuple =
      std::tuple<typename DoFHandler<dim>::active_cell_iterator,
                 typename DoFHandler<dim>::active_cell_iterator>;
 
    using IteratorPair = SynchronousIterators<IteratorTuple>;

    ...

        WorkStream::run(
          IteratorPair(IteratorTuple(dof_handler_velocity.begin_active(),
                                                        dof_handler_pressure.begin_active())),
          IteratorPair(IteratorTuple(dof_handler_velocity.end(),
                                                        dof_handler_pressure.end())),
          *this,
          &NavierStokesProjection<dim>::assemble_one_cell_of_gradient,
          &NavierStokesProjection<dim>::copy_gradient_local_to_global,
          scratch_data,
          per_task_data);

"""

Best,

Tao