Account for metric term in fluxes due to inhomogeneous Neumann boundary conditions (#4579)

This fix targets the calculation of fluxes following a linear solve
using an MLCellABecLapT linear operator with a metric term. Previously,
there was an inconsistency between fluxes within the domain and at the
domain boundaries - the presence of the metric term was only accounted
for in the calculation of fluxes within the domain but not in the
calculation of fluxes at the domain boundary. By changing the order of
operations in the MLCellABecLapT::getFluxes method, this fix solves that
issue and fluxes are calculated consistently across the domain.

Author: jjschoser

Src/LinearSolvers/MLMG/AMReX_MLCellABecLap.H

@@ -285,12 +285,14 @@ MLCellABecLapT<MF>::getFluxes (const Vector<Array<MF*,AMREX_SPACEDIM> >& a_flux,
     for (int alev = 0; alev < nlevs; ++alev) {
         this->compFlux(alev, a_flux[alev], *a_sol[alev], a_loc);
         for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
-            this->unapplyMetricTerm(alev, 0, *a_flux[alev][idim]);
             if (betainv != RT(1.0)) {
                 a_flux[alev][idim]->mult(betainv, 0, ncomp);
             }
         }
         this->addInhomogNeumannFlux(alev, a_flux[alev], *a_sol[alev], true);
+        for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
+            this->unapplyMetricTerm(alev, 0, *a_flux[alev][idim]);
+        }
     }
 }