MatSetValues_MPIAIJ() triggers malloc error at line 565 during libMesh termination

[Ricardotu]

MatSetValues_MPIAIJ() triggers malloc error at line 565 during libMesh termination

[Ricardotu]

Thank you. The issue has been resolved by inserting zeros. In some time steps, even though it's known that certain entries are zero, inserting zeros was necessary to maintain consistency in the distribution of non-zero throughout all time step.

[roystgnr]

You've got something like a Discontinuous Galerkin method, which creates a coupling between neighboring elements' DoFs with non-overlapping support, because of jump terms on the interface between the neighbors?

We normally try to detect this by looking for systems with all discontinuous variables (in which case you must be doing interface coupling, right? apologies to anyone doing a global solve for a bunch of local L2 problems who got an overallocated matrix to boot...) - but do you have any non-discontinuous variables in the same system? In that case we fear you might just be doing something like a continuous-velocity+discontinuous-pressure formulation and we don't turn on that neighbor coupling for fear of grossly overallocating.

If you don't have any non-discontinuous variables in the same system, could you set up a small test code for us to replicate? That would probably be a libMesh bug, but I couldn't guess exactly how without getting into gdb.

If you do have non-discontinuous variables (or discontinuous-in-special-ways, like HDiv or HCurl vector variables) that are foiling the automatic implicit coupling, or if you don't but you want to try to workaround whatever bug might be hitting you, you can manually turn on implicit coupling between neighbors by using --implicit-neighbor-dofs on the command line (the easiest thing to do to check that that's the issue), or by calling set_implicit_neighbor_dofs() on your DofMap (the most robust way to ensure it, if you want everything coupled between neighbors), or by adding a custom GhostingFunctor subclass object as a DofMap coupling_functor (the most complicated but flexible way to do it, if you want some variables coupled between neighbors but not every variable coupled to every neighbor variable).

[roystgnr]

Hmm... looking at your description more closely, I'm wondering if there's something foiling us at the PETSc level here instead? You've only got discontinuous variables, you're getting successful time steps so you must be getting the implicit coupling to start with, and unless you're changing the topology of your problem as time marches forward (contact? an immersed body mesh?) there shouldn't be any issues going forward to the next step ... except that IIRC PETSc has a trick sometimes where it assumes you've probably overallocated your Jacobian, so when you finish your first assembly it gets rid of any entries that are exactly zero, but then if it turns out you didn't overallocate and one of those entries needs to become non-zero later it'll scream at you just as if you'd failed to preallocate properly. This could happen in an nonlinear transient problem where the Jacobian only starts to become "interesting" after the evolving solution breaks some initial symmetry.

I only vaguely remember others running into this, though; I don't remember the exact details or the workaround.

[lindsayad]

By default PETSc keeps zero entries during assembly. This is the canonical way to avoid the issue you are experiencing.

something that can cause your issue is if you have something like an upwind advection scheme and you don’t insert a 0 for dependence on a downstream cell but in a later time step the wind switches direction and now you try to add that dependency into the matrix

[Ricardotu]

Based on some constraints, I am unable to provide a small test case (I won't say I don't know how to do it). Here, I'll provide more details about my code. The systems involves only two Discontinuous Galerkin (DG) variables with L2 Lagrange basis functions.The problem itself is tackled using the penalty DG method, and there is also a need to handle upwind terms. I intend to solve what I'll refer to as "system a" (the main solving system, but it exhibits strong local nonlinearity leading to not nonlinear convergence). Therefore, I want to solve another system "b" before each time step, using the solution of "system b" as the initial guess for the solution of "system a". This is what is known as nonlinear preconditioning. The problem arises during the solution of "system b" and a crucial point is that the assembly of "system b" divides the computational domain into two parts (strongly and weakly nonlinear) based on the situation of "system a". just like that:

void b_residual(const NumericVector<Number> &soln,
                  NumericVector<Number> &residual,
                  NonlinearImplicitSystem &sys)
{
   // ...
  for (const auto &elem : mesh.active_local_element_ptr_range())
  {
    if (this elem is strongly nonlinear)
    {
      // ...
    }
    else
    {
        // ...
    }
  }

}

but ,in every time step , the division changes with each time step. The distribution of non-zero entries in the Jacobian matrix is entirely different between the two regions. In the strongly nonlinear region, it is similar to "a" (including dofs from current and neighboring elements). However, in the weak region, there are nonzero entries only along the main diagonal. This may be a contributing factor to the issue I am facing.

[lindsayad]

Great to hear! Can we regard this issue as closed?