Inverting a matrix

[colmenaresj]

Hello,

For an explicit solver, would it make sense to invert the mass matrix once, and then just multiply the RHS by it, instead of having to do a "Solve" each time? This is assuming the mass matrix doesn't change over time. Thanks.

 // RHS for Explicit Solver
   void Mult(const Vector &u, Vector &du_dt) const override
   {
      K->Mult(u, t1);
      t1.Neg();
      t1 += (*b);
      M_solver.Mult(t1, du_dt); // Change for something like Minv.Mult( )...
      du_dt.SetSubVector(ess_tdofs_, 0.0);
   }

[colmenaresj]

If it makes sense to do that, how would I do it in MFEM?

[artv3]

I think it would depending on the space. If it is L^2, then certainly. One way is to use the MassIntegrator and call the element assembly (https://github.com/mfem/mfem/blob/master/fem/integ/bilininteg_mass_ea.cpp#L259) method to store the element matrices in an mfem::Vector, then use methods from the MFEM linear algebra files to invert the matrices.

If this is H1, then I think it gets more complicated since you may have an MPI distributed data structure.

[colmenaresj]

Thank you @artv3 . I am indeed working with H1 and was hoping to extend it later to DG. It sounds like it might be cumbersome.

[colmenaresj]

The reason why I am interested in this is because when you compare the speed of MFEM vs a finite difference code with explicit solver with similar resolution, the FD code is much faster. My thought that this was due to having to invert the mass matrix every timestep. Is that a fair assumption? Or what else could help get some speedup?
I am running on CPU only for now, with 8 MPI ranks.

[artv3]

In general FD methods will have substantially less flops, and for cartesian geometries may be hard to beat. Finite elements are more flexible and will be able to handle more complex geometries easier than FD methods. One suggestion I have is to make sure you are using the partial assembly options in the ParBilinearForms,

a.SetAssemblyLevel(AssemblyLevel::PARTIAL);

to avoid costly matrix assembly. You can also pair your solver with a matrix-free pre-conditioner such as:

Solver *prec = new OperatorJacobiSmoother(a, ess_tdof_list);

Some of these concepts are captures in example ex1p.cpp for the diffusion integrator but is same concepts can be applied to the mass integrator.

[colmenaresj]

excellent, thanks again for response