From 3c02cade2a09565bcf90af00246a87d0acc265e8 Mon Sep 17 00:00:00 2001 From: Candace Gilet Date: Tue, 2 Apr 2024 13:38:12 -0400 Subject: [PATCH] Add nonconservative tracer (#94) --- src/CMakeLists.txt | 1 + src/Make.package | 1 + .../incflo_compute_advection_term.cpp | 96 ++++++--- src/diffusion/DiffusionScalarOp.H | 2 +- src/diffusion/DiffusionScalarOp.cpp | 12 +- src/incflo.H | 11 +- src/incflo.cpp | 2 +- src/incflo_apply_corrector.cpp | 80 +------- src/incflo_apply_predictor.cpp | 78 +------ src/incflo_compute_forces.cpp | 8 +- src/incflo_explicit_update.cpp | 192 ++++++++++++++++++ src/setup/init.cpp | 2 +- 12 files changed, 290 insertions(+), 195 deletions(-) create mode 100644 src/incflo_explicit_update.cpp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 09533db0..a696aac4 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -10,6 +10,7 @@ target_sources(incflo incflo_compute_forces.cpp incflo_correct_small_cells.cpp incflo_redistribute.cpp + incflo_explicit_update.cpp incflo_tagging.cpp incflo_regrid.cpp main.cpp diff --git a/src/Make.package b/src/Make.package index 9038fd57..e98e6a5a 100644 --- a/src/Make.package +++ b/src/Make.package @@ -4,6 +4,7 @@ CEXE_sources += incflo_apply_predictor.cpp CEXE_sources += incflo_apply_corrector.cpp CEXE_sources += incflo_compute_dt.cpp CEXE_sources += incflo_compute_forces.cpp +CEXE_sources += incflo_explicit_update.cpp CEXE_sources += incflo_tagging.cpp CEXE_sources += incflo_regrid.cpp CEXE_sources += main.cpp diff --git a/src/convection/incflo_compute_advection_term.cpp b/src/convection/incflo_compute_advection_term.cpp index ef29eeed..16d2df1b 100644 --- a/src/convection/incflo_compute_advection_term.cpp +++ b/src/convection/incflo_compute_advection_term.cpp @@ -28,9 +28,19 @@ void incflo::init_advection () m_iconserv_density.resize(1, 1); m_iconserv_density_d.resize(1, 1); - // We update (rho * tracer), not tracer itself, hence we update conservatively + // Advect scalars conservatively? m_iconserv_tracer.resize(m_ntrac, 1); - m_iconserv_tracer_d.resize(m_ntrac, 1); + ParmParse pp("incflo"); + pp.queryarr("trac_is_conservative", m_iconserv_tracer, 0, m_ntrac ); + m_iconserv_tracer_d.resize(m_ntrac); + // copy +#ifdef AMREX_USE_GPU + Gpu::htod_memcpy +#else + std::memcpy +#endif + (m_iconserv_tracer_d.data(), m_iconserv_tracer.data(), sizeof(int)*m_ntrac); + } void @@ -80,6 +90,21 @@ incflo::compute_convective_term (Vector const& conv_u, Vector > fluxes(finest_level+1); Vector > faces(finest_level+1); + bool any_conserv_trac = false; + for (auto& i : m_iconserv_tracer){ + if ( i == 1 ) { + any_conserv_trac = true; + break; + } + } + bool any_convective_trac = false; + for (auto& i : m_iconserv_tracer){ + if ( i == 0 ) { + any_convective_trac = true; + break; + } + } + for (int lev = 0; lev <= finest_level; ++lev) { AMREX_D_TERM( face_x[lev].define(u_mac[lev]->boxArray(),dmap[lev],n_flux_comp,0,MFInfo(),Factory(lev));, @@ -94,7 +119,7 @@ incflo::compute_convective_term (Vector const& conv_u, if (m_advect_momentum) rhovel[lev].define(vel[lev]->boxArray(),dmap[lev],AMREX_SPACEDIM, vel[lev]->nGrow(),MFInfo(),Factory(lev)); - if (m_advect_tracer && m_ntrac > 0) + if (m_advect_tracer && m_ntrac > 0 && any_conserv_trac) rhotrac[lev].define(vel[lev]->boxArray(),dmap[lev],tracer[lev]->nComp(), tracer[lev]->nGrow(),MFInfo(),Factory(lev)); @@ -151,7 +176,8 @@ incflo::compute_convective_term (Vector const& conv_u, if (nghost_force() > 0) fillpatch_force(m_cur_time, vel_forces, nghost_force()); - // Note this is forcing for (rho s), not for s + // Note that for conservative tracers, this is forcing for (rho s) + // and for non-conservative, this is forcing for s if (m_advect_tracer) { compute_tra_forces(tra_forces, get_density_old_const()); @@ -409,7 +435,7 @@ incflo::compute_convective_term (Vector const& conv_u, } // ************************************************************************ - // (Rho*Tracer) + // Tracer // ************************************************************************ // Make a FAB holding (rho * tracer) that is the same size as the original tracer FAB if (m_advect_tracer && (m_ntrac>0)) { @@ -421,13 +447,22 @@ incflo::compute_convective_term (Vector const& conv_u, Array4 tra = tracer[lev]->const_array(mfi); Array4 rho = density[lev]->const_array(mfi); - Array4 ro_trac = rhotrac[lev].array(mfi); + Array4 trac_tmp; - amrex::ParallelFor(bxg, m_ntrac, - [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept - { - ro_trac(i,j,k,n) = rho(i,j,k) * tra(i,j,k,n); - }); + auto const* iconserv = get_tracer_iconserv_device_ptr(); + if ( any_conserv_trac ) { + trac_tmp = rhotrac[lev].array(mfi); + + amrex::ParallelFor(bxg, m_ntrac, + [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept + { + if ( iconserv[n] ){ + trac_tmp(i,j,k,n) = rho(i,j,k) * tra(i,j,k,n); + } else { + trac_tmp(i,j,k,n) = tra(i,j,k,n); + } + }); + } if (m_constant_density) face_comp = AMREX_SPACEDIM; @@ -437,7 +472,8 @@ incflo::compute_convective_term (Vector const& conv_u, is_velocity = false; Array4 const& tracbc_arr = tracBC_MF ? (*tracBC_MF).const_array(mfi) : Array4{}; - HydroUtils::ComputeFluxesOnBoxFromState( bx, ncomp, mfi, ro_trac, + HydroUtils::ComputeFluxesOnBoxFromState( bx, ncomp, mfi, + any_conserv_trac ? trac_tmp : tracer[lev]->const_array(mfi), AMREX_D_DECL(flux_x[lev].array(mfi,face_comp), flux_y[lev].array(mfi,face_comp), flux_z[lev].array(mfi,face_comp)), @@ -616,13 +652,10 @@ incflo::compute_convective_term (Vector const& conv_u, } // mfi } // not constant density - // Note: (rho*trac) is always updated conservatively -- we do not provide an option for - // updating (rho*trac) convectively if (m_advect_tracer && m_ntrac > 0) { int flux_comp = (m_constant_density) ? AMREX_SPACEDIM : AMREX_SPACEDIM+1; -//Was this OMP intentionally left off? #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif @@ -650,15 +683,31 @@ incflo::compute_convective_term (Vector const& conv_u, (flagfab.getType(bx) != FabType::regular) ? ebfact->getBndryNormal().const_array(mfi) : Array4{}); #else - auto const& update_arr = conv_t[lev]->array(mfi); - HydroUtils::ComputeDivergence(bx, update_arr, - AMREX_D_DECL(flux_x[lev].const_array(mfi,flux_comp), - flux_y[lev].const_array(mfi,flux_comp), - flux_z[lev].const_array(mfi,flux_comp)), - m_ntrac, geom[lev], mult, - fluxes_are_area_weighted); + auto const& update_arr = conv_t[lev]->array(mfi); + HydroUtils::ComputeDivergence(bx, update_arr, + AMREX_D_DECL(flux_x[lev].const_array(mfi,flux_comp), + flux_y[lev].const_array(mfi,flux_comp), + flux_z[lev].const_array(mfi,flux_comp)), + m_ntrac, geom[lev], mult, + fluxes_are_area_weighted); #endif + if ( any_convective_trac ) + { + // If convective, we define u dot grad trac = div (u trac) - trac div(u) + HydroUtils::ComputeConvectiveTerm(bx, m_ntrac, mfi, + tracer[lev]->array(mfi,0), + AMREX_D_DECL(face_x[lev].array(mfi), + face_y[lev].array(mfi), + face_z[lev].array(mfi)), + divu[lev].array(mfi), + update_arr, + get_tracer_iconserv_device_ptr(), +#ifdef AMREX_USE_EB + *ebfact, +#endif + m_advection_type); + } } // mfi } // advect tracer @@ -699,7 +748,8 @@ incflo::compute_convective_term (Vector const& conv_u, if (m_advect_tracer) { auto const& bc_tra = get_tracer_bcrec_device_ptr(); redistribute_term(mfi, *conv_t[lev], dtdt_tmp, - rhotrac[lev], bc_tra, lev); + any_conserv_trac ? rhotrac[lev] : *tracer[lev], + bc_tra, lev); } } // mfi #endif diff --git a/src/diffusion/DiffusionScalarOp.H b/src/diffusion/DiffusionScalarOp.H index 826462c2..a77ac786 100644 --- a/src/diffusion/DiffusionScalarOp.H +++ b/src/diffusion/DiffusionScalarOp.H @@ -28,7 +28,7 @@ public: void compute_laps (amrex::Vector const& laps, amrex::Vector const& a_scalar, - amrex::Vector const& density, + amrex::Vector const& /*density*/, amrex::Vector const& eta); void compute_divtau (amrex::Vector const& a_divtau, diff --git a/src/diffusion/DiffusionScalarOp.cpp b/src/diffusion/DiffusionScalarOp.cpp index f0517f72..76553317 100644 --- a/src/diffusion/DiffusionScalarOp.cpp +++ b/src/diffusion/DiffusionScalarOp.cpp @@ -434,7 +434,7 @@ DiffusionScalarOp::diffuse_vel_components (Vector const& vel, void DiffusionScalarOp::compute_laps (Vector const& a_laps, Vector const& a_scalar, - Vector const& a_density, + Vector const& /*a_density*/, Vector const& a_eta) { BL_PROFILE("DiffusionScalarOp::compute_laps"); @@ -524,11 +524,6 @@ void DiffusionScalarOp::compute_laps (Vector const& a_laps, // We want to return div (mu grad)) phi m_reg_scal_apply_op->setScalars(0.0, -1.0); - // This should have no effect since the first scalar is 0 - for (int lev = 0; lev <= finest_level; ++lev) { - m_reg_scal_apply_op->setACoeffs(lev, *a_density[lev]); - } - for (int comp = 0; comp < m_incflo->m_ntrac; ++comp) { int eta_comp = comp; @@ -651,11 +646,6 @@ void DiffusionScalarOp::compute_divtau (Vector const& a_divtau, // We want to return div (mu grad)) phi m_reg_vel_apply_op->setScalars(0.0, -1.0); - // This should have no effect since the first scalar is 0 - for (int lev = 0; lev <= finest_level; ++lev) { - m_reg_vel_apply_op->setACoeffs(lev, *a_density[lev]); - } - int eta_comp = 0; Vector divtau_single; Vector vel_single; diff --git a/src/incflo.H b/src/incflo.H index ddd9934e..ac042e7c 100644 --- a/src/incflo.H +++ b/src/incflo.H @@ -81,8 +81,6 @@ public: // Delete level data void ClearLevel (int lev) override; -// for cuda - void ComputeDt (int initialization, bool explicit_diffusion); amrex::Real vol_wgt_sum (amrex::Vector const& mf, int icomp); @@ -115,10 +113,6 @@ public: std::string const& field); amrex::iMultiFab make_nodalBC_mask (int lev); amrex::Vector make_robinBC_MFs(int lev, amrex::MultiFab* state = nullptr); - // void make_ccBC_mask (int lev, const amrex::BoxArray& ba, - // const amrex::DistributionMapping& dm); - // void make_nodalBC_mask (int lev, const amrex::BoxArray& ba, - // const amrex::DistributionMapping& dm); #ifdef AMREX_USE_EB void set_eb_velocity (int lev, amrex::Real time, amrex::MultiFab& eb_vel, int nghost); @@ -156,6 +150,9 @@ public: amrex::Vector const& vel_forces, amrex::Real time); + void tracer_explicit_update(amrex::Vector const& tra_forces); + void tracer_explicit_update_corrector(amrex::Vector const& tra_forces); + /////////////////////////////////////////////////////////////////////////// // // derive @@ -199,6 +196,8 @@ public: [[nodiscard]] amrex::Array average_scalar_eta_to_faces (int lev, int comp, amrex::MultiFab const& cc_eta) const; + // FIXME? compute_laps never used density in any meaningful way, but + // then removing density would break anything downstream.... void compute_laps (amrex::Vector const& laps, amrex::Vector const& scalar, amrex::Vector const& density, diff --git a/src/incflo.cpp b/src/incflo.cpp index 9d7e72e6..15ec2f75 100644 --- a/src/incflo.cpp +++ b/src/incflo.cpp @@ -59,7 +59,7 @@ void incflo::InitData () InitialRedistribution(); #endif - if (m_do_initial_proj) { + if (m_do_initial_proj) { InitialProjection(); } diff --git a/src/incflo_apply_corrector.cpp b/src/incflo_apply_corrector.cpp index af22d130..4b148306 100644 --- a/src/incflo_apply_corrector.cpp +++ b/src/incflo_apply_corrector.cpp @@ -166,13 +166,11 @@ void incflo::ApplyCorrector() // Define local variables for lambda to capture. // ************************************************************************************* Real l_dt = m_dt; - bool l_constant_density = m_constant_density; - int l_ntrac = (m_advect_tracer) ? m_ntrac : 0; // ************************************************************************************* // Update density first // ************************************************************************************* - if (l_constant_density) + if (m_constant_density) { for (int lev = 0; lev <= finest_level; lev++) MultiFab::Copy(density_nph[lev], m_leveldata[lev]->density_o, 0, 0, 1, 0); @@ -229,80 +227,8 @@ void incflo::ApplyCorrector() // ************************************************************************************* if (m_advect_tracer) { - for (int lev = 0; lev <= finest_level; lev++) - { - auto& ld = *m_leveldata[lev]; -#ifdef _OPENMP -#pragma omp parallel if (Gpu::notInLaunchRegion()) -#endif - for (MFIter mfi(ld.tracer,TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - Box const& bx = mfi.tilebox(); - Array4 const& tra_o = ld.tracer_o.const_array(mfi); - Array4 const& rho_o = ld.density_o.const_array(mfi); - Array4 const& tra = ld.tracer.array(mfi); - Array4 const& rho = ld.density.const_array(mfi); - Array4 const& dtdt_o = ld.conv_tracer_o.const_array(mfi); - Array4 const& dtdt = ld.conv_tracer.const_array(mfi); - Array4 const& tra_f = (l_ntrac > 0) ? tra_forces[lev].const_array(mfi) - : Array4{}; - - if (m_diff_type == DiffusionType::Explicit) - { - Array4 const& laps_o = (l_ntrac > 0) ? ld.laps_o.const_array(mfi) - : Array4{}; - Array4 const& laps = (l_ntrac > 0) ? ld.laps.const_array(mfi) - : Array4{}; - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - for (int n = 0; n < l_ntrac; ++n) - { - Real tra_new = rho_o(i,j,k)*tra_o(i,j,k,n) + l_dt * ( - m_half*( dtdt(i,j,k,n) + dtdt_o(i,j,k,n)) - +m_half*(laps_o(i,j,k,n) + laps(i,j,k,n)) - + tra_f(i,j,k,n) ); - - tra_new /= rho(i,j,k); - tra(i,j,k,n) = tra_new; - } - }); - } - else if (m_diff_type == DiffusionType::Crank_Nicolson) - { - Array4 const& laps_o = (l_ntrac > 0) ? ld.laps_o.const_array(mfi) - : Array4{}; - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - for (int n = 0; n < l_ntrac; ++n) - { - Real tra_new = rho_o(i,j,k)*tra_o(i,j,k,n) + l_dt * ( - m_half*( dtdt(i,j,k,n) + dtdt_o(i,j,k,n)) - +m_half*(laps_o(i,j,k,n) ) - + tra_f(i,j,k,n) ); - - tra_new /= rho(i,j,k); - tra(i,j,k,n) = tra_new; - } - }); - } - else if (m_diff_type == DiffusionType::Implicit) - { - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - for (int n = 0; n < l_ntrac; ++n) - { - Real tra_new = rho_o(i,j,k)*tra_o(i,j,k,n) + l_dt * ( - m_half*( dtdt(i,j,k,n)+dtdt_o(i,j,k,n)) - + tra_f(i,j,k,n) ); - - tra_new /= rho(i,j,k); - tra(i,j,k,n) = tra_new; - } - }); - } - } // mfi - } // lev - } // if (m_advect_tracer) + tracer_explicit_update_corrector(tra_forces); + } // ************************************************************************************* // Solve diffusion equation for tracer diff --git a/src/incflo_apply_predictor.cpp b/src/incflo_apply_predictor.cpp index 0583e9a5..f353cadf 100644 --- a/src/incflo_apply_predictor.cpp +++ b/src/incflo_apply_predictor.cpp @@ -181,12 +181,11 @@ void incflo::ApplyPredictor (bool incremental_projection) // Define local variables for lambda to capture. // ************************************************************************************* Real l_dt = m_dt; - bool l_constant_density = m_constant_density; // ************************************************************************************* // Update density first // ************************************************************************************* - if (l_constant_density) + if (m_constant_density) { for (int lev = 0; lev <= finest_level; lev++) MultiFab::Copy(density_nph[lev], m_leveldata[lev]->density_o, 0, 0, 1, 1); @@ -244,7 +243,8 @@ void incflo::ApplyPredictor (bool incremental_projection) // ************************************************************************************* - // Compute (or if Godunov, re-compute) the tracer forcing terms (forcing for (rho s), not for s) + // Compute (or if Godunov, re-compute) the tracer forcing terms + // (forcing for (rho s) if conservative) // ************************************************************************************* if (m_advect_tracer) compute_tra_forces(GetVecOfPtrs(tra_forces), GetVecOfConstPtrs(density_nph)); @@ -252,78 +252,10 @@ void incflo::ApplyPredictor (bool incremental_projection) // ************************************************************************************* // Update the tracer next // ************************************************************************************* - int l_ntrac = (m_advect_tracer) ? m_ntrac : 0; - if (m_advect_tracer) { - for (int lev = 0; lev <= finest_level; lev++) - { - auto& ld = *m_leveldata[lev]; -#ifdef _OPENMP -#pragma omp parallel if (Gpu::notInLaunchRegion()) -#endif - for (MFIter mfi(ld.tracer,TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - Box const& bx = mfi.tilebox(); - Array4 const& tra_o = ld.tracer_o.const_array(mfi); - Array4 const& rho_o = ld.density_o.const_array(mfi); - Array4 const& tra = ld.tracer.array(mfi); - Array4 const& rho = ld.density.const_array(mfi); - Array4 const& dtdt_o = ld.conv_tracer_o.const_array(mfi); - Array4 const& tra_f = (l_ntrac > 0) ? tra_forces[lev].const_array(mfi) - : Array4{}; - - if (m_diff_type == DiffusionType::Explicit) - { - Array4 const& laps_o = (l_ntrac > 0) ? ld.laps_o.const_array(mfi) - : Array4{}; - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - // (rho trac)^new = (rho trac)^old + dt * ( - // div(rho trac u) + div (mu grad trac) + rho * f_t - for (int n = 0; n < l_ntrac; ++n) - { - Real tra_new = rho_o(i,j,k)*tra_o(i,j,k,n) + l_dt * - ( dtdt_o(i,j,k,n) + tra_f(i,j,k,n) + laps_o(i,j,k,n) ); - - tra_new /= rho(i,j,k); - tra(i,j,k,n) = tra_new; - } - }); - } - else if (m_diff_type == DiffusionType::Crank_Nicolson) - { - Array4 const& laps_o = (l_ntrac > 0) ? ld.laps_o.const_array(mfi) - : Array4{}; - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - for (int n = 0; n < l_ntrac; ++n) - { - Real tra_new = rho_o(i,j,k)*tra_o(i,j,k,n) + l_dt * - ( dtdt_o(i,j,k,n) + tra_f(i,j,k,n) + m_half * laps_o(i,j,k,n) ); - - tra_new /= rho(i,j,k); - tra(i,j,k,n) = tra_new; - } - }); - } - else if (m_diff_type == DiffusionType::Implicit) - { - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - for (int n = 0; n < l_ntrac; ++n) - { - Real tra_new = rho_o(i,j,k)*tra_o(i,j,k,n) + l_dt * - ( dtdt_o(i,j,k,n) + tra_f(i,j,k,n) ); - - tra_new /= rho(i,j,k); - tra(i,j,k,n) = tra_new; - } - }); - } - } // mfi - } // lev - } // if (m_advect_tracer) + tracer_explicit_update(tra_forces); + } // ************************************************************************************* // Solve diffusion equation for tracer diff --git a/src/incflo_compute_forces.cpp b/src/incflo_compute_forces.cpp index ff2a92de..4df2a15c 100644 --- a/src/incflo_compute_forces.cpp +++ b/src/incflo_compute_forces.cpp @@ -7,6 +7,8 @@ void incflo::compute_tra_forces (Vector const& tra_forces, { // NOTE: this routine must return the force term for the update of (rho s), NOT just s. if (m_advect_tracer) { + + auto const* iconserv = get_tracer_iconserv_device_ptr(); #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif @@ -23,8 +25,10 @@ void incflo::compute_tra_forces (Vector const& tra_forces, // For now we don't have any external forces on the scalars tra_f(i,j,k,n) = 0.0; - // Return the force term for the update of (rho s), NOT just s. - tra_f(i,j,k,n) *= rho(i,j,k); + if (iconserv[n]){ + // Return the force term for the update of (rho s), NOT just s. + tra_f(i,j,k,n) *= rho(i,j,k); + } }); } } diff --git a/src/incflo_explicit_update.cpp b/src/incflo_explicit_update.cpp new file mode 100644 index 00000000..ea6ca31e --- /dev/null +++ b/src/incflo_explicit_update.cpp @@ -0,0 +1,192 @@ +#include + +using namespace amrex; + +void incflo::tracer_explicit_update (Vector const& tra_forces) +{ + if (m_advect_tracer == 0) { return; } + + constexpr Real m_half = Real(0.5); + Real l_dt = m_dt; + int l_ntrac = m_ntrac; + for (int lev = 0; lev <= finest_level; lev++) + { + auto& ld = *m_leveldata[lev]; +#ifdef _OPENMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(ld.tracer,TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + Box const& bx = mfi.tilebox(); + Array4 const& tra_o = ld.tracer_o.const_array(mfi); + Array4 const& rho_o = ld.density_o.const_array(mfi); + Array4 const& tra = ld.tracer.array(mfi); + Array4 const& rho = ld.density.const_array(mfi); + Array4 const& dtdt_o = ld.conv_tracer_o.const_array(mfi); + Array4 const& tra_f = tra_forces[lev].const_array(mfi); + + auto const* iconserv = get_tracer_iconserv_device_ptr(); + + if (m_diff_type == DiffusionType::Explicit) + { + Array4 const& laps_o = ld.laps_o.const_array(mfi); + + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + // If conservative + // (rho trac)^new = (rho trac)^old + dt * ( + // div(rho trac u) + div (mu grad trac) + rho * f_t ) + // else non-conservative + // (trac)^new = (trac)^old + dt * ( + // u dot grad trac + div (mu grad trac) + f_t ) + for (int n = 0; n < l_ntrac; ++n) + { + if ( iconserv[n] ) { + Real tra_new = rho_o(i,j,k)*tra_o(i,j,k,n) + l_dt * + ( dtdt_o(i,j,k,n) + tra_f(i,j,k,n) + laps_o(i,j,k,n) ); + tra(i,j,k,n) = tra_new/rho(i,j,k); + } else { + tra(i,j,k,n) = tra_o(i,j,k,n) + l_dt * + ( dtdt_o(i,j,k,n) + tra_f(i,j,k,n) + laps_o(i,j,k,n) ); + } + } + }); + } + else if (m_diff_type == DiffusionType::Crank_Nicolson) + { + Array4 const& laps_o = ld.laps_o.const_array(mfi); + + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + for (int n = 0; n < l_ntrac; ++n) + { + if ( iconserv[n] ) { + Real tra_new = rho_o(i,j,k)*tra_o(i,j,k,n) + l_dt * + ( dtdt_o(i,j,k,n) + tra_f(i,j,k,n) + m_half * laps_o(i,j,k,n) ); + tra(i,j,k,n) = tra_new/rho(i,j,k); + } else { + tra(i,j,k,n) = tra_o(i,j,k,n) + l_dt * + ( dtdt_o(i,j,k,n) + tra_f(i,j,k,n) + m_half * laps_o(i,j,k,n) ); + } + } + }); + } + else if (m_diff_type == DiffusionType::Implicit) + { + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + for (int n = 0; n < l_ntrac; ++n) + { + if ( iconserv[n] ) { + Real tra_new = rho_o(i,j,k)*tra_o(i,j,k,n) + l_dt * + ( dtdt_o(i,j,k,n) + tra_f(i,j,k,n) ); + tra(i,j,k,n) = tra_new/rho(i,j,k); + } else { + tra(i,j,k,n) = tra_o(i,j,k,n) + l_dt * + ( dtdt_o(i,j,k,n) + tra_f(i,j,k,n) ); + } + + } + }); + } + } // mfi + } // lev +} + +void incflo::tracer_explicit_update_corrector (Vector const& tra_forces) +{ + if (m_advect_tracer == 0) { return; } + + constexpr Real m_half = Real(0.5); + Real l_dt = m_dt; + int l_ntrac = m_ntrac; + for (int lev = 0; lev <= finest_level; lev++) + { + auto& ld = *m_leveldata[lev]; +#ifdef _OPENMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(ld.tracer,TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + Box const& bx = mfi.tilebox(); + Array4 const& tra_o = ld.tracer_o.const_array(mfi); + Array4 const& rho_o = ld.density_o.const_array(mfi); + Array4 const& tra = ld.tracer.array(mfi); + Array4 const& rho = ld.density.const_array(mfi); + Array4 const& dtdt_o = ld.conv_tracer_o.const_array(mfi); + Array4 const& dtdt = ld.conv_tracer.const_array(mfi); + Array4 const& tra_f = tra_forces[lev].const_array(mfi); + auto const* iconserv = get_tracer_iconserv_device_ptr(); + + if (m_diff_type == DiffusionType::Explicit) + { + Array4 const& laps_o = ld.laps_o.const_array(mfi); + Array4 const& laps = ld.laps.const_array(mfi); + + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + for (int n = 0; n < l_ntrac; ++n) + { + if ( iconserv[n] ) { + Real tra_new = rho_o(i,j,k)*tra_o(i,j,k,n) + l_dt * ( + m_half*( dtdt(i,j,k,n) + dtdt_o(i,j,k,n)) + +m_half*(laps_o(i,j,k,n) + laps(i,j,k,n)) + + tra_f(i,j,k,n) ); + + tra(i,j,k,n) = tra_new / rho(i,j,k); + } else { + tra(i,j,k,n) = tra_o(i,j,k,n) + l_dt * ( + m_half*( dtdt(i,j,k,n) + dtdt_o(i,j,k,n)) + +m_half*(laps_o(i,j,k,n) + laps(i,j,k,n)) + + tra_f(i,j,k,n) ); + } + } + }); + } + else if (m_diff_type == DiffusionType::Crank_Nicolson) + { + Array4 const& laps_o = ld.laps_o.const_array(mfi); + + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + for (int n = 0; n < l_ntrac; ++n) + { + if ( iconserv[n] ) { + Real tra_new = rho_o(i,j,k)*tra_o(i,j,k,n) + l_dt * ( + m_half*( dtdt(i,j,k,n) + dtdt_o(i,j,k,n)) + +m_half*(laps_o(i,j,k,n) ) + + tra_f(i,j,k,n) ); + + tra(i,j,k,n) = tra_new / rho(i,j,k); + } else { + tra(i,j,k,n) = tra_o(i,j,k,n) + l_dt * ( + m_half*( dtdt(i,j,k,n) + dtdt_o(i,j,k,n)) + +m_half*(laps_o(i,j,k,n) ) + + tra_f(i,j,k,n) ); + } + } + }); + } + else if (m_diff_type == DiffusionType::Implicit) + { + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + for (int n = 0; n < l_ntrac; ++n) + { + if ( iconserv[n] ) { + Real tra_new = rho_o(i,j,k)*tra_o(i,j,k,n) + l_dt * ( + m_half*( dtdt(i,j,k,n)+dtdt_o(i,j,k,n)) + + tra_f(i,j,k,n) ); + + tra(i,j,k,n) = tra_new / rho(i,j,k); + } else { + tra(i,j,k,n) = tra_o(i,j,k,n) + l_dt * ( + m_half*( dtdt(i,j,k,n)+dtdt_o(i,j,k,n)) + + tra_f(i,j,k,n) ); + } + } + }); + } + } // mfi + } // lev +} diff --git a/src/setup/init.cpp b/src/setup/init.cpp index d0761ca8..ceaffc11 100644 --- a/src/setup/init.cpp +++ b/src/setup/init.cpp @@ -28,7 +28,7 @@ void incflo::ReadParameters () } // end prefix amr { // Prefix incflo - ParmParse pp("incflo"); + ParmParse pp("incflo"); pp.query("verbose", m_verbose);