Deterministic YAFluxRegister on GPU

[BenWibking]

On GPU, YAFluxRegister is nondeterministic due to the use of atomics in yafluxreg_fineadd:

> Src/Boundary/AMReX_YAFluxRegister_3D_K.H:70:

      case 0 :
      {
          for (int n = 0; n < nc; ++n) {
              for     (int k = lo.z; k <= hi.z; ++k) {
                  for (int j = lo.y; j <= hi.y; ++j) {
                      const int i = lo.x;
                      const int ii = (i+1)*rr.x;
                      T* AMREX_RESTRICT dp = &(d(i,j,k,n));
                      for (    int koff = 0; koff < rr.z; ++koff) {
                          const     int kk =      k*rr.z  + koff;
                          for (int joff = 0; joff < rr.y; ++joff) {
                              const int jj =      j*rr.y  + joff;
                              T tmp = -dtdx*f(ii,jj,kk,n);
                              HostDevice::Atomic::Add(dp, tmp);

>  Src/Boundary/AMReX_YAFluxRegister_3D_K.H:92:

      case 1 :
      {
          for (int n = 0; n < nc; ++n) {
              for     (int k = lo.z; k <= hi.z; ++k) {
                  for (int j = lo.y; j <= hi.y; ++j) {
                      const int i = lo.x;
                      const int ii = i*rr.x;
                      T* AMREX_RESTRICT dp = &(d(i,j,k,n));
                      for (int     koff = 0; koff < rr.z; ++koff) {
                          const     int kk =      k*rr.z  + koff;
                          for (int joff = 0; joff < rr.y; ++joff) {
                              const int jj =      j*rr.y  + joff;
                              T tmp = dtdx*f(ii,jj,kk,n);
                              HostDevice::Atomic::Add(dp, tmp);

This means that bitwise reproducibility is currently not possible for our code when using AMR.

[BenWibking]

It is also nondeterministic because it uses ParallelCopy(..., FabArrayBase::ADD), which uses atomics. That should be fixed by #4658.

[WeiqunZhang]

You could use amrex::FluxRegister instead. I believe that it is deterministic on GPUs. You need to have a MultiFab storing the fluxes. Then if it's a coarse level, you call FluxRegister::CrseInit. If it's a fine level, you call FluxRegister::FineAdd. Finally, you call FluxRegister::Reflux. See amrex/Tests/Amr/Advection_AmrLevel and amrex/Tests/Amr/Advection_AmrCore an examples.

https://github.com/AMReX-Codes/amrex/tree/development/Tests/Amr/Advection_AmrLevel
https://github.com/AMReX-Codes/amrex/tree/development/Tests/Amr/Advection_AmrCore

[BenWibking]

@WeiqunZhang Does this line cause an issue? It looks like ParallelCopy with ADD is used underneath:

amrex/Src/AmrCore/AMReX_FluxRegister.cpp

Line 334 in df20009

bndry[face].plusFrom(mf,0,0,destcomp,numcomp,geom.periodicity());

[BenWibking]

If we use amrex::FluxRegister, we need a way to save and restore the fine fluxes for timestep retries. Or we need to apply the reflux after each dt/2 substep on the fine level.

Edit: Actually, we can just avoid filling the flux register if the step fails, since we have a global MultiFab for the whole step's fluxes now. So no changes should be necessary.

[BenWibking]

@WeiqunZhang Does this line cause an issue? It looks like ParallelCopy with ADD is used underneath:

amrex/Src/AmrCore/AMReX_FluxRegister.cpp

Line 334 in df20009

bndry[face].plusFrom(mf,0,0,destcomp,numcomp,geom.periodicity());

Ah, nevermind. I see this is avoided if I call CrseInit instead of CrseAdd.

[BenWibking]

Ok, I've switched Quokka to use amrex::FluxRegister, and everything works reproducibly now. I'll close this issue.