From b0bd7caa81d97fc2c35676a932ea8ae7836d0e12 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Mon, 25 Sep 2023 18:42:12 +0200 Subject: [PATCH] Update coords and driver --- .github/workflows/ci.yml | 1 + src/hydro/diffusion/conduction.cpp | 53 ++++----- src/hydro/hydro_driver.cpp | 176 +++++++++-------------------- 3 files changed, 82 insertions(+), 148 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index e24e197b..05049ba7 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -53,6 +53,7 @@ jobs: build/tst/regression/outputs/cluster_hse/analytic_comparison.png build/tst/regression/outputs/cluster_tabular_cooling/convergence.png build/tst/regression/outputs/aniso_therm_cond_ring_conv/ring_convergence.png + build/tst/regression/outputs/aniso_therm_cond_gauss_conv/cond.png build/tst/regression/outputs/field_loop/field_loop.png build/tst/regression/outputs/riemann_hydro/shock_tube.png build/tst/regression/outputs/turbulence/parthenon.hst diff --git a/src/hydro/diffusion/conduction.cpp b/src/hydro/diffusion/conduction.cpp index 0103ed1b..d6e57d30 100644 --- a/src/hydro/diffusion/conduction.cpp +++ b/src/hydro/diffusion/conduction.cpp @@ -78,14 +78,14 @@ Real EstimateConductionTimestep(MeshData *md) { {1, 1, 1, ib.e + 1 - ib.s}), KOKKOS_LAMBDA(const int b, const int k, const int j, const int i, Real &min_dt) { const auto &coords = prim_pack.GetCoords(b); - min_dt = fmin(min_dt, SQR(coords.Dx(parthenon::X1DIR, k, j, i)) / - (thermal_diff_coeff + TINY_NUMBER)); + min_dt = fmin(min_dt, + SQR(coords.Dxc<1>(k, j, i)) / (thermal_diff_coeff + TINY_NUMBER)); if (ndim >= 2) { - min_dt = fmin(min_dt, SQR(coords.Dx(parthenon::X2DIR, k, j, i)) / + min_dt = fmin(min_dt, SQR(coords.Dxc<2>(k, j, i)) / (thermal_diff_coeff + TINY_NUMBER)); } if (ndim >= 3) { - min_dt = fmin(min_dt, SQR(coords.Dx(parthenon::X3DIR, k, j, i)) / + min_dt = fmin(min_dt, SQR(coords.Dxc<3>(k, j, i)) / (thermal_diff_coeff + TINY_NUMBER)); } }, @@ -106,20 +106,20 @@ Real EstimateConductionTimestep(MeshData *md) { const auto dTdx = 0.5 * (prim(IPR, k, j, i + 1) / prim(IDN, k, j, i + 1) - prim(IPR, k, j, i - 1) / prim(IDN, k, j, i - 1)) / - coords.dx1v(i); + coords.Dxc<1>(i); const auto dTdy = ndim >= 2 ? 0.5 * (prim(IPR, k, j + 1, i) / prim(IDN, k, j + 1, i) - prim(IPR, k, j - 1, i) / prim(IDN, k, j - 1, i)) / - coords.dx2v(j) + coords.Dxc<2>(j) : 0.0; const auto dTdz = ndim >= 3 ? 0.5 * (prim(IPR, k + 1, j, i) / prim(IDN, k + 1, j, i) - prim(IPR, k - 1, j, i) / prim(IDN, k - 1, j, i)) / - coords.dx3v(k) + coords.Dxc<3>(k) : 0.0; const auto gradTmag = sqrt(SQR(dTdx) + SQR(dTdy) + SQR(dTdz)); @@ -130,15 +130,12 @@ Real EstimateConductionTimestep(MeshData *md) { auto thermal_diff_coeff = thermal_diff.Get(p, rho); if (thermal_diff.GetType() == Conduction::isotropic) { - min_dt = fmin(min_dt, - SQR(coords.Dx(parthenon::X1DIR, k, j, i)) / thermal_diff_coeff); + min_dt = fmin(min_dt, SQR(coords.Dxc<1>(k, j, i)) / thermal_diff_coeff); if (ndim >= 2) { - min_dt = fmin(min_dt, SQR(coords.Dx(parthenon::X2DIR, k, j, i)) / - thermal_diff_coeff); + min_dt = fmin(min_dt, SQR(coords.Dxc<2>(k, j, i)) / thermal_diff_coeff); } if (ndim >= 3) { - min_dt = fmin(min_dt, SQR(coords.Dx(parthenon::X3DIR, k, j, i)) / - thermal_diff_coeff); + min_dt = fmin(min_dt, SQR(coords.Dxc<3>(k, j, i)) / thermal_diff_coeff); } return; } @@ -166,16 +163,16 @@ Real EstimateConductionTimestep(MeshData *md) { const auto costheta = fabs(Bx * dTdx + By * dTdy + Bz * dTdz) / (Bmag * gradTmag); - min_dt = fmin(min_dt, SQR(coords.Dx(parthenon::X1DIR, k, j, i)) / + min_dt = fmin(min_dt, SQR(coords.Dxc<1>(k, j, i)) / (thermal_diff_coeff * fabs(Bx) / Bmag * costheta + TINY_NUMBER)); if (ndim >= 2) { - min_dt = fmin(min_dt, SQR(coords.Dx(parthenon::X2DIR, k, j, i)) / + min_dt = fmin(min_dt, SQR(coords.Dxc<2>(k, j, i)) / (thermal_diff_coeff * fabs(By) / Bmag * costheta + TINY_NUMBER)); } if (ndim >= 3) { - min_dt = fmin(min_dt, SQR(coords.Dx(parthenon::X3DIR, k, j, i)) / + min_dt = fmin(min_dt, SQR(coords.Dxc<3>(k, j, i)) / (thermal_diff_coeff * fabs(Bz) / Bmag * costheta + TINY_NUMBER)); } @@ -217,7 +214,7 @@ void ThermalFluxIsoFixed(MeshData *md) { const auto &prim = prim_pack(b); const auto T_i = prim(IPR, k, j, i) / prim(IDN, k, j, i); const auto T_im1 = prim(IPR, k, j, i - 1) / prim(IDN, k, j, i - 1); - const auto dTdx = (T_i - T_im1) / coords.Dx(parthenon::X1DIR, k, j, i); + const auto dTdx = (T_i - T_im1) / coords.Dxc<1>(k, j, i); const auto denf = 0.5 * (prim(IDN, k, j, i) + prim(IDN, k, j, i - 1)); cons.flux(X1DIR, IEN, k, j, i) -= thermal_diff_coeff * denf * dTdx; }); @@ -236,7 +233,7 @@ void ThermalFluxIsoFixed(MeshData *md) { const auto T_j = prim(IPR, k, j, i) / prim(IDN, k, j, i); const auto T_jm1 = prim(IPR, k, j - 1, i) / prim(IDN, k, j - 1, i); - const auto dTdy = (T_j - T_jm1) / coords.Dx(parthenon::X2DIR, k, j, i); + const auto dTdy = (T_j - T_jm1) / coords.Dxc<2>(k, j, i); const auto denf = 0.5 * (prim(IDN, k, j, i) + prim(IDN, k, j - 1, i)); cons.flux(X2DIR, IEN, k, j, i) -= thermal_diff_coeff * denf * dTdy; }); @@ -255,7 +252,7 @@ void ThermalFluxIsoFixed(MeshData *md) { const auto T_k = prim(IPR, k, j, i) / prim(IDN, k, j, i); const auto T_km1 = prim(IPR, k - 1, j, i) / prim(IDN, k - 1, j, i); - const auto dTdz = (T_k - T_km1) / coords.Dx(parthenon::X3DIR, k, j, i); + const auto dTdz = (T_k - T_km1) / coords.Dxc<3>(k, j, i); const auto denf = 0.5 * (prim(IDN, k, j, i) + prim(IDN, k - 1, j, i)); cons.flux(X3DIR, IEN, k, j, i) -= thermal_diff_coeff * denf * dTdz; }); @@ -304,7 +301,7 @@ void ThermalFluxGeneral(MeshData *md) { prim(IPR, k, j , i - 1) / prim(IDN, k, j , i - 1), prim(IPR, k, j , i - 1) / prim(IDN, k, j , i - 1) - prim(IPR, k, j - 1, i - 1) / prim(IDN, k, j - 1, i - 1)) / - coords.Dx(parthenon::X2DIR, k, j, i); + coords.Dxc<2>( k, j, i); if (ndim >= 3) { /* Monotonized temperature difference dT/dz, 3D problem ONLY */ @@ -316,13 +313,13 @@ void ThermalFluxGeneral(MeshData *md) { prim(IPR, k , j, i - 1) / prim(IDN, k , j, i - 1), prim(IPR, k , j, i - 1) / prim(IDN, k , j, i - 1) - prim(IPR, k - 1, j, i - 1) / prim(IDN, k - 1, j, i - 1)) / - coords.Dx(parthenon::X3DIR, k, j, i); + coords.Dxc<3>( k, j, i); } // clang-format on const auto T_i = prim(IPR, k, j, i) / prim(IDN, k, j, i); const auto T_im1 = prim(IPR, k, j, i - 1) / prim(IDN, k, j, i - 1); - const auto dTdx = (T_i - T_im1) / coords.Dx(parthenon::X1DIR, k, j, i); + const auto dTdx = (T_i - T_im1) / coords.Dxc<1>(k, j, i); const auto denf = 0.5 * (prim(IDN, k, j, i) + prim(IDN, k, j, i - 1)); const auto thermal_diff_f = @@ -399,7 +396,7 @@ void ThermalFluxGeneral(MeshData *md) { prim(IPR, k, j - 1, i ) / prim(IDN, k, j - 1, i ), prim(IPR, k, j - 1, i ) / prim(IDN, k, j - 1, i ) - prim(IPR, k, j - 1, i - 1) / prim(IDN, k, j - 1, i - 1)) / - coords.Dx(parthenon::X1DIR, k, j, i); + coords.Dxc<1>(k, j, i); if (ndim >= 3) { /* Monotonized temperature difference dT/dz, 3D problem ONLY */ @@ -411,14 +408,14 @@ void ThermalFluxGeneral(MeshData *md) { prim(IPR, k , j - 1, i) / prim(IDN, k , j - 1, i), prim(IPR, k , j - 1, i) / prim(IDN, k , j - 1, i) - prim(IPR, k - 1, j - 1, i) / prim(IDN, k - 1, j - 1, i)) / - coords.Dx(parthenon::X3DIR, k, j, i); + coords.Dxc<3>(k, j, i); } // clang-format on const auto T_j = prim(IPR, k, j, i) / prim(IDN, k, j, i); const auto T_jm1 = prim(IPR, k, j - 1, i) / prim(IDN, k, j - 1, i); - const auto dTdy = (T_j - T_jm1) / coords.Dx(parthenon::X2DIR, k, j, i); + const auto dTdy = (T_j - T_jm1) / coords.Dxc<2>(k, j, i); const auto denf = 0.5 * (prim(IDN, k, j, i) + prim(IDN, k, j - 1, i)); const auto gradTmag = sqrt(SQR(dTdx) + SQR(dTdy) + SQR(dTdz)); @@ -489,7 +486,7 @@ void ThermalFluxGeneral(MeshData *md) { prim(IPR, k - 1, j, i ) / prim(IDN, k - 1, j, i ), prim(IPR, k - 1, j, i ) / prim(IDN, k - 1, j, i ) - prim(IPR, k - 1, j, i - 1) / prim(IDN, k - 1, j, i - 1)) / - coords.Dx(parthenon::X1DIR, k, j, i); + coords.Dxc<1>(k, j, i); /* Monotonized temperature difference dT/dy */ const auto dTdy = @@ -501,12 +498,12 @@ void ThermalFluxGeneral(MeshData *md) { prim(IPR, k - 1, j , i) / prim(IDN, k - 1, j , i), prim(IPR, k - 1, j , i) / prim(IDN, k - 1, j , i) - prim(IPR, k - 1, j - 1, i) / prim(IDN, k - 1, j - 1, i)) / - coords.Dx(parthenon::X2DIR, k, j, i); + coords.Dxc<2>(k, j, i); // clang-format on const auto T_k = prim(IPR, k, j, i) / prim(IDN, k, j, i); const auto T_km1 = prim(IPR, k - 1, j, i) / prim(IDN, k - 1, j, i); - const auto dTdz = (T_k - T_km1) / coords.Dx(parthenon::X3DIR, k, j, i); + const auto dTdz = (T_k - T_km1) / coords.Dxc<3>(k, j, i); const auto denf = 0.5 * (prim(IDN, k, j, i) + prim(IDN, k - 1, j, i)); const auto gradTmag = sqrt(SQR(dTdx) + SQR(dTdy) + SQR(dTdz)); diff --git a/src/hydro/hydro_driver.cpp b/src/hydro/hydro_driver.cpp index 60295eec..bbc16e9d 100644 --- a/src/hydro/hydro_driver.cpp +++ b/src/hydro/hydro_driver.cpp @@ -250,8 +250,6 @@ void AddSTSTasks(TaskCollection *ptask_coll, Mesh *pmesh, BlockList_t &blocks, auto &pmb = blocks[i]; auto &tl = region_init[i]; auto &base = pmb->meshblock_data.Get(); - auto start_recv = tl.AddTask(none, &MeshBlockData::StartReceiving, base.get(), - BoundaryCommSubset::all); // Add extra registers. No-op for existing variables so it's safe to call every // time. @@ -266,6 +264,11 @@ void AddSTSTasks(TaskCollection *ptask_coll, Mesh *pmesh, BlockList_t &blocks, for (int i = 0; i < num_partitions; i++) { auto &tl = region_calc_fluxes_step_init[i]; auto &base = pmesh->mesh_data.GetOrAdd("base", i); + const auto any = parthenon::BoundaryType::any; + auto start_bnd = tl.AddTask(none, parthenon::StartReceiveBoundBufs, base); + auto start_flxcor_recv = + tl.AddTask(none, parthenon::StartReceiveFluxCorrections, base); + // Reset flux arrays (not guaranteed to be zero) auto reset_fluxes = tl.AddTask(none, ResetFluxes, base.get()); @@ -274,27 +277,19 @@ void AddSTSTasks(TaskCollection *ptask_coll, Mesh *pmesh, BlockList_t &blocks, // (in every subsetp). auto hydro_diff_fluxes = tl.AddTask(reset_fluxes, CalcDiffFluxes, hydro_pkg.get(), base.get()); - } - TaskRegion ®ion_flux_correct_step_init = ptask_coll->AddRegion(blocks.size()); - for (int i = 0; i < blocks.size(); i++) { - auto &tl = region_flux_correct_step_init[i]; - auto &base = blocks[i]->meshblock_data.Get("base"); - auto send_flux = - tl.AddTask(none, &MeshBlockData::SendFluxCorrection, base.get()); - auto recv_flux = - tl.AddTask(none, &MeshBlockData::ReceiveFluxCorrection, base.get()); - } + auto send_flx = + tl.AddTask(hydro_diff_fluxes, parthenon::LoadAndSendFluxCorrections, base); + auto recv_flx = + tl.AddTask(start_flxcor_recv, parthenon::ReceiveFluxCorrections, base); + auto set_flx = + tl.AddTask(recv_flx | hydro_diff_fluxes, parthenon::SetFluxCorrections, base); - TaskRegion ®ion_rkl2_step_init = ptask_coll->AddRegion(num_partitions); - for (int i = 0; i < num_partitions; i++) { - auto &tl = region_rkl2_step_init[i]; auto &Y0 = pmesh->mesh_data.GetOrAdd("u1", i); auto &MY0 = pmesh->mesh_data.GetOrAdd("MY0", i); - auto &base = pmesh->mesh_data.GetOrAdd("base", i); auto &Yjm2 = pmesh->mesh_data.GetOrAdd("Yjm2", i); - auto init_MY0 = tl.AddTask(none, parthenon::Update::FluxDivergence>, + auto init_MY0 = tl.AddTask(set_flx, parthenon::Update::FluxDivergence>, base.get(), MY0.get()); // Initialize Y0 and Y1 and the recursion relation starting with j = 2 needs data from @@ -302,36 +297,18 @@ void AddSTSTasks(TaskCollection *ptask_coll, Mesh *pmesh, BlockList_t &blocks, auto rkl2_step_first = tl.AddTask(init_MY0, RKL2StepFirst, Y0.get(), base.get(), Yjm2.get(), MY0.get(), s_rkl, tau); - // Update ghost cells of Y1 (as MY1 is calculated for each Y_j) + // Update ghost cells of Y1 (as MY1 is calculated for each Y_j). + // Y1 stored in "base", see rkl2_step_first task. + // Update ghost cells (local and non local), prolongate and apply bound cond. + // TODO(someone) experiment with split (local/nonlocal) comms with respect to + // performance for various tests (static, amr, block sizes) and then decide on the + // best impl. Go with default call (split local/nonlocal) for now. // TODO(pgrete) optimize (in parthenon) to only send subset of updated vars - auto send = tl.AddTask(rkl2_step_first, - parthenon::cell_centered_bvars::SendBoundaryBuffers, base); - auto recv = - tl.AddTask(send, parthenon::cell_centered_bvars::ReceiveBoundaryBuffers, base); - auto fill_from_bufs = - tl.AddTask(recv, parthenon::cell_centered_bvars::SetBoundaries, base); - } + auto bounds_exchange = parthenon::AddBoundaryExchangeTasks( + rkl2_step_first | start_bnd, tl, base, pmesh->multilevel); - TaskRegion ®ion_clear_bnd = ptask_coll->AddRegion(blocks.size()); - for (int i = 0; i < blocks.size(); i++) { - auto &tl = region_clear_bnd[i]; - auto &base = blocks[i]->meshblock_data.Get(); - auto clear_comm_flags = tl.AddTask(none, &MeshBlockData::ClearBoundary, - base.get(), BoundaryCommSubset::all); - auto prolongBound = none; - if (pmesh->multilevel) { - prolongBound = tl.AddTask(none, parthenon::ProlongateBoundaries, base); - } - - // set physical boundaries - auto set_bc = tl.AddTask(prolongBound, parthenon::ApplyBoundaryConditions, base); - } - TaskRegion ®ion_cons_to_prim = ptask_coll->AddRegion(num_partitions); - for (int i = 0; i < num_partitions; i++) { - auto &tl = region_cons_to_prim[i]; - auto &base = pmesh->mesh_data.GetOrAdd("base", i); - auto fill_derived = - tl.AddTask(none, parthenon::Update::FillDerived>, base.get()); + tl.AddTask(bounds_exchange, parthenon::Update::FillDerived>, + base.get()); } // Compute coefficients. Meyer+2012 eq. (16) @@ -351,82 +328,52 @@ void AddSTSTasks(TaskCollection *ptask_coll, Mesh *pmesh, BlockList_t &blocks, mu_tilde_j = mu_j * w1; gamma_tilde_j = -(1.0 - b_jm1) * mu_tilde_j; // -a_jm1*mu_tilde_j - TaskRegion ®ion_init_other = ptask_coll->AddRegion(blocks.size()); - for (int i = 0; i < blocks.size(); i++) { - auto &pmb = blocks[i]; - auto &tl = region_init_other[i]; - auto &base = pmb->meshblock_data.Get(); - // Only need boundaries for base as it's the only "active" container exchanging - // data/fluxes with neighbors. All other containers are passive (i.e., data is only - // used but not exchanged). - auto start_recv = tl.AddTask(none, &MeshBlockData::StartReceiving, base.get(), - BoundaryCommSubset::all); - } - TaskRegion ®ion_calc_fluxes_step_other = ptask_coll->AddRegion(num_partitions); for (int i = 0; i < num_partitions; i++) { auto &tl = region_calc_fluxes_step_other[i]; auto &base = pmesh->mesh_data.GetOrAdd("base", i); + // Only need boundaries for base as it's the only "active" container exchanging + // data/fluxes with neighbors. All other containers are passive (i.e., data is only + // used but not exchanged). + const auto any = parthenon::BoundaryType::any; + auto start_bnd = tl.AddTask(none, parthenon::StartReceiveBoundBufs, base); + auto start_flxcor_recv = + tl.AddTask(none, parthenon::StartReceiveFluxCorrections, base); + // Reset flux arrays (not guaranteed to be zero) auto reset_fluxes = tl.AddTask(none, ResetFluxes, base.get()); // Calculate the diffusive fluxes for Yjm1 (here u1) auto hydro_diff_fluxes = tl.AddTask(reset_fluxes, CalcDiffFluxes, hydro_pkg.get(), base.get()); - } - TaskRegion ®ion_flux_correct_step_other = ptask_coll->AddRegion(blocks.size()); - for (int i = 0; i < blocks.size(); i++) { - auto &tl = region_flux_correct_step_other[i]; - auto &base = blocks[i]->meshblock_data.Get(); - auto send_flux = - tl.AddTask(none, &MeshBlockData::SendFluxCorrection, base.get()); - auto recv_flux = - tl.AddTask(none, &MeshBlockData::ReceiveFluxCorrection, base.get()); - } + auto send_flx = + tl.AddTask(hydro_diff_fluxes, parthenon::LoadAndSendFluxCorrections, base); + auto recv_flx = + tl.AddTask(start_flxcor_recv, parthenon::ReceiveFluxCorrections, base); + auto set_flx = + tl.AddTask(recv_flx | hydro_diff_fluxes, parthenon::SetFluxCorrections, base); - TaskRegion ®ion_rkl2_step_other = ptask_coll->AddRegion(num_partitions); - for (int i = 0; i < num_partitions; i++) { - auto &tl = region_rkl2_step_other[i]; auto &Y0 = pmesh->mesh_data.GetOrAdd("u1", i); auto &MY0 = pmesh->mesh_data.GetOrAdd("MY0", i); - auto &base = pmesh->mesh_data.GetOrAdd("base", i); auto &Yjm2 = pmesh->mesh_data.GetOrAdd("Yjm2", i); auto rkl2_step_other = - tl.AddTask(none, RKL2StepOther, Y0.get(), base.get(), Yjm2.get(), MY0.get(), + tl.AddTask(set_flx, RKL2StepOther, Y0.get(), base.get(), Yjm2.get(), MY0.get(), mu_j, nu_j, mu_tilde_j, gamma_tilde_j, tau); // update ghost cells of base (currently storing Yj) + // Update ghost cells (local and non local), prolongate and apply bound cond. + // TODO(someone) experiment with split (local/nonlocal) comms with respect to + // performance for various tests (static, amr, block sizes) and then decide on the + // best impl. Go with default call (split local/nonlocal) for now. // TODO(pgrete) optimize (in parthenon) to only send subset of updated vars - auto send = tl.AddTask(rkl2_step_other, - parthenon::cell_centered_bvars::SendBoundaryBuffers, base); - auto recv = - tl.AddTask(send, parthenon::cell_centered_bvars::ReceiveBoundaryBuffers, base); - auto fill_from_bufs = - tl.AddTask(recv, parthenon::cell_centered_bvars::SetBoundaries, base); - } - TaskRegion ®ion_clear_bnd_other = ptask_coll->AddRegion(blocks.size()); - for (int i = 0; i < blocks.size(); i++) { - auto &tl = region_clear_bnd_other[i]; - auto &base = blocks[i]->meshblock_data.Get(); - auto clear_comm_flags = tl.AddTask(none, &MeshBlockData::ClearBoundary, - base.get(), BoundaryCommSubset::all); - auto prolongBound = none; - if (pmesh->multilevel) { - prolongBound = tl.AddTask(none, parthenon::ProlongateBoundaries, base); - } - - // set physical boundaries - auto set_bc = tl.AddTask(prolongBound, parthenon::ApplyBoundaryConditions, base); - } - TaskRegion ®ion_cons_to_prim_other = ptask_coll->AddRegion(num_partitions); - for (int i = 0; i < num_partitions; i++) { - auto &tl = region_cons_to_prim_other[i]; - auto &base = pmesh->mesh_data.GetOrAdd("base", i); - auto fill_derived = - tl.AddTask(none, parthenon::Update::FillDerived>, base.get()); + auto bounds_exchange = parthenon::AddBoundaryExchangeTasks( + rkl2_step_other | start_bnd, tl, base, pmesh->multilevel); + + tl.AddTask(bounds_exchange, parthenon::Update::FillDerived>, + base.get()); } b_jm2 = b_jm1; @@ -639,7 +586,11 @@ TaskCollection HydroDriver::MakeTaskCollection(BlockList_t &blocks, int stage) { auto &tl = single_tasklist_per_pack_region[i]; auto &mu0 = pmesh->mesh_data.GetOrAdd("base", i); auto &mu1 = pmesh->mesh_data.GetOrAdd("u1", i); - tl.AddTask(none, parthenon::StartReceiveFluxCorrections, mu0); + + const auto any = parthenon::BoundaryType::any; + auto start_bnd = tl.AddTask(none, parthenon::StartReceiveBoundBufs, mu0); + auto start_flxcor_recv = + tl.AddTask(none, parthenon::StartReceiveFluxCorrections, mu0); const auto flux_str = (stage == 1) ? "flux_first_stage" : "flux_other_stage"; FluxFun_t *calc_flux_fun = hydro_pkg->Param(flux_str); @@ -660,9 +611,9 @@ TaskCollection HydroDriver::MakeTaskCollection(BlockList_t &blocks, int stage) { auto send_flx = tl.AddTask(first_order_flux_correct, parthenon::LoadAndSendFluxCorrections, mu0); - auto recv_flx = - tl.AddTask(first_order_flux_correct, parthenon::ReceiveFluxCorrections, mu0); - auto set_flx = tl.AddTask(recv_flx, parthenon::SetFluxCorrections, mu0); + auto recv_flx = tl.AddTask(start_flxcor_recv, parthenon::ReceiveFluxCorrections, mu0); + auto set_flx = tl.AddTask(recv_flx | first_order_flux_correct, + parthenon::SetFluxCorrections, mu0); // compute the divergence of fluxes of conserved variables auto update = tl.AddTask( @@ -694,29 +645,14 @@ TaskCollection HydroDriver::MakeTaskCollection(BlockList_t &blocks, int stage) { tl.AddTask(source_split_strang_final, AddSplitSourcesFirstOrder, mu0.get(), tm); } - // Update ghost cells (local and non local) + // Update ghost cells (local and non local), prolongate and apply bound cond. // TODO(someone) experiment with split (local/nonlocal) comms with respect to // performance for various tests (static, amr, block sizes) and then decide on the // best impl. Go with default call (split local/nonlocal) for now. - parthenon::AddBoundaryExchangeTasks(source_split_first_order, tl, mu0, + parthenon::AddBoundaryExchangeTasks(source_split_first_order | start_bnd, tl, mu0, pmesh->multilevel); } - TaskRegion &async_region_3 = tc.AddRegion(num_task_lists_executed_independently); - for (int i = 0; i < blocks.size(); i++) { - auto &tl = async_region_3[i]; - auto &u0 = blocks[i]->meshblock_data.Get("base"); - auto prolongBound = none; - // Currently taken care of by AddBoundaryExchangeTasks above. - // Needs to be reintroduced once we reintroduce split (local/nonlocal) communication. - // if (pmesh->multilevel) { - // prolongBound = tl.AddTask(none, parthenon::ProlongateBoundaries, u0); - //} - - // set physical boundaries - auto set_bc = tl.AddTask(prolongBound, parthenon::ApplyBoundaryConditions, u0); - } - // Single task in single (serial) region to reset global vars used in reductions in the // first stage. if (stage == integrator->nstages && hydro_pkg->Param("calc_c_h")) {