Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Implicit add filtering #5086

Merged
Merged
Show file tree
Hide file tree
Changes from 3 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion Docs/source/developers/fields.rst
Original file line number Diff line number Diff line change
Expand Up @@ -112,7 +112,7 @@ Bilinear filter

The multi-pass bilinear filter (applied on the current density) is implemented in ``Source/Filter/``, and class ``WarpX`` holds an instance of this class in member variable ``WarpX::bilinear_filter``. For performance reasons (to avoid creating too many guard cells), this filter is directly applied in communication routines, see ``WarpX::AddCurrentFromFineLevelandSumBoundary`` above and

.. doxygenfunction:: WarpX::ApplyFilterJ(const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>, 3>> &current, int lev, int idim)
.. doxygenfunction:: WarpX::ApplyFilterMF(const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>, 3>> &mfvec, int lev, int idim)

.. doxygenfunction:: WarpX::SumBoundaryJ(const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>, 3>> &current, int lev, int idim, const amrex::Periodicity &period)

Expand Down
Original file line number Diff line number Diff line change
@@ -0,0 +1,31 @@
{
"lev=0": {
"Bx": 60986.14027968379,
"By": 81807.16925197575,
"Bz": 51620.23830412359,
"Ex": 56569963234064.06,
"Ey": 15735696324387.898,
"Ez": 53301895885240.7,
"divE": 5.10172977366128e+22,
"jx": 1.7511741282649858e+19,
"jy": 2.1441037234086957e+19,
"jz": 1.8956446467052007e+19,
"rho": 451716735861.51776
},
"electrons": {
"particle_momentum_x": 4.864934380680235e-19,
"particle_momentum_y": 4.862649855848387e-19,
"particle_momentum_z": 4.884107133942756e-19,
"particle_position_x": 0.004250258285605058,
"particle_position_y": 0.004250550955854736,
"particle_weight": 2823958719279159.5
},
"protons": {
"particle_momentum_x": 2.088034298114596e-17,
"particle_momentum_y": 2.0869487719950244e-17,
"particle_momentum_z": 2.094045881084313e-17,
"particle_position_x": 0.004251274413462257,
"particle_position_y": 0.004251274300177191,
"particle_weight": 2823958719279159.5
}
}
14 changes: 14 additions & 0 deletions Regression/WarpX-tests.ini
Original file line number Diff line number Diff line change
Expand Up @@ -3894,6 +3894,20 @@ useOMP = 0
numthreads = 1
analysisRoutine = Examples/Tests/Implicit/analysis_vandb_jfnk_2d.py

[ThetaImplicitJFNK_VandB_2d_filtered]
buildDir = .
inputFile = Examples/Tests/Implicit/inputs_vandb_jfnk_2d
runtime_params = warpx.abort_on_warning_threshold=high warpx.use_filter=1
dim = 2
addToCompileString =
cmakeSetupOpts = -DWarpX_DIMS=2
restartTest = 0
useMPI = 1
numprocs = 2
useOMP = 0
numthreads = 1
analysisRoutine = Examples/Tests/Implicit/analysis_vandb_jfnk_2d.py

[SemiImplicitPicard_1d]
buildDir = .
inputFile = Examples/Tests/Implicit/inputs_1d_semiimplicit
Expand Down
6 changes: 3 additions & 3 deletions Source/Evolve/WarpXEvolve.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -569,7 +569,7 @@ void WarpX::SyncCurrentAndRho ()
{
// TODO This works only without mesh refinement
const int lev = 0;
if (use_filter) { ApplyFilterJ(current_fp_vay, lev); }
if (use_filter) { ApplyFilterMF(current_fp_vay, lev); }
}
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should ApplyFilterJ() be different than applying a filter for E or B when using curvilinear geometries?

I believe the filter should be applied to the current before dividing by the volume element to convert to current density. Perhaps this is being done somewhere else?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You are probably correct, but filtering doesn't work for RZ as it is now (at least for FDTD). As it is now, the inverse volume scaling is done before the filtering. That order could be changed but that would need to be a separate PR.

}
}
Expand Down Expand Up @@ -803,7 +803,7 @@ WarpX::OneStep_sub1 (Real cur_time)
PushParticlesandDeposit(fine_lev, cur_time, DtType::FirstHalf);
RestrictCurrentFromFineToCoarsePatch(current_fp, current_cp, fine_lev);
RestrictRhoFromFineToCoarsePatch(rho_fp, rho_cp, fine_lev);
if (use_filter) { ApplyFilterJ(current_fp, fine_lev); }
if (use_filter) { ApplyFilterMF(current_fp, fine_lev); }
SumBoundaryJ(current_fp, fine_lev, Geom(fine_lev).periodicity());
ApplyFilterandSumBoundaryRho(rho_fp, rho_cp, fine_lev, PatchType::fine, 0, 2*ncomps);

Expand Down Expand Up @@ -864,7 +864,7 @@ WarpX::OneStep_sub1 (Real cur_time)
PushParticlesandDeposit(fine_lev, cur_time + dt[fine_lev], DtType::SecondHalf);
RestrictCurrentFromFineToCoarsePatch(current_fp, current_cp, fine_lev);
RestrictRhoFromFineToCoarsePatch(rho_fp, rho_cp, fine_lev);
if (use_filter) { ApplyFilterJ(current_fp, fine_lev); }
if (use_filter) { ApplyFilterMF(current_fp, fine_lev); }
SumBoundaryJ(current_fp, fine_lev, Geom(fine_lev).periodicity());
ApplyFilterandSumBoundaryRho(rho_fp, rho_cp, fine_lev, PatchType::fine, 0, ncomps);

Expand Down
2 changes: 2 additions & 0 deletions Source/FieldSolver/ImplicitSolvers/WarpXImplicitOps.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -60,6 +60,8 @@ WarpX::ImplicitPreRHSOp ( amrex::Real a_cur_time,
using namespace amrex::literals;
amrex::ignore_unused( a_full_dt, a_nl_iter, a_from_jacobian );

if (use_filter) { ApplyFilterMF(Efield_fp, 0); }

// Advance the particle positions by 1/2 dt,
// particle velocities by dt, then take average of old and new v,
// deposit currents, giving J at n+1/2
Expand Down
34 changes: 17 additions & 17 deletions Source/Parallelization/WarpXComm.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -1030,7 +1030,7 @@ WarpX::SyncCurrent (
// Now it's safe to apply filter and sumboundary on J_cp
if (use_filter)
{
ApplyFilterJ(J_cp, lev+1, idim);
ApplyFilterMF(J_cp, lev+1, idim);
}
SumBoundaryJ(J_cp, lev+1, idim, period);
}
Expand Down Expand Up @@ -1064,7 +1064,7 @@ WarpX::SyncCurrent (

if (use_filter)
{
ApplyFilterJ(J_fp, lev, idim);
ApplyFilterMF(J_fp, lev, idim);
}
SumBoundaryJ(J_fp, lev, idim, period);
}
Expand Down Expand Up @@ -1176,30 +1176,30 @@ void WarpX::RestrictCurrentFromFineToCoarsePatch (
ablastr::coarsen::average::Coarsen(*crse[2], *fine[2], refinement_ratio );
}

void WarpX::ApplyFilterJ (
const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>,3>>& current,
void WarpX::ApplyFilterMF (
const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>,3>>& mfvec,
const int lev,
const int idim)
{
amrex::MultiFab& J = *current[lev][idim];
amrex::MultiFab& mf = *mfvec[lev][idim];

const int ncomp = J.nComp();
const amrex::IntVect ngrow = J.nGrowVect();
amrex::MultiFab Jf(J.boxArray(), J.DistributionMap(), ncomp, ngrow);
bilinear_filter.ApplyStencil(Jf, J, lev);
const int ncomp = mf.nComp();
const amrex::IntVect ngrow = mf.nGrowVect();
amrex::MultiFab mf_filtered(mf.boxArray(), mf.DistributionMap(), ncomp, ngrow);
bilinear_filter.ApplyStencil(mf_filtered, mf, lev);

const int srccomp = 0;
const int dstcomp = 0;
amrex::MultiFab::Copy(J, Jf, srccomp, dstcomp, ncomp, ngrow);
amrex::MultiFab::Copy(mf, mf_filtered, srccomp, dstcomp, ncomp, ngrow);
}

void WarpX::ApplyFilterJ (
const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>,3>>& current,
void WarpX::ApplyFilterMF (
const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>,3>>& mfvec,
const int lev)
{
for (int idim=0; idim<3; ++idim)
{
ApplyFilterJ(current, lev, idim);
ApplyFilterMF(mfvec, lev, idim);
}
}

Expand Down Expand Up @@ -1275,7 +1275,7 @@ void WarpX::AddCurrentFromFineLevelandSumBoundary (

if (use_filter)
{
ApplyFilterJ(J_fp, lev);
ApplyFilterMF(J_fp, lev);
}
SumBoundaryJ(J_fp, lev, period);

Expand All @@ -1294,8 +1294,8 @@ void WarpX::AddCurrentFromFineLevelandSumBoundary (

if (use_filter && J_buffer[lev+1][idim])
{
ApplyFilterJ(J_cp, lev+1, idim);
ApplyFilterJ(J_buffer, lev+1, idim);
ApplyFilterMF(J_cp, lev+1, idim);
ApplyFilterMF(J_buffer, lev+1, idim);

MultiFab::Add(
*J_buffer[lev+1][idim], *J_cp[lev+1][idim],
Expand All @@ -1309,7 +1309,7 @@ void WarpX::AddCurrentFromFineLevelandSumBoundary (
}
else if (use_filter) // but no buffer
{
ApplyFilterJ(J_cp, lev+1, idim);
ApplyFilterMF(J_cp, lev+1, idim);

ablastr::utils::communication::ParallelAdd(
mf, *J_cp[lev+1][idim], 0, 0,
Expand Down
8 changes: 4 additions & 4 deletions Source/WarpX.H
Original file line number Diff line number Diff line change
Expand Up @@ -1312,12 +1312,12 @@ private:
int lev);
void StoreCurrent (int lev);
void RestoreCurrent (int lev);
void ApplyFilterJ (
const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>,3>>& current,
void ApplyFilterMF (
const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>,3>>& mfvec,
int lev,
int idim);
void ApplyFilterJ (
const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>,3>>& current,
void ApplyFilterMF (
const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>,3>>& mfvec,
int lev);
void SumBoundaryJ (
const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>,3>>& current,
Expand Down
Loading