Source term cannot be applied when non_equilibrium_initial_residuum is activated.
If the function of
compensate_non_equilibrium_initial_residuum is enabled, an initial non equilibrium initial residuum vector is computed by taking account of the contribution of the nodal source terms, which are computed from Neumann BCs or other source terms. However the contribution of the nodal source terms to the residuum vector vector disappears after the global assembly of the linearization step with the following computation
if (_r_neq != nullptr) LinAlg::axpy(res, -1, *_r_neq);
Therefore the prescribed Neumann BCs or source terms take no effect with
I will start to fix this issue by dropping the assembly of the source term in the calculation of non equilibrium initial residuum.