aboutsummaryrefslogtreecommitdiff
path: root/Source
diff options
context:
space:
mode:
Diffstat (limited to 'Source')
-rw-r--r--Source/Diagnostics/FieldIO.H5
-rw-r--r--Source/Diagnostics/FieldIO.cpp31
-rw-r--r--Source/Diagnostics/WarpXIO.cpp23
-rw-r--r--Source/Filter/BilinearFilter.H26
-rw-r--r--Source/Filter/BilinearFilter.cpp156
-rw-r--r--Source/Filter/Make.package3
-rw-r--r--Source/Filter/filter_module.F9087
-rw-r--r--Source/FortranInterface/WarpX_f.H6
-rw-r--r--Source/Initialization/WarpXInitData.cpp15
-rw-r--r--Source/Parallelization/WarpXComm.cpp69
-rw-r--r--Source/WarpX.H10
-rw-r--r--Source/WarpX.cpp153
12 files changed, 372 insertions, 212 deletions
diff --git a/Source/Diagnostics/FieldIO.H b/Source/Diagnostics/FieldIO.H
index fc3aa9672..121fb9374 100644
--- a/Source/Diagnostics/FieldIO.H
+++ b/Source/Diagnostics/FieldIO.H
@@ -75,4 +75,9 @@ getInterpolatedVector(
const int r_ratio, const Real* dx,
const int ngrow );
+coarsenCellCenteredFields(
+ Vector<MultiFab>& coarse_mf, Vector<Geometry>& coarse_geom,
+ const Vector<MultiFab>& source_mf, const Vector<Geometry>& source_geom,
+ int coarse_ratio, int finest_level );
+
#endif
diff --git a/Source/Diagnostics/FieldIO.cpp b/Source/Diagnostics/FieldIO.cpp
index 03754404b..965324d02 100644
--- a/Source/Diagnostics/FieldIO.cpp
+++ b/Source/Diagnostics/FieldIO.cpp
@@ -298,7 +298,36 @@ WarpX::AverageAndPackFields ( Vector<std::string>& varnames,
}
BL_ASSERT(dcomp == ncomp);
- } // end loop over levels of refinement
+
+ } // end loop over levels of refinement
+
+};
+
+/** \brief Reduce the size of all the fields in `source_mf`
+ * by `coarse_ratio` and store the results in `coarse_mf`.
+ * Calculate the corresponding coarse Geometry from `source_geom`
+ * and store the results in `coarse_geom` */
+void
+coarsenCellCenteredFields(
+ Vector<MultiFab>& coarse_mf, Vector<Geometry>& coarse_geom,
+ const Vector<MultiFab>& source_mf, const Vector<Geometry>& source_geom,
+ int coarse_ratio, int finest_level )
+{
+ // Check that the Vectors to be filled have an initial size of 0
+ AMREX_ALWAYS_ASSERT( coarse_mf.size()==0 );
+ AMREX_ALWAYS_ASSERT( coarse_geom.size()==0 );
+
+ // Fill the vectors with one element per level
+ int ncomp = source_mf[0].nComp();
+ for (int lev=0; lev<=finest_level; lev++) {
+ AMREX_ALWAYS_ASSERT( source_mf[lev].is_cell_centered() );
+
+ coarse_geom.push_back(amrex::coarsen( source_geom[lev], IntVect(coarse_ratio)));
+
+ BoxArray small_ba = amrex::coarsen(source_mf[lev].boxArray(), coarse_ratio);
+ coarse_mf.push_back( MultiFab(small_ba, source_mf[lev].DistributionMap(), ncomp, 0) );
+ average_down(source_mf[lev], coarse_mf[lev], 0, ncomp, IntVect(coarse_ratio));
+ }
};
diff --git a/Source/Diagnostics/WarpXIO.cpp b/Source/Diagnostics/WarpXIO.cpp
index 1551a7df6..ee1991c21 100644
--- a/Source/Diagnostics/WarpXIO.cpp
+++ b/Source/Diagnostics/WarpXIO.cpp
@@ -470,13 +470,26 @@ WarpX::WritePlotFile () const
const std::string& plotfilename = amrex::Concatenate(plot_file,istep[0]);
amrex::Print() << " Writing plotfile " << plotfilename << "\n";
- // Average the fields from the simulation to the cell centers
+ // Average the fields from the simulation grid to the cell centers
const int ngrow = 0;
Vector<std::string> varnames; // Name of the written fields
// mf_avg will contain the averaged, cell-centered fields
Vector<MultiFab> mf_avg;
WarpX::AverageAndPackFields( varnames, mf_avg, ngrow );
+ // Coarsen the fields, if requested by the user
+ Vector<const MultiFab*> output_mf; // will point to the data to be written
+ Vector<MultiFab> coarse_mf; // will remain empty if there is no coarsening
+ Vector<Geometry> output_geom;
+ if (plot_coarsening_ratio != 1) {
+ coarsenCellCenteredFields( coarse_mf, output_geom, mf_avg, Geom(),
+ plot_coarsening_ratio, finest_level );
+ output_mf = amrex::GetVecOfConstPtrs(coarse_mf);
+ } else { // No averaging necessary, simply point to mf_avg
+ output_mf = amrex::GetVecOfConstPtrs(mf_avg);
+ output_geom = Geom();
+ }
+
// Write the fields contained in `mf_avg`, and corresponding to the
// names `varnames`, into a plotfile.
// Prepare extra directory (filled later), for the raw fields
@@ -484,11 +497,9 @@ WarpX::WritePlotFile () const
VisMF::Header::Version current_version = VisMF::GetHeaderVersion();
VisMF::SetHeaderVersion(plotfile_headerversion);
if (plot_raw_fields) rfs.emplace_back("raw_fields");
- amrex::WriteMultiLevelPlotfile(
- plotfilename, finest_level+1,
- amrex::GetVecOfConstPtrs(mf_avg),
- varnames, Geom(), t_new[0], istep,
- refRatio(),
+ amrex::WriteMultiLevelPlotfile(plotfilename, finest_level+1,
+ output_mf, varnames, output_geom,
+ t_new[0], istep, refRatio(),
"HyperCLaw-V1.1",
"Level_",
"Cell",
diff --git a/Source/Filter/BilinearFilter.H b/Source/Filter/BilinearFilter.H
new file mode 100644
index 000000000..f94f4467d
--- /dev/null
+++ b/Source/Filter/BilinearFilter.H
@@ -0,0 +1,26 @@
+#include <AMReX_MultiFab.H>
+#include <AMReX_Vector.H>
+
+#ifndef WARPX_FILTER_H_
+#define WARPX_FILTER_H_
+
+class BilinearFilter
+{
+public:
+ BilinearFilter () = default;
+
+ void ComputeStencils();
+ void ApplyStencil(amrex::MultiFab& dstmf, const amrex::MultiFab& srcmf, int scomp=0, int dcomp=0, int ncomp=10000);
+
+ amrex::IntVect npass_each_dir;
+ amrex::IntVect stencil_length_each_dir;
+
+private:
+
+ void Filter(const amrex::Box& tbx, amrex::FArrayBox const& tmpfab, amrex::FArrayBox &dstfab, int scomp, int dcomp, int ncomp);
+
+ amrex::Vector<amrex::Real> stencil_x, stencil_y, stencil_z;
+
+ amrex::Dim3 slen;
+};
+#endif
diff --git a/Source/Filter/BilinearFilter.cpp b/Source/Filter/BilinearFilter.cpp
new file mode 100644
index 000000000..17857a416
--- /dev/null
+++ b/Source/Filter/BilinearFilter.cpp
@@ -0,0 +1,156 @@
+#include <WarpX.H>
+#include <BilinearFilter.H>
+#include <WarpX_f.H>
+
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+
+using namespace amrex;
+
+namespace {
+void compute_stencil(Vector<Real> &stencil, int npass){
+ Vector<Real> old_s(1+npass,0.);
+ Vector<Real> new_s(1+npass,0.);
+
+ old_s[0] = 1.;
+ int jmax = 1;
+ amrex::Real loc;
+ // Convolve the filter with itself npass times
+ for(int ipass=1; ipass<npass+1; ipass++){
+ // element 0 has to be treated in its own way
+ new_s[0] = 0.5 * old_s[0];
+ if (1<jmax) new_s[0] += 0.5 * old_s[1];
+ loc = 0.;
+ // For each element j, apply the filter to
+ // old_s to get new_s[j]. loc stores the tmp
+ // filtered value.
+ for(int j=1; j<jmax+1; j++){
+ loc = 0.5 * old_s[j];
+ loc += 0.25 * old_s[j-1];
+ if (j<jmax) loc += 0.25 * old_s[j+1];
+ new_s[j] = loc;
+ }
+ // copy new_s into old_s
+ old_s = new_s;
+ // extend the stencil length for next iteration
+ jmax += 1;
+ }
+ // we use old_s here to make sure the stencil
+ // is corrent even when npass = 0
+ stencil = old_s;
+ stencil[0] *= 0.5; // because we will use it twice
+}
+}
+
+void BilinearFilter::ComputeStencils(){
+ BL_PROFILE("BilinearFilter::ComputeStencils()");
+ stencil_length_each_dir = npass_each_dir;
+ stencil_length_each_dir += 1.;
+#if (AMREX_SPACEDIM == 3)
+ // npass_each_dir = npass_x npass_y npass_z
+ stencil_x.resize( 1 + npass_each_dir[0] );
+ stencil_y.resize( 1 + npass_each_dir[1] );
+ stencil_z.resize( 1 + npass_each_dir[2] );
+ compute_stencil(stencil_x, npass_each_dir[0]);
+ compute_stencil(stencil_y, npass_each_dir[1]);
+ compute_stencil(stencil_z, npass_each_dir[2]);
+#elif (AMREX_SPACEDIM == 2)
+ // npass_each_dir = npass_x npass_z
+ stencil_x.resize( 1 + npass_each_dir[0] );
+ stencil_z.resize( 1 + npass_each_dir[1] );
+ compute_stencil(stencil_x, npass_each_dir[0]);
+ compute_stencil(stencil_z, npass_each_dir[1]);
+#endif
+ slen = stencil_length_each_dir.dim3();
+#if (AMREX_SPACEDIM == 2)
+ slen.z = 1;
+#endif
+}
+
+
+void
+BilinearFilter::ApplyStencil (MultiFab& dstmf, const MultiFab& srcmf, int scomp, int dcomp, int ncomp)
+{
+ BL_PROFILE("BilinearFilter::ApplyStencil()");
+ ncomp = std::min(ncomp, srcmf.nComp());
+#ifdef _OPENMP
+#pragma omp parallel
+#endif
+ {
+ FArrayBox tmpfab;
+ for (MFIter mfi(dstmf,true); mfi.isValid(); ++mfi){
+ const auto& srcfab = srcmf[mfi];
+ auto& dstfab = dstmf[mfi];
+ const Box& tbx = mfi.growntilebox();
+ const Box& gbx = amrex::grow(tbx,stencil_length_each_dir-1);
+ // tmpfab has enough ghost cells for the stencil
+ tmpfab.resize(gbx,ncomp);
+ tmpfab.setVal(0.0, gbx, 0, ncomp);
+ // Copy values in srcfab into tmpfab
+ const Box& ibx = gbx & srcfab.box();
+ tmpfab.copy(srcfab, ibx, scomp, ibx, 0, ncomp);
+ // Apply filter
+ Filter(tbx, tmpfab, dstfab, 0, dcomp, ncomp);
+ }
+ }
+}
+
+void BilinearFilter::Filter (const Box& tbx, FArrayBox const& tmpfab, FArrayBox &dstfab,
+ int scomp, int dcomp, int ncomp)
+{
+ const auto lo = amrex::lbound(tbx);
+ const auto hi = amrex::ubound(tbx);
+ const auto tmp = tmpfab.array();
+ const auto dst = dstfab.array();
+ // tmp and dst are of type Array4 (Fortran ordering)
+ amrex::Real const* AMREX_RESTRICT sx = stencil_x.data();
+ amrex::Real const* AMREX_RESTRICT sy = stencil_y.data();
+ amrex::Real const* AMREX_RESTRICT sz = stencil_z.data();
+ for (int n = 0; n < ncomp; ++n) {
+ // Set dst value to 0.
+ for (int k = lo.z; k <= hi.z; ++k) {
+ for (int j = lo.y; j <= hi.y; ++j) {
+ for (int i = lo.x; i <= hi.x; ++i) {
+ dst(i,j,k,dcomp+n) = 0.0;
+ }
+ }
+ }
+ // 3 nested loop on 3D stencil
+ for (int iz=0; iz < slen.z; ++iz){
+ for (int iy=0; iy < slen.y; ++iy){
+ for (int ix=0; ix < slen.x; ++ix){
+#if (AMREX_SPACEDIM == 3)
+ Real sss = sx[ix]*sy[iy]*sz[iz];
+#else
+ Real sss = sx[ix]*sz[iy];
+#endif
+ // 3 nested loop on 3D array
+ for (int k = lo.z; k <= hi.z; ++k) {
+ for (int j = lo.y; j <= hi.y; ++j) {
+ AMREX_PRAGMA_SIMD
+ for (int i = lo.x; i <= hi.x; ++i) {
+#if (AMREX_SPACEDIM == 3)
+ dst(i,j,k,dcomp+n) += sss*(tmp(i-ix,j-iy,k-iz,scomp+n)
+ +tmp(i+ix,j-iy,k-iz,scomp+n)
+ +tmp(i-ix,j+iy,k-iz,scomp+n)
+ +tmp(i+ix,j+iy,k-iz,scomp+n)
+ +tmp(i-ix,j-iy,k+iz,scomp+n)
+ +tmp(i+ix,j-iy,k+iz,scomp+n)
+ +tmp(i-ix,j+iy,k+iz,scomp+n)
+ +tmp(i+ix,j+iy,k+iz,scomp+n));
+#else
+ dst(i,j,k,dcomp+n) += sss*(tmp(i-ix,j-iy,k,scomp+n)
+ +tmp(i+ix,j-iy,k,scomp+n)
+ +tmp(i-ix,j+iy,k,scomp+n)
+ +tmp(i+ix,j+iy,k,scomp+n));
+#endif
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+}
+
diff --git a/Source/Filter/Make.package b/Source/Filter/Make.package
index 56c0e7852..36e74bfb0 100644
--- a/Source/Filter/Make.package
+++ b/Source/Filter/Make.package
@@ -1,4 +1,5 @@
-F90EXE_sources += filter_module.F90
+CEXE_sources += BilinearFilter.cpp
+CEXE_headers += BilinearFilter.H
INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Filter
VPATH_LOCATIONS += $(WARPX_HOME)/Source/Filter
diff --git a/Source/Filter/filter_module.F90 b/Source/Filter/filter_module.F90
deleted file mode 100644
index ea87be7ca..000000000
--- a/Source/Filter/filter_module.F90
+++ /dev/null
@@ -1,87 +0,0 @@
-module warpx_filter_module
-
- use iso_c_binding
- use amrex_fort_module, only : rt => amrex_real
- use amrex_constants_module, only : fourth, eighth
-
- implicit none
-
-contains
-
- subroutine warpx_filter_3d(lo, hi, src, slo, shi, dst, dlo, dhi, nc) &
- bind(c, name='warpx_filter_3d')
- integer, intent(in), value :: nc
- integer, dimension(3), intent(in) :: lo, hi, dlo, dhi, slo, shi
- real(rt), intent(inout) :: dst(dlo(1):dhi(1), dlo(2):dhi(2), dlo(3):dhi(3),nc)
- real(rt), intent(in ) :: src(slo(1):shi(1), slo(2):shi(2), slo(3):shi(3),nc)
-
- real(rt), parameter :: sixteenth = 1.d0/16.d0
- real(rt), parameter :: thirtysecond = 1.d0/32.d0
- real(rt), parameter :: sixtyfourth = 1.d0/64.d0
- integer :: i,j,k,c
-
- do c = 1, nc
- do k = lo(3), hi(3)
- do j = lo(2), hi(2)
- do i = lo(1), hi(1)
- dst(i,j,k,c) = sixtyfourth*(src(i-1,j-1,k-1,c) &
- & +src(i+1,j-1,k-1,c) &
- & +src(i-1,j+1,k-1,c) &
- & +src(i+1,j+1,k-1,c) &
- & +src(i-1,j-1,k+1,c) &
- & +src(i+1,j-1,k+1,c) &
- & +src(i-1,j+1,k+1,c) &
- & +src(i+1,j+1,k+1,c)) &
- + thirtysecond*(src(i-1,j-1,k ,c) &
- & +src(i+1,j-1,k ,c) &
- & +src(i-1,j+1,k ,c) &
- & +src(i+1,j+1,k ,c) &
- & +src(i-1,j ,k-1,c) &
- & +src(i+1,j ,k-1,c) &
- & +src(i-1,j ,k+1,c) &
- & +src(i+1,j ,k+1,c) &
- & +src(i ,j-1,k-1,c) &
- & +src(i ,j+1,k-1,c) &
- & +src(i ,j-1,k+1,c) &
- & +src(i ,j+1,k+1,c)) &
- + sixteenth *(src(i-1,j ,k ,c) &
- & +src(i+1,j ,k ,c) &
- & +src(i ,j-1,k ,c) &
- & +src(i ,j+1,k ,c) &
- & +src(i ,j ,k-1,c) &
- & +src(i ,j ,k+1,c)) &
- + eighth * src(i ,j ,k ,c)
- end do
- end do
- end do
- end do
- end subroutine warpx_filter_3d
-
- subroutine warpx_filter_2d(lo, hi, src, slo, shi, dst, dlo, dhi, nc) &
- bind(c, name='warpx_filter_2d')
- integer, intent(in), value :: nc
- integer, dimension(2), intent(in) :: lo, hi, dlo, dhi, slo, shi
- real(rt), intent(inout) :: dst(dlo(1):dhi(1), dlo(2):dhi(2), nc)
- real(rt), intent(in ) :: src(slo(1):shi(1), slo(2):shi(2), nc)
-
- real(rt), parameter :: sixteenth = 1.d0/16.d0
- integer :: i,j,comp
-
- do comp = 1, nc
- do j = lo(2), hi(2)
- do i = lo(1), hi(1)
- dst(i, j, comp) = sixteenth*(src(i-1,j-1,comp) &
- & +src(i+1,j-1,comp) &
- & +src(i-1,j+1,comp) &
- & +src(i+1,j+1,comp)) &
- + eighth *(src(i ,j-1,comp) &
- & +src(i ,j+1,comp) &
- & +src(i-1,j ,comp) &
- & +src(i+1,j ,comp)) &
- + fourth * src(i ,j ,comp)
- end do
- end do
- end do
- end subroutine warpx_filter_2d
-
-end module warpx_filter_module
diff --git a/Source/FortranInterface/WarpX_f.H b/Source/FortranInterface/WarpX_f.H
index ba9c4b044..1ccf82f87 100644
--- a/Source/FortranInterface/WarpX_f.H
+++ b/Source/FortranInterface/WarpX_f.H
@@ -40,7 +40,6 @@
#define WRPX_PUSH_LEAPFROG_POSITIONS warpx_push_leapfrog_positions_3d
#define WRPX_LORENTZ_TRANSFORM_Z warpx_lorentz_transform_z
-#define WRPX_FILTER warpx_filter_3d
#define WRPX_COPY_SLICE warpx_copy_slice_3d
#define WRPX_PXR_NCI_CORR_INIT init_godfrey_filter_coeffs
#define WRPX_PXR_GODFREY_FILTER apply_godfrey_filter_z_3d
@@ -71,7 +70,6 @@
#define WRPX_PUSH_LEAPFROG_POSITIONS warpx_push_leapfrog_positions_2d
#define WRPX_LORENTZ_TRANSFORM_Z warpx_lorentz_transform_z
-#define WRPX_FILTER warpx_filter_2d
#define WRPX_COPY_SLICE warpx_copy_slice_2d
#define WRPX_PXR_NCI_CORR_INIT init_godfrey_filter_coeffs
#define WRPX_PXR_GODFREY_FILTER apply_godfrey_filter_z_2d
@@ -450,10 +448,6 @@ extern "C"
const BL_FORT_FAB_ARG_ANYD(fine),
const int* ncomp);
- void WRPX_FILTER (const int* lo, const int* hi,
- const amrex_real*, const int*, const int*,
- amrex_real*, const int*, const int*, int);
-
void WRPX_PXR_NCI_CORR_INIT(amrex::Real*, amrex::Real*, const int,
const amrex::Real, const int);
diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp
index ff5442b00..9190b5a61 100644
--- a/Source/Initialization/WarpXInitData.cpp
+++ b/Source/Initialization/WarpXInitData.cpp
@@ -6,6 +6,7 @@
#include <WarpX.H>
#include <WarpX_f.H>
+#include <BilinearFilter.H>
#ifdef BL_USE_SENSEI_INSITU
#include <AMReX_AmrMeshInSituBridge.H>
@@ -21,7 +22,7 @@ WarpX::InitData ()
if (restart_chkfile.empty())
{
ComputeDt();
- InitFromScratch();
+ InitFromScratch();
}
else
{
@@ -38,6 +39,10 @@ WarpX::InitData ()
WarpX::InitNCICorrector();
}
+ if (WarpX::use_filter) {
+ WarpX::InitFilter();
+ }
+
BuildBufferMasks();
InitDiagnostics();
@@ -178,6 +183,14 @@ WarpX::InitNCICorrector ()
}
void
+WarpX::InitFilter (){
+ if (WarpX::use_filter){
+ WarpX::bilinear_filter.npass_each_dir = WarpX::filter_npass_each_dir;
+ WarpX::bilinear_filter.ComputeStencils();
+ }
+}
+
+void
WarpX::PostRestart ()
{
#ifdef WARPX_USE_PSATD
diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp
index 09767c265..28971eb0c 100644
--- a/Source/Parallelization/WarpXComm.cpp
+++ b/Source/Parallelization/WarpXComm.cpp
@@ -1,4 +1,3 @@
-
#include <WarpX.H>
#include <WarpX_f.H>
@@ -381,35 +380,42 @@ WarpX::SyncCurrent ()
if (WarpX::use_filter) {
for (int lev = 0; lev <= finest_level; ++lev) {
IntVect ng = current_fp[lev][0]->nGrowVect();
- ng += 1;
+ ng += bilinear_filter.stencil_length_each_dir-1;
for (int idim = 0; idim < 3; ++idim) {
+ // Create new MultiFab j_jp with enough guard cells for the
+ // (potentially large) stencil of the multi-pass bilinear filter.
j_fp[lev][idim].reset(new MultiFab(current_fp[lev][idim]->boxArray(),
current_fp[lev][idim]->DistributionMap(),
1, ng));
- applyFilter(*j_fp[lev][idim], *current_fp[lev][idim]);
+ // Apply the filter to current_fp, store the result in j_fp.
+ bilinear_filter.ApplyStencil(*j_fp[lev][idim], *current_fp[lev][idim]);
+ // Then swap j_fp and current_fp
std::swap(j_fp[lev][idim], current_fp[lev][idim]);
+ // At this point, current_fp may have false values close to the
+ // edges of each FAB. This will be solved with a SumBoundary later.
+ // j_fp contains the exact MultiFab current_fp before this step.
}
}
for (int lev = 1; lev <= finest_level; ++lev) {
IntVect ng = current_cp[lev][0]->nGrowVect();
- ng += 1;
+ ng += bilinear_filter.stencil_length_each_dir-1;
for (int idim = 0; idim < 3; ++idim) {
j_cp[lev][idim].reset(new MultiFab(current_cp[lev][idim]->boxArray(),
current_cp[lev][idim]->DistributionMap(),
1, ng));
- applyFilter(*j_cp[lev][idim], *current_cp[lev][idim]);
+ bilinear_filter.ApplyStencil(*j_cp[lev][idim], *current_cp[lev][idim]);
std::swap(j_cp[lev][idim], current_cp[lev][idim]);
}
}
for (int lev = 1; lev <= finest_level; ++lev) {
if (current_buf[lev][0]) {
IntVect ng = current_buf[lev][0]->nGrowVect();
- ng += 1;
+ ng += bilinear_filter.stencil_length_each_dir-1;
for (int idim = 0; idim < 3; ++idim) {
j_buf[lev][idim].reset(new MultiFab(current_buf[lev][idim]->boxArray(),
current_buf[lev][idim]->DistributionMap(),
1, ng));
- applyFilter(*j_buf[lev][idim], *current_buf[lev][idim]);
+ bilinear_filter.ApplyStencil(*j_buf[lev][idim], *current_buf[lev][idim]);
std::swap(*j_buf[lev][idim], *current_buf[lev][idim]);
}
}
@@ -420,6 +426,9 @@ WarpX::SyncCurrent ()
for (int lev = 0; lev <= finest_level; ++lev)
{
const auto& period = Geom(lev).periodicity();
+ // When using a bilinear filter with many passes, current_fp has
+ // temporarily more ghost cells here, so that its value inside
+ // the domain is correct at the end of this stage.
current_fp[lev][0]->SumBoundary(period);
current_fp[lev][1]->SumBoundary(period);
current_fp[lev][2]->SumBoundary(period);
@@ -461,8 +470,14 @@ WarpX::SyncCurrent ()
for (int lev = 0; lev <= finest_level; ++lev)
{
for (int idim = 0; idim < 3; ++idim) {
+ // swap j_fp and current_fp so that j_fp has correct values inside
+ // the domain and wrong number of ghost cells.
+ // current_fp has right number of ghost cells.
std::swap(j_fp[lev][idim], current_fp[lev][idim]);
+ // Then copy the interior of j_fp to current_fp.
MultiFab::Copy(*current_fp[lev][idim], *j_fp[lev][idim], 0, 0, 1, 0);
+ // current_fp has right number of ghost cells and
+ // correct filtered values here.
}
}
for (int lev = 1; lev <= finest_level; ++lev)
@@ -557,32 +572,32 @@ WarpX::SyncRho (amrex::Vector<std::unique_ptr<amrex::MultiFab> >& rhof,
for (int lev = 0; lev <= finest_level; ++lev) {
const int ncomp = rhof[lev]->nComp();
IntVect ng = rhof[lev]->nGrowVect();
- ng += 1;
+ ng += bilinear_filter.stencil_length_each_dir-1;
rho_f_g[lev].reset(new MultiFab(rhof[lev]->boxArray(),
rhof[lev]->DistributionMap(),
ncomp, ng));
- applyFilter(*rho_f_g[lev], *rhof[lev]);
+ bilinear_filter.ApplyStencil(*rho_f_g[lev], *rhof[lev]);
std::swap(rho_f_g[lev], rhof[lev]);
}
for (int lev = 1; lev <= finest_level; ++lev) {
const int ncomp = rhoc[lev]->nComp();
IntVect ng = rhoc[lev]->nGrowVect();
- ng += 1;
+ ng += bilinear_filter.stencil_length_each_dir-1;
rho_c_g[lev].reset(new MultiFab(rhoc[lev]->boxArray(),
rhoc[lev]->DistributionMap(),
ncomp, ng));
- applyFilter(*rho_c_g[lev], *rhoc[lev]);
+ bilinear_filter.ApplyStencil(*rho_c_g[lev], *rhoc[lev]);
std::swap(rho_c_g[lev], rhoc[lev]);
}
for (int lev = 1; lev <= finest_level; ++lev) {
if (charge_buf[lev]) {
const int ncomp = charge_buf[lev]->nComp();
IntVect ng = charge_buf[lev]->nGrowVect();
- ng += 1;
+ ng += bilinear_filter.stencil_length_each_dir-1;
rho_buf_g[lev].reset(new MultiFab(charge_buf[lev]->boxArray(),
charge_buf[lev]->DistributionMap(),
ncomp, ng));
- applyFilter(*rho_buf_g[lev], *charge_buf[lev]);
+ bilinear_filter.ApplyStencil(*rho_buf_g[lev], *charge_buf[lev]);
std::swap(*rho_buf_g[lev], *charge_buf[lev]);
}
}
@@ -711,9 +726,9 @@ WarpX::ApplyFilterandSumBoundaryJ (int lev, PatchType patch_type)
for (int idim = 0; idim < 3; ++idim) {
if (use_filter) {
IntVect ng = j[idim]->nGrowVect();
- ng += 1;
+ ng += bilinear_filter.stencil_length_each_dir-1;
MultiFab jf(j[idim]->boxArray(), j[idim]->DistributionMap(), 1, ng);
- applyFilter(jf, *j[idim]);
+ bilinear_filter.ApplyStencil(jf, *j[idim]);
jf.SumBoundary(period);
MultiFab::Copy(*j[idim], jf, 0, 0, 1, 0);
} else {
@@ -752,15 +767,15 @@ WarpX::AddCurrentFromFineLevelandSumBoundary (int lev)
{
// coarse patch of fine level
IntVect ng = current_cp[lev+1][idim]->nGrowVect();
- ng += 1;
+ ng += bilinear_filter.stencil_length_each_dir-1;
MultiFab jfc(current_cp[lev+1][idim]->boxArray(),
current_cp[lev+1][idim]->DistributionMap(), 1, ng);
- applyFilter(jfc, *current_cp[lev+1][idim]);
+ bilinear_filter.ApplyStencil(jfc, *current_cp[lev+1][idim]);
// buffer patch of fine level
MultiFab jfb(current_buf[lev+1][idim]->boxArray(),
current_buf[lev+1][idim]->DistributionMap(), 1, ng);
- applyFilter(jfb, *current_buf[lev+1][idim]);
+ bilinear_filter.ApplyStencil(jfb, *current_buf[lev+1][idim]);
MultiFab::Add(jfb, jfc, 0, 0, 1, ng);
mf.ParallelAdd(jfb, 0, 0, 1, ng, IntVect::TheZeroVector(), period);
@@ -772,10 +787,10 @@ WarpX::AddCurrentFromFineLevelandSumBoundary (int lev)
{
// coarse patch of fine level
IntVect ng = current_cp[lev+1][idim]->nGrowVect();
- ng += 1;
+ ng += bilinear_filter.stencil_length_each_dir-1;
MultiFab jf(current_cp[lev+1][idim]->boxArray(),
current_cp[lev+1][idim]->DistributionMap(), 1, ng);
- applyFilter(jf, *current_cp[lev+1][idim]);
+ bilinear_filter.ApplyStencil(jf, *current_cp[lev+1][idim]);
mf.ParallelAdd(jf, 0, 0, 1, ng, IntVect::TheZeroVector(), period);
jf.SumBoundary(period);
MultiFab::Copy(*current_cp[lev+1][idim], jf, 0, 0, 1, 0);
@@ -822,9 +837,9 @@ WarpX::ApplyFilterandSumBoundaryRho (int lev, PatchType patch_type, int icomp, i
if (r == nullptr) return;
if (use_filter) {
IntVect ng = r->nGrowVect();
- ng += 1;
+ ng += bilinear_filter.stencil_length_each_dir-1;
MultiFab rf(r->boxArray(), r->DistributionMap(), ncomp, ng);
- applyFilter(rf, *r, icomp, 0, ncomp);
+ bilinear_filter.ApplyStencil(rf, *r, icomp, 0, ncomp);
rf.SumBoundary(period);
MultiFab::Copy(*r, rf, 0, icomp, ncomp, 0);
} else {
@@ -858,15 +873,15 @@ WarpX::AddRhoFromFineLevelandSumBoundary(int lev, int icomp, int ncomp)
{
// coarse patch of fine level
IntVect ng = rho_cp[lev+1]->nGrowVect();
- ng += 1;
+ ng += bilinear_filter.stencil_length_each_dir-1;
MultiFab rhofc(rho_cp[lev+1]->boxArray(),
rho_cp[lev+1]->DistributionMap(), ncomp, ng);
- applyFilter(rhofc, *rho_cp[lev+1], icomp, 0, ncomp);
+ bilinear_filter.ApplyStencil(rhofc, *rho_cp[lev+1], icomp, 0, ncomp);
// buffer patch of fine level
MultiFab rhofb(charge_buf[lev+1]->boxArray(),
charge_buf[lev+1]->DistributionMap(), ncomp, ng);
- applyFilter(rhofb, *charge_buf[lev+1], icomp, 0, ncomp);
+ bilinear_filter.ApplyStencil(rhofb, *charge_buf[lev+1], icomp, 0, ncomp);
MultiFab::Add(rhofb, rhofc, 0, 0, ncomp, ng);
mf.ParallelAdd(rhofb, 0, 0, ncomp, ng, IntVect::TheZeroVector(), period);
@@ -877,9 +892,9 @@ WarpX::AddRhoFromFineLevelandSumBoundary(int lev, int icomp, int ncomp)
else if (use_filter) // but no buffer
{
IntVect ng = rho_cp[lev+1]->nGrowVect();
- ng += 1;
+ ng += bilinear_filter.stencil_length_each_dir-1;
MultiFab rf(rho_cp[lev+1]->boxArray(), rho_cp[lev+1]->DistributionMap(), ncomp, ng);
- applyFilter(rf, *rho_cp[lev+1], icomp, 0, ncomp);
+ bilinear_filter.ApplyStencil(rf, *rho_cp[lev+1], icomp, 0, ncomp);
mf.ParallelAdd(rf, 0, 0, ncomp, ng, IntVect::TheZeroVector(), period);
rf.SumBoundary(0, ncomp, period);
MultiFab::Copy(*rho_cp[lev+1], rf, 0, icomp, ncomp, 0);
diff --git a/Source/WarpX.H b/Source/WarpX.H
index 5ee2b518c..79488661b 100644
--- a/Source/WarpX.H
+++ b/Source/WarpX.H
@@ -22,6 +22,7 @@
#include <MultiParticleContainer.H>
#include <PML.H>
#include <BoostedFrameDiagnostic.H>
+#include <BilinearFilter.H>
#ifdef WARPX_USE_PSATD
#include <fftw3.h>
@@ -140,6 +141,9 @@ public:
}
}
+ static amrex::IntVect filter_npass_each_dir;
+ BilinearFilter bilinear_filter;
+
void ComputeDt ();
int MoveWindow (bool move_j);
void UpdatePlasmaInjectionPosition (amrex::Real dt);
@@ -342,6 +346,8 @@ private:
void InitPML ();
void ComputePMLFactors ();
+ void InitFilter ();
+
void InitDiagnostics ();
void InitNCICorrector ();
@@ -383,9 +389,6 @@ private:
void LoadBalance ();
- static void applyFilter (amrex::MultiFab& dstmf, const amrex::MultiFab& srcmf,
- int scomp = 0, int dcomp = 0, int ncomp = 10000);
-
void BuildBufferMasks ();
const amrex::iMultiFab* getCurrentBufferMasks (int lev) const {
return current_buffer_masks[lev].get();
@@ -510,6 +513,7 @@ private:
bool plot_crsepatch = false;
bool plot_raw_fields = false;
bool plot_raw_fields_guards = false;
+ int plot_coarsening_ratio = 1;
amrex::VisMF::Header::Version checkpoint_headerversion = amrex::VisMF::Header::NoFabHeader_v1;
// amrex::VisMF::Header::Version plotfile_headerversion = amrex::VisMF::Header::NoFabHeader_v1;
diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp
index cc721ea03..dd6c94945 100644
--- a/Source/WarpX.cpp
+++ b/Source/WarpX.cpp
@@ -91,6 +91,8 @@ IntVect WarpX::jy_nodal_flag(1,1); // y is the missing dimension to 2D AMReX
IntVect WarpX::jz_nodal_flag(1,0); // z is the second dimension to 2D AMReX
#endif
+IntVect WarpX::filter_npass_each_dir(1);
+
int WarpX::n_field_gather_buffer = 0;
int WarpX::n_current_deposition_buffer = -1;
@@ -98,8 +100,6 @@ int WarpX::do_nodal = false;
WarpX* WarpX::m_instance = nullptr;
-
-
WarpX&
WarpX::GetInstance ()
{
@@ -260,14 +260,14 @@ WarpX::ReadParameters ()
pp.query("cfl", cfl);
pp.query("verbose", verbose);
pp.query("regrid_int", regrid_int);
- pp.query("do_subcycling", do_subcycling);
-
- AMREX_ALWAYS_ASSERT_WITH_MESSAGE(do_subcycling != 1 || max_level <= 1,
- "Subcycling method 1 only works for 2 levels.");
-
- ReadBoostedFrameParameters(gamma_boost, beta_boost, boost_direction);
-
- pp.queryarr("B_external", B_external);
+ pp.query("do_subcycling", do_subcycling);
+
+ AMREX_ALWAYS_ASSERT_WITH_MESSAGE(do_subcycling != 1 || max_level <= 1,
+ "Subcycling method 1 only works for 2 levels.");
+
+ ReadBoostedFrameParameters(gamma_boost, beta_boost, boost_direction);
+
+ pp.queryarr("B_external", B_external);
pp.query("do_moving_window", do_moving_window);
if (do_moving_window)
@@ -298,57 +298,67 @@ WarpX::ReadParameters ()
pp.query("do_plasma_injection", do_plasma_injection);
if (do_plasma_injection) {
- pp.get("num_injected_species", num_injected_species);
- injected_plasma_species.resize(num_injected_species);
- pp.getarr("injected_plasma_species", injected_plasma_species,
- 0, num_injected_species);
- if (moving_window_v >= 0){
- // Inject particles continuously from the right end of the box
- current_injection_position = geom[0].ProbHi(moving_window_dir);
- } else {
- // Inject particles continuously from the left end of the box
- current_injection_position = geom[0].ProbLo(moving_window_dir);
- }
+ pp.get("num_injected_species", num_injected_species);
+ injected_plasma_species.resize(num_injected_species);
+ pp.getarr("injected_plasma_species", injected_plasma_species,
+ 0, num_injected_species);
+ if (moving_window_v >= 0){
+ // Inject particles continuously from the right end of the box
+ current_injection_position = geom[0].ProbHi(moving_window_dir);
+ } else {
+ // Inject particles continuously from the left end of the box
+ current_injection_position = geom[0].ProbLo(moving_window_dir);
+ }
}
- pp.query("do_boosted_frame_diagnostic", do_boosted_frame_diagnostic);
- if (do_boosted_frame_diagnostic) {
-
- AMREX_ALWAYS_ASSERT_WITH_MESSAGE(gamma_boost > 1.0,
- "gamma_boost must be > 1 to use the boosted frame diagnostic.");
-
- std::string s;
- pp.get("boost_direction", s);
- AMREX_ALWAYS_ASSERT_WITH_MESSAGE( (s == "z" || s == "Z"),
- "The boosted frame diagnostic currently only works if the boost is in the z direction.");
-
- pp.get("num_snapshots_lab", num_snapshots_lab);
- pp.get("dt_snapshots_lab", dt_snapshots_lab);
- pp.get("gamma_boost", gamma_boost);
-
- pp.query("do_boosted_frame_fields", do_boosted_frame_fields);
- pp.query("do_boosted_frame_particles", do_boosted_frame_particles);
-
-
- AMREX_ALWAYS_ASSERT_WITH_MESSAGE(do_moving_window,
- "The moving window should be on if using the boosted frame diagnostic.");
-
- pp.get("moving_window_dir", s);
- AMREX_ALWAYS_ASSERT_WITH_MESSAGE( (s == "z" || s == "Z"),
- "The boosted frame diagnostic currently only works if the moving window is in the z direction.");
- }
+ pp.query("do_boosted_frame_diagnostic", do_boosted_frame_diagnostic);
+ if (do_boosted_frame_diagnostic) {
+
+ AMREX_ALWAYS_ASSERT_WITH_MESSAGE(gamma_boost > 1.0,
+ "gamma_boost must be > 1 to use the boosted frame diagnostic.");
+
+ std::string s;
+ pp.get("boost_direction", s);
+ AMREX_ALWAYS_ASSERT_WITH_MESSAGE( (s == "z" || s == "Z"),
+ "The boosted frame diagnostic currently only works if the boost is in the z direction.");
+
+ pp.get("num_snapshots_lab", num_snapshots_lab);
+ pp.get("dt_snapshots_lab", dt_snapshots_lab);
+ pp.get("gamma_boost", gamma_boost);
+
+ pp.query("do_boosted_frame_fields", do_boosted_frame_fields);
+ pp.query("do_boosted_frame_particles", do_boosted_frame_particles);
+
+
+ AMREX_ALWAYS_ASSERT_WITH_MESSAGE(do_moving_window,
+ "The moving window should be on if using the boosted frame diagnostic.");
+
+ pp.get("moving_window_dir", s);
+ AMREX_ALWAYS_ASSERT_WITH_MESSAGE( (s == "z" || s == "Z"),
+ "The boosted frame diagnostic currently only works if the moving window is in the z direction.");
+ }
- pp.query("do_electrostatic", do_electrostatic);
- pp.query("n_buffer", n_buffer);
- pp.query("const_dt", const_dt);
+ pp.query("do_electrostatic", do_electrostatic);
+ pp.query("n_buffer", n_buffer);
+ pp.query("const_dt", const_dt);
pp.query("use_laser", use_laser);
+ // Read filter and fill IntVect filter_npass_each_dir with
+ // proper size for AMREX_SPACEDIM
pp.query("use_filter", use_filter);
+ Vector<int> parse_filter_npass_each_dir(AMREX_SPACEDIM,1);
+ pp.queryarr("filter_npass_each_dir", parse_filter_npass_each_dir);
+ filter_npass_each_dir[0] = parse_filter_npass_each_dir[0];
+ filter_npass_each_dir[1] = parse_filter_npass_each_dir[1];
+#if (AMREX_SPACEDIM == 3)
+ filter_npass_each_dir[2] = parse_filter_npass_each_dir[2];
+#endif
+
pp.query("serialize_ics", serialize_ics);
pp.query("refine_plasma", refine_plasma);
- pp.query("do_dive_cleaning", do_dive_cleaning);
- pp.query("n_field_gather_buffer", n_field_gather_buffer);
- pp.query("n_current_deposition_buffer", n_current_deposition_buffer);
+ pp.query("do_dive_cleaning", do_dive_cleaning);
+ pp.query("n_field_gather_buffer", n_field_gather_buffer);
+ pp.query("n_current_deposition_buffer", n_current_deposition_buffer);
pp.query("sort_int", sort_int);
pp.query("do_pml", do_pml);
@@ -368,6 +378,16 @@ WarpX::ReadParameters ()
pp.query("plot_divb" , plot_divb);
pp.query("plot_rho" , plot_rho);
pp.query("plot_F" , plot_F);
+ pp.query("plot_coarsening_ratio", plot_coarsening_ratio);
+ // Check that the coarsening_ratio can divide the blocking factor
+ for (int lev=0; lev<maxLevel(); lev++){
+ for (int comp=0; comp<AMREX_SPACEDIM; comp++){
+ if ( blockingFactor(lev)[comp] % plot_coarsening_ratio != 0 ){
+ amrex::Abort("plot_coarsening_ratio should be an integer divisor of the blocking factor.");
+ }
+ }
+ }
+
if (plot_F){
AMREX_ALWAYS_ASSERT_WITH_MESSAGE(do_dive_cleaning,
"plot_F only works if warpx.do_dive_cleaning = 1");
@@ -946,33 +966,6 @@ WarpX::ComputeDivE (MultiFab& divE, int dcomp,
}
void
-WarpX::applyFilter (MultiFab& dstmf, const MultiFab& srcmf, int scomp, int dcomp, int ncomp)
-{
- ncomp = std::min(ncomp, srcmf.nComp());
-#ifdef _OPENMP
-#pragma omp parallel
-#endif
- {
- FArrayBox tmpfab;
- for (MFIter mfi(dstmf,true); mfi.isValid(); ++mfi)
- {
- const auto& srcfab = srcmf[mfi];
- auto& dstfab = dstmf[mfi];
- const Box& tbx = mfi.growntilebox();
- const Box& gbx = amrex::grow(tbx,1);
- tmpfab.resize(gbx,ncomp);
- tmpfab.setVal(0.0, gbx, 0, ncomp);
- const Box& ibx = gbx & srcfab.box();
- tmpfab.copy(srcfab, ibx, scomp, ibx, 0, ncomp);
- WRPX_FILTER(BL_TO_FORTRAN_BOX(tbx),
- BL_TO_FORTRAN_ANYD(tmpfab),
- BL_TO_FORTRAN_N_ANYD(dstfab,dcomp),
- ncomp);
- }
- }
-}
-
-void
WarpX::BuildBufferMasks ()
{
int ngbuffer = std::max(n_field_gather_buffer, n_current_deposition_buffer);
generated by cgit v1.2.3 (git 2.46.0) at 2025-08-19 10:11:03 +0000