Fix the extrapolation for guard cells in the insulator region (#5499)

For the nodal guard cells in the insulator region, the extrapolation was
being done relative to cells only inside the domain, not using the value
on the boundary. This fixes this to do extrapolation from the value on
the boundary and the next inward cell.

Source/BoundaryConditions/PEC_Insulator.cpp

@@ -62,8 +62,9 @@ namespace
                                   amrex::GpuArray<FieldBoundaryType, 3> const fbndry_hi)
     {
         using namespace amrex::literals;
+        amrex::IntVect ijk_next = ijk_vec;
+        amrex::IntVect ijk_nextp1 = ijk_vec;
         amrex::IntVect ijk_mirror = ijk_vec;
-        amrex::IntVect ijk_mirrorp1 = ijk_vec;
         bool OnBoundary = false;
         bool GuardCell = false;
         bool isInsulatorBoundary = false;
@@ -114,12 +115,14 @@ namespace
                     } else if (ig > 0) {
                         GuardCell = true;
 
+                        // Location of the next cells inward
+                        ijk_next[idim] = ijk_vec[idim] - ig*iside;
+                        ijk_nextp1[idim] = ijk_next[idim] - ig*iside;
+
                         // Mirror location inside the domain by "ig" number of cells
                         ijk_mirror[idim] = ( (iside == -1)
                                         ? (dom_lo[idim] + ig - (1 - is_nodal[idim]))
-                                        : (dom_hi[idim] + 1 - ig));
-                        // Location twice as far in, for extrapolation
-                        ijk_mirrorp1[idim] = 2*ijk_mirror[idim] - ijk_vec[idim];
+                                        : (dom_hi[idim] - ig + 1));
 
                         // Check for components with even symmetry.
                         // True for E_like and tangential, and B_like and normal
@@ -156,12 +159,12 @@ namespace
                 // The value on the boundary is left unmodified
                 // The values in the guard cells are extrapolated
                 if (GuardCell) {
-                    field(ijk_vec, n) = 2._rt*field(ijk_mirror, n) - field(ijk_mirrorp1, n);
+                    field(ijk_vec, n) = 2._rt*field(ijk_next, n) - field(ijk_nextp1, n);
                 }
             } else if ((OnBoundary || GuardCell) && set_field) {
                 field(ijk_vec, n) = field_value;
             } else if (GuardCell) {
-                field(ijk_vec, n) = 2._rt*field(ijk_mirror, n) - field(ijk_mirrorp1, n);
+                field(ijk_vec, n) = 2._rt*field(ijk_next, n) - field(ijk_nextp1, n);
             }
         } else {
             if (OnBoundary && (E_like ^ is_normal_to_boundary)) {