Using the DG limiter on composition in an iterative Advection scheme halts nonlinear convergence

[anne-glerum]

As can be seen in the screen output of the two tests in #5531, switching on the DG limiter when using an iterative Advection scheme halts the nonlinear solver convergence. The relative residual of the advection does not improve after the ~second nonlinear iteration and therefore stalls the total residual. Note that this only occurs when the global limits actually limit values of the field; if the limits are set wider/higher than the actual values of the fields, convergence is the same as when no limit is set. Also, the limiting works as intended: it keeps values between the specified global min and max, avoiding under/overshoots.

I've seen this behaviour in several setups while working on #4370. I have not tested the limiter with DG temperature. The current_residual that is returned by solve_advection() reduces much less with nonlinear iterations with the limiter on than with the limiter off. The initial residual at time zero (t0) is the same with and without the limiter, as expected.

The limiter is applied to the solution after each advection solve. I wonder whether at t0 it actually needs to be applied during the first nonlinear iteration. The reaction terms and rates are zero at t0 and in the first nonlinear iteration, velocity is also zero. (Switching it off doesn't fix the issue though).

[gassmoeller]

I looked a bit into this, but dont have a satisfactory solution yet. I think I have narrowed down the issue though.
Here is the current sequence of operations:

0. Before solving the advection equations we compute the current residual using the current_linearization_point (which for nonlinear_iterations > 0 already contains the limited solution)
1. We solve the advection equation, this obviously satisfies the equation and significantly reduces the nonlinear residual (but we dont compute it on this solution).
2. Then we apply the limiter. This modifies the solution (removing overshoots), but now the solution no longer satisfies the equation and the residual stalls. This solution is what we copy into the current_linearization_point for the next iteration.

So the issue seems to be that we do not compute the residual after solving a nonlinear iteration with the actual solution (before limiting), but before the next one (with the limited solution). I do not remember the rationale for this in detail, but I think it was related to the fact that this way we automatically compute the initial nonlinear residual in the first nonlinear iteration. Maybe it is time to rethink this approach? This could also allow us to reduce the number of nonlinear iterations by one (currently we always do at least 2 nonlinear iterations if we have to solve at least once, so one iteration too many.

[anne-glerum]

Thanks for looking in to this Rene. In trying to understand what's going on again, I've listed pseudo-code that summarizes when what residual is being computed:

1. solve_iterated_advection_and_stokes()
   initial_composition_residual = 0
   relative_residual = max
   relative_composition_residual = assemble_and_solve_composition(initial_composition_residual, NI == 0 ?
   &initial_composition_residual : nullptr);
   relative_residual = max(relative_temperature_residual, relative_composition_residual)
   relative_residual = max(relative_residual, relative_nonlinear_stokes_residual)
   check relative_residual whether to stop iterating
2. assemble_and_solve_composition()
   current_residual = 0
   if (initial_composition_residual)
   (*initial_composition_residual)[c] = system_rhs.block(introspection.block_indices.compositional_fields[c]).l2_norm();
   current_residual = solve_advection()
   current_residual /= initial_composition_residual
   return current_residual
3. solve_advection()
   distributed_solution = current_linearization_point
   const double initial_residual = system_matrix.residual(temp,
   distributed_solution,
   system_rhs);
   solve()
   apply_limiter_to_dg_solutions()
   return initial_residual

So if I understand this correctly, the initial_composition_residual is calculated as the rhs l2_norm in the first nonlinear iteration (NI==0). This is separate from the current_residual for composition, which is computed in the call to solve_advection() before solving and before applying the limiter. In point 3, could we move the computation of initial_residual (i.e. current_residual) to after the solve(), but before the apply_limiter_to_dg_solutions() call? The initial_composition_residual would not be affected, but the current_residual would be computed based on the solution of the current iteration, not the limited solution of the previous iteration.

[anne-glerum]

Talking to Juliane yesterday, we do the calculation of the residual before solving so that we can use the already assembled matrix and rhs based on the previous solution. We could replace the current_linearization_point in the computation of the current residual before the solve with the unlimited solution. This unlimited solution would have to be available then though; at the moment it is not stored.