diff --git a/Docs/sphinx_documentation/source/FluidEquations.rst b/Docs/sphinx_documentation/source/FluidEquations.rst index d6e74a2d..4a317f69 100644 --- a/Docs/sphinx_documentation/source/FluidEquations.rst +++ b/Docs/sphinx_documentation/source/FluidEquations.rst @@ -12,12 +12,16 @@ Fluid Variables +-----------------------+--------------------------------------------------+ | :math:`\mu_s` | scalar diffusivity | +-----------------------+--------------------------------------------------+ + | :math:`\mu_T` | thermal conductivity | + +-----------------------+--------------------------------------------------+ | :math:`{\bf g}` | Gravitational acceleration | +-----------------------+--------------------------------------------------+ | :math:`{\bf H}_U` | :math:`= (H_x , H_y , H_z )`, External Forces | +-----------------------+--------------------------------------------------+ | :math:`H_s` | External sources | +-----------------------+--------------------------------------------------+ + | :math:`H_T` | External heat sources | + +-----------------------+--------------------------------------------------+ Fluid Equations =============== @@ -43,5 +47,9 @@ or, for non-conservative, .. math:: \frac{\partial s}{\partial t} + U \cdot \nabla s = \nabla \cdot \mu_s \nabla s + H_s -By default, :math:`H_s = 0` and :math:`{\bf H}_U = {\bf 0}`. +Optional temperature: + +.. math:: \rho c_p \left( \frac{\partial T}{\partial t} + U \cdot \nabla T \right) = \nabla \cdot \mu_T \nabla T + H_T + +By default, :math:`H_s = 0`, :math:`H_T = 0`, and :math:`{\bf H}_U = {\bf 0}`. If gravity is set during runtime, then :math:`{\bf H}_U` defaults to :math:`\rho {\bf g}`. diff --git a/Docs/sphinx_documentation/source/InputsAlgorithm.rst b/Docs/sphinx_documentation/source/InputsAlgorithm.rst index ae3ba821..ce7814f8 100644 --- a/Docs/sphinx_documentation/source/InputsAlgorithm.rst +++ b/Docs/sphinx_documentation/source/InputsAlgorithm.rst @@ -18,6 +18,8 @@ The following inputs must be preceded by "incflo." +----------------------+-----------------------------------------------------------------------+-------------+--------------+ | ntrac | number of tracers | int | 1 | +----------------------+-----------------------------------------------------------------------+-------------+--------------+ +| use_temperature | include temperature equation? | bool | false | ++----------------------+-----------------------------------------------------------------------+-------------+--------------+ | constant_density | Only evolve the continuity equation if false | bool | true | +----------------------+-----------------------------------------------------------------------+-------------+--------------+ | rho_0 | density (if constant) | Real | 1.0 | @@ -28,6 +30,8 @@ The following inputs must be preceded by "incflo." +----------------------+-----------------------------------------------------------------------+-------------+--------------+ | mu_s | scalar diffusivity | Real(s) | 0.0 | +----------------------+-----------------------------------------------------------------------+-------------+--------------+ +| mu_T | thermal conductivity | Real(s) | 0.0 | ++----------------------+-----------------------------------------------------------------------+-------------+--------------+ | use_cc_proj | Use cell-centered rather than nodal pressure; this changes the | bool | false | | | which approximate projection we use | | | +----------------------+-----------------------------------------------------------------------+-------------+--------------+ diff --git a/Docs/sphinx_documentation/source/InputsProblemDefinition.rst b/Docs/sphinx_documentation/source/InputsProblemDefinition.rst index 4434800e..0f61c62b 100644 --- a/Docs/sphinx_documentation/source/InputsProblemDefinition.rst +++ b/Docs/sphinx_documentation/source/InputsProblemDefinition.rst @@ -50,15 +50,18 @@ Setting basic boundary conditions can be specified by inputs preceded by "xlo", | | * 'po' or 'pressure_outflow' | | | | | * 'mi' or 'mass_inflow' | | | | | * 'nsw' or 'no_slip_wall' | | | +| | * 'sw' or 'slip_wall' | | | | | * 'mixed' | | | +--------------------+---------------------------------------------------------------------------+-------------+-----------+ | pressure | Sets boundary pressure for pressure inflows, outflows and mass inflows | Real | None | +--------------------+---------------------------------------------------------------------------+-------------+-----------+ -| velocity | Sets boundary velocity for mass inflows | Reals | (0, 0, 0) | +| velocity | Sets Dirichlet boundary condition for velocity | Reals | (0, 0, 0) | +--------------------+---------------------------------------------------------------------------+-------------+-----------+ -| density | Sets boundary density for mass inflows | Real | 1.0 | +| density | Sets Dirichlet boundary condition for density | Real | 1.0 | +--------------------+---------------------------------------------------------------------------+-------------+-----------+ -| tracer | Sets boundary tracer for mass inflows | Real | 0.0 | +| tracer | Sets Dirichlet boundary condition for tracer | Real | 0.0 | ++--------------------+---------------------------------------------------------------------------+-------------+-----------+ +| temperature | Sets Dirichlet boundary condition for temperature | Real | 1.0 | +--------------------+---------------------------------------------------------------------------+-------------+-----------+ The 'mixed' boundary type allows for inflow and outflow on the same domain face. diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index f6602817..b7d673ec 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -16,6 +16,8 @@ target_sources(incflo incflo_update_density.cpp incflo_update_tracer.cpp incflo_update_velocity.cpp + incflo_update_temperature.cpp + incflo_specific_heat.cpp incflo_utils.cpp main.cpp ) diff --git a/src/Make.package b/src/Make.package index 135785eb..86112b3f 100644 --- a/src/Make.package +++ b/src/Make.package @@ -10,6 +10,8 @@ CEXE_sources += incflo_tagging.cpp CEXE_sources += incflo_update_density.cpp CEXE_sources += incflo_update_tracer.cpp CEXE_sources += incflo_update_velocity.cpp +CEXE_sources += incflo_update_temperature.cpp +CEXE_sources += incflo_specific_heat.cpp CEXE_sources += incflo_utils.cpp CEXE_sources += main.cpp diff --git a/src/boundary_conditions/README.org b/src/boundary_conditions/README.org index 16147d81..6008a5b9 100644 --- a/src/boundary_conditions/README.org +++ b/src/boundary_conditions/README.org @@ -1,17 +1,19 @@ * fillpatch -| | pi | po | mi | nsw | sw | dir_dep | -|--------+----------+----------+----------+----------------------+---------------------------------------|----------------------- -| v_n | foextrap | foextrap | ext_dir | ext_dir (0) | ext_dir (0) | ext_dir if inflowing | -| | | | | | | foextrap otherwise | -| v_t | foextrap | foextrap | ext_dir | ext_dir (0) | hoextrap | ext_dir if inflowing | -| | | | | | | foextrap otherwise | -| rho | foextrap | foextrap | ext_dir | foextrap | if (advection_type == "BDS") foextrap | ext_dir if inflowing | -| | | | | | else hoextrap | foextrap otherwise | -| tracer | foextrap | foextrap | ext_dir | ext_dir if specified | ext_dir if specified, otherwise | ext_dir if inflowing | -| | | | | foextrap otherwise | if (advection_type == "BDS") foextrap | foextrap otherwise | -| | | | | | else hoextrap | | -| force | foextrap | foextrap | foextrap | foextrap | foextrap | foextrap | +| | pi | po | mi | nsw | sw | dir_dep | +|-------------+----------+----------+----------+------------------------+---------------------------------------|----------------------- +| v_n | foextrap | foextrap | ext_dir | ext_dir (0) | ext_dir (0) | ext_dir if inflowing | +| | | | | | | foextrap otherwise | +| v_t | foextrap | foextrap | ext_dir | ext_dir (0) | hoextrap | ext_dir if inflowing | +| | | | | | | foextrap otherwise | +| rho | foextrap | foextrap | ext_dir | foextrap | if (advection_type == "BDS") foextrap | ext_dir if inflowing | +| | | | | | else hoextrap | foextrap otherwise | +| tracer | foextrap | foextrap | ext_dir | ext_dir if specified | ext_dir if specified, otherwise | ext_dir if inflowing | +| | | | | foextrap otherwise | if (advection_type == "BDS") foextrap | foextrap otherwise | +| | | | | | else hoextrap | | +| temperature | foextrap | foextrap | ext_dir | ext_dir if specified | ext_dir if specified, otherwise | ext_dir if inflowing | +| | | | | reflect_even otherwise | reflect_even | foextrap otherwise | +| force | foextrap | foextrap | foextrap | foextrap | foextrap | foextrap | * projection @@ -32,5 +34,5 @@ |---------+---------+---------+-----------+------------------------+------------------------|------------ | v_n | Neumann | Neumann | Dirichlet | Dirichlet | Dirichlet | Dirichlet | | v_t | Neumann | Neumann | Dirichlet | Dirichlet | Neumann | Dirichlet | -| tracer | Neumann | Neumann | Dirichlet | Dirchelet if specified | Dirichlet if specified | Dirichlet | +| scalars | Neumann | Neumann | Dirichlet | Dirchelet if specified | Dirichlet if specified | Dirichlet | | | | | | Neumann otherwise | Neumann otherwise | | diff --git a/src/boundary_conditions/boundary_conditions.cpp b/src/boundary_conditions/boundary_conditions.cpp index 8cdbb262..f4afacb7 100644 --- a/src/boundary_conditions/boundary_conditions.cpp +++ b/src/boundary_conditions/boundary_conditions.cpp @@ -13,6 +13,7 @@ void incflo::init_bcs () m_bcrec_velocity.resize(AMREX_SPACEDIM); m_bcrec_density.resize(1); if (m_ntrac > 0) { m_bcrec_tracer.resize(m_ntrac); } + m_bcrec_temperature.resize(1); auto f = [this] (std::string const& bcid, Orientation ori) { @@ -21,6 +22,7 @@ void incflo::init_bcs () m_bc_velocity[ori][1] = 0.0;, m_bc_velocity[ori][2] = 0.0;); m_bc_tracer[ori].resize(m_ntrac,0.0); + m_bc_temperature[ori] = 1.0; ParmParse pp(bcid); std::string bc_type_in = "null"; @@ -41,6 +43,7 @@ void incflo::init_bcs () m_bcrec_velocity[2].set(ori, BCType::foextrap);); m_bcrec_density[0].set(ori, BCType::foextrap); for (auto& b : m_bcrec_tracer) { b.set(ori, BCType::foextrap); } + m_bcrec_temperature[0].set(ori, BCType::foextrap); } else if (bc_type == "pressure_outflow" || bc_type == "po") { @@ -56,6 +59,7 @@ void incflo::init_bcs () m_bcrec_velocity[2].set(ori, BCType::foextrap);); m_bcrec_density[0].set(ori, BCType::foextrap); for (auto& b : m_bcrec_tracer) { b.set(ori, BCType::foextrap); } + m_bcrec_temperature[0].set(ori, BCType::foextrap); } else if (bc_type == "mass_inflow" || bc_type == "mi") { @@ -72,6 +76,7 @@ void incflo::init_bcs () pp.query("density", m_bc_density[ori]); pp.queryarr("tracer", m_bc_tracer[ori], 0, m_ntrac); + pp.query("temperature", m_bc_temperature[ori]); // Set mathematical BCs AMREX_D_TERM(m_bcrec_velocity[0].set(ori, BCType::ext_dir);, @@ -79,6 +84,7 @@ void incflo::init_bcs () m_bcrec_velocity[2].set(ori, BCType::ext_dir);); m_bcrec_density[0].set(ori, BCType::ext_dir); for (auto& b : m_bcrec_tracer) { b.set(ori, BCType::ext_dir); } + m_bcrec_temperature[0].set(ori, BCType::ext_dir); } else if (bc_type == "direction_dependent" || bc_type == "dd" ) { @@ -97,12 +103,14 @@ void incflo::init_bcs () pp.query("density", m_bc_density[ori]); pp.queryarr("tracer", m_bc_tracer[ori], 0, m_ntrac); + pp.query("temperature", m_bc_temperature[ori]); AMREX_D_TERM(m_bcrec_velocity[0].set(ori, BCType::direction_dependent);, m_bcrec_velocity[1].set(ori, BCType::direction_dependent);, m_bcrec_velocity[2].set(ori, BCType::direction_dependent);); m_bcrec_density[0].set(ori, BCType::direction_dependent); for (auto& b : m_bcrec_tracer) { b.set(ori, BCType::direction_dependent); } + m_bcrec_temperature[0].set(ori, BCType::direction_dependent); } else if (bc_type == "no_slip_wall" || bc_type == "nsw") { @@ -126,6 +134,7 @@ void incflo::init_bcs () // We potentially read in values at no-slip walls in the event that the // tracer has Dirichlet bcs pp.queryarr("tracer", m_bc_tracer[ori], 0, m_ntrac); + pp.query("temperature", m_bc_temperature[ori]); // Set mathematical BCs AMREX_D_TERM(m_bcrec_velocity[0].set(ori, BCType::ext_dir);, @@ -137,6 +146,11 @@ void incflo::init_bcs () } else { for (auto& b : m_bcrec_tracer) { b.set(ori, BCType::foextrap); } } + if ( pp.contains("temperature") ) { + for (auto& b : m_bcrec_temperature) { b.set(ori, BCType::ext_dir); } + } else { + for (auto& b : m_bcrec_temperature) { b.set(ori, BCType::reflect_even); } + } } else if (bc_type == "slip_wall" || bc_type == "sw") { @@ -151,6 +165,7 @@ void incflo::init_bcs () // We potentially read in values at slip walls in the event that the // tracer has Dirichlet bcs pp.queryarr("tracer", m_bc_tracer[ori], 0, m_ntrac); + pp.query("temperature", m_bc_temperature[ori]); // Tangential directions have hoextrap AMREX_D_TERM(m_bcrec_velocity[0].set(ori, BCType::hoextrap);, @@ -175,6 +190,11 @@ void incflo::init_bcs () } } } + if ( pp.contains("temperature") ) { + for (auto& b : m_bcrec_temperature) { b.set(ori, BCType::ext_dir); } + } else { + for (auto& b : m_bcrec_temperature) { b.set(ori, BCType::reflect_even); } + } } else if (bc_type == "mixed" ) { @@ -208,6 +228,7 @@ void incflo::init_bcs () pp.query("density", m_bc_density[ori]); pp.queryarr("tracer", m_bc_tracer[ori], 0, m_ntrac); + pp.query("temperature", m_bc_temperature[ori]); // Set mathematical BCs. BC_mask will handle Dirichlet part. AMREX_D_TERM(m_bcrec_velocity[0].set(ori, BCType::foextrap);, @@ -215,6 +236,7 @@ void incflo::init_bcs () m_bcrec_velocity[2].set(ori, BCType::foextrap);); m_bcrec_density[0].set(ori, BCType::foextrap); for (auto& b : m_bcrec_tracer) { b.set(ori, BCType::foextrap); } + m_bcrec_temperature[0].set(ori, BCType::foextrap); } else { @@ -231,6 +253,7 @@ void incflo::init_bcs () m_bcrec_velocity[2].set(ori, BCType::int_dir);); m_bcrec_density[0].set(ori, BCType::int_dir); for (auto& b : m_bcrec_tracer) { b.set(ori, BCType::int_dir); } + m_bcrec_temperature[0].set(ori, BCType::int_dir); } else { amrex::Abort("Wrong BC type for periodic boundary"); } @@ -298,6 +321,16 @@ void incflo::init_bcs () (m_bcrec_tracer_d.data(), m_bcrec_tracer.data(), sizeof(BCRec)*m_ntrac); } + if (m_use_temperature) { + m_bcrec_temperature_d.resize(1); +#ifdef AMREX_USE_GPU + Gpu::htod_memcpy +#else + std::memcpy +#endif + (m_bcrec_temperature_d.data(), m_bcrec_temperature.data(), sizeof(BCRec)); + } + // force { const int ncomp = std::max(m_ntrac, AMREX_SPACEDIM); diff --git a/src/boundary_conditions/incflo_fillpatch.cpp b/src/boundary_conditions/incflo_fillpatch.cpp index a8a9e692..b3a02632 100644 --- a/src/boundary_conditions/incflo_fillpatch.cpp +++ b/src/boundary_conditions/incflo_fillpatch.cpp @@ -114,6 +114,43 @@ void incflo::fillpatch_tracer (int lev, Real time, MultiFab& tracer, int ng) } } +void incflo::fillpatch_temperature (int lev, Real time, MultiFab& temperature, int ng) +{ + if (!m_use_temperature) return; + + if (lev == 0) { + PhysBCFunct > physbc + (geom[lev], get_temperature_bcrec(), IncfloTempFill{m_probtype, m_bc_temperature, m_bc_velocity}); + FillPatchSingleLevel(temperature, IntVect(ng), time, + {&(m_leveldata[lev]->temperature_o), + &(m_leveldata[lev]->temperature)}, + {m_t_old[lev], m_t_new[lev]}, 0, 0, 1, geom[lev], + physbc, 0); + } else { + const auto& bcrec = get_temperature_bcrec(); + PhysBCFunct > cphysbc + (geom[lev-1], bcrec, IncfloTempFill{m_probtype, m_bc_temperature, m_bc_velocity}); + PhysBCFunct > fphysbc + (geom[lev], bcrec, IncfloTempFill{m_probtype, m_bc_temperature, m_bc_velocity}); +#ifdef AMREX_USE_EB + Interpolater* mapper = (EBFactory(0).isAllRegular()) ? + (Interpolater*)(&cell_cons_interp) : (Interpolater*)(&eb_cell_cons_interp); +#else + Interpolater* mapper = &cell_cons_interp; +#endif + FillPatchTwoLevels(temperature, IntVect(ng), time, + {&(m_leveldata[lev-1]->temperature_o), + &(m_leveldata[lev-1]->temperature)}, + {m_t_old[lev-1], m_t_new[lev-1]}, + {&(m_leveldata[lev]->temperature_o), + &(m_leveldata[lev]->temperature)}, + {m_t_old[lev], m_t_new[lev]}, + 0, 0, 1, geom[lev-1], geom[lev], + cphysbc, 0, fphysbc, 0, + refRatio(lev-1), mapper, bcrec, 0); + } +} + void incflo::fillpatch_gradp (int lev, Real time, MultiFab& gp, int ng) { if (lev == 0) { @@ -234,6 +271,28 @@ void incflo::fillcoarsepatch_tracer (int lev, Real time, MultiFab& tracer, int n refRatio(lev-1), mapper, bcrec, 0); } +void incflo::fillcoarsepatch_temperature (int lev, Real time, MultiFab& temperature, int ng) +{ + if (m_use_temperature) return; + + const auto& bcrec = get_temperature_bcrec(); + PhysBCFunct > cphysbc + (geom[lev-1], bcrec, IncfloTempFill{m_probtype, m_bc_temperature, m_bc_velocity}); + PhysBCFunct > fphysbc + (geom[lev], bcrec, IncfloTempFill{m_probtype, m_bc_temperature, m_bc_velocity}); +#ifdef AMREX_USE_EB + Interpolater* mapper = (EBFactory(0).isAllRegular()) ? + (Interpolater*)(&cell_cons_interp) : (Interpolater*)(&eb_cell_cons_interp); +#else + Interpolater* mapper = &cell_cons_interp; +#endif + amrex::InterpFromCoarseLevel(temperature, IntVect(ng), time, + m_leveldata[lev-1]->temperature, 0, 0, 1, + geom[lev-1], geom[lev], + cphysbc, 0, fphysbc, 0, + refRatio(lev-1), mapper, bcrec, 0); +} + void incflo::fillcoarsepatch_gradp (int lev, Real time, MultiFab& gp, int ng) { const auto& bcrec = get_force_bcrec(); diff --git a/src/boundary_conditions/incflo_fillphysbc.cpp b/src/boundary_conditions/incflo_fillphysbc.cpp index f233b97f..d0fd830a 100644 --- a/src/boundary_conditions/incflo_fillphysbc.cpp +++ b/src/boundary_conditions/incflo_fillphysbc.cpp @@ -26,3 +26,14 @@ void incflo::fillphysbc_tracer (int lev, Real time, MultiFab& tracer, int ng) physbc.FillBoundary(tracer, 0, m_ntrac, IntVect(ng), time, 0); } } + +void incflo::fillphysbc_temperature (int lev, Real time, MultiFab& temperature, int ng) +{ + if (m_use_temperature) + { + PhysBCFunct > physbc + (geom[lev], get_temperature_bcrec(), + IncfloTempFill{m_probtype, m_bc_temperature, m_bc_velocity}); + physbc.FillBoundary(temperature, 0, 1, IntVect(ng), time, 0); + } +} diff --git a/src/boundary_conditions/incflo_set_bcs.cpp b/src/boundary_conditions/incflo_set_bcs.cpp index b9f16aa8..dde7a5c6 100644 --- a/src/boundary_conditions/incflo_set_bcs.cpp +++ b/src/boundary_conditions/incflo_set_bcs.cpp @@ -429,4 +429,73 @@ incflo::set_eb_tracer (int lev, Real /*time*/, MultiFab& eb_tracer, int nghost) IntVect ng_vect(AMREX_D_DECL(nghost,nghost,nghost)); eb_tracer.FillBoundary(0,m_ntrac,ng_vect,gm.periodicity()); } + +void +incflo::set_eb_temperature (int lev, Real /*time*/, MultiFab& eb_temperature, int nghost) +{ + Geometry const& gm = Geom(lev); + eb_temperature.setVal(0.); + + const auto& factory = + dynamic_cast(eb_temperature.Factory()); + +#ifdef _OPENMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(eb_temperature, TilingIfNotGPU()); mfi.isValid(); ++mfi) { + const Box& bx = mfi.tilebox(); + const auto& flagfab = factory.getMultiEBCellFlagFab()[mfi]; + + if (flagfab.getType(bx) == FabType::singlevalued) { + const auto& flags_arr = flagfab.const_array(); + const auto& eb_temperature_arr = eb_temperature[mfi].array(); + const auto& norm_arr = factory.getBndryNormal()[mfi].const_array(); + + bool has_normal = m_eb_flow.has_normal; + GpuArray normal{0.}; + if (has_normal) { + AMREX_D_TERM( + normal[0] = m_eb_flow.normal[0];, + normal[1] = m_eb_flow.normal[1];, + normal[2] = m_eb_flow.normal[2]); + } + Real pad = std::numeric_limits::epsilon(); + Real normal_tol = m_eb_flow.normal_tol; + Real norm_tol_lo = Real(-1.) - (normal_tol + pad); + Real norm_tol_hi = Real(-1.) + (normal_tol + pad); + + Real eb_flow_temperature = m_eb_flow.temperature[0]; + + ParallelFor(bx, [flags_arr,eb_temperature_arr,norm_arr,has_normal,normal, + norm_tol_lo, norm_tol_hi, eb_flow_temperature] + AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + if (flags_arr(i,j,k).isSingleValued()) { + Real mask = Real(1.0); + + if(has_normal) { +#if (AMREX_SPACEDIM==3) + Real dotprod = norm_arr(i,j,k,0)*normal[0] + + norm_arr(i,j,k,1)*normal[1] + + norm_arr(i,j,k,2)*normal[2]; +#else + Real dotprod = norm_arr(i,j,k,0)*normal[0] + + norm_arr(i,j,k,1)*normal[1]; +#endif + + mask = ((norm_tol_lo <= dotprod) && + (dotprod <= norm_tol_hi)) ? Real(1.0) : Real(0.0); + } + + eb_temperature_arr(i,j,k) = mask*eb_flow_temperature; + } + }); + } + } + + // We make sure to only fill "nghost" ghost cells so we don't accidentally + // over-write good ghost cell values with unfilled ghost cell values + IntVect ng_vect(AMREX_D_DECL(nghost,nghost,nghost)); + eb_temperature.FillBoundary(0,1,ng_vect,gm.periodicity()); +} #endif diff --git a/src/convection/incflo_compute_advection_term.cpp b/src/convection/incflo_compute_advection_term.cpp index 47614d7c..7c71a485 100644 --- a/src/convection/incflo_compute_advection_term.cpp +++ b/src/convection/incflo_compute_advection_term.cpp @@ -28,6 +28,10 @@ void incflo::init_advection () m_iconserv_density.resize(1, 1); m_iconserv_density_d.resize(1, 1); + // Temperature is always updated non-conservatively + m_iconserv_temperature.resize(1, 0); + m_iconserv_temperature_d.resize(1, 0); + // Advect scalars conservatively? m_iconserv_tracer.resize(m_ntrac, 1); ParmParse pp("incflo"); @@ -46,15 +50,18 @@ void incflo::init_advection () void incflo::compute_convective_term (Vector const& conv_u, Vector const& conv_r, - Vector const& conv_t, + Vector const& conv_tra, + Vector const& conv_tem, Vector const& vel, Vector const& density, Vector const& tracer, + Vector const& temperature, AMREX_D_DECL(Vector const& u_mac, Vector const& v_mac, Vector const& w_mac), Vector const& vel_forces, Vector const& tra_forces, + Vector const& tem_forces, Real /*time*/) { bool fluxes_are_area_weighted = false; @@ -68,6 +75,7 @@ incflo::compute_convective_term (Vector const& conv_u, int n_flux_comp = AMREX_SPACEDIM; if (!m_constant_density) n_flux_comp += 1; if ( m_advect_tracer) n_flux_comp += m_ntrac; + if ( m_use_temperature) n_flux_comp += 1; // This will hold state on faces Vector face_x(finest_level+1); @@ -191,6 +199,41 @@ incflo::compute_convective_term (Vector const& conv_u, fillpatch_force(m_cur_time, tra_forces, nghost_force()); } + if (m_use_temperature) + { + compute_tem_forces(m_cur_time, tem_forces); + 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(*density[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) { + Box const& bx = mfi.tilebox(); + FArrayBox cp_fab(bx, 1, The_Async_Arena()); + Array4 const& cp = cp_fab.array(); + Array4 const& rho = density[lev]->array(mfi); + Array4 const& tem_f = tem_forces[lev]->array(mfi); + + compute_cp(lev, mfi, cp_fab); + if (m_godunov_include_diff_in_forcing) { + Array4 const& laps = ld.laps_tem_o.const_array(mfi); + ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + tem_f(i,j,k) = (tem_f(i,j,k) + laps(i,j,k)) / (rho(i,j,k)*cp(i,j,k)); + }); + } else { + ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + tem_f(i,j,k) /= ( rho(i,j,k)*cp(i,j,k) ); + }); + + } + } + } + if (nghost_force() > 0) { + fillpatch_force(m_cur_time, tem_forces, nghost_force()); } + } + } // end m_advection_type for (int lev = 0; lev <= finest_level; ++lev) @@ -298,10 +341,17 @@ incflo::compute_convective_term (Vector const& conv_u, // ************************************************************************************* // Define domain boundary conditions at half-time to be used for fluxes if using Godunov // ************************************************************************************* - // + // For time-dependent Dirichlet BCs MultiFab vel_nph( vel[lev]->boxArray(), vel[lev]->DistributionMap(),AMREX_SPACEDIM,1); MultiFab rho_nph(density[lev]->boxArray(),density[lev]->DistributionMap(),1,1); - MultiFab trac_nph( tracer[lev]->boxArray(), tracer[lev]->DistributionMap(),m_ntrac,1); + MultiFab trac_nph; + if (m_advect_tracer && (m_ntrac>0)) { + trac_nph.define( tracer[lev]->boxArray(), tracer[lev]->DistributionMap(),m_ntrac,1); + } + MultiFab temp_nph; + if (m_use_temperature) { + temp_nph.define(temperature[lev]->boxArray(),temperature[lev]->DistributionMap(),1,1); + } if (m_advection_type != "MOL") { vel_nph.setVal(0.); @@ -330,6 +380,11 @@ incflo::compute_convective_term (Vector const& conv_u, } } } + + if (m_use_temperature) { + temp_nph.setVal(0.); + fillphysbc_temperature(lev, time_nph, temp_nph, 1); + } } // ************************************************************************ @@ -352,6 +407,12 @@ incflo::compute_convective_term (Vector const& conv_u, tracBC_MF = make_BC_MF(lev, m_bcrec_tracer_d, "tracer"); } } + // FIXME? this may just work as long as we add temp in make_BC_MF() + std::unique_ptr tempBC_MF; + if (m_use_temperature && m_has_mixedBC) { + Abort("Temperature equation with mixed BC not completed yet"); + tempBC_MF = make_BC_MF(lev, m_bcrec_temperature_d, "temperature"); + } // ************************************************************************ // Compute advective fluxes @@ -549,6 +610,48 @@ incflo::compute_convective_term (Vector const& conv_u, m_advection_type, PPM::default_limiter, allow_inflow_on_outflow, tracbc_arr); } + + // ************************************************************************ + // Temperature + // ************************************************************************ + if (m_use_temperature) { + // Temperature adveciton is always non-conservative + + face_comp = (m_advect_tracer && (m_ntrac>0)) ? m_ntrac : 0; + face_comp += (m_constant_density) ? AMREX_SPACEDIM : AMREX_SPACEDIM+1; + ncomp = 1; + is_velocity = false; + allow_inflow_on_outflow = false; + Array4 const& tempbc_arr = tempBC_MF ? (*tempBC_MF).const_array(mfi) + : Array4{}; + HydroUtils::ComputeFluxesOnBoxFromState( bx, ncomp, mfi, + temperature[lev]->const_array(mfi), + temp_nph.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)), + AMREX_D_DECL(face_x[lev].array(mfi,face_comp), + face_y[lev].array(mfi,face_comp), + face_z[lev].array(mfi,face_comp)), + knownFaceStates, + AMREX_D_DECL(u_mac[lev]->const_array(mfi), + v_mac[lev]->const_array(mfi), + w_mac[lev]->const_array(mfi)), + divu_arr, + (!tem_forces.empty()) ? tem_forces[lev]->const_array(mfi) : Array4{}, + geom[lev], m_dt, + get_temperature_bcrec(), + get_temperature_bcrec_device_ptr(), + m_iconserv_temperature_d.data(), +#ifdef AMREX_USE_EB + ebfact, + m_eb_flow.enabled ? get_temperature_eb()[lev]->const_array(mfi) : Array4{}, +#endif + m_godunov_ppm, m_godunov_use_forces_in_trans, + is_velocity, fluxes_are_area_weighted, + m_advection_type, PPM::default_limiter, + allow_inflow_on_outflow, tempbc_arr); + } } // mfi } // lev @@ -574,11 +677,13 @@ incflo::compute_convective_term (Vector const& conv_u, MultiFab dvdt_tmp(vel[lev]->boxArray(),dmap[lev],AMREX_SPACEDIM,3,MFInfo(),Factory(lev)); MultiFab drdt_tmp(vel[lev]->boxArray(),dmap[lev],1 ,3,MFInfo(),Factory(lev)); MultiFab dtdt_tmp(vel[lev]->boxArray(),dmap[lev],m_ntrac ,3,MFInfo(),Factory(lev)); + MultiFab dtemdt_tmp(vel[lev]->boxArray(),dmap[lev],1 ,3,MFInfo(),Factory(lev)); // Must initialize to zero because not all values may be set, e.g. outside the domain. dvdt_tmp.setVal(0.); drdt_tmp.setVal(0.); dtdt_tmp.setVal(0.); + dtemdt_tmp.setVal(0.); const EBFArrayBoxFactory* ebfact = &EBFactory(lev); auto const& vfrac = ebfact->getVolFrac(); @@ -666,7 +771,7 @@ incflo::compute_convective_term (Vector const& conv_u, if (!m_constant_density) { int flux_comp = AMREX_SPACEDIM; -//Was this OMP intentionally left off? + #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif @@ -710,7 +815,7 @@ incflo::compute_convective_term (Vector const& conv_u, #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif - for (MFIter mfi(*conv_t[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) + for (MFIter mfi(*conv_tra[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) { Box const& bx = mfi.tilebox(); @@ -734,7 +839,7 @@ 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); + auto const& update_arr = conv_tra[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), @@ -762,12 +867,70 @@ incflo::compute_convective_term (Vector const& conv_u, } // mfi } // advect tracer + if (m_use_temperature) + { + int flux_comp = (m_advect_tracer && (m_ntrac>0)) ? m_ntrac : 0; + flux_comp += (m_constant_density) ? AMREX_SPACEDIM : AMREX_SPACEDIM+1; + +#ifdef _OPENMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(*conv_tem[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + Box const& bx = mfi.tilebox(); + +#ifdef AMREX_USE_EB + EBCellFlagFab const& flagfab = ebfact->getMultiEBCellFlagFab()[mfi]; + auto const& update_arr = dtemdt_tmp.array(mfi); + if (flagfab.getType(bx) != FabType::covered) + HydroUtils::EB_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)), + vfrac.const_array(mfi), 1, geom[lev], mult, + fluxes_are_area_weighted, + m_eb_flow.enabled ? + get_velocity_eb()[lev]->const_array(mfi) : Array4{}, + m_eb_flow.enabled ? + get_temperature_eb()[lev]->const_array(mfi) : Array4{}, + flagfab.const_array(), + (flagfab.getType(bx) != FabType::regular) ? + ebfact->getBndryArea().const_array(mfi) : Array4{}, + (flagfab.getType(bx) != FabType::regular) ? + ebfact->getBndryNormal().const_array(mfi) : Array4{}); +#else + auto const& update_arr = conv_tem[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)), + 1, geom[lev], mult, + fluxes_are_area_weighted); +#endif + + // For convective, we define u dot grad trac = div (u trac) - trac div(u) + HydroUtils::ComputeConvectiveTerm(bx, 1, mfi, + temperature[lev]->array(mfi,0), + AMREX_D_DECL(face_x[lev].array(mfi,flux_comp), + face_y[lev].array(mfi,flux_comp), + face_z[lev].array(mfi,flux_comp)), + divu[lev].array(mfi), + update_arr, + m_iconserv_temperature_d.data(), +#ifdef AMREX_USE_EB + *ebfact, +#endif + m_advection_type); + } // mfi + } // advect temperature + #ifdef AMREX_USE_EB // We only filled these on the valid cells so we fill same-level interior ghost cells here. // (We don't need values outside the domain or at a coarser level so we can call just FillBoundary) dvdt_tmp.FillBoundary(geom[lev].periodicity()); drdt_tmp.FillBoundary(geom[lev].periodicity()); dtdt_tmp.FillBoundary(geom[lev].periodicity()); + dtemdt_tmp.FillBoundary(geom[lev].periodicity()); #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) @@ -798,10 +961,17 @@ 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, + redistribute_term(mfi, *conv_tra[lev], dtdt_tmp, any_conserv_trac ? rhotrac[lev] : *tracer[lev], bc_tra, lev); } + + if (m_use_temperature) { + auto const& bc_tem = get_tracer_bcrec_device_ptr(); + redistribute_term(mfi, *conv_tem[lev], dtemdt_tmp,//fixme + *temperature[lev], + bc_tem, lev); + } } // mfi #endif } // lev diff --git a/src/diffusion/DiffusionScalarOp.H b/src/diffusion/DiffusionScalarOp.H index df74be14..ef843825 100644 --- a/src/diffusion/DiffusionScalarOp.H +++ b/src/diffusion/DiffusionScalarOp.H @@ -15,9 +15,12 @@ class DiffusionScalarOp public: DiffusionScalarOp (incflo* a_incflo); - void diffuse_scalar (amrex::Vector const& tracer, - amrex::Vector const& density, + void diffuse_scalar (amrex::Vector const& a_scalar, + amrex::Vector const& chi, amrex::Vector const& eta, + amrex::Vector const& eb_dirichlet, + amrex::Vector const& use_chi, + amrex::Vector bcrec, amrex::Real dt); void diffuse_vel_components ( @@ -28,7 +31,8 @@ public: void compute_laps (amrex::Vector const& laps, amrex::Vector const& a_scalar, - amrex::Vector const& eta); + amrex::Vector const& eta, + amrex::Vector bcrec); void compute_divtau (amrex::Vector const& a_divtau, amrex::Vector const& a_vel, diff --git a/src/diffusion/DiffusionScalarOp.cpp b/src/diffusion/DiffusionScalarOp.cpp index e8a9b2d4..b277c1d0 100644 --- a/src/diffusion/DiffusionScalarOp.cpp +++ b/src/diffusion/DiffusionScalarOp.cpp @@ -32,11 +32,6 @@ DiffusionScalarOp::DiffusionScalarOp (incflo* a_incflo) m_incflo->DistributionMap(0,finest_level), info_solve, ebfact); m_eb_scal_solve_op->setMaxOrder(m_mg_maxorder); - // For now, code requires (in more than 1 place) that m_ntrac>=1 and all the tracers have the same BCs - m_eb_scal_solve_op->setDomainBC(m_incflo->get_diffuse_scalar_bc(Orientation::low, - m_incflo->m_bcrec_tracer[0].lo()), - m_incflo->get_diffuse_scalar_bc(Orientation::high, - m_incflo->m_bcrec_tracer[0].hi())); if (!m_incflo->useTensorSolve()) { @@ -56,11 +51,6 @@ DiffusionScalarOp::DiffusionScalarOp (incflo* a_incflo) m_incflo->DistributionMap(0,finest_level), info_apply, ebfact); m_eb_scal_apply_op->setMaxOrder(m_mg_maxorder); - - m_eb_scal_apply_op->setDomainBC(m_incflo->get_diffuse_scalar_bc(Orientation::low, - m_incflo->m_bcrec_tracer[0].lo()), - m_incflo->get_diffuse_scalar_bc(Orientation::high, - m_incflo->m_bcrec_tracer[0].hi())); } if ( (m_incflo->need_divtau() && !m_incflo->useTensorSolve()) || @@ -83,10 +73,6 @@ DiffusionScalarOp::DiffusionScalarOp (incflo* a_incflo) m_incflo->DistributionMap(0,m_incflo->finestLevel()), info_solve); m_reg_scal_solve_op->setMaxOrder(m_mg_maxorder); - m_reg_scal_solve_op->setDomainBC(m_incflo->get_diffuse_scalar_bc(Orientation::low, - m_incflo->m_bcrec_tracer[0].lo()), - m_incflo->get_diffuse_scalar_bc(Orientation::high, - m_incflo->m_bcrec_tracer[0].hi())); if (!m_incflo->useTensorSolve()) { @@ -104,10 +90,6 @@ DiffusionScalarOp::DiffusionScalarOp (incflo* a_incflo) m_incflo->DistributionMap(0,m_incflo->finestLevel()), info_apply); m_reg_scal_apply_op->setMaxOrder(m_mg_maxorder); - m_reg_scal_apply_op->setDomainBC(m_incflo->get_diffuse_scalar_bc(Orientation::low, - m_incflo->m_bcrec_tracer[0].lo()), - m_incflo->get_diffuse_scalar_bc(Orientation::high, - m_incflo->m_bcrec_tracer[0].hi())); } if ( (m_incflo->need_divtau() && !m_incflo->useTensorSolve()) || @@ -147,30 +129,38 @@ DiffusionScalarOp::readParameters () } void -DiffusionScalarOp::diffuse_scalar (Vector const& tracer, - Vector const& density, +DiffusionScalarOp::diffuse_scalar (Vector const& a_scalar, + Vector const& chi, Vector const& eta, + Vector const& eb_dirichlet, + amrex::Vector const& use_chi, + amrex::Vector bcrec, Real dt) { // // Solves - // alpha a - beta div ( b grad ) + // [alpha a - beta div ( b grad )] sca = RHS // - // So for conservative: d(rho tra) / dt - div mu grad tra = -div(U rho tra) + rho H --> - // ( rho - dt div mu grad ) tra^(n+1) = rho tra^(*,n+1) - // alpha: 1 - // a: rho - // beta: dt - // b: mu - // RHS: density * tracer + // If use_chi, solve + // ( chi - dt div eta grad ) sca^(n+1) = chi sca^(*,n+1) + // alpha: 1 + // a: chi + // beta: dt + // b: eta + // RHS: chi * a_scalar + // This corresponds to the conservative scalar equation with chi -> rho: + // d(rho sca) / dt - div mu grad sca = -div(U rho sca) + rho H + // And also to the (non-conservative) temperature equation with chi -> rhoCp // - // So for non- conservative: d tra / dt - div mu grad tra = -U dot grad tra + H --> - // ( 1 - dt div mu grad ) tra^(n+1) = tra^(*,n+1) - // alpha: 1 - // a: 1 - // beta: dt - // b: mu - // RHS: tracer + // If !use_chi, solve + // ( 1 - dt div eta grad ) sca^(n+1) = sca^(*,n+1) + // alpha: 1 + // a: 1 + // beta: dt + // b: eta + // RHS: a_scalar + // Which corresponds to the non-conservative scalar equation: + // d(sca) / dt - div mu grad sca = -U dot grad sca + H if (m_verbose > 0) { amrex::Print() << "Diffusing scalars one at a time ..." << std::endl; @@ -181,17 +171,16 @@ DiffusionScalarOp::diffuse_scalar (Vector const& tracer, Vector rhs_c(finest_level+1); // Note only conservative uses this rhs_c container for (int lev = 0; lev <= finest_level; ++lev) { - rhs_c[lev].define(tracer[lev]->boxArray(), tracer[lev]->DistributionMap(), 1, 0); + rhs_c[lev].define(a_scalar[lev]->boxArray(), a_scalar[lev]->DistributionMap(), 1, 0); } - auto iconserv = m_incflo->get_tracer_iconserv(); #ifdef AMREX_USE_EB if (m_eb_scal_solve_op) { m_eb_scal_solve_op->setScalars(1.0, dt); for (int lev = 0; lev <= finest_level; ++lev) { - if ( iconserv[0] ) { - m_eb_scal_solve_op->setACoeffs(lev, *density[lev]); + if ( use_chi[0] ) { + m_eb_scal_solve_op->setACoeffs(lev, *chi[lev]); } else { m_eb_scal_solve_op->setACoeffs(lev, 1.0); } @@ -202,19 +191,24 @@ DiffusionScalarOp::diffuse_scalar (Vector const& tracer, { m_reg_scal_solve_op->setScalars(1.0, dt); for (int lev = 0; lev <= finest_level; ++lev) { - if ( iconserv[0] ) { - m_reg_scal_solve_op->setACoeffs(lev, *density[lev]); + if ( use_chi[0] ) { + m_reg_scal_solve_op->setACoeffs(lev, *chi[lev]); } else { m_reg_scal_solve_op->setACoeffs(lev, 1.0); } } } - for (int comp = 0; comp < tracer[0]->nComp(); ++comp) + for (int comp = 0; comp < a_scalar[0]->nComp(); ++comp) { #ifdef AMREX_USE_EB if (m_eb_scal_solve_op) { + m_eb_scal_solve_op->setDomainBC(m_incflo->get_diffuse_scalar_bc(Orientation::low, + bcrec[0].lo()), + m_incflo->get_diffuse_scalar_bc(Orientation::high, + bcrec[0].hi())); + if ( m_incflo->m_has_mixedBC && comp>0 ) { // Must reset scalars (and Acoef, done below) to reuse solver with Robin BC m_eb_scal_solve_op->setScalars(1.0, dt); @@ -223,16 +217,16 @@ DiffusionScalarOp::diffuse_scalar (Vector const& tracer, for (int lev = 0; lev <= finest_level; ++lev) { // Only reset Acoeff if necessary. if (comp > 0 && - ( m_incflo->m_has_mixedBC || (iconserv[comp] != iconserv[comp-1]) ) ) { - if ( iconserv[comp] ) { - m_eb_scal_solve_op->setACoeffs(lev, *density[lev]); + ( m_incflo->m_has_mixedBC || (use_chi[comp] != use_chi[comp-1]) ) ) { + if ( use_chi[comp] ) { + m_eb_scal_solve_op->setACoeffs(lev, *chi[lev]); } else { m_eb_scal_solve_op->setACoeffs(lev, 1.0); } } - if (m_incflo->hasEBFlow()) { - MultiFab phi(*m_incflo->get_tracer_eb()[lev], amrex::make_alias, comp, 1); + if (!eb_dirichlet[lev]->empty()) { + MultiFab phi(*eb_dirichlet[lev], amrex::make_alias, comp, 1); m_eb_scal_solve_op->setEBDirichlet(lev, phi, *eta[lev]); } // else use default homogeneous Neumann on EB @@ -244,10 +238,16 @@ DiffusionScalarOp::diffuse_scalar (Vector const& tracer, else #endif { + amrex::ignore_unused(eb_dirichlet); + m_reg_scal_solve_op->setDomainBC(m_incflo->get_diffuse_scalar_bc(Orientation::low, + bcrec[comp].lo()), + m_incflo->get_diffuse_scalar_bc(Orientation::high, + bcrec[comp].hi())); + for (int lev = 0; lev <= finest_level; ++lev) { - if ( comp > 0 && (iconserv[comp] != iconserv[comp-1]) ) { - if ( iconserv[comp] ) { - m_reg_scal_solve_op->setACoeffs(lev, *density[lev]); + if ( comp > 0 && (use_chi[comp] != use_chi[comp-1]) ) { + if ( use_chi[comp] ) { + m_reg_scal_solve_op->setACoeffs(lev, *chi[lev]); } else { m_reg_scal_solve_op->setACoeffs(lev, 1.0); } @@ -261,10 +261,10 @@ DiffusionScalarOp::diffuse_scalar (Vector const& tracer, Vector phi; Vector rhs; for (int lev = 0; lev <= finest_level; ++lev) { - phi.emplace_back(*tracer[lev], amrex::make_alias, comp, 1); + phi.emplace_back(*a_scalar[lev], amrex::make_alias, comp, 1); - if ( !iconserv[comp] ) { - rhs.emplace_back(*tracer[lev], amrex::make_alias, comp, 1); + if ( !use_chi[comp] ) { + rhs.emplace_back(*a_scalar[lev], amrex::make_alias, comp, 1); } else { rhs.emplace_back(rhs_c[lev], amrex::make_alias, 0, 1); #ifdef _OPENMP @@ -273,11 +273,11 @@ DiffusionScalarOp::diffuse_scalar (Vector const& tracer, for (MFIter mfi(rhs[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) { Box const& bx = mfi.tilebox(); Array4 const& rhs_a = rhs[lev].array(mfi); - Array4 const& tra_a = tracer[lev]->const_array(mfi,comp); - Array4 const& rho_a = density[lev]->const_array(mfi); + Array4 const& sca_a = a_scalar[lev]->const_array(mfi,comp); + Array4 const& chi_a = chi[lev]->const_array(mfi); ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - rhs_a(i,j,k) = rho_a(i,j,k) * tra_a(i,j,k); + rhs_a(i,j,k) = chi_a(i,j,k) * sca_a(i,j,k); }); } } @@ -487,11 +487,13 @@ DiffusionScalarOp::diffuse_vel_components (Vector const& vel, void DiffusionScalarOp::compute_laps (Vector const& a_laps, Vector const& a_scalar, - Vector const& a_eta) + Vector const& a_eta, + amrex::Vector bcrec) { BL_PROFILE("DiffusionScalarOp::compute_laps"); int finest_level = m_incflo->finestLevel(); + int n_comp = a_laps[0]->nComp(); #ifdef AMREX_USE_EB if (m_eb_scal_apply_op) @@ -501,7 +503,7 @@ void DiffusionScalarOp::compute_laps (Vector const& a_laps, for (int lev = 0; lev <= finest_level; ++lev) { laps_tmp[lev].define(a_laps[lev]->boxArray(), a_laps[lev]->DistributionMap(), - m_incflo->m_ntrac, tmp_comp, MFInfo(), + n_comp, tmp_comp, MFInfo(), a_laps[lev]->Factory()); laps_tmp[lev].setVal(0.0); } @@ -513,8 +515,12 @@ void DiffusionScalarOp::compute_laps (Vector const& a_laps, // For when we use the stencil for centroid values // m_eb_scal_apply_op->setPhiOnCentroid(); - // FIXME? Can we do the solve together now? - for (int comp = 0; comp < m_incflo->m_ntrac; ++comp) { + for (int comp = 0; comp < n_comp; ++comp) { + m_eb_scal_apply_op->setDomainBC(m_incflo->get_diffuse_scalar_bc(Orientation::low, + bcrec[comp].lo()), + m_incflo->get_diffuse_scalar_bc(Orientation::high, + bcrec[comp].hi())); + int eta_comp = comp; if ( m_incflo->m_has_mixedBC && comp>0 ){ @@ -551,9 +557,13 @@ void DiffusionScalarOp::compute_laps (Vector const& a_laps, for(int lev = 0; lev <= finest_level; lev++) { + // Note that this is strictly flux redistribution + // May also now be written in a way that it's safe to have laps_tmp = a_laps amrex::single_level_redistribute(laps_tmp[lev], - *a_laps[lev], 0, m_incflo->m_ntrac, + *a_laps[lev], 0, n_comp, m_incflo->Geom(lev)); + // FIXME? This will use SRD if that's what's set for m_redistribution_type + // Need to think about how to make this work with temperature b/c BC is different // auto const& bc = m_incflo->get_tracer_bcrec_device_ptr(); // m_incflo->redistribute_term(*a_laps[lev], laps_tmp[lev], *a_scalar[lev], // bc, lev); @@ -565,7 +575,11 @@ 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); - for (int comp = 0; comp < m_incflo->m_ntrac; ++comp) { + for (int comp = 0; comp < n_comp; ++comp) { + m_reg_scal_apply_op->setDomainBC(m_incflo->get_diffuse_scalar_bc(Orientation::low, + bcrec[comp].lo()), + m_incflo->get_diffuse_scalar_bc(Orientation::high, + bcrec[comp].hi())); int eta_comp = comp; diff --git a/src/diffusion/incflo_diffusion.cpp b/src/diffusion/incflo_diffusion.cpp index e869692c..248bd8c8 100644 --- a/src/diffusion/incflo_diffusion.cpp +++ b/src/diffusion/incflo_diffusion.cpp @@ -59,7 +59,18 @@ incflo::compute_laps(Vector const& laps, Vector const& scalar, Vector const& eta) { - get_diffusion_scalar_op()->compute_laps(laps, scalar, eta); + get_diffusion_scalar_op()->compute_laps(laps, scalar, eta, + get_tracer_bcrec()); + +} + +void +incflo::compute_laps_T(Vector const& laps, + Vector const& scalar, + Vector const& eta) +{ + get_diffusion_scalar_op()->compute_laps(laps, scalar, eta, + get_tracer_bcrec()); } void @@ -68,9 +79,22 @@ incflo::diffuse_scalar(Vector const& scalar, Vector const& eta, Real dt_diff) { - get_diffusion_scalar_op()->diffuse_scalar(scalar, density, eta, dt_diff); + get_diffusion_scalar_op()->diffuse_scalar(scalar, density, eta, get_tracer_eb(), + get_tracer_iconserv(), + get_tracer_bcrec(), dt_diff); } +void +incflo::diffuse_temperature(Vector const& temperature, + Vector const& rhocp, + Vector const& eta, + Real dt_diff) +{ + get_diffusion_scalar_op()->diffuse_scalar(temperature, rhocp, eta, + get_temperature_eb(), + {1} /* use rhocp */, + get_temperature_bcrec(), dt_diff); +} void incflo::diffuse_velocity(Vector const& vel, diff --git a/src/incflo.H b/src/incflo.H index c39940e9..c852567f 100644 --- a/src/incflo.H +++ b/src/incflo.H @@ -96,6 +96,8 @@ public: void compute_tra_forces (amrex::Vector const& tra_forces, amrex::Vector const& density); + void compute_tem_forces (amrex::Real time, + amrex::Vector const& tem_forces); void compute_vel_forces (amrex::Vector const& vel_forces, amrex::Vector const& velocity, amrex::Vector const& density, @@ -129,6 +131,7 @@ public: void set_eb_velocity (int lev, amrex::Real time, amrex::MultiFab& eb_vel, int nghost); void set_eb_density (int lev, amrex::Real time, amrex::MultiFab& eb_density, int nghost); void set_eb_tracer (int lev, amrex::Real time, amrex::MultiFab& eb_tracer, int nghost); + void set_eb_temperature (int lev, amrex::Real time, amrex::MultiFab& eb_temperature, int nghost); #endif /////////////////////////////////////////////////////////////////////////// @@ -141,17 +144,41 @@ public: void ApplyCorrector (); void compute_convective_term (amrex::Vector const& conv_u, amrex::Vector const& conv_r, - amrex::Vector const& conv_t, + amrex::Vector const& conv_tra, + amrex::Vector const& conv_tem, amrex::Vector const& vel, amrex::Vector const& density, amrex::Vector const& tracer, + amrex::Vector const& temperature, AMREX_D_DECL(amrex::Vector const& u_mac, amrex::Vector const& v_mac, amrex::Vector const& w_mac), amrex::Vector const& vel_forces, amrex::Vector const& tra_forces, + amrex::Vector const& tem_forces, amrex::Real time); + void compute_convective_term (amrex::Vector const& conv_u, + amrex::Vector const& conv_r, + amrex::Vector const& conv_tra, + amrex::Vector const& vel, + amrex::Vector const& density, + amrex::Vector const& tracer, + AMREX_D_DECL(amrex::Vector const& u_mac, + amrex::Vector const& v_mac, + amrex::Vector const& w_mac), + amrex::Vector const& vel_forces, + amrex::Vector const& tra_forces, + amrex::Real time) + { + if (m_use_temperature) { amrex::Abort("Calling wrong compute_convective_term()"); } + compute_convective_term(conv_u, conv_r, conv_tra, {}, + vel, density, tracer, {}, + AMREX_D_DECL(u_mac, v_mac, w_mac), + vel_forces, tra_forces, {}, time); + + } + void compute_MAC_projected_velocities ( amrex::Vector const& vel, amrex::Vector const& density, @@ -167,6 +194,8 @@ public: void update_density (StepType step_type); void update_tracer (StepType step_type, amrex::Vector& tra_eta, amrex::Vector& tra_forces); + void update_temperature (StepType step_type, amrex::Vector& tem_eta, + amrex::Vector& scratch); void update_velocity (StepType step_type, amrex::Vector& vel_eta, amrex::Vector& vel_forces); @@ -220,11 +249,20 @@ public: amrex::Vector const& scalar, amrex::Vector const& eta); + void compute_laps_T (amrex::Vector const& laps, + amrex::Vector const& scalar, + amrex::Vector const& eta); + void diffuse_scalar (amrex::Vector const& scalar, amrex::Vector const& density, amrex::Vector const& eta, amrex::Real dt_diff); + void diffuse_temperature (amrex::Vector const& temperature, + amrex::Vector const& rhocp, + amrex::Vector const& eta, + amrex::Real dt_diff); + void fixup_eta_on_domain_faces (int lev, amrex::Array& fc, amrex::MultiFab const& cc) const; @@ -305,7 +343,11 @@ public: amrex::MultiFab* vel, amrex::Geometry& lev_geom, amrex::Real time, int nghost); - void compute_tracer_diff_coeff (amrex::Vector const& tra_eta, int nghost); + void compute_tracer_diff_coeff (amrex::Vector const& tra_eta, int nghost) const; + void compute_temperature_diff_coeff (amrex::Real time, + amrex::Vector const& eta) const; + + void compute_cp (int lev, amrex::MFIter& mfi, amrex::FArrayBox& cp) const; #ifdef AMREX_USE_EB /////////////////////////////////////////////////////////////////////////// @@ -358,6 +400,7 @@ private: amrex::Real m_ic_w = amrex::Real(0.0); amrex::Real m_ic_p = amrex::Real(0.0); amrex::Vector m_ic_t = amrex::Vector (m_ntrac, amrex::Real(0.0)); + amrex::Real m_ic_tem = amrex::Real(1.0); amrex::Vector m_t_old; amrex::Vector m_t_new; @@ -467,6 +510,10 @@ private: // Scalar diffusive coefficient amrex::Vector m_mu_s; + // Thermal diffusivity + amrex::Real m_mu_T = 0.; + // specific heat capacity + amrex::Real m_cp = 1.; // Density (if constant) amrex::Real m_ro_0 = 1.0; @@ -503,6 +550,7 @@ private: // scalars amrex::Real density{1.}; amrex::Vector tracer; + amrex::Vector temperature; EBFlow_t () @@ -632,6 +680,9 @@ private: amrex::MultiFab tracer; amrex::MultiFab tracer_eb; amrex::MultiFab tracer_o; + amrex::MultiFab temperature; + amrex::MultiFab temperature_eb; + amrex::MultiFab temperature_o; amrex::MultiFab mac_phi; // cell-centered pressure used in MAC projection @@ -648,13 +699,19 @@ private: amrex::MultiFab conv_density_o; amrex::MultiFab conv_tracer; amrex::MultiFab conv_tracer_o; + amrex::MultiFab conv_temperature; + amrex::MultiFab conv_temperature_o; amrex::MultiFab divtau; amrex::MultiFab divtau_o; amrex::MultiFab laps; amrex::MultiFab laps_o; + amrex::MultiFab laps_tem; + amrex::MultiFab laps_tem_o; }; + bool m_use_temperature = false; + amrex::Vector > m_leveldata; amrex::Vector > > m_factory; @@ -672,6 +729,7 @@ private: amrex::GpuArray , AMREX_SPACEDIM*2> m_bc_velocity; amrex::GpuArray , AMREX_SPACEDIM*2> m_bc_tracer; + amrex::GpuArray m_bc_temperature; amrex::Vector m_bc_eb_velocity; // amrex::Vector cannot be used on gpu, so ... @@ -684,6 +742,8 @@ private: amrex::Gpu::DeviceVector m_bcrec_density_d; amrex::Vector m_bcrec_tracer; amrex::Gpu::DeviceVector m_bcrec_tracer_d; + amrex::Vector m_bcrec_temperature; + amrex::Gpu::DeviceVector m_bcrec_temperature_d; amrex::Vector m_bcrec_force; amrex::Gpu::DeviceVector m_bcrec_force_d; @@ -693,6 +753,8 @@ private: amrex::Gpu::DeviceVector m_iconserv_density_d; amrex::Vector m_iconserv_tracer; amrex::Gpu::DeviceVector m_iconserv_tracer_d; + amrex::Vector m_iconserv_temperature; + amrex::Gpu::DeviceVector m_iconserv_temperature_d; std::unique_ptr m_diffusion_tensor_op; std::unique_ptr m_diffusion_scalar_op; @@ -782,24 +844,32 @@ private: amrex::Vector get_velocity_eb () noexcept; amrex::Vector get_density_eb () noexcept; amrex::Vector get_tracer_eb () noexcept; + amrex::Vector get_temperature_eb () noexcept; amrex::Vector get_density_old () noexcept; amrex::Vector get_density_new () noexcept; amrex::Vector get_density_nph () noexcept; amrex::Vector get_tracer_old () noexcept; amrex::Vector get_tracer_new () noexcept; + amrex::Vector get_temperature_old () noexcept; + amrex::Vector get_temperature_new () noexcept; amrex::Vector get_mac_phi () noexcept; amrex::Vector get_vel_forces () noexcept; amrex::Vector get_tra_forces () noexcept; + amrex::Vector get_tem_forces () noexcept; amrex::Vector get_conv_velocity_old () noexcept; amrex::Vector get_conv_velocity_new () noexcept; amrex::Vector get_conv_density_old () noexcept; amrex::Vector get_conv_density_new () noexcept; amrex::Vector get_conv_tracer_old () noexcept; amrex::Vector get_conv_tracer_new () noexcept; + amrex::Vector get_conv_temperature_old () noexcept; + amrex::Vector get_conv_temperature_new () noexcept; amrex::Vector get_divtau_old () noexcept; amrex::Vector get_divtau_new () noexcept; amrex::Vector get_laps_old () noexcept; amrex::Vector get_laps_new () noexcept; + amrex::Vector get_laps_tem_old () noexcept; + amrex::Vector get_laps_tem_new () noexcept; // [[nodiscard]] amrex::Vector get_velocity_old_const () const noexcept; [[nodiscard]] amrex::Vector get_velocity_new_const () const noexcept; @@ -811,8 +881,11 @@ private: [[nodiscard]] amrex::Vector get_density_nph_const () const noexcept; [[nodiscard]] amrex::Vector get_tracer_old_const () const noexcept; [[nodiscard]] amrex::Vector get_tracer_new_const () const noexcept; + [[nodiscard]] amrex::Vector get_temperature_old_const () const noexcept; + [[nodiscard]] amrex::Vector get_temperature_new_const () const noexcept; [[nodiscard]] amrex::Vector get_vel_forces_const () const noexcept; [[nodiscard]] amrex::Vector get_tra_forces_const () const noexcept; + [[nodiscard]] amrex::Vector get_tem_forces_const () const noexcept; [[nodiscard]] amrex::Vector const& get_velocity_iconserv () const noexcept { return m_iconserv_velocity; } [[nodiscard]] amrex::Vector const& get_density_iconserv () const noexcept { return m_iconserv_density; } @@ -828,6 +901,7 @@ private: [[nodiscard]] amrex::Vector const& get_velocity_bcrec () const noexcept { return m_bcrec_velocity; } [[nodiscard]] amrex::Vector const& get_density_bcrec () const noexcept { return m_bcrec_density; } [[nodiscard]] amrex::Vector const& get_tracer_bcrec () const noexcept { return m_bcrec_tracer; } + [[nodiscard]] amrex::Vector const& get_temperature_bcrec () const noexcept { return m_bcrec_temperature; } [[nodiscard]] amrex::Vector const& get_force_bcrec () const noexcept { return m_bcrec_force; } // [[nodiscard]] amrex::BCRec const* get_velocity_bcrec_device_ptr () const noexcept { @@ -836,6 +910,8 @@ private: return m_bcrec_density_d.data(); } [[nodiscard]] amrex::BCRec const* get_tracer_bcrec_device_ptr () const noexcept { return m_bcrec_tracer_d.data(); } + [[nodiscard]] amrex::BCRec const* get_temperature_bcrec_device_ptr () const noexcept { + return m_bcrec_temperature_d.data(); } [[nodiscard]] amrex::BCRec const* get_force_bcrec_device_ptr () const noexcept { return m_bcrec_force_d.data(); } @@ -856,17 +932,20 @@ private: void fillpatch_velocity (int lev, amrex::Real time, amrex::MultiFab& vel, int ng); void fillpatch_density (int lev, amrex::Real time, amrex::MultiFab& density, int ng); void fillpatch_tracer (int lev, amrex::Real time, amrex::MultiFab& tracer, int ng); + void fillpatch_temperature (int lev, amrex::Real time, amrex::MultiFab& temperature, int ng); void fillpatch_gradp (int lev, amrex::Real time, amrex::MultiFab& gp, int ng); void fillpatch_force (amrex::Real time, amrex::Vector const& force, int ng); void fillcoarsepatch_velocity (int lev, amrex::Real time, amrex::MultiFab& vel, int ng); void fillcoarsepatch_density (int lev, amrex::Real time, amrex::MultiFab& density, int ng); void fillcoarsepatch_tracer (int lev, amrex::Real time, amrex::MultiFab& tracer, int ng); + void fillcoarsepatch_temperature (int lev, amrex::Real time, amrex::MultiFab& temperature, int ng); void fillcoarsepatch_gradp (int lev, amrex::Real time, amrex::MultiFab& gp, int ng); void fillphysbc_velocity (int lev, amrex::Real time, amrex::MultiFab& vel, int ng); void fillphysbc_density (int lev, amrex::Real time, amrex::MultiFab& density, int ng); void fillphysbc_tracer (int lev, amrex::Real time, amrex::MultiFab& tracer, int ng); + void fillphysbc_temperature (int lev, amrex::Real time, amrex::MultiFab& temperature, int ng); void copy_from_new_to_old_velocity ( amrex::IntVect const& ng = amrex::IntVect{0}); void copy_from_new_to_old_velocity (int lev, amrex::IntVect const& ng = amrex::IntVect{0}); @@ -874,6 +953,7 @@ private: void copy_from_new_to_old_density (int lev, amrex::IntVect const& ng = amrex::IntVect{0}); void copy_from_new_to_old_tracer ( amrex::IntVect const& ng = amrex::IntVect{0}); void copy_from_new_to_old_tracer (int lev, amrex::IntVect const& ng = amrex::IntVect{0}); + void copy_from_new_to_old_temperature ( amrex::IntVect const& ng = amrex::IntVect{0}); // void copy_from_old_to_new_velocity ( amrex::IntVect const& ng = amrex::IntVect{0}); void copy_from_old_to_new_velocity (int lev, amrex::IntVect const& ng = amrex::IntVect{0}); @@ -881,6 +961,7 @@ private: void copy_from_old_to_new_density (int lev, amrex::IntVect const& ng = amrex::IntVect{0}); void copy_from_old_to_new_tracer ( amrex::IntVect const& ng = amrex::IntVect{0}); void copy_from_old_to_new_tracer (int lev, amrex::IntVect const& ng = amrex::IntVect{0}); + void copy_from_old_to_new_temperature ( amrex::IntVect const& ng = amrex::IntVect{0}); void Advance (); bool writeNow () { return writeNow(m_plot_int, m_plot_per_approx, m_plot_per_exact); } diff --git a/src/incflo_advance.cpp b/src/incflo_advance.cpp index 9b387433..9dee90ba 100644 --- a/src/incflo_advance.cpp +++ b/src/incflo_advance.cpp @@ -32,6 +32,7 @@ void incflo::Advance() copy_from_new_to_old_velocity(); copy_from_new_to_old_density(); copy_from_new_to_old_tracer(); + copy_from_new_to_old_temperature(); int ng = nghost_state(); for (int lev = 0; lev <= finest_level; ++lev) { @@ -40,6 +41,9 @@ void incflo::Advance() if (m_advect_tracer) { fillpatch_tracer(lev, m_t_old[lev], m_leveldata[lev]->tracer_o, ng); } + if (m_use_temperature) { + fillpatch_temperature(lev, m_t_old[lev], m_leveldata[lev]->temperature_o, ng); + } } #ifdef AMREX_USE_EB @@ -47,9 +51,18 @@ void incflo::Advance() for (int lev = 0; lev <= finest_level; ++lev) { set_eb_velocity(lev, m_t_old[lev], *get_velocity_eb()[lev], 1); set_eb_density(lev, m_t_old[lev], *get_density_eb()[lev], 1); - set_eb_tracer(lev, m_t_old[lev], *get_tracer_eb()[lev], 1); } } + if (m_advect_tracer && !m_eb_flow.tracer.empty()) { + for (int lev = 0; lev <= finest_level; ++lev) { + set_eb_tracer(lev, m_t_old[lev], *get_tracer_eb()[lev], 1); + } + } + if (m_use_temperature && !m_eb_flow.temperature.empty()) { + for (int lev = 0; lev <= finest_level; ++lev) { + set_eb_temperature(lev, m_t_old[lev], *get_temperature_eb()[lev], 1); + } + } #endif ApplyPredictor(); @@ -61,6 +74,9 @@ void incflo::Advance() if (m_advect_tracer) { fillpatch_tracer(lev, m_t_new[lev], m_leveldata[lev]->tracer, ng); } + if (m_use_temperature) { + fillpatch_temperature(lev, m_t_new[lev], m_leveldata[lev]->temperature, ng); + } } ApplyCorrector(); diff --git a/src/incflo_apply_predictor.cpp b/src/incflo_apply_predictor.cpp index 58296b86..114e3a7e 100644 --- a/src/incflo_apply_predictor.cpp +++ b/src/incflo_apply_predictor.cpp @@ -6,15 +6,17 @@ using namespace amrex; // // 1. Use u = vel_old to compute // -// if (!advect_momentum) then -// conv_u = - u grad u -// else -// conv_u = - del dot (rho u u) +// +// conv_u = - u grad u , if (!advect_momentum) +// = - del dot (rho u u), otherwise +// // conv_r = - div( u rho ) -// if (m_iconserv_tracer) then -// conv_t = - div( u trac ) -// else -// conv_t = - u dot grad trac +// +// conv_tra = - div( u trac ) , if (m_iconserv_tracer) +// = - u dot grad trac, otherwise +// +// conv_tem = - u dot grad temperature +// // eta_old = visosity at m_cur_time // if (m_diff_type == DiffusionType::Explicit) // divtau _old = div( eta ( (grad u) + (grad u)^T ) ) / rho^n @@ -23,24 +25,30 @@ using namespace amrex; // divtau_old = 0.0 // rhs = u + dt * conv // -// eta = eta at new_time -// // 2. Add explicit forcing term i.e. gravity + lagged pressure gradient // // rhs += dt * ( g - grad(p + p0) / rho^nph ) // // 3. A. If (m_diff_type == DiffusionType::Implicit) -// solve implicit diffusion equation for u* +// solve implicit diffusion equation for u*, tracers, temperature +// +// ( 1 - dt / rho^nph * div ( eta grad ) ) u* = u^n + dt * conv_u +// + dt * ( g - grad(p + p0) / rho^nph ) // -// ( 1 - dt / rho^nph * div ( eta grad ) ) u* = u^n + dt * conv_u -// + dt * ( g - grad(p + p0) / rho^nph ) +// Tracer, for conservative: +// ( rho - dt div mu grad ) tra^(n+1) = rho tra^(n) - dt * div(U rho tra) + dt * rho H +// or for non-conservative: +// ( 1 - dt div mu grad ) tra^(n+1) = tra^(n) - dt * U dot grad tra + dt * H +// +// Temperature: +// ( rho*cp - dt div mu_T grad ) T^(n+1) = rho*cp T^n - rho*cp dt U dot grad T + H_T // // B. If (m_diff_type == DiffusionType::Crank-Nicolson) -// solve semi-implicit diffusion equation for u* +// solve semi-implicit diffusion equation for u*, tracers, temperature // -// ( 1 - (dt/2) / rho^nph * div ( eta_old grad ) ) u* = u^n + -// dt * conv_u + (dt/2) / rho * div (eta_old grad) u^n -// + dt * ( g - grad(p + p0) / rho^nph ) +// ( 1 - (dt/2) / rho^nph * div ( eta_old grad ) ) u* = u^n + +// dt * conv_u + (dt/2) / rho * div (eta_old grad) u^n +// + dt * ( g - grad(p + p0) / rho^nph ) // // 4. Apply projection // @@ -97,25 +105,30 @@ void incflo::ApplyPredictor (bool incremental_projection) // ************************************************************************************* // Allocate space for half-time density // ************************************************************************************* - // Forcing terms - Vector vel_forces, tra_forces; + // Forcing terms for velocity, tracers, temperature + Vector vel_forces, tra_forces, tem_forces; - Vector vel_eta, tra_eta; + Vector vel_eta, tra_eta, tem_eta; // ************************************************************************************* - // Allocate space for the forcing terms + // Allocate space for the forcing terms and viscosity / diffusive coefficients // ************************************************************************************* + int nghost_eta = 1; for (int lev = 0; lev <= finest_level; ++lev) { vel_forces.emplace_back(grids[lev], dmap[lev], AMREX_SPACEDIM, nghost_force(), MFInfo(), Factory(lev)); + vel_eta.emplace_back(grids[lev], dmap[lev], 1, nghost_eta, MFInfo(), Factory(lev)); if (m_advect_tracer) { tra_forces.emplace_back(grids[lev], dmap[lev], m_ntrac, nghost_force(), MFInfo(), Factory(lev)); + tra_eta.emplace_back(grids[lev], dmap[lev], m_ntrac, nghost_eta, + MFInfo(), Factory(lev)); } - vel_eta.emplace_back(grids[lev], dmap[lev], 1, 1, MFInfo(), Factory(lev)); - if (m_advect_tracer) { - tra_eta.emplace_back(grids[lev], dmap[lev], m_ntrac, 1, MFInfo(), Factory(lev)); + if (m_use_temperature) { + tem_forces.emplace_back(grids[lev], dmap[lev], 1, nghost_force(), + MFInfo(), Factory(lev)); + tem_eta.emplace_back(grids[lev], dmap[lev], 1, nghost_eta, MFInfo(), Factory(lev)); } } @@ -124,7 +137,7 @@ void incflo::ApplyPredictor (bool incremental_projection) // ************************************************************************************* compute_viscosity(GetVecOfPtrs(vel_eta), get_density_old(), get_velocity_old(), - m_cur_time, 1); + m_cur_time, nghost_eta); // ************************************************************************************* // Compute explicit viscous term @@ -141,11 +154,18 @@ void incflo::ApplyPredictor (bool incremental_projection) // ************************************************************************************* if (m_advect_tracer) { - compute_tracer_diff_coeff(GetVecOfPtrs(tra_eta),1); + compute_tracer_diff_coeff(GetVecOfPtrs(tra_eta), nghost_eta); if (need_divtau()) { compute_laps(get_laps_old(), get_tracer_old_const(), GetVecOfConstPtrs(tra_eta)); } } + if (m_use_temperature) + { + compute_temperature_diff_coeff(m_cur_time, GetVecOfPtrs(tem_eta)); + if (need_divtau()) { + compute_laps_T(get_laps_tem_old(), get_temperature_old_const(), GetVecOfConstPtrs(tem_eta)); + } + } // ********************************************************************************************** // Compute the forcing terms @@ -170,11 +190,13 @@ void incflo::ApplyPredictor (bool incremental_projection) // Note that if advection_type != "MOL" then we call compute_tra_forces inside this routine // ************************************************************************************* compute_convective_term(get_conv_velocity_old(), get_conv_density_old(), get_conv_tracer_old(), + get_conv_temperature_old(), get_velocity_old_const(), get_density_old_const(), get_tracer_old_const(), + get_temperature_old_const(), AMREX_D_DECL(GetVecOfPtrs(u_mac), GetVecOfPtrs(v_mac), - GetVecOfPtrs(w_mac)), + GetVecOfPtrs(w_mac)), GetVecOfPtrs(vel_forces), GetVecOfPtrs(tra_forces), - m_cur_time); + GetVecOfPtrs(tem_forces), m_cur_time); // ************************************************************************************* // Update density @@ -186,6 +208,11 @@ void incflo::ApplyPredictor (bool incremental_projection) // ********************************************************************************************** update_tracer(StepType::Predictor, tra_eta, tra_forces); + // ********************************************************************************************** + // Update temperature + // ********************************************************************************************** + update_temperature(StepType::Predictor, tem_eta, tem_forces); + // ********************************************************************************************** // Update velocity // ********************************************************************************************** diff --git a/src/incflo_compute_forces.cpp b/src/incflo_compute_forces.cpp index 8ad6af7b..003fe183 100644 --- a/src/incflo_compute_forces.cpp +++ b/src/incflo_compute_forces.cpp @@ -2,10 +2,36 @@ using namespace amrex; +// Compute temperature forcing terms. +void incflo::compute_tem_forces (Real /*time*/, Vector const& tem_forces) +{ + if (m_use_temperature) { + +#ifdef _OPENMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (int lev = 0; lev <= finest_level; ++lev) { + for (MFIter mfi(*tem_forces[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + Box const& bx = mfi.tilebox(); + Array4 const& tem_f = tem_forces[lev]->array(mfi); + + ParallelFor(bx, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + // Just H_T here, as in + // rho*cp ( dT/dt + U dot grad T) = div mu_T grad T + H_T + // For now we don't have any external forces on temperature + tem_f(i,j,k) = 0.0; + }); + } + } + } +} + void incflo::compute_tra_forces (Vector const& tra_forces, Vector const& density) { - // 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(); diff --git a/src/incflo_regrid.cpp b/src/incflo_regrid.cpp index 87c507d0..2663d924 100644 --- a/src/incflo_regrid.cpp +++ b/src/incflo_regrid.cpp @@ -33,6 +33,9 @@ void incflo::MakeNewLevelFromCoarse (int lev, if (m_ntrac > 0) { fillcoarsepatch_tracer(lev, time, new_leveldata->tracer, 0); } + if (m_use_temperature) { + fillcoarsepatch_temperature(lev, time, new_leveldata->temperature, 0); + } fillcoarsepatch_gradp(lev, time, new_leveldata->gp, 0); if (m_use_cc_proj) { @@ -86,6 +89,9 @@ void incflo::RemakeLevel (int lev, Real time, const BoxArray& ba, if (m_ntrac > 0) { fillpatch_tracer(lev, time, new_leveldata->tracer, 0); } + if (m_use_temperature) { + fillpatch_temperature(lev, time, new_leveldata->temperature, 0); + } fillpatch_gradp(lev, time, new_leveldata->gp, 0); if (m_use_cc_proj) { diff --git a/src/incflo_specific_heat.cpp b/src/incflo_specific_heat.cpp new file mode 100644 index 00000000..2c43d6f8 --- /dev/null +++ b/src/incflo_specific_heat.cpp @@ -0,0 +1,17 @@ +#include + +using namespace amrex; + +void incflo::compute_cp (int /*lev*/, MFIter& /*mfi*/, FArrayBox& cp) const +{ + // Get leveldata if desired, e.g. + // Array4 const& rho = m_leveldata[lev]->density.const_array(mfi); + + Box const& bx = cp.box(); + Array4 const& cp_a = cp.array(); + Real l_cp = m_cp; + ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + cp_a(i,j,k) = l_cp; + }); +} diff --git a/src/incflo_update_temperature.cpp b/src/incflo_update_temperature.cpp new file mode 100644 index 00000000..f2d7c05f --- /dev/null +++ b/src/incflo_update_temperature.cpp @@ -0,0 +1,126 @@ +#include + +using namespace amrex; + +void incflo::update_temperature (StepType step_type, Vector& tem_eta, Vector& scratch) +{ + BL_PROFILE("incflo::update_temperature"); + + if (!m_use_temperature) { return; } + + Real const new_time = m_cur_time + m_dt; + Real const half_time = m_cur_time + m_dt/2.; + + // ************************************************************************************* + // Compute the temperature forcing terms + // ************************************************************************************* + compute_tem_forces(half_time, GetVecOfPtrs(scratch)); + + // ************************************************************************************* + // Compute explicit diffusive term (if corrector) + // ************************************************************************************* + if (step_type == StepType::Corrector) + { + compute_temperature_diff_coeff(new_time, GetVecOfPtrs(tem_eta)); + if (m_diff_type == DiffusionType::Explicit) { + compute_laps_T(get_laps_new(), get_temperature_new_const(), GetVecOfConstPtrs(tem_eta)); + } + } + + // ************************************************************************************* + // Update the temperature with time-explicit terms + // ************************************************************************************* + if (step_type == StepType::Predictor) { + constexpr Real m_half = Real(0.5); + Real l_dt = m_dt; + + 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& tem_o = ld.temperature_o.const_array(mfi); + Array4 const& tem = ld.temperature.array(mfi); + Array4 const& rho_h = ld.density_nph.const_array(mfi); + Array4 const& dtdt_o = ld.conv_temperature_o.const_array(mfi); + // temperature forcing term (Q) is in scratch + Array4 const& tem_f = scratch[lev].array(mfi); + + FArrayBox cp_fab(bx, 1, The_Async_Arena()); + compute_cp(lev, mfi, cp_fab); + Array4 const& cp = cp_fab.array(); + + if (m_diff_type == DiffusionType::Explicit) + { + Array4 const& laps_o = ld.laps_tem_o.const_array(mfi); + + ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + tem(i,j,k) = tem_o(i,j,k) + l_dt * + ( dtdt_o(i,j,k) + (tem_f(i,j,k) + laps_o(i,j,k))/(rho_h(i,j,k) * cp(i,j,k)) ); + }); + } + else if (m_diff_type == DiffusionType::Crank_Nicolson) + { + Array4 const& laps_o = ld.laps_o.const_array(mfi); + + ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + tem(i,j,k) = tem_o(i,j,k) + l_dt * + ( dtdt_o(i,j,k) + (tem_f(i,j,k) + m_half*laps_o(i,j,k))/(rho_h(i,j,k) * cp(i,j,k)) ); + // Save rhoCp for use in implicit solve. + // Reuse scratch space since we are done with forcing now. + tem_f(i,j,k) = rho_h(i,j,k) * cp(i,j,k); + }); + } + else if (m_diff_type == DiffusionType::Implicit) + { + ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + tem(i,j,k) = tem_o(i,j,k) + l_dt * + (dtdt_o(i,j,k) + tem_f(i,j,k)) / (rho_h(i,j,k) * cp(i,j,k)); + // Save rhoCp for use in implicit solve. + // Reuse scratch space since we are done with forcing now. + tem_f(i,j,k) = rho_h(i,j,k) * cp(i,j,k); + }); + } + } // mfi + } // lev + + } else if (step_type == StepType::Corrector) { + Abort("incflo::update_temperature does not yet work with the corrector"); + } + + // ************************************************************************************* + // Solve implicit diffusion equation for temperature + // ************************************************************************************* + if (m_diff_type == DiffusionType::Crank_Nicolson || m_diff_type == DiffusionType::Implicit) + { + const int ng_diffusion = 1; + for (int lev = 0; lev <= finest_level; ++lev) { + fillphysbc_temperature(lev, new_time, m_leveldata[lev]->temperature, ng_diffusion); + } + Real dt_diff = (m_diff_type == DiffusionType::Implicit) ? m_dt : Real(0.5)*m_dt; + // scratch holds rhoCp + diffuse_temperature(get_temperature_new(), GetVecOfPtrs(scratch), GetVecOfConstPtrs(tem_eta), + dt_diff); + } + else + { + // Need to average down temperature since the diffusion solver didn't do it for us. + for (int lev = finest_level-1; lev >= 0; --lev) { +#ifdef AMREX_USE_EB + amrex::EB_average_down(m_leveldata[lev+1]->temperature, m_leveldata[lev]->temperature, + 0, m_ntrac, refRatio(lev)); +#else + amrex::average_down(m_leveldata[lev+1]->temperature, m_leveldata[lev]->temperature, + 0, m_ntrac, refRatio(lev)); +#endif + } + } +} diff --git a/src/incflo_utils.cpp b/src/incflo_utils.cpp index 3e155e1a..c3cf793d 100644 --- a/src/incflo_utils.cpp +++ b/src/incflo_utils.cpp @@ -32,6 +32,18 @@ Vector incflo::get_tracer_eb () noexcept return r; } +Vector incflo::get_temperature_eb () noexcept +{ + Vector r; + if (m_use_temperature) { + r.reserve(finest_level+1); + for (int lev = 0; lev <= finest_level; ++lev) { + r.push_back(&(m_leveldata[lev]->temperature_eb)); + } + } + return r; +} + Vector incflo::get_velocity_old () noexcept { Vector r; @@ -102,6 +114,30 @@ Vector incflo::get_tracer_new () noexcept return r; } +Vector incflo::get_temperature_old () noexcept +{ + Vector r; + if (m_use_temperature) { + r.reserve(finest_level+1); + for (int lev = 0; lev <= finest_level; ++lev) { + r.push_back(&(m_leveldata[lev]->temperature_o)); + } + } + return r; +} + +Vector incflo::get_temperature_new () noexcept +{ + Vector r; + if (m_use_temperature) { + r.reserve(finest_level+1); + for (int lev = 0; lev <= finest_level; ++lev) { + r.push_back(&(m_leveldata[lev]->temperature)); + } + } + return r; +} + Vector incflo::get_mac_phi() noexcept { Vector r; @@ -172,6 +208,30 @@ Vector incflo::get_conv_tracer_new () noexcept return r; } +Vector incflo::get_conv_temperature_old () noexcept +{ + Vector r; + if (m_use_temperature) { + r.reserve(finest_level+1); + for (int lev = 0; lev <= finest_level; ++lev) { + r.push_back(&(m_leveldata[lev]->conv_temperature_o)); + } + } + return r; +} + +Vector incflo::get_conv_temperature_new () noexcept +{ + Vector r; + if (m_use_temperature) { + r.reserve(finest_level+1); + for (int lev = 0; lev <= finest_level; ++lev) { + r.push_back(&(m_leveldata[lev]->conv_temperature)); + } + } + return r; +} + Vector incflo::get_divtau_old () noexcept { Vector r; @@ -212,6 +272,30 @@ Vector incflo::get_laps_new () noexcept return r; } +Vector incflo::get_laps_tem_old () noexcept +{ + Vector r; + if (m_use_temperature) { + r.reserve(finest_level+1); + for (int lev = 0; lev <= finest_level; ++lev) { + r.push_back(&(m_leveldata[lev]->laps_tem_o)); + } + } + return r; +} + +Vector incflo::get_laps_tem_new () noexcept +{ + Vector r; + if (m_use_temperature) { + r.reserve(finest_level+1); + for (int lev = 0; lev <= finest_level; ++lev) { + r.push_back(&(m_leveldata[lev]->laps_tem)); + } + } + return r; +} + Vector incflo::get_velocity_old_const () const noexcept { Vector r; @@ -282,6 +366,30 @@ Vector incflo::get_tracer_new_const () const noexcept return r; } +Vector incflo::get_temperature_old_const () const noexcept +{ + Vector r; + if (m_use_temperature) { + r.reserve(finest_level+1); + for (int lev = 0; lev <= finest_level; ++lev) { + r.push_back(&(m_leveldata[lev]->temperature_o)); + } + } + return r; +} + +Vector incflo::get_temperature_new_const () const noexcept +{ + Vector r; + if (m_use_temperature) { + r.reserve(finest_level+1); + for (int lev = 0; lev <= finest_level; ++lev) { + r.push_back(&(m_leveldata[lev]->temperature)); + } + } + return r; +} + void incflo::copy_from_new_to_old_velocity (IntVect const& ng) { for (int lev = 0; lev <= finest_level; ++lev) { @@ -363,3 +471,19 @@ void incflo::copy_from_old_to_new_tracer (int lev, IntVect const& ng) m_leveldata[lev]->tracer_o, 0, 0, m_ntrac, ng); } } + +void incflo::copy_from_new_to_old_temperature (IntVect const& ng) +{ + for (int lev = 0; lev <= finest_level; ++lev) { + MultiFab::Copy(m_leveldata[lev]->temperature_o, + m_leveldata[lev]->temperature, 0, 0, 1, ng); + } +} + +void incflo::copy_from_old_to_new_temperature (IntVect const& ng) +{ + for (int lev = 0; lev <= finest_level; ++lev) { + MultiFab::Copy(m_leveldata[lev]->temperature, + m_leveldata[lev]->temperature_o, 0, 0, 1, ng); + } +} diff --git a/src/prob/prob_bc.H b/src/prob/prob_bc.H index f1cf47ad..8f60e800 100644 --- a/src/prob/prob_bc.H +++ b/src/prob/prob_bc.H @@ -502,6 +502,115 @@ struct IncfloTracFill } }; +struct IncfloTempFill +{ + int probtype; + amrex::GpuArray bcv_tem; + amrex::GpuArray, AMREX_SPACEDIM*2> bcv_vel; + + AMREX_GPU_HOST + constexpr IncfloTempFill (int a_probtype, + amrex::GpuArray const& a_bcv_tem, + amrex::GpuArray,AMREX_SPACEDIM*2> const& a_bcv_vel) + : probtype(a_probtype), bcv_tem(a_bcv_tem), bcv_vel(a_bcv_vel) {} + + AMREX_GPU_DEVICE + void operator() (const amrex::IntVect& iv, amrex::Array4 const& temperature, + const int /*dcomp*/, const int /*numcomp*/, + amrex::GeometryData const& geom, const amrex::Real /*time*/, + const amrex::BCRec* bcr, const int bcomp, + const int /*orig_comp*/) const + { + using namespace amrex; + + // do something for external Dirichlet (amrex::BCType::ext_dir) + const int i = iv[0]; + const int j = iv[1]; +#if (AMREX_SPACEDIM == 3) + const int k = iv[2]; +#else + const int k = 0; +#endif + + const Box& domain_box = geom.Domain(); + + const BCRec& bc = bcr[bcomp]; + + if (1101 == probtype && i < domain_box.smallEnd(0)) + { + int direction = 1; + int half_num_cells = domain_box.length(direction) / 2; + if (j > half_num_cells) { + // Here we take the inflow BC specified in inputs file + temperature(i,j,k) = bcv_tem[amrex::Orientation(amrex::Direction::x,amrex::Orientation::low)]; + } + } + else if (1101 == probtype && i > domain_box.bigEnd(0)) + { + int direction = 1; + int half_num_cells = domain_box.length(direction) / 2; + if (j <= half_num_cells) { + temperature(i,j,k) = bcv_tem[amrex::Orientation(amrex::Direction::x,amrex::Orientation::high)]; + } + } +#if (AMREX_SPACEDIM == 3) + else if (1102 == probtype && j > domain_box.bigEnd(1)) + { + int direction = 2; + int half_num_cells = domain_box.length(direction) / 2; + if (k <= half_num_cells) { + temperature(i,j,k) = bcv_tem[amrex::Orientation(amrex::Direction::y,amrex::Orientation::high)]; + } + } +#endif + else if ( (i < domain_box.smallEnd(0)) && + ( (bc.lo(0) == amrex::BCType::ext_dir) || + (bc.lo(0) == amrex::BCType::direction_dependent && + bcv_vel[amrex::Orientation(amrex::Direction::x,amrex::Orientation::low)][0] >= 0.) ) ) + { + temperature(i,j,k) = bcv_tem[amrex::Orientation(amrex::Direction::x,amrex::Orientation::low)]; + } + else if ( (i > domain_box.bigEnd(0)) && + ( (bc.hi(0) == amrex::BCType::ext_dir) || + (bc.hi(0) == amrex::BCType::direction_dependent && + bcv_vel[amrex::Orientation(amrex::Direction::x,amrex::Orientation::high)][0] <= 0.) ) ) + { + temperature(i,j,k) = bcv_tem[amrex::Orientation(amrex::Direction::x,amrex::Orientation::high)]; + } + + if ( (j < domain_box.smallEnd(1)) && + ( (bc.lo(1) == amrex::BCType::ext_dir) || + (bc.lo(1) == amrex::BCType::direction_dependent && + bcv_vel[amrex::Orientation(amrex::Direction::y,amrex::Orientation::low)][1] >= 0.) ) ) + { + temperature(i,j,k) = bcv_tem[amrex::Orientation(amrex::Direction::y,amrex::Orientation::low)]; + } + else if ( (j > domain_box.bigEnd(1)) && + ( (bc.hi(1) == amrex::BCType::ext_dir) || + (bc.hi(1) == amrex::BCType::direction_dependent && + bcv_vel[amrex::Orientation(amrex::Direction::y,amrex::Orientation::high)][1] <= 0.) ) ) + { + temperature(i,j,k) = bcv_tem[amrex::Orientation(amrex::Direction::y,amrex::Orientation::high)]; + } +#if (AMREX_SPACEDIM == 3) + if ( (k < domain_box.smallEnd(2)) && + ( (bc.lo(2) == amrex::BCType::ext_dir) || + (bc.lo(2) == amrex::BCType::direction_dependent && + bcv_vel[amrex::Orientation(amrex::Direction::z,amrex::Orientation::low)][2] >= 0.) ) ) + { + temperature(i,j,k) = bcv_tem[amrex::Orientation(amrex::Direction::z,amrex::Orientation::low)]; + } + else if ( (k > domain_box.bigEnd(2)) && + ( (bc.hi(2) == amrex::BCType::ext_dir) || + (bc.hi(2) == amrex::BCType::direction_dependent && + bcv_vel[amrex::Orientation(amrex::Direction::z,amrex::Orientation::high)][2] <= 0.) ) ) + { + temperature(i,j,k) = bcv_tem[amrex::Orientation(amrex::Direction::z,amrex::Orientation::high)]; + } +#endif + } +}; + struct IncfloForFill { int probtype; diff --git a/src/prob/prob_init_fluid.cpp b/src/prob/prob_init_fluid.cpp index e92288f2..8f3b20e2 100644 --- a/src/prob/prob_init_fluid.cpp +++ b/src/prob/prob_init_fluid.cpp @@ -29,6 +29,10 @@ void incflo::prob_init_fluid (int lev) ld.tracer.setVal(m_ic_t[comp], comp, 1); } + if (m_use_temperature) { + ld.temperature.setVal(m_ic_tem); + } + for (MFIter mfi(ld.density); mfi.isValid(); ++mfi) { const Box& vbx = mfi.validbox(); diff --git a/src/projection/incflo_apply_nodal_projection.cpp b/src/projection/incflo_apply_nodal_projection.cpp index dad0239f..fafd193f 100644 --- a/src/projection/incflo_apply_nodal_projection.cpp +++ b/src/projection/incflo_apply_nodal_projection.cpp @@ -130,8 +130,6 @@ void incflo::ApplyNodalProjection (Vector const& density, #ifdef AMREX_USE_EB if (m_eb_flow.enabled) { set_eb_velocity(lev, time, *get_velocity_eb()[lev], 1); - set_eb_density(lev, time, *get_density_eb()[lev], 1); - set_eb_tracer(lev, time, *get_tracer_eb()[lev], 1); } #endif vel[lev]->setBndry(0.0); diff --git a/src/rheology/incflo_rheology.cpp b/src/rheology/incflo_rheology.cpp index ca2d4b23..88d3dea6 100644 --- a/src/rheology/incflo_rheology.cpp +++ b/src/rheology/incflo_rheology.cpp @@ -134,7 +134,7 @@ void incflo::compute_viscosity_at_level (int /*lev*/, } } -void incflo::compute_tracer_diff_coeff (Vector const& tra_eta, int nghost) +void incflo::compute_tracer_diff_coeff (Vector const& tra_eta, int nghost) const { for (auto *mf : tra_eta) { for (int n = 0; n < m_ntrac; ++n) { @@ -142,3 +142,11 @@ void incflo::compute_tracer_diff_coeff (Vector const& tra_eta, int ng } } } + +void incflo::compute_temperature_diff_coeff (Real /*time*/, Vector const& tem_eta) const +{ + for (auto *mf : tem_eta) { // loop over levels + mf->setVal(m_mu_T); + } +} + diff --git a/src/setup/incflo_arrays.cpp b/src/setup/incflo_arrays.cpp index dc035226..77af5aae 100644 --- a/src/setup/incflo_arrays.cpp +++ b/src/setup/incflo_arrays.cpp @@ -28,33 +28,58 @@ incflo::LevelData::LevelData (amrex::BoxArray const& ba, } else { p_nd.define(convert(ba,IntVect::TheNodeVector()), dm, 1, 0, MFInfo(), fact); } + if (my_incflo->m_use_temperature) { + temperature.define (ba, dm, 1, my_incflo->nghost_state(), MFInfo(), fact); + temperature_o.define (ba, dm, 1, my_incflo->nghost_state(), MFInfo(), fact); + + conv_temperature_o.define(ba, dm, 1, 0, MFInfo(), fact); + } #ifdef AMREX_USE_EB if (my_incflo->hasEBFlow()) { velocity_eb.define(ba, dm, AMREX_SPACEDIM, my_incflo->nghost_state(), MFInfo(), fact); density_eb.define (ba, dm, 1 , my_incflo->nghost_state(), MFInfo(), fact); + } + // Allow for Dirichlet BC on EB even if there's no flow through the EB + if (my_incflo->m_advect_tracer && !my_incflo->m_eb_flow.tracer.empty()) { tracer_eb.define (ba, dm, my_incflo->m_ntrac, my_incflo->nghost_state(), MFInfo(), fact); } + if (my_incflo->m_use_temperature && !my_incflo->m_eb_flow.temperature.empty()) { + temperature_eb.define(ba, dm, 1, my_incflo->nghost_state(), MFInfo(), fact); + } #endif if (my_incflo->m_advection_type != "MOL") { divtau_o.define(ba, dm, AMREX_SPACEDIM, 0, MFInfo(), fact); if (my_incflo->m_advect_tracer) { laps_o.define(ba, dm, my_incflo->m_ntrac, 0, MFInfo(), fact); } + if (my_incflo->m_use_temperature) { + laps_tem_o.define(ba, dm, 1, 0, MFInfo(), fact); + } } else { conv_velocity.define(ba, dm, AMREX_SPACEDIM , 0, MFInfo(), fact); conv_density.define (ba, dm, 1 , 0, MFInfo(), fact); conv_tracer.define (ba, dm, my_incflo->m_ntrac, 0, MFInfo(), fact); + if (my_incflo->m_use_temperature) { + conv_temperature.define(ba, dm, 1, 0, MFInfo(), fact); + } + bool implicit_diffusion = my_incflo->m_diff_type == DiffusionType::Implicit; if (!implicit_diffusion || my_incflo->use_tensor_correction) { divtau.define (ba, dm, AMREX_SPACEDIM, 0, MFInfo(), fact); divtau_o.define(ba, dm, AMREX_SPACEDIM, 0, MFInfo(), fact); } - if (!implicit_diffusion && my_incflo->m_advect_tracer) + if (!implicit_diffusion) { - laps.define (ba, dm, my_incflo->m_ntrac, 0, MFInfo(), fact); - laps_o.define(ba, dm, my_incflo->m_ntrac, 0, MFInfo(), fact); + if ( my_incflo->m_advect_tracer) { + laps.define (ba, dm, my_incflo->m_ntrac, 0, MFInfo(), fact); + laps_o.define(ba, dm, my_incflo->m_ntrac, 0, MFInfo(), fact); + } + if (my_incflo->m_use_temperature) { + laps_tem.define (ba, dm, 1, 0, MFInfo(), fact); + laps_tem_o.define(ba, dm, 1, 0, MFInfo(), fact); + } } } } diff --git a/src/setup/init.cpp b/src/setup/init.cpp index bc9da8da..75b6b613 100644 --- a/src/setup/init.cpp +++ b/src/setup/init.cpp @@ -141,6 +141,7 @@ void incflo::ReadParameters () if ( !pp.queryarr("ic_t", m_ic_t, 0, m_ntrac) ) { m_ic_t.resize(m_ntrac, 0.); } + pp.query("ic_tem", m_ic_tem); // Viscosity (if constant) pp.query("mu", m_mu); @@ -157,6 +158,22 @@ void incflo::ReadParameters () for (int i = 0; i < m_ntrac; i++) { amrex::Print() << "Tracer diffusion coeff: " << i << ":" << m_mu_s[i] << std::endl; } + + pp.query("use_temperature", m_use_temperature); + // Checks for things not yet implemented/checked + if (m_use_temperature && m_advection_type == "MOL") { + // temperature equation not added to the corrector + amrex::Abort("Temperature equation not yet implemented with MOL option"); + } +#ifdef AMREX_USE_EB + std::string geom_type = "all_regular"; + pp.query("geometry", geom_type); +#endif + + // Thermal diffusivity (if constant) + pp.query("mu_T", m_mu_T); + pp.query("cp", m_cp); + } // end prefix incflo ReadIOParameters(); @@ -182,8 +199,13 @@ void incflo::ReadParameters () pp_eb_flow.query("density", m_eb_flow.density); - m_eb_flow.tracer.resize(m_ntrac, 0.0); - pp_eb_flow.queryarr("tracer", m_eb_flow.tracer, 0, m_ntrac); + if (m_advect_tracer) { + pp_eb_flow.queryarr("tracer", m_eb_flow.tracer); + } + + if (m_use_temperature) { + pp_eb_flow.queryarr("temperature", m_eb_flow.temperature); + } if (pp_eb_flow.contains("vel_mag")) { m_eb_flow.enabled = true; @@ -202,6 +224,13 @@ void incflo::ReadParameters () pp_eb_flow.query("normal_tol", tol_deg); m_eb_flow.normal_tol = tol_deg*M_PI/amrex::Real(180.); } + + if (m_advect_tracer && m_eb_flow.enabled && m_eb_flow.tracer.empty()) { + Abort("Must specify tracer EB value for flow through EB"); + } + if (m_use_temperature && m_eb_flow.enabled && m_eb_flow.temperature.empty()) { + Abort("Must specify temperature EB value for flow through EB"); + } } // end prefix eb_flow #endif @@ -256,6 +285,11 @@ void incflo::ReadIOParameters() #endif } + if (m_use_temperature) { + // Add temperature to the default list + m_plotVars.push_back("temperature"); + } + // Helper function to update m_plotVars according to m_plt_* flags auto update_plotVars = [this] (std::string const& a_name, bool plot_flag) { for (int n = 0; n < m_plotVars.size(); n++) { @@ -406,6 +440,7 @@ void incflo::InitialIterations () copy_from_new_to_old_velocity(); copy_from_new_to_old_density(); copy_from_new_to_old_tracer(); + copy_from_new_to_old_temperature(); int initialisation = 1; bool explicit_diffusion = (m_diff_type == DiffusionType::Explicit); @@ -428,6 +463,9 @@ void incflo::InitialIterations () if (m_advect_tracer) { fillpatch_tracer(lev, m_t_old[lev], m_leveldata[lev]->tracer_o, ng); } + if (m_use_temperature) { + fillpatch_temperature(lev, m_t_old[lev], m_leveldata[lev]->temperature_o, ng); + } } for (int iter = 0; iter < m_initial_iterations; ++iter) @@ -439,6 +477,7 @@ void incflo::InitialIterations () copy_from_old_to_new_velocity(); copy_from_old_to_new_density(); copy_from_old_to_new_tracer(); + copy_from_old_to_new_temperature(); } // Reset dt to get initial step as specified, otherwise we can see increase to dt @@ -606,6 +645,12 @@ incflo::InitialRedistribution () MultiFab::Copy(ld.tracer_o, ld.tracer, 0, 0, m_ntrac, ld.tracer.nGrow()); fillpatch_tracer(lev, m_t_new[lev], ld.tracer_o, 3); } + if (m_use_temperature) + { + ld.temperature.FillBoundary(geom[lev].periodicity()); + MultiFab::Copy(ld.temperature_o, ld.temperature, 0, 0, 1, ld.temperature.nGrow()); + fillpatch_temperature(lev, m_t_new[lev], ld.temperature_o, 3); + } for (MFIter mfi(ld.density,TilingIfNotGPU()); mfi.isValid(); ++mfi) { @@ -656,6 +701,16 @@ incflo::InitialRedistribution () AMREX_D_DECL(fcx, fcy, fcz), ccc, bc_tra, geom[lev], m_redistribution_type); } + if (m_use_temperature) + { + ncomp = 1; + auto const& bc_T = get_temperature_bcrec_device_ptr(); + ApplyInitialRedistribution( bx,ncomp, + ld.temperature.array(mfi), ld.temperature_o.array(mfi), + flag, AMREX_D_DECL(apx, apy, apz), vfrac, + AMREX_D_DECL(fcx, fcy, fcz), ccc, + bc_T, geom[lev], m_redistribution_type); + } } } @@ -663,6 +718,9 @@ incflo::InitialRedistribution () ld.velocity.FillBoundary(geom[lev].periodicity()); ld.density.FillBoundary(geom[lev].periodicity()); ld.tracer.FillBoundary(geom[lev].periodicity()); + if (m_use_temperature) { + ld.temperature.FillBoundary(geom[lev].periodicity()); + } } } } diff --git a/src/utilities/io.cpp b/src/utilities/io.cpp index 7883a93b..a1e9f186 100644 --- a/src/utilities/io.cpp +++ b/src/utilities/io.cpp @@ -93,6 +93,11 @@ void incflo::WriteCheckPointFile() const amrex::MultiFabFileFullPrefix(lev, checkpointname, level_prefix, "tracer")); } + if (m_use_temperature) { + VisMF::Write(m_leveldata[lev]->temperature, + amrex::MultiFabFileFullPrefix(lev, checkpointname, level_prefix, "temperature")); + } + VisMF::Write(m_leveldata[lev]->gp, amrex::MultiFabFileFullPrefix(lev, checkpointname, level_prefix, "gradp")); @@ -224,6 +229,11 @@ void incflo::ReadCheckpointFile() amrex::MultiFabFileFullPrefix(lev, m_restart_file, level_prefix, "tracer")); } + if (m_use_temperature) { + VisMF::Read(m_leveldata[lev]->temperature, + amrex::MultiFabFileFullPrefix(lev, m_restart_file, level_prefix, "temperature")); + } + VisMF::Read(m_leveldata[lev]->gp, amrex::MultiFabFileFullPrefix(lev, m_restart_file, level_prefix, "gradp")); @@ -384,6 +394,8 @@ void incflo::WritePlotVariables(Vector vars, const std::string& plo fillpatch_velocity(lev, m_cur_time, m_leveldata[lev]->velocity, ng); fillpatch_density(lev, m_cur_time, m_leveldata[lev]->density, ng); fillpatch_tracer(lev, m_cur_time, m_leveldata[lev]->tracer, ng); + // Whether temperature fillpatch is needed depends on form of forcing term + // fillpatch_temperature(lev, m_cur_time, m_leveldata[lev]->temperature, ng); } break; } @@ -493,6 +505,13 @@ void incflo::WritePlotVariables(Vector vars, const std::string& plo } icomp += m_ntrac; } + else if (vars[n] == "temperature") { + for (int lev = 0; lev <= finest_level; ++lev) { + MultiFab::Copy(mf[lev], m_leveldata[lev]->temperature, 0, icomp, 1, 0); + } + pltscaVarsName.push_back("temperature"); + ++icomp; + } else if (vars[n] == "p") { if (m_use_cc_proj) { for (int lev = 0; lev <= finest_level; ++lev) { diff --git a/test_3d/benchmark.channel_spherecube b/test_3d/benchmark.channel_spherecube index 26025ed6..1b8feeb6 100644 --- a/test_3d/benchmark.channel_spherecube +++ b/test_3d/benchmark.channel_spherecube @@ -29,6 +29,9 @@ incflo.mu = 1000.0 # Dynamic viscosity coefficient incflo.tau_0 = 340700.0 # Yield stress incflo.papa_reg = 1.0e-1 # Regularisation parameter +incflo.use_temperature = true +incflo.mu_T = 1900.0 + #¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# # ADAPTIVE MESH REFINEMENT # #.......................................# @@ -49,12 +52,16 @@ xhi.type = "nsw" xhi.velocity = 0. 0. 1. ylo.type = "nsw" ylo.velocity = 0. 0. 1. +ylo.temperature = 3. yhi.type = "nsw" yhi.velocity = 0. 0. 1. # Add sphere incflo.geometry = "spherecube" +# Dirichlet BC on the EB for temperature +eb_flow.temperature = 1. + #¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# # INITIAL CONDITIONS # #.......................................# @@ -62,6 +69,7 @@ incflo.ic_u = 0.0 # incflo.ic_v = 0.0 # incflo.ic_w = 1.0 # incflo.ic_p = 0.0 # +incflo.ic_tem = 2.0 # #¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# # VERBOSITY #