diff --git a/src/boundary_conditions/README.org b/src/boundary_conditions/README.org index 16147d817..6008a5b95 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 8cdbb2624..9d7f56348 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,8 @@ 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); +// FIXME? do i want something reasonable here, or force users to set inputs? + m_bc_temperature[ori] = 1.0; ParmParse pp(bcid); std::string bc_type_in = "null"; @@ -41,6 +44,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 +60,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 +77,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 +85,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 +104,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 +135,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], 0, 1); // Set mathematical BCs AMREX_D_TERM(m_bcrec_velocity[0].set(ori, BCType::ext_dir);, @@ -137,6 +147,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 +166,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], 0, 1); // Tangential directions have hoextrap AMREX_D_TERM(m_bcrec_velocity[0].set(ori, BCType::hoextrap);, @@ -175,6 +191,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 +229,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 +237,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 +254,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"); } @@ -285,7 +309,8 @@ void incflo::init_bcs () #else std::memcpy #endif - (m_bcrec_density_d.data(), m_bcrec_density.data(), sizeof(BCRec)); + (m_bcrec_density_d.data(), m_bcrec_density.data(), sizeof(BCRec)); + } if (m_ntrac > 0) { @@ -298,6 +323,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/convection/incflo_compute_advection_term.cpp b/src/convection/incflo_compute_advection_term.cpp index 07a64715b..f46653650 100644 --- a/src/convection/incflo_compute_advection_term.cpp +++ b/src/convection/incflo_compute_advection_term.cpp @@ -46,15 +46,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; @@ -191,6 +194,17 @@ incflo::compute_convective_term (Vector const& conv_u, fillpatch_force(m_cur_time, tra_forces, nghost_force()); } + if (m_use_temperature) + { + compute_T_forces(tem_forces, get_density_old_const()); + if (m_godunov_include_diff_in_forcing) + for (int lev = 0; lev <= finest_level; ++lev) + // FIXME? Saxpy? with 1/rhocp - could be as low as point function or by level? + MultiFab::Add(*tem_forces[lev], m_leveldata[lev]->laps_T_o, 0, 0, 1, 0); + 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) @@ -300,7 +314,7 @@ 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); @@ -354,6 +368,13 @@ incflo::compute_convective_term (Vector const& conv_u, tracBC_MF = make_BC_MF(lev, m_bcrec_tracer_d, "tracer"); } } + // FIXME? do all the scalars really need individual BC MFs? + // Need to think about general mixed BC vs inflow-outflow only? + 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 @@ -551,6 +572,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 + + // FIXME- check this!!! + face_comp += m_ntrac; + 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(), + get_temperature_iconserv_device_ptr(), +#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 diff --git a/src/diffusion/DiffusionScalarOp.H b/src/diffusion/DiffusionScalarOp.H index df74be14c..fbe1020b6 100644 --- a/src/diffusion/DiffusionScalarOp.H +++ b/src/diffusion/DiffusionScalarOp.H @@ -15,9 +15,10 @@ class DiffusionScalarOp public: DiffusionScalarOp (incflo* a_incflo); - void diffuse_scalar (amrex::Vector const& tracer, + void diffuse_scalar (amrex::Vector const& a_scalar, amrex::Vector const& density, amrex::Vector const& eta, + amrex::Vector const& iconserv, amrex::Real dt); void diffuse_vel_components ( diff --git a/src/diffusion/DiffusionScalarOp.cpp b/src/diffusion/DiffusionScalarOp.cpp index 5ccd8af6f..083a665ee 100644 --- a/src/diffusion/DiffusionScalarOp.cpp +++ b/src/diffusion/DiffusionScalarOp.cpp @@ -147,30 +147,31 @@ DiffusionScalarOp::readParameters () } void -DiffusionScalarOp::diffuse_scalar (Vector const& tracer, +DiffusionScalarOp::diffuse_scalar (Vector const& a_scalar, Vector const& density, Vector const& eta, + amrex::Vector const& iconserv, Real dt) { // // Solves // alpha a - beta div ( b grad ) // - // 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) + // So for conservative: d(rho sca) / dt - div mu grad sca = -div(U rho sca) + rho H --> + // ( rho - dt div mu grad ) sca^(n+1) = rho sca^(*,n+1) // alpha: 1 // a: rho // beta: dt // b: mu - // RHS: density * tracer + // RHS: density * a_scalar // - // 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) + // So for non- conservative: d sca / dt - div mu grad sca = -U dot grad sca + H --> + // ( 1 - dt div mu grad ) sca^(n+1) = sca^(*,n+1) // alpha: 1 // a: 1 // beta: dt // b: mu - // RHS: tracer + // RHS: a_scalar if (m_verbose > 0) { amrex::Print() << "Diffusing scalars one at a time ..." << std::endl; @@ -181,10 +182,9 @@ 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) { @@ -210,7 +210,7 @@ DiffusionScalarOp::diffuse_scalar (Vector const& tracer, } } - 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) @@ -255,10 +255,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); + 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 @@ -267,11 +267,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& sca_a = a_scalar[lev]->const_array(mfi,comp); Array4 const& rho_a = density[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) = rho_a(i,j,k) * sca_a(i,j,k); }); } } @@ -486,6 +486,7 @@ void DiffusionScalarOp::compute_laps (Vector const& a_laps, 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) @@ -495,7 +496,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); } @@ -508,7 +509,7 @@ void DiffusionScalarOp::compute_laps (Vector const& a_laps, // 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) { int eta_comp = comp; if ( m_incflo->m_has_mixedBC && comp>0 ){ @@ -545,9 +546,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)); + // 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); @@ -559,7 +564,7 @@ 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) { int eta_comp = comp; diff --git a/src/incflo.H b/src/incflo.H index 08670e23a..dbd771da2 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_T_forces (amrex::Vector const& T_forces, + amrex::Vector const& density); void compute_vel_forces (amrex::Vector const& vel_forces, amrex::Vector const& velocity, amrex::Vector const& density, @@ -139,15 +141,18 @@ 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_MAC_projected_velocities ( @@ -165,6 +170,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& T_eta, + amrex::Vector& T_forces); void update_velocity (StepType step_type, amrex::Vector& vel_eta, amrex::Vector& vel_forces); @@ -285,6 +292,7 @@ public: amrex::Geometry& lev_geom, amrex::Real time, int nghost); void compute_tracer_diff_coeff (amrex::Vector const& tra_eta, int nghost); + void compute_temperature_diff_coeff (amrex::Vector const& eta, int nghost); #ifdef AMREX_USE_EB /////////////////////////////////////////////////////////////////////////// @@ -612,6 +620,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 @@ -627,13 +638,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_T; + amrex::MultiFab laps_T_o; }; + bool m_use_temperature = false; + amrex::Vector > m_leveldata; amrex::Vector > > m_factory; @@ -651,6 +668,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 ... @@ -663,6 +681,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; @@ -766,6 +786,8 @@ private: 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; @@ -775,10 +797,14 @@ private: 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_T_old () noexcept; + amrex::Vector get_laps_T_new () noexcept; // [[nodiscard]] amrex::Vector get_velocity_old_const () const noexcept; [[nodiscard]] amrex::Vector get_velocity_new_const () const noexcept; @@ -790,8 +816,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_T_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; } @@ -807,6 +836,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 { @@ -815,6 +845,10 @@ 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_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(); } @@ -835,17 +869,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}); @@ -853,6 +890,8 @@ 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_new_to_old_temperature (int lev, 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}); @@ -860,6 +899,8 @@ 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 copy_from_old_to_new_temperature (int lev, amrex::IntVect const& ng = amrex::IntVect{0}); void Advance (); bool writeNow (); diff --git a/src/incflo_apply_predictor.cpp b/src/incflo_apply_predictor.cpp index 8039574cc..1cf2be9b5 100644 --- a/src/incflo_apply_predictor.cpp +++ b/src/incflo_apply_predictor.cpp @@ -97,10 +97,10 @@ 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, T_forces; - Vector vel_eta, tra_eta; + Vector vel_eta, tra_eta, T_eta; // ************************************************************************************* // Allocate space for the forcing terms @@ -108,15 +108,18 @@ void incflo::ApplyPredictor (bool incremental_projection) 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, 1, MFInfo(), Factory(lev)); + if (m_advect_tracer) { tra_forces.emplace_back(grids[lev], dmap[lev], m_ntrac, nghost_force(), 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) { + T_forces.emplace_back(grids[lev], dmap[lev], 1, nghost_force(), + MFInfo(), Factory(lev)); + T_eta.emplace_back(grids[lev], dmap[lev], 1, 1, MFInfo(), Factory(lev)); + } } // ************************************************************************************* @@ -141,11 +144,19 @@ 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),tra_eta[0].nComp()); if (need_divtau()) { compute_laps(get_laps_old(), get_tracer_old_const(), GetVecOfConstPtrs(tra_eta)); } } + if (m_use_temperature) + { + compute_temperature_diff_coeff(GetVecOfPtrs(T_eta),T_eta[0].nComp()); + if (need_divtau()) { + compute_laps(get_laps_T_old(), get_temperature_old_const(), GetVecOfConstPtrs(T_eta)); + // Mulitply by rho cp + } + } // ********************************************************************************************** // Compute the forcing terms diff --git a/src/setup/incflo_arrays.cpp b/src/setup/incflo_arrays.cpp index 031f295a4..97b7c5b49 100644 --- a/src/setup/incflo_arrays.cpp +++ b/src/setup/incflo_arrays.cpp @@ -28,26 +28,47 @@ incflo::LevelData::LevelData (amrex::BoxArray const& ba, conv_density_o (ba, dm, 1 , 0, MFInfo(), fact), conv_tracer_o (ba, dm, my_incflo->m_ntrac, 0, MFInfo(), fact) { + if (my_incflo->m_use_temperature) { + temperature (ba, dm, 1, my_incflo->nghost_state(), MFInfo(), fact); + temperature_eb (ba, dm, 1, my_incflo->nghost_state(), MFInfo(), fact); + temperature_o (ba, dm, 1, my_incflo->nghost_state(), MFInfo(), fact); + + conv_temperature_o (ba, dm, 1, 0, MFInfo(), fact); + } + 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_T_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 (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_T.define (ba, dm, 1, 0, MFInfo(), fact); + laps_T_o.define(ba, dm, 1, 0, MFInfo(), fact); + } } } } diff --git a/src/setup/init.cpp b/src/setup/init.cpp index 33a62b0da..c3ce88477 100644 --- a/src/setup/init.cpp +++ b/src/setup/init.cpp @@ -154,6 +154,19 @@ 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 + if (m_use_temperature) { + amrex::Abort("Temperature equation not yet tested with EB"); + // Maybe will want to disallow EB flow with T at first... + } +#endif } // end prefix incflo ReadIOParameters(); @@ -502,6 +515,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); + } } } @@ -509,6 +532,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/test_2d/GNUmakefile b/test_2d/GNUmakefile index 67768a03a..7633efd25 100644 --- a/test_2d/GNUmakefile +++ b/test_2d/GNUmakefile @@ -1,4 +1,4 @@ -AMREX_HOME ?= ../../amrex +AMREX_HOME ?= ../../amrex_fork AMREX_HYDRO_HOME ?= ../../AMReX-Hydro TOP = .. @@ -19,7 +19,7 @@ COMP = gnu USE_MPI = TRUE USE_OMP = FALSE -USE_CUDA = FALSE +USE_CUDA = TRUE USE_HYPRE = FALSE diff --git a/test_3d/GNUmakefile b/test_3d/GNUmakefile index df7e424ea..3aafc56bb 100644 --- a/test_3d/GNUmakefile +++ b/test_3d/GNUmakefile @@ -1,4 +1,4 @@ -AMREX_HOME ?= ../../amrex +AMREX_HOME ?= ../../amrex_fork AMREX_HYDRO_HOME ?= ../../AMReX-Hydro TOP = .. diff --git a/test_no_eb_2d/GNUmakefile b/test_no_eb_2d/GNUmakefile index ef9332d27..b77936f03 100644 --- a/test_no_eb_2d/GNUmakefile +++ b/test_no_eb_2d/GNUmakefile @@ -1,4 +1,4 @@ -AMREX_HOME ?= ../../amrex +AMREX_HOME ?= ../../amrex_fork AMREX_HYDRO_HOME ?= ../../AMReX-Hydro TOP = .. @@ -10,12 +10,12 @@ PRECISION = DOUBLE TINY_PROFILE = FALSE -DEBUG = TRUE DEBUG = FALSE +DEBUG = TRUE COMP = gnu -USE_MPI = TRUE +USE_MPI = FALSE USE_OMP = FALSE USE_CUDA = FALSE diff --git a/test_no_eb_3d/GNUmakefile b/test_no_eb_3d/GNUmakefile index 0299d3eba..db969760f 100644 --- a/test_no_eb_3d/GNUmakefile +++ b/test_no_eb_3d/GNUmakefile @@ -1,4 +1,4 @@ -AMREX_HOME ?= ../../amrex +AMREX_HOME ?= ../../amrex_fork AMREX_HYDRO_HOME ?= ../../AMReX-Hydro TOP = .. @@ -10,12 +10,12 @@ PRECISION = DOUBLE TINY_PROFILE = FALSE -DEBUG = TRUE DEBUG = FALSE +DEBUG = TRUE COMP = gnu -USE_MPI = TRUE +USE_MPI = FALSE USE_OMP = FALSE USE_CUDA = FALSE