amrex-astro / emu Goto Github PK

Currently the code crashes with an arithmetic error if everything is not set to zero first (msw_test branch). Should figure out what needs to be initialized and explicitly initialize it in main.cpp.

Add a visualization example using yt

Andrew has suggested an optimization where instead of redistributing in the post_update step, we add an additional ghost zone and redistribute at the end of every timestep.

This requires that particles in that extra ghost zone deposit & interpolate.

WarpX's shape factors by default won't allow this, as they use integer casting to get the leftmost cell where a particle deposits. Because C++ integer casting rounds towards 0, this will say that for particles in grid cell -1, the leftmost cell is 0 instead of -1.

We can fix this by modifying the shape factor structs like this:

const T x = xmid + T(0.5);
const auto i = static_cast<int>(x);
const auto j = std::signbit(x) ? i-1 : i;

And while we're at it, remove the Esirkepov shape factors we're not using.

Right now the implementation of the "length" of the particle needs to be generalized to three flavors. Right now it only works for two. In the case where there are initially no heavy lepton neutrinos it doesn't matter (the evolution equations are the same), but should be fixed for the future.

For the Cori GPU test profiling, particle to mesh is 65%, and mesh to particle is 17%.

If we did a sort of tree reduction for the atomics, by reducing in warps first, then to global memory, maybe we could greatly reduce this disparity.

Save particles to plotfiles (optional)

Save particles to checkpoint files

We can report the L2 and max norm of the error in the trace of the particles' density matrices.

Currently, for each particle the flavor vector length L is calculated from N_ab, the matrix storing the number of neutrinos in each state. We also have the particle weight N_p which is just the number of neutrinos in each particle.

Error in Tr(rho) = 1 can significantly affect the number of neutrinos that transform if L ~ N_p, since the size of L controls the maximum amount of flavor transformation. We can define an error E = |L(rho) * N_p - L_0| where L_0 is the initial length stored in the particles and L(rho) * N_p is the flavor vector length calculated from rho, scaled by the particle weight. Then take the L2 and max norms of E by computing Ep^2 for each particle and depositing it onto the grid E, then use sum and max to get L2 and max norms on the grid.

We can then add an input parameter for the error reporting frequency, like report_error_every.

The idea is that this will tell us when our simulations would benefit from symbolically making the assumption that Tr(rho_ab) = 1 => Tr(N_ab) = N_p, eliminating one of the diagonal components algebraically, and evolving a traceless matrix N_ab_p for particles.

This will make the code more efficient (requiring fewer ghost zones) and enable us to use third order particle shapes and to design better time integrators.

At the end of each timestep, go through the particles and make sure the sum of the diagonals is 1

Requires creating a copy of the particles.

Add a job_info file like Castro's where we write out the git hashes for Emu and amrex we're currently using, along with the compiler details.

Adapt the timestep to the current error estimate, perhaps taking a min of the adaptive timestep and the timestep determined from the neutrino potential as in #57

We could use, e.g. RK45 or Richardson extrapolation.

We currently store number of neutrinos in a grid cell Nab, same for flux, because we divide by cell volume in the interpolation step.

We can move the cell volume division into the deposition step, so our grid quantities are more typical number density & flux density.

Quietly get wrong results but no crashing for 1D simulations with shape factor 3 (e.g. bipolar test). I don't know where the assert went...

Removed gitpod to minimize maintenance area since I am the only one maintaining for now. Would be great to update it and put it back.

To get the code running on GPUs we made the redistribute global. With local redistribute, it was asserting that the number of particles out of bounds must be zero. We should return to local redistribute for efficiency.

The timestep is reported wrongly when recovering from output. If a checkpoint is written at timestep 1000, when restarting from that checkpoint the output says "restarting from checkpoint 999".

The changes to FlavoredNeutrinoContainer.cpp_Renormalize_fill do not allow for compiling on gpus because of the use of std::string

Sometimes we use Abort() which just generates a backtrace.

It would be much easier to understand why the code stopped if we replace all uses of amrex::Abort() with amrex::Error("An error message.") -- it will behave the same way as Abort() but print a helpful message :)

When compiling in parallel, we need to be able to serially (re) generate the code first, then actually compile it all.

It would make accessing fields easier if we could read the particle fields from the Header and provide a dictionary access mode for the raw particle data.

the timestep uses the largest diagonal component. Rather, we should calculate the length of the Hamiltonian isospin vector for each grid cell. This could lead to significant speedup at late times when decoherence shrinks the grid isospin vector length.

Before releasing we should write a readme that describes what the code does, basic configuration and running, and have an open-source license. The header of each file should also contain a brief copyright.

Make timestep limiter based on the grid potentials

Compute runtime diagnostics every so many steps.

To start, can simply implement the sums & averages from reduce_data.py using AMReX reductions.

We'll have to generate code for the reduction kernels & reduction variables, ofc.

Save the reduced data into an HDF5 file with an entry per diagnostic timestep.

Also, when we restart from a checkpoint, if that HDF5 file is in the current directory, we should delete entries at timesteps later than the restart timestep before resuming the simulation.