It cames out it is not so easy to provide simplest working example. I've
tried to remake step-37, but surprisingly I was not able to reproduce the
problem. I'll take a closer look what is causing a problem.
W dniu wtorek, 10 marca 2020 07:46:11 UTC+1 użytkownik peterrum napisał:
>
> Hi Michal,
>
> any chance that you post or send me a small runnable test program.
>
> By the way, there is an open PR (
> https://github.com/dealii/dealii/pull/9343) for the computation of the
> diagonal in a matrix-free manner. Once this is merged, I will work the
> matrix-free assembly of sparse matrices.
>
> Thanks,
> Peter
>
> On Monday, 9 March 2020 15:26:10 UTC+1, Michał Wichrowski wrote:
>>
>> Dear all,
>> I've got matrix-free multigrid solver for Stokes problem. The main
>> bottleneck is solution of coarse problem, so I tried to assemble the
>> regular sparse matrix and use direct solver. Since the coarse problem is
>> (relatively) small, I used vmults by unit vector to obtain columns of the
>> matrix. This is my code:
>>
>>
>> setup_matrix(); // generates sparsity patters and reinits matrices
>>
>> LinearAlgebra::distributed::BlockVector dst;
>> LinearAlgebra::distributed::BlockVector src;
>>
>> src.reinit(2);
>> for (unsigned int b = 0; b < 2; ++b)
>> stokes_matrix.get_matrix_free()->initialize_dof_vector(
>> src.block(b), b);
>> src.collect_sizes();
>> src =0;
>>
>>
>> dst.reinit(2);
>> dst.block(0).reinit(owned_dofs_u, relevant_dofs_u, MPI_COMM_WORLD);
>> dst.block(1).reinit(owned_dofs_p, relevant_dofs_p, MPI_COMM_WORLD);
>> dst.collect_sizes();
>>
>> for(types::global_dof_index i =0; i< owned_dofs_u.size(); ++i){
>> src=0;
>> dst=0;
>> if(owned_dofs_u.is_element(i) )
>> src.block(0)(i)=1;
>>
>> src.compress(VectorOperation::insert);
>> stokes_matrix.vmult(dst, src);
>>
>> dst.update_ghost_values();
>>
>> for(IndexSet::ElementIterator index_iter =owned_dofs_u.begin();
>> index_iter != owned_dofs_u.end();
>> ++index_iter){
>> if(dst.block(0)(*index_iter)!=0 )
>> block_A->set(*index_iter,i, dst.block(0)(*index_iter) );
>> }
>>
>> }
>>
>> block_A->compress(VectorOperation::unknown);
>>
>> Without constrains the matrix matches the matrix-free operator, but with
>> constrains present it does not. What is the proper way to assemble the
>> matrix with vmult?
>>
>> Best,
>> Michał Wichrowski
>>
>>
>>
>>
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/3944b062-fc59-4894-bc89-c43e017c16e9%40googlegroups.com.
