Solution mismatch (FullMatrix vs. PETScWrappers::MPI::SparseMatrix)

[Hermes Sampedro]

Dear all, 

I am experiencing strange behaviour in what I am trying to do and I would like to ask for a bit of support.

I have a running solver done with a distributed MPI implementation. The solver is similar to step-29, where real and imaginary parts are used and a 2N system is solved.

Now I want to apply reduced basis methods to it and to do that I will solve a FullMatrix system Lx=b by x=L^1*b.

For now, and to have a sanity check,  what I am doing is to convert the operators (PETScWrappers::MPI::SparseMatrix) to FullMatrix after assemble_system(),  and solve the system by inverting the operator and using vmult(). However, the results are different when I solve with sparse matrix compared to FullMatrix.

The solver code is bellow, where locally_relevant_solution is the solution for the distributed system and solution_rb is thee

for the FullMatrix system.

__________________________________________________________________________

void Problem<dim>::solve()

{

//Distributed solver

PETScWrappers::MPI::Vector completely_distributed_solution(locally_owned_dofs,mpi_communicator);

SolverControl cn;

PETScWrappers::SparseDirectMUMPS solver(cn, mpi_communicator);

solver.solve(system_matrix, completely_distributed_solution, system_rhs);

constraints.distribute(completely_distributed_solution);

locally_relevant_solution = completely_distributed_solution;

// Full matrix solver

FullMatrix<double> Linv(dof_handler.n_dofs(),dof_handler.n_dofs());

Vector<double> solution_rb(dof_handler.n_dofs());

Linv.invert(L);

Linv.vmult(solution_rb,b,false);

}

__________________________________________________________________________

I checked that in both cases, the operators/system_matrix have the same values just in the assemble_system() function. I was wondering if there is something internally that change somehow the systme_matrix. I can see that after checking the values of the operators I have the following lines which I am not sure if there is something happening:

_________________

...

constraints.distribute_local_to_global(cell_matrix,

cell_rhs,

local_dof_indices,

system_matrix,

system_rhs);

}

system_matrix.compress(VectorOperation::add);

system_rhs.compress(VectorOperation::add);

_________________

An example of how the results mismatch is given in the following figure where the dotted line corresponds to the solution with FullMatrix and the line corresponds to the distributed solution.

I have to add that the distributed solver that I have is validated and thus, I trust the solutions. It looks like there is a multiplication factor somewhere, that I am missing where using the FullMatrix.

Thank you very much for your help. Regards,

[Wolfgang Bangerth]

On 2/8/22 09:06, Hermes Sampedro wrote:
>
> I checked that in both cases, the operators/system_matrix have the same
> values just in the assemble_system() function. I was wondering if there
> is something internally that change somehow the systme_matrix. I can see
> that after checking the values of the operators I have the following
> lines which I am not sure if there is something happening:

No functions in deal.II touch your linear system unless you specifically
call a function that takes a matrix or vector as argument. So, if you
have the distribute_local_to_global() call for one type of matrix, you
need to call it for the other kind of matrix as well. But these
functions do not internally store some kind of magic pointer to the
matrix and right hand side so that at some later point, when you're not
watching, they do something nefarious to your matrix :-)

My suggestion would be to make sure that before the solve, the two
matrices and corresponding right hand sides are the same. This is often
awkward if the objects are large, but you can test this by outputting
the norm of these objects.

Then check that the solution is also the same. My guess is that for some
reason, either the matrix or the rhs are different.

Best
W.

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

[Hermes Sampedro]

Dear Wolfgang, 

Thank you for the clarification. In that case, I am missing some operation somewhere...

Let me clarify what I am doing. Apart from building the system_matrix, I save the operator L in the assembling process. But I did not do any distribute_local_to_global nor compress to that operator, as I do for the system matrix. Is that wrong?

________________________________________________________________________________________

for (const auto &cell : dof_handler.active_cell_iterators())

if (cell->is_locally_owned())

{

for (unsigned int i = 0; i < dofs_per_cell; ++i)

{ ......

for (unsigned int j = 0; j < dofs_per_cell; ++j)

{. .....

cell_matrix(i, j)= ....

....

L_operator(local_dof_indices[i],local_dof_indices[j]) = cell_matrix(i,j);

}

}

cell->get_dof_indices(local_dof_indices);

constraints.distribute_local_to_global(cell_matrix,

cell_rhs,

local_dof_indices,

system_matrix,

system_rhs);

}

system_matrix.compress(VectorOperation::add);

system_rhs.compress(VectorOperation::add);

________________________________________________________________________________________

I checked the content of the matrics by solving a very small case. I found that indeed system_matrix is different to the saved operator L.  Do you have any suggestion where the different is coming from? Both are taken from the same assembling, the difference is that system_matrix applies distribute_local_to_global() and compress().

Thank you again.

Regards

[Timo Heister]

> L_operator(local_dof_indices[i],local_dof_indices[j]) = cell_matrix(i,j);

This looks incorrect, because global pairs of indices will happen more than once in a normal assembly (with a continuous finite element). You have to add the contribution using +=.

--
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+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/0b63fba4-0c3f-4727-9ce2-18ddc897aba9n%40googlegroups.com.

[Hermes Sampedro]

Thank you very much

Regards