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.

Sign up for GitHub

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

Remove particles that are initialized in the EB #4585

Merged
merged 1 commit into from
Jan 29, 2024
Merged

Remove particles that are initialized in the EB #4585

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
16 changes: 16 additions & 0 deletions Source/Particles/PhysicalParticleContainer.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,6 +39,10 @@
#include "Utils/WarpXAlgorithmSelection.H"
#include "Utils/WarpXConst.H"
#include "Utils/WarpXProfilerWrapper.H"
#ifdef AMREX_USE_EB
# include "EmbeddedBoundary/ParticleBoundaryProcess.H"
# include "EmbeddedBoundary/ParticleScraper.H"
#endif
#include "WarpX.H"

#include <ablastr/warn_manager/WarnManager.H>
Expand Down Expand Up @@ -1430,6 +1434,12 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int
}
}

// Remove particles that are inside the embedded boundaries
#ifdef AMREX_USE_EB
auto & distance_to_eb = WarpX::GetInstance().GetDistanceToEB();
scrapeParticles( *this, amrex::GetVecOfConstPtrs(distance_to_eb), ParticleBoundaryProcess::Absorb());
#endif

// The function that calls this is responsible for redistributing particles.
}

Expand Down Expand Up @@ -1923,6 +1933,12 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector,
}
}

// Remove particles that are inside the embedded boundaries
#ifdef AMREX_USE_EB
auto & distance_to_eb = WarpX::GetInstance().GetDistanceToEB();
scrapeParticles(tmp_pc, amrex::GetVecOfConstPtrs(distance_to_eb), ParticleBoundaryProcess::Absorb());
#endif

// Redistribute the new particles that were added to the temporary container.
// (This eliminates invalid particles, and makes sure that particles
// are in the right tile.)
Expand Down
11 changes: 10 additions & 1 deletion Source/Particles/WarpXParticleContainer.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -60,7 +60,10 @@
#include <AMReX_ParticleUtil.H>
#include <AMReX_Random.H>
#include <AMReX_Utility.H>

#ifdef AMREX_USE_EB
# include "EmbeddedBoundary/ParticleBoundaryProcess.H"
# include "EmbeddedBoundary/ParticleScraper.H"
#endif

#ifdef AMREX_USE_OMP
# include <omp.h>
Expand Down Expand Up @@ -291,6 +294,12 @@ WarpXParticleContainer::AddNParticles (int /*lev*/, long n,
);
}

// Remove particles that are inside the embedded boundaries
Copy link
Member

Choose a reason for hiding this comment

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

This will work but will result in extra work setting up the particles that will be then removed.

Looking at what scrapeParticles is doing, the check for whether a particle is in the EB is a small bit of code and could easily be added in the loop adding the particles so any particles in the EB are immediately rejected. Would this be better?

                int ieb, jeb, keb;
                amrex::Real Web[AMREX_SPACEDIM][2];
                ablastr::particles::compute_weights_nodal(ppos.x, ppos.y, ppos.z, plo, dxi, ieb, jeb, keb, Web);
                amrex::Real phi_value = ablastr::particles::interp_field_nodal(ieb, jeb, keb, Web, phi);
                if (phi_value < 0.0) {
                    p.id() = -1;
                    continue;
                }

Though I realize that it wouldn't save much since it would have to happen after the UpdatePosition anyway.

Copy link
Member

@roelof-groenewald roelof-groenewald Jan 8, 2024 •
edited
Loading

Choose a reason for hiding this comment

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

Your comment about removing particles after the UpdatePosition call makes me think there is actually a bug in the code presently (besides giving confusing diagnostics on the first step). If I initialize a thermal plasma everywhere in the domain (including over the EB), but we don't remove those immediately, some particles (inside the EB) with velocities directed across the boundary will make it into the valid domain on the first push. Those particles will therefore lead to an undesired preferential velocity (away from the EB) along the edge of the boundary. I guess all this is to say that particles should be removed before the initial push (and also before UpdatePosition, right?).

Copy link
Member

Choose a reason for hiding this comment

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

My comment about UpdatePosition only applies to the flux injection, where particles are loaded onto a plane and UpdatePosition is used to advance the particles a fraction of a time step off of the plane. But your comment would still apply, so it would be better to do the check before the position update, which means doing the check as I suggested. This would also mean that any particles that started from a valid location but then crossed into an EB on its initial particle step would be included in the lost particle diagnostic.

Copy link
Member Author

Choose a reason for hiding this comment

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

Good points ; we could create a small device function isInsideEB to do this.

#ifdef AMREX_USE_EB
auto & distance_to_eb = WarpX::GetInstance().GetDistanceToEB();
scrapeParticles( *this, amrex::GetVecOfConstPtrs(distance_to_eb), ParticleBoundaryProcess::Absorb());
#endif

Redistribute();
}

Expand Down
4 changes: 3 additions & 1 deletion Source/WarpX.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -132,7 +132,9 @@ public:
MacroscopicProperties& GetMacroscopicProperties () { return *m_macroscopic_properties; }
HybridPICModel& GetHybridPICModel () { return *m_hybrid_pic_model; }
MultiDiagnostics& GetMultiDiags () {return *multi_diags;}

#ifdef AMREX_USE_EB
amrex::Vector<std::unique_ptr<amrex::MultiFab> >& GetDistanceToEB () {return m_distance_to_eb;}
#endif
ParticleBoundaryBuffer& GetParticleBoundaryBuffer () { return *m_particle_boundary_buffer; }

static void shiftMF (amrex::MultiFab& mf, const amrex::Geometry& geom,
Expand Down
Loading