Add a few additional diags (divE etc.) (#844)

* Start implementation of new averaging with staggering:

- face-to-cell-center and edge-to-cell-center replaced so far;
- TODO: node-to-cell-center and 1D behavior (AMREX_SPACEDIM=1).

* first implementation of Diags base classes

* add example, temporarily

* Continue implementation of new averaging with staggering:

- new function takes reference to single MultiFab (no vector);
- TODO: node-to-cell-center still in progress.

* Fix small bug and clean up

* Fix bug in loop over n=0,...,ncomp-1 and clean up

* add more functions

* Add Doxygen documentation and clean up

* Small clean-up in Doxygen documentation

* Compile in single precision: add _rt suffix to avoid unnecessary conversions

* Avoid accessing staggering index directly from IntVect in innermost loops

* Replace do-while loop with for loop (default ncomp=1)

* Remove temporary pointer and pass reference to MultiFab (instead of MultiFab*)

* Replace AMREX_LAUNCH_HOST_DEVICE_LAMBDA with ParallelFor

* cleaning and initialize output mf

* use general average routine

* move flush in new class, and implemented the Plotfile derived class

* add comments

* eol

* free memory in destructor

* typo

* typo

* no need to clear MF pointers there

* though shalt not break existing tests

* FlushRaw doesnt have to be virtual for now

* The importance of being constant

* Capability to select fields in output files

* EOL

* revert to old inputs

* const in right place

* avoid brace initializer there

* oops, fix logic error in is_in

* user can choose flush interval, same behavior as plot_int

* Add option to plot raw fields

* eol

* replace ter flush with dump to avoid confusion

* add options

* Diagnostics stores a vector of functors to compute diags on the fly

* eol

* Field gather from diags to sync particle quantities

* New diagnostics handle RZ with same behavior as old ones

* cleaning and doc

* const ref for string

* smarter for loop from Axel and typo fix from Reva

* Functors to compute some fields

* simplify code following Dave's comments

* Create subfolders and add more output options (divE etc.)

* eol

* rename mode_avg to convertRZmodes2cartesian

* Update CellCenterFunctor.H

* fill varnames and vector of functors at the same time

* output rho_new, not rho_old

* WarpX instance not needed here

* add const

* little bit more of reorganization

* Apply suggestions from code review

Co-Authored-By: Axel Huebl <axel.huebl@plasma.ninja>

* add a bunch of const

* make derived classes final

Co-authored-by: Edoardo Zoni <ezoni@lbl.gov>
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>

1 files changed, 148 insertions, 0 deletions

diff --git a/Source/Utils/Interpolate.cpp b/Source/Utils/Interpolate.cpp
new file mode 100644
index 000000000..2dbfddb96
--- /dev/null
+++ b/Source/Utils/Interpolate.cpp
@@ -0,0 +1,148 @@
+#include "Interpolate.H"
+#include <AMReX_FillPatchUtil_F.H>
+
+namespace Interpolate
+{
+    using namespace amrex;
+
+    std::unique_ptr<MultiFab>
+    getInterpolatedScalar(
+        const MultiFab& F_cp, const MultiFab& F_fp,
+        const DistributionMapping& dm, const int r_ratio,
+        const Real* /*dx*/, const int ngrow )
+    {
+        // Prepare the structure that will contain the returned fields
+        std::unique_ptr<MultiFab> interpolated_F;
+        interpolated_F.reset( new MultiFab(F_fp.boxArray(), dm, 1, ngrow) );
+        interpolated_F->setVal(0.);
+
+        // Loop through the boxes and interpolate the values from the _cp data
+#ifdef _OPEMP
+#pragma omp parallel
+#endif
+        {
+            FArrayBox ffab; // Temporary array ; contains interpolated fields
+            for (MFIter mfi(*interpolated_F); mfi.isValid(); ++mfi)
+            {
+                Box finebx = mfi.fabbox();
+                finebx.coarsen(r_ratio).refine(r_ratio); // so that finebx is coarsenable
+
+                const FArrayBox& cfab = (F_cp)[mfi];
+                ffab.resize(finebx);
+
+                // - Fully nodal
+                if ( F_fp.is_nodal() ){
+                    IntVect refinement_vector{AMREX_D_DECL(r_ratio, r_ratio, r_ratio)};
+                    node_bilinear_interp.interp(cfab, 0, ffab, 0, 1,
+                                                finebx, refinement_vector, {}, {}, {}, 0, 0, RunOn::Cpu);
+                } else {
+                    amrex::Abort("Unknown field staggering.");
+                }
+
+                // Add temporary array to the returned structure
+                const Box& bx = (*interpolated_F)[mfi].box();
+                (*interpolated_F)[mfi].plus<RunOn::Host>(ffab, bx, bx, 0, 0, 1);
+            }
+        }
+        return interpolated_F;
+    }
+
+    std::array<std::unique_ptr<MultiFab>, 3>
+    getInterpolatedVector(
+        const MultiFab* Fx_cp,
+        const MultiFab* Fy_cp,
+        const MultiFab* Fz_cp,
+        const MultiFab* Fx_fp,
+        const MultiFab* Fy_fp,
+        const MultiFab* Fz_fp,
+        const DistributionMapping& dm, const int r_ratio,
+        const Real* dx, const int ngrow )
+    {
+
+        // Prepare the structure that will contain the returned fields
+        std::array<std::unique_ptr<MultiFab>, 3> interpolated_F;
+        interpolated_F[0].reset( new MultiFab(Fx_fp->boxArray(), dm, 1, ngrow) );
+        interpolated_F[1].reset( new MultiFab(Fy_fp->boxArray(), dm, 1, ngrow) );
+        interpolated_F[2].reset( new MultiFab(Fz_fp->boxArray(), dm, 1, ngrow) );
+        for (int i=0; i<3; i++) interpolated_F[i]->setVal(0.);
+
+        // Loop through the boxes and interpolate the values from the _cp data
+        const int use_limiter = 0;
+#ifdef _OPEMP
+#pragma omp parallel
+#endif
+        {
+            std::array<FArrayBox,3> ffab; // Temporary array ; contains interpolated fields
+            for (MFIter mfi(*interpolated_F[0]); mfi.isValid(); ++mfi)
+            {
+                Box ccbx = mfi.fabbox();
+                ccbx.enclosedCells();
+                ccbx.coarsen(r_ratio).refine(r_ratio); // so that ccbx is coarsenable
+
+                const FArrayBox& cxfab = (*Fx_cp)[mfi];
+                const FArrayBox& cyfab = (*Fy_cp)[mfi];
+                const FArrayBox& czfab = (*Fz_cp)[mfi];
+                ffab[0].resize(amrex::convert(ccbx,(*Fx_fp)[mfi].box().type()));
+                ffab[1].resize(amrex::convert(ccbx,(*Fy_fp)[mfi].box().type()));
+                ffab[2].resize(amrex::convert(ccbx,(*Fz_fp)[mfi].box().type()));
+
+                // - Face centered, in the same way as B on a Yee grid
+                if ( (*Fx_fp)[mfi].box().type() == IntVect{AMREX_D_DECL(1,0,0)} ){
+#if (AMREX_SPACEDIM == 3)
+                    amrex_interp_div_free_bfield(ccbx.loVect(), ccbx.hiVect(),
+                                                 BL_TO_FORTRAN_ANYD(ffab[0]),
+                                                 BL_TO_FORTRAN_ANYD(ffab[1]),
+                                                 BL_TO_FORTRAN_ANYD(ffab[2]),
+                                                 BL_TO_FORTRAN_ANYD(cxfab),
+                                                 BL_TO_FORTRAN_ANYD(cyfab),
+                                                 BL_TO_FORTRAN_ANYD(czfab),
+                                                 dx, &r_ratio, &use_limiter);
+#else
+                    amrex_interp_div_free_bfield(ccbx.loVect(), ccbx.hiVect(),
+                                                 BL_TO_FORTRAN_ANYD(ffab[0]),
+                                                 BL_TO_FORTRAN_ANYD(ffab[2]),
+                                                 BL_TO_FORTRAN_ANYD(cxfab),
+                                                 BL_TO_FORTRAN_ANYD(czfab),
+                                                 dx, &r_ratio, &use_limiter);
+                    amrex_interp_cc_bfield(ccbx.loVect(), ccbx.hiVect(),
+                                           BL_TO_FORTRAN_ANYD(ffab[1]),
+                                           BL_TO_FORTRAN_ANYD(cyfab),
+                                           &r_ratio, &use_limiter);
+#endif
+                    // - Edge centered, in the same way as E on a Yee grid
+                } else if ( (*Fx_fp)[mfi].box().type() == IntVect{AMREX_D_DECL(0,1,1)} ){
+#if (AMREX_SPACEDIM == 3)
+                    amrex_interp_efield(ccbx.loVect(), ccbx.hiVect(),
+                                        BL_TO_FORTRAN_ANYD(ffab[0]),
+                                        BL_TO_FORTRAN_ANYD(ffab[1]),
+                                        BL_TO_FORTRAN_ANYD(ffab[2]),
+                                        BL_TO_FORTRAN_ANYD(cxfab),
+                                        BL_TO_FORTRAN_ANYD(cyfab),
+                                        BL_TO_FORTRAN_ANYD(czfab),
+                                        &r_ratio, &use_limiter);
+#else
+                    amrex_interp_efield(ccbx.loVect(), ccbx.hiVect(),
+                                        BL_TO_FORTRAN_ANYD(ffab[0]),
+                                        BL_TO_FORTRAN_ANYD(ffab[2]),
+                                        BL_TO_FORTRAN_ANYD(cxfab),
+                                        BL_TO_FORTRAN_ANYD(czfab),
+                                        &r_ratio,&use_limiter);
+                    amrex_interp_nd_efield(ccbx.loVect(), ccbx.hiVect(),
+                                           BL_TO_FORTRAN_ANYD(ffab[1]),
+                                           BL_TO_FORTRAN_ANYD(cyfab),
+                                           &r_ratio);
+#endif
+                } else {
+                    amrex::Abort("Unknown field staggering.");
+                }
+
+                // Add temporary array to the returned structure
+                for (int i = 0; i < 3; ++i) {
+                    const Box& bx = (*interpolated_F[i])[mfi].box();
+                    (*interpolated_F[i])[mfi].plus<RunOn::Host>(ffab[i], bx, bx, 0, 0, 1);
+                }
+            }
+        }
+        return interpolated_F;
+    }
+}