Hydrodynamic profile is showing small periodic unphjysical kinks, can't figure out what is manifesting this?

[moisesae]

I have taken the <u_x> moment using the particle_fields.do_average feature in the input file to get the hydrodynamic flow profile that WarpX interpolates onto the grid which I can then analyze. I have been taking numerical gradients of this profile and have noticed that I consistently get periodic spikes showing up. I thought that if I did a time-average of the velocity profile (~1000 time-steps), I would be able to see that the spikes go away when I apply the gradient. However, they persist even with time-averaging and I am having a difficult time trying to figure out why. They seem to occur every ~1mm but there isn't a length-scale that would correspond to this such as the debye length, gyro-radius, etc, which are sufficiently resolved as well (time-scales too). Pre-gradient, the spikes look small but once you apply the gradient, they get significantly amplified. I can see this in all of the species being simulated, with and without the B-field.

I'm not sure if it is an issue of how the moments are getting interpolated onto the grid, but I would appreciate if anyone could point me to some sanity-checks I can try to pinpoint this down.

(Each curve is normalized by their maximum value so they can be plotted on the same scale and be dimensionless)

The following is an example where I used significant time-averaging but the periodic spikes persist.

[roelof-groenewald]

Hi @moisesae. If you look at the WarpX output you should see a printout saying how many "grids" were used and the size of each. Does it look like the spikes correspond to the edges of those grids?

[roelof-groenewald]

Thanks for checking this. This kind of issue sometimes come up when the guard cells for fields are not properly updated.

Could you please add the following line:

red_mf.FillBoundary(warpx.Geom(m_lev).periodicity());

at line 82 of the file WarpX/Source/Diagnostics/ComputeDiagFunctors/ParticleReductionFunctor.cpp, then rebuild WarpX and see if the problem persists?

So just to be clear, after adding the line, that part of the code should look like this:

...
                // Fix dimensions since parser assumes u = gamma * v / c
                const amrex::ParticleReal ux = p.rdata(PIdx::ux) / PhysConst::c;
                const amrex::ParticleReal uy = p.rdata(PIdx::uy) / PhysConst::c;
                const amrex::ParticleReal uz = p.rdata(PIdx::uz) / PhysConst::c;
                const bool filtered_out_flag = ((do_filter) && (filter_fn(xw, yw, zw, ux, uy, uz) == 0.0_prt));
                const amrex::Real value = (filtered_out_flag) ? (0._rt):(map_fn(xw, yw, zw, ux, uy, uz));
                amrex::Gpu::Atomic::AddNoRet(&out_array(ii, jj, kk, 0), (amrex::Real)(p.rdata(PIdx::w) * value));
            });
    red_mf.FillBoundary(warpx.Geom(m_lev).periodicity());
    if (m_do_average) {
        amrex::MultiFab ppc_mf(warpx.boxArray(m_lev), warpx.DistributionMap(m_lev), 1, ng);
        // Add the weight for each particle -- total number of particles of this species
        ParticleToMesh(pc, ppc_mf, m_lev,
                [=] AMREX_GPU_DEVICE (const WarpXParticleContainer::SuperParticleType& p,
                    amrex::Array4<amrex::Real> const& out_array,
                    amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& plo,
                    amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& dxi)
                {
                    // Get position in WarpX convention to use in parser. Will be different from
                    // p.pos() for 1D and 2D simulations.
                    amrex::ParticleReal xw = 0._rt, yw = 0._rt, zw = 0._rt;
                    get_particle_position(p, xw, yw, zw);
...

[moisesae]

new_input_file.txt

@roelof-groenewald Thank you for the help. I took a look at the source code and added this correction and recompiled. Unfortunately, I am still seeing the spikes. I do want to mention that I am using an older sheath branch version of warpX but I was able to find where this line fits in that .cpp file despite some minor differences. I am sharing my input file so you can see what I'm working with.

Would it mess up the calculation or make the simulation slower if the block size is set to be the domain size? I am going to try. your correction with the main branch and see how it turns out.

[roelof-groenewald]

Thanks for trying it! I'm sorry it didn't fix the problem though.

WarpX does not have separate particle & grid distribution, in other words, if you set the block size equal to the domain size you would only be able to run on 1 processor. But could you please check if you run with half the number of processors (and setting the block size to be twice as big) if you see the distance between the spikes double as well? This would be a good check to confirm that the problem really is with guard cell values.

[moisesae]

@roelof-groenewald

Yes, I tried two cases using 16 and 32 processors. You can clearly see the spacing between the spikes double in length for the 16 processor case.

[roelof-groenewald]

Thanks for confirming that, so we at least know it is a problem with the guard cells. Could you please try moving the new line we added to right before the calculated particle function is interpolated onto the diagnostic grid (at the end of the same file I asked you to modify before)?
The end of that file (ParticleReductionFunctor.cpp) should then be:

    red_mf.FillBoundary(warpx.Geom(m_lev).periodicity());
    // Coarsen and interpolate from ppc_mf to the output diagnostic MultiFab, mf_dst.
    ablastr::coarsen::sample::Coarsen(mf_dst, red_mf, dcomp, 0, nComp(), 0, m_crse_ratio);

The reason I'm asking this is that perhaps the issue was with the logic that averages the particle function value over the number of particles. By doing the guard cell exchange after all the diagnostic calculations are completed, we can rule out the above statement. If this doesn't work either I'll block off some time for detailed debugging of this issue.