Refactor AmrMesh::MakeNewGrids

The refactor addresses two issues.

First, there was a bug when projecting grids during regrid. TagBoxArray
might not have enough ghost cells to guarantee proper nesting when the new
fine level grows by a lot.

Second, it could be very inefficient to use the previous approach of adding
extra ghost cells to TagBoxArray for proper nesting. In some users' cases,
it may need more than 30 ghost cells.

In the new approach, instead of using extra ghost cells, we explicitly add
the projected tags outside the TagBoxArray's domain.

Author: WeiqunZhang

Src/AmrCore/AMReX_AmrMesh.cpp

@@ -645,28 +645,8 @@ AmrMesh::MakeNewGrids (int lbase, Real time, int& new_finest, Vector<BoxArray>&
     for (int levc = max_crse; levc >= lbase; levc--)
     {
         int levf = levc+1;
-        //
-        // Construct TagBoxArray with sufficient grow factor to contain
-        // new levels projected down to this level.
-        //
-        IntVect ngt = n_error_buf[levc];
-        BoxArray ba_proj;
-        if (levf < new_finest)
-        {
-            ba_proj = new_grids[levf+1].simplified();
-            ba_proj.coarsen(ref_ratio[levf]);
-            ba_proj.growcoarsen(n_proper, ref_ratio[levc]);
 
-            BoxArray levcBA = grids[levc].simplified();
-            int ngrow = 0;
-            while (!levcBA.contains(ba_proj))
-            {
-                levcBA.grow(1);
-                ++ngrow;
-            }
-            ngt.max(IntVect(ngrow));
-        }
-        TagBoxArray tags(grids[levc],dmap[levc],ngt);
+        TagBoxArray tags(grids[levc],dmap[levc],n_error_buf[levc]);
 
         //
         // Only use error estimation to tag cells for the creation of new grids
@@ -713,12 +693,92 @@ AmrMesh::MakeNewGrids (int lbase, Real time, int& new_finest, Vector<BoxArray>&
         ManualTagsPlacement(levc, tags, bf_lev);
         //
         // If new grids have been constructed above this level, project
-        // those grids down and tag cells on intersections to ensure
-        // proper nesting.
+        // those grids down and tag cells on intersections to ensure proper
+        // nesting. Note that the projected BoxArray may contain cells
+        // outside the TagBoxArray's domain. For those, we collect them in
+        // Vector tag_proj.
         //
+        Vector<IntVect> tags_proj;
         if (levf < new_finest) {
-            ba_proj.coarsen(bf_lev[levc]);
+            BoxArray ba_tags = tags.boxArray();
+            BoxList bl_proj = new_grids[levf+1].simplified_list();
+            auto& bxs = bl_proj.data();
+            Long nbxs = bl_proj.size();
+            Vector<Vector<IntVect>> tags_proj_priv(OpenMP::get_max_threads());
+            Box domain = Geom(levc).Domain();
+            domain.coarsen(bf_lev[levc]);
+            auto const& domain_lo = domain.smallEnd();
+            auto const& domain_hi = domain.bigEnd();
+            auto const& domain_len = domain.length();
+            auto const& is_periodic = Geom(levc).isPeriodicArray();
+#ifdef AMREX_USE_OMP
+#pragma omp parallel
+#endif
+            {
+                auto& tv = tags_proj_priv[OpenMP::get_thread_num()];
+                std::vector<std::pair<int,Box>> isects;
+                BoxList bl(ba_tags.ixType());
+                BoxList bl2(ba_tags.ixType());
+                BoxList bltmp;
+#ifdef AMREX_USE_OMP
+#pragma omp for
+#endif
+                for (Long ibx = 0; ibx < nbxs; ++ibx) {
+                    Box& b = bxs[ibx];
+                    b.coarsen(ref_ratio[levf]);
+                    b.grow(bf_lev[levf]*n_proper);
+                    b.coarsen(ref_ratio[levc]);
+                    b.coarsen(bf_lev[levc]);
+
+                    ba_tags.intersections(b, isects, false, tags.nGrowVect());
+                    bl.clear();
+                    bl.push_back(b);
+                    for (auto const& kv : isects) {
+                        bl2.clear();
+                        for (auto & btmp : bl) {
+                            amrex::boxDiff(bltmp, btmp, kv.second);
+                            bl2.join(bltmp);
+                        }
+                        std::swap(bl,bl2);
+                    }
+                    for (auto const& bleft : bl) {
+                        amrex::LoopOnCpu(bleft, [&] (IntVect iv)
+                        {
+                            for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
+                                if (iv[idim] < domain_lo[idim]) {
+                                    if (is_periodic[idim]) {
+                                        iv[idim] += domain_len[idim];
+                                    } else {
+                                        return;
+                                    }
+                                } else if (iv[idim] > domain_hi[idim]) {
+                                    if (is_periodic[idim]) {
+                                        iv[idim] -= domain_len[idim];
+                                    } else {
+                                        return;
+                                    }
+                                }
+                            }
+                            tv.push_back(iv);
+                        });
+                    }
+                }
+            }
+
+            BoxArray ba_proj(std::move(bl_proj));
             tags.setVal(ba_proj,TagBox::SET);
+
+            Long nextra = 0;
+            for (auto const& tv : tags_proj_priv) {
+                nextra += tv.size();
+            }
+            if (nextra > 0) {
+                tags_proj.reserve(nextra);
+                for (auto const& tv : tags_proj_priv) {
+                    tags_proj.insert(std::end(tags_proj), std::begin(tv), std::end(tv));
+                }
+                amrex::RemoveDuplicates(tags_proj);
+            }
         }
         //
         // Map tagged points through periodic boundaries, if any.
@@ -739,6 +799,10 @@ AmrMesh::MakeNewGrids (int lbase, Real time, int& new_finest, Vector<BoxArray>&
         tags.collate(tagvec);
         tags.clear();
 
+        if (!tags_proj.empty()) {
+            tagvec.insert(tagvec.end(), tags_proj.data(), tags_proj.data()+tags_proj.size());
+        }
+
         if (!tagvec.empty())
         {
             //