Update convergence criterion for MLMG (#6408)

In #2410, the convergence
criterion of the MLMG solver was modified so as to use the norm of `b`
(whenever `|b|` is not zero), to evaluate the relative error in the
Poisson equation (written as `Ax = b`) and decide when to stop the
iterations.

However, there are cases where this it not suitable. A typical case is
when the user fixes a large value of the potential at the boundaries (or
on the embedded boundary), and also has a very low-charge beam in the
simulation. In this case, the residual is dominated by the potential
created by the boundaries, but the norm of `b` (which depends on the
charge density of the beam) is very low, resulting in a large
`resid/bnorm`. In some cases, this large `resid/bnorm` can even trigger
the error `MLMG failing so lets stop here` because of this check in
AMREX:

https://github.com/AMReX-Codes/amrex/blob/development/Src/LinearSolvers/MLMG/AMReX_MLMG.H#L539

Instead, with this PR, we always use the maximum of the norm of `b` and
of the initial residual `resid0` to evaluate the relative error in the
Poisson equation - instead of using directly the maximum of the norm of
`b`.

I also checked that this fixes the issue mentioned here:
#6161 (comment)

The changes in the checksum is because of the change in convergence
criterion.

Author: RemiLehe

Regression/Checksum/benchmarks_json/test_2d_background_mcc.json

@@ -1,22 +1,22 @@
 {
   "lev=0": {
-    "rho_electrons": 0.03558889420026204,
-    "rho_he_ions": 0.041762340951146074
+    "rho_electrons": 0.03558889420351498,
+    "rho_he_ions": 0.04176234095095327
   },
   "he_ions": {
-    "particle_momentum_x": 2.8840921556995102e-18,
-    "particle_momentum_y": 2.1957048736351022e-18,
-    "particle_momentum_z": 2.198216674541737e-18,
-    "particle_position_x": 17607.425457505866,
-    "particle_position_y": 1100.0247860590894,
+    "particle_momentum_x": 2.8840921570976e-18,
+    "particle_momentum_y": 2.195704873634808e-18,
+    "particle_momentum_z": 2.1982166745272215e-18,
+    "particle_position_x": 17607.425457559715,
+    "particle_position_y": 1100.024786059876,
     "particle_weight": 71976747650.1465
   },
   "electrons": {
-    "particle_momentum_x": 1.0105526130672624e-18,
-    "particle_momentum_y": 2.8197429889743996e-19,
-    "particle_momentum_z": 2.8091808266114265e-19,
-    "particle_position_x": 17136.01865030665,
-    "particle_position_y": 936.3651769884851,
+    "particle_momentum_x": 1.010552606677337e-18,
+    "particle_momentum_y": 2.81974298750407e-19,
+    "particle_momentum_z": 2.8091808348750305e-19,
+    "particle_position_x": 17136.0186559433,
+    "particle_position_y": 936.3651746710019,
     "particle_weight": 61113170379.63868
   }
-}
+}

Regression/Checksum/benchmarks_json/test_2d_embedded_boundary_cube.json

@@ -3,8 +3,8 @@
     "Bx": 9.2636945391314e-05,
     "By": 0.00031905198911870116,
     "Bz": 7.328424778796838e-05,
-    "Ex": 8553.90669231672,
+    "Ex": 8553.906692257373,
     "Ey": 60867.0482641368,
-    "Ez": 4.223901049143958e-06
+    "Ez": 3.35209259952839e-08
   }
 }

Regression/Checksum/benchmarks_json/test_2d_particle_reflection_picmi.json

@@ -1,7 +1,7 @@
 {
   "lev=0": {
-    "Ex": 4.865922376234882e-11,
+    "Ex": 1.5548810315827429e-18,
     "Ey": 0.0,
-    "Ez": 2.3293326580399806e-10
+    "Ez": 9.69073951936753e-18
   }
-}
+}

Source/ablastr/fields/EffectivePotentialPoissonSolver.H

@@ -150,10 +150,6 @@ computeEffectivePotentialPhi (
 
     auto const finest_level = static_cast<int>(rho.size() - 1);
 
-    // determine if rho is zero everywhere
-    const amrex::Real max_norm_b = getMaxNormRho(
-        amrex::GetVecOfConstPtrs(rho), finest_level, absolute_tolerance);
-
     const amrex::LPInfo info;
 
     for (int lev=0; lev<=finest_level; lev++) {
@@ -231,7 +227,7 @@ computeEffectivePotentialPhi (
         amrex::MLMG mlmg(*linop); // actual solver defined here
         mlmg.setVerbose(verbosity);
         mlmg.setMaxIter(max_iters);
-        mlmg.setConvergenceNormType((max_norm_b > 0) ? amrex::MLMGNormType::bnorm : amrex::MLMGNormType::greater);
+        mlmg.setConvergenceNormType(amrex::MLMGNormType::greater);
 
         const int ng = int(grid_type == utils::enums::GridType::Collocated); // ghost cells
         if (ng) {

Source/ablastr/fields/PoissonSolver.H

@@ -261,10 +261,6 @@ computePhi (
             [](const auto b){return (b<1.0_rt);}),
         "Components of beta_solver must be < 1.");
 
-    // determine if rho is zero everywhere
-    const amrex::Real max_norm_b = getMaxNormRho(
-        amrex::GetVecOfConstPtrs(rho), finest_level, absolute_tolerance);
-
     amrex::LPInfo info;
 
     for (int lev=0; lev<=finest_level; lev++) {
@@ -274,6 +270,9 @@ computePhi (
                               geom[lev].CellSize(2)/std::sqrt(1._rt-beta_solver[2]*beta_solver[2]))};
 
 #if (defined(ABLASTR_USE_FFT) && defined(WARPX_DIM_3D))
+        // determine if rho is zero everywhere
+        const amrex::Real max_norm_b = getMaxNormRho(
+            amrex::GetVecOfConstPtrs(rho), finest_level, absolute_tolerance);
         // Use the Integrated Green Function solver (FFT) on the coarsest level if it was selected
         if(is_solver_igf_on_lev0 && lev==0){
             if ( max_norm_b == 0 ) {
@@ -404,7 +403,7 @@ computePhi (
         amrex::MLMG mlmg(*linop); // actual solver defined here
         mlmg.setVerbose(verbosity);
         mlmg.setMaxIter(max_iters);
-        mlmg.setConvergenceNormType((max_norm_b > 0) ? amrex::MLMGNormType::bnorm : amrex::MLMGNormType::greater);
+        mlmg.setConvergenceNormType(amrex::MLMGNormType::greater);
 
         const int ng = int(grid_type == utils::enums::GridType::Collocated); // ghost cells
         if (ng) {