search

mfem / mfem

Lightweight, general, scalable C++ library for finite element methods
http://mfem.org
BSD 3-Clause "New" or "Revised" License
1.73k stars 502 forks source link

Some quesitons about the ImplicitSolve() from ex9p #4491

Closed Wayne901 closed 2 months ago

Wayne901 commented 2 months ago

Hi there, The DG_Solver in ex9p, as I understood, solve the time-dependent PDE with backward Euler, (M-dt*K) du_dt = K*u + b. However, I can not see the form of M-dt*K in DG_Solver class's SetTimeStep(). We get a matrix A=-dt*K and a diagnoal matrix A_diag=M_diag-dt*K_diag, where A_diag was not used latter. Is the diagonal of A linked to A_diag through A->GetDiag(A_diag)? Why there's no form of M-dt*K? I believe mass lumping was not applied in this example so M shouldn't be a mere diagonal ...

thartland commented 2 months ago

The diagonal of A is linked to A_diag through A->GetDiag(A_diag). See the comment hypre.hpp, which shows that A_diag is not a deep copy of the local block diagonal data of A. So, A = Add(-dt, K, 0.0, K) initially sets A = -dt*K, but then by changing A_diag here you are also changing A. This can be verified by calling A->Print("filename"); before and after altering A_diag with different filenames and you will see that the matrix entries of A differ.

You are right that M is not a mere diagonal, it seems to comes from this MassIntegrator. It seems to me that a BilinearForm whose assembled matrix is M - dt*K is another approach. I am not sure why it wasn't used.

Wayne901 commented 2 months ago

Thank you @thartland for your explanation. It's very clear. I'll leave the issue opened so maybe someone could explain why only the diagonal of M is used.

cjvogl commented 2 months ago

My understanding is that HypreParMatrix::GetDiag(SparseMatrix &) sets a reference to the local diagonal block of the HypreParMatrix (as opposed to strictly the diagonal entries). Because a DG finite element space is used, the mass matrix M is strictly block diagonal, so M_diag contains the entire M matrix (or at least the part that is local to the given MPI rank).

Wayne901 commented 2 months ago

So when diag is a SparseMatrix, HypreParMatrix::GetDiag(diag) returns a block diagonal. Thank you @cjvogl.