From 3fb5c8e156a2432e197acc956bcbca60d184d761 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Sun, 23 Aug 2026 13:03:01 -0700 Subject: [PATCH] fixes issue 91 92 93 94 95 --- Source/CAMR.H | 15 +++++--- Source/CAMR_advance.cpp | 11 +++++- Source/Hydro/CAMR_construct_hydro_source.cpp | 39 ++++++++++++++------ Source/Hydro/Godunov/Godunov_2D_eb.cpp | 4 +- Source/Hydro/Hydro_umdrv_eb.cpp | 9 ++++- 5 files changed, 57 insertions(+), 21 deletions(-) diff --git a/Source/CAMR.H b/Source/CAMR.H index d549478..c510d4a 100644 --- a/Source/CAMR.H +++ b/Source/CAMR.H @@ -479,11 +479,16 @@ CAMR::numGrow() { int ng = 4; #ifdef AMREX_USE_EB - // IF MOL + FRD then 5 - // IF Godunov + FRD then 6 - // IF MOL + SRD then 6 - // IF Godunov + SRD then 7 - if (redistribution_type == "FluxRedist") + // IF MOL + FRD then 5 + // IF Godunov + FRD then 6 + // IF MOL + SRD then 6 + // IF Godunov + SRD then 7 + // IF MOL + none then 5 + // IF Godunov + none then 6 + // + // NoRedist takes the same ngrow_bx = 2 as FluxRedist in + // construct_hydro_source, so it needs the same number of ghost cells + if (redistribution_type == "FluxRedist" || redistribution_type == "NoRedist") { if (do_mol) { ng = 5; diff --git a/Source/CAMR_advance.cpp b/Source/CAMR_advance.cpp index 8f8125e..15f9431 100644 --- a/Source/CAMR_advance.cpp +++ b/Source/CAMR_advance.cpp @@ -123,7 +123,6 @@ CAMR::CAMR_advance (Real time, MultiFab::Saxpy(S_new, dt, *old_sources[src_list[n]], 0, 0, NVAR, 0); } - sources_for_hydro.setVal(0.0); // // Now build and add the hydro source term(s) to S_new // @@ -140,6 +139,16 @@ CAMR::CAMR_advance (Real time, MultiFab::LinComb(S_new, 0.5, Sborder, 0, 0.5, S_old, 0, 0, NVAR, 0); MultiFab::Saxpy (S_new, 0.5*dt, new_hydro_source, 0, 0, NVAR, 0); + // The LinComb above averaged S^{n+1,*} -- which already held + // dt * old_sources -- against S^n, which does not, so only half of the + // old-time sources survived. Put the other half back, so that once the + // new-source correction (which is only 0.5*(Src_new - Src_old)) is + // added below the state carries the time-centered + // 0.5 * dt * (Src_old + Src_new), just as the Godunov branch does. + for (int n = 0; n < src_list.size(); ++n) { + MultiFab::Saxpy(S_new, 0.5*dt, *old_sources[src_list[n]], 0, 0, NVAR, 0); + } + } else { construct_hydro_source(Sborder, hydro_source, time, dt); diff --git a/Source/Hydro/CAMR_construct_hydro_source.cpp b/Source/Hydro/CAMR_construct_hydro_source.cpp index 74cb35e..9ac05f6 100644 --- a/Source/Hydro/CAMR_construct_hydro_source.cpp +++ b/Source/Hydro/CAMR_construct_hydro_source.cpp @@ -33,6 +33,10 @@ CAMR::construct_hydro_source (const MultiFab& S, // at a coarse-fine boundary. (The MOL path never builds srcQ.) int ng = (do_mol) ? 0 : numGrow(); + // Note that we must zero this here rather than once per advance because + // the Saxpy below accumulates, and the MOL path calls this routine twice. + sources_for_hydro.setVal(0.0); + for (int n = 0; n < src_list.size(); ++n) { MultiFab::Saxpy(sources_for_hydro, 1.0, *old_sources[src_list[n]], 0, 0, NVAR, ng); } @@ -59,6 +63,16 @@ CAMR::construct_hydro_source (const MultiFab& S, #ifdef AMREX_USE_EB const auto& ebfact = dynamic_cast(Factory()); + + // The boundary conditions are loop-invariant, so we copy them to the device + // once here rather than once per box. Doing it inside the loop would also + // free this memory while the redistribution kernels that read it are still + // in flight, and hand the same block to the next iteration. + const amrex::StateDescriptor* desc = state[State_Type].descriptor(); + const auto& bcs = desc->getBCs(); + amrex::Gpu::DeviceVector bcs_d(desc->nComp()); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, bcs.begin(), bcs.end(), bcs_d.begin()); #endif #ifdef _OPENMP @@ -67,9 +81,6 @@ CAMR::construct_hydro_source (const MultiFab& S, { #ifdef AMREX_USE_EB int ncomp = src_to_fill.nComp(); - FArrayBox dm_as_fine(Box::TheUnitBox(),ncomp); - FArrayBox fab_drho_as_crse(Box::TheUnitBox(),ncomp); - IArrayBox fab_rrflag_as_crse(Box::TheUnitBox()); #endif amrex::MFItInfo tiling = amrex::TilingIfNotGPU() ? amrex::MFItInfo().EnableTiling(hydro_tile_size) : amrex::MFItInfo(); @@ -171,24 +182,30 @@ CAMR::construct_hydro_source (const MultiFab& S, int as_crse = (fr_as_crse != nullptr); int as_fine = (fr_as_fine != nullptr); + // These are only placeholders when this box is not itself a + // coarse-fine interface; the flux register owns the data otherwise. + FArrayBox fab_drho_as_crse(Box::TheUnitBox(),ncomp,The_Async_Arena()); + IArrayBox fab_rrflag_as_crse(Box::TheUnitBox(),1,The_Async_Arena()); + FArrayBox* p_drho_as_crse = (fr_as_crse) ? fr_as_crse->getCrseData(mfi) : &fab_drho_as_crse; const IArrayBox* p_rrflag_as_crse = (fr_as_crse) ? fr_as_crse->getCrseFlag(mfi) : &fab_rrflag_as_crse; + // dm_as_fine is written by the redistribution and then read by the + // asynchronous FineAdd below, so it must be a fresh allocation from + // the async arena on every box -- reusing one fab across boxes + // would zero it while the previous box's FineAdd is still reading. + Box bx_for_dm(Box::TheUnitBox()); if (fr_as_fine) { const Box dbox1 = geom.growPeriodicDomain(1); - Box bx_for_dm(amrex::grow(bx,1) & dbox1); - dm_as_fine.resize(bx_for_dm,ncomp); + bx_for_dm = amrex::grow(bx,1) & dbox1; + } + FArrayBox dm_as_fine(bx_for_dm,ncomp,The_Async_Arena()); + if (fr_as_fine) { dm_as_fine.setVal(0.0); } - const amrex::StateDescriptor* desc = state[State_Type].descriptor(); - const auto& bcs = desc->getBCs(); - amrex::Gpu::DeviceVector bcs_d(desc->nComp()); - amrex::Gpu::copy( - amrex::Gpu::hostToDevice, bcs.begin(), bcs.end(), bcs_d.begin()); - const auto& dxInv = geom.InvCellSizeArray(); // Return hyd_src - centered at half-time if using Godunov method diff --git a/Source/Hydro/Godunov/Godunov_2D_eb.cpp b/Source/Hydro/Godunov/Godunov_2D_eb.cpp index eee3ac6..b6b9d12 100644 --- a/Source/Hydro/Godunov/Godunov_2D_eb.cpp +++ b/Source/Hydro/Godunov/Godunov_2D_eb.cpp @@ -182,7 +182,7 @@ Godunov_umeth_eb ( // This box must be grown by one in transverse direction so we can // do tangential interpolation when taking divergence later - const Box& xfxbx = surroundingNodes( grow(bx_to_fill, 1, 1-cdir), cdir); + const Box& xfxbx = surroundingNodes( grow(bx_to_fill, 1-cdir, 1), cdir); // Final Riemann problem X ParallelFor(xfxbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept @@ -209,7 +209,7 @@ Godunov_umeth_eb ( // This box must be grown by one in transverse direction so we can // do tangential interpolation when taking divergence later - const Box& yfxbx = surroundingNodes( grow(bx_to_fill, 1, 1-cdir), cdir); + const Box& yfxbx = surroundingNodes( grow(bx_to_fill, 1-cdir, 1), cdir); // Final Riemann problem Y ParallelFor(yfxbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept diff --git a/Source/Hydro/Hydro_umdrv_eb.cpp b/Source/Hydro/Hydro_umdrv_eb.cpp index 95ee059..7412de2 100644 --- a/Source/Hydro/Hydro_umdrv_eb.cpp +++ b/Source/Hydro/Hydro_umdrv_eb.cpp @@ -79,7 +79,11 @@ hydro_umdrv_eb( const bool do_mol, Box const& bx, // **************************************************************** // Quantities for redistribution // **************************************************************** - FArrayBox divc,redistwgt; + // These are put on the async arena (like qec above) because the + // redistribution kernels that read them are still in flight when this + // function returns -- only the FluxRedist branch of ApplyMLRedistribution + // ends with a stream synchronize. + FArrayBox divc(amrex::The_Async_Arena()), redistwgt(amrex::The_Async_Arena()); if (l_redistribution_type == "StateRedist") { divc.resize(bxg_i,NVAR); // This will hold "dUdt" before redistribution @@ -103,7 +107,8 @@ hydro_umdrv_eb( const bool do_mol, Box const& bx, // **************************************************************** FArrayBox flux_tmp[AMREX_SPACEDIM]; for (int idim=0; idim < AMREX_SPACEDIM; ++idim) { - flux_tmp[idim].resize(amrex::surroundingNodes(bxg_ii,idim),NVAR); + flux_tmp[idim].resize(amrex::surroundingNodes(bxg_ii,idim),NVAR, + amrex::The_Async_Arena()); flux_tmp[idim].setVal(0.); }