diff --git a/amr-wind/CFDSim.H b/amr-wind/CFDSim.H index 6356e5e96b..6909367965 100644 --- a/amr-wind/CFDSim.H +++ b/amr-wind/CFDSim.H @@ -117,6 +117,8 @@ public: bool has_mesh_mapping() const { return m_mesh_mapping; } + bool is_anelastic() const; + private: amrex::AmrCore& m_mesh; diff --git a/amr-wind/CFDSim.cpp b/amr-wind/CFDSim.cpp index b02e05f78a..f1b786c8d9 100644 --- a/amr-wind/CFDSim.cpp +++ b/amr-wind/CFDSim.cpp @@ -4,6 +4,7 @@ #include "amr-wind/utilities/PostProcessing.H" #include "amr-wind/overset/OversetManager.H" #include "amr-wind/core/ExtSolver.H" +#include "amr-wind/wind_energy/ABL.H" #include "AMReX_ParmParse.H" @@ -72,4 +73,12 @@ void CFDSim::activate_mesh_map() } } +bool CFDSim::is_anelastic() const +{ + if (m_physics_mgr.contains("ABL")) { + const auto& abl = m_physics_mgr.template get(); + return abl.anelastic().is_anelastic(); + } + return false; +} } // namespace amr_wind diff --git a/amr-wind/equation_systems/AdvOp_Godunov.H b/amr-wind/equation_systems/AdvOp_Godunov.H index 2a93b8e6ea..c27a80fc21 100644 --- a/amr-wind/equation_systems/AdvOp_Godunov.H +++ b/amr-wind/equation_systems/AdvOp_Godunov.H @@ -28,7 +28,8 @@ struct AdvectionOp< PDEFields& fields_in, bool /* has_overset */, bool /* variable density */, - bool /* mesh mapping */) + bool /* mesh mapping */, + bool /* is_anelastic */) : fields(fields_in) , density(fields_in.repo.get_field("density")) , u_mac(fields_in.repo.get_field("u_mac")) diff --git a/amr-wind/equation_systems/AdvOp_MOL.H b/amr-wind/equation_systems/AdvOp_MOL.H index 689df94721..8803ab624b 100644 --- a/amr-wind/equation_systems/AdvOp_MOL.H +++ b/amr-wind/equation_systems/AdvOp_MOL.H @@ -27,7 +27,8 @@ struct AdvectionOp< PDEFields& fields_in, bool /* has_overset */, bool /* variable density */, - bool /* mesh mapping */) + bool /* mesh mapping */, + bool /* is_anelastic */) : fields(fields_in) , density(fields_in.repo.get_field("density")) , u_mac(fields_in.repo.get_field("u_mac")) diff --git a/amr-wind/equation_systems/PDE.H b/amr-wind/equation_systems/PDE.H index 57217ecec3..974c256cdd 100644 --- a/amr-wind/equation_systems/PDE.H +++ b/amr-wind/equation_systems/PDE.H @@ -63,13 +63,13 @@ public: new TurbulenceOp(m_sim.turbulence_model(), m_fields)); } - bool variable_density = + const bool variable_density = (!m_sim.pde_manager().constant_density() || m_sim.physics_manager().contains("MultiPhase")); m_adv_op.reset(new AdvectionOp( m_fields, m_sim.has_overset(), variable_density, - m_sim.has_mesh_mapping())); + m_sim.has_mesh_mapping(), m_sim.is_anelastic())); m_src_op.init_source_terms(m_sim); // Post-solve operations should also apply after initialization @@ -86,13 +86,13 @@ public: m_fields, m_sim.has_overset(), m_sim.has_mesh_mapping())); } - bool variable_density = + const bool variable_density = (!m_sim.pde_manager().constant_density() || m_sim.physics_manager().contains("MultiPhase")); m_adv_op.reset(new AdvectionOp( m_fields, m_sim.has_overset(), variable_density, - m_sim.has_mesh_mapping())); + m_sim.has_mesh_mapping(), m_sim.is_anelastic())); // Post-solve operations should also apply after a regrid m_post_solve_op(m_time.current_time()); diff --git a/amr-wind/equation_systems/PDEBase.cpp b/amr-wind/equation_systems/PDEBase.cpp index c6dfb556d7..10775cbb61 100644 --- a/amr-wind/equation_systems/PDEBase.cpp +++ b/amr-wind/equation_systems/PDEBase.cpp @@ -24,6 +24,21 @@ PDEMgr::PDEMgr(CFDSim& sim) : m_sim(sim) pp.query("use_godunov", m_use_godunov); pp.query("constant_density", m_constant_density); + bool is_anelastic = false; + { + amrex::ParmParse pp_abl("ABL"); + pp_abl.query("anelastic", is_anelastic); + } + + if ((is_anelastic) && (m_constant_density)) { + amrex::Print() + << "WARNING: Anelastic implies variable density. Setting " + "contant_density to false. Add incflo.constant_density=false to " + "your input file to remove this warning." + << std::endl; + m_constant_density = false; + } + m_scheme = m_use_godunov ? fvm::Godunov::scheme_name() : fvm::MOL::scheme_name(); } diff --git a/amr-wind/equation_systems/icns/icns_advection.H b/amr-wind/equation_systems/icns/icns_advection.H index e618c3817a..be5218929c 100644 --- a/amr-wind/equation_systems/icns/icns_advection.H +++ b/amr-wind/equation_systems/icns/icns_advection.H @@ -25,7 +25,8 @@ public: FieldRepo& /*repo*/, bool /*has_overset*/, bool /*variable_density*/, - bool /*mesh_mapping*/); + bool /*mesh_mapping*/, + bool /*is_anelastic*/); void operator()(const FieldState fstate, const amrex::Real dt); @@ -51,6 +52,7 @@ private: bool m_need_init{true}; bool m_variable_density{false}; bool m_mesh_mapping{false}; + bool m_is_anelastic{false}; amrex::Real m_rho_0{1.0}; }; @@ -64,12 +66,18 @@ struct AdvectionOp PDEFields& fields_in, bool has_overset, bool variable_density, - bool mesh_mapping) + bool mesh_mapping, + bool is_anelastic) : fields(fields_in) , u_mac(fields_in.repo.get_field("u_mac")) , v_mac(fields_in.repo.get_field("v_mac")) , w_mac(fields_in.repo.get_field("w_mac")) - , m_macproj_op(fields.repo, has_overset, variable_density, mesh_mapping) + , m_macproj_op( + fields.repo, + has_overset, + variable_density, + mesh_mapping, + is_anelastic) { amrex::ParmParse pp("incflo"); @@ -533,14 +541,19 @@ struct AdvectionOp PDEFields& fields_in, bool has_overset, bool variable_density, - bool mesh_mapping) + bool mesh_mapping, + bool is_anelastic) : fields(fields_in) , u_mac(fields_in.repo.get_field("u_mac")) , v_mac(fields_in.repo.get_field("v_mac")) , w_mac(fields_in.repo.get_field("w_mac")) , m_mesh_mapping(mesh_mapping) , m_macproj_op( - fields.repo, has_overset, variable_density, m_mesh_mapping) + fields.repo, + has_overset, + variable_density, + m_mesh_mapping, + is_anelastic) {} void preadvect( diff --git a/amr-wind/equation_systems/icns/icns_advection.cpp b/amr-wind/equation_systems/icns/icns_advection.cpp index 25348f3066..6bffd88035 100644 --- a/amr-wind/equation_systems/icns/icns_advection.cpp +++ b/amr-wind/equation_systems/icns/icns_advection.cpp @@ -41,12 +41,17 @@ amrex::Array get_projection_bc( } // namespace MacProjOp::MacProjOp( - FieldRepo& repo, bool has_overset, bool variable_density, bool mesh_mapping) + FieldRepo& repo, + bool has_overset, + bool variable_density, + bool mesh_mapping, + bool is_anelastic) : m_repo(repo) , m_options("mac_proj") , m_has_overset(has_overset) , m_variable_density(variable_density) , m_mesh_mapping(mesh_mapping) + , m_is_anelastic(is_anelastic) { amrex::ParmParse pp("incflo"); pp.query("density", m_rho_0); @@ -136,11 +141,24 @@ void MacProjOp::operator()(const FieldState fstate, const amrex::Real dt) amrex::Real factor = m_has_overset ? 0.5 * dt : 1.0; + std::unique_ptr ref_rho_xf, ref_rho_yf, ref_rho_zf; + if (m_is_anelastic) { + ref_rho_xf = + m_repo.create_scratch_field(1, 0, amr_wind::FieldLoc::XFACE); + ref_rho_yf = + m_repo.create_scratch_field(1, 0, amr_wind::FieldLoc::YFACE); + ref_rho_zf = + m_repo.create_scratch_field(1, 0, amr_wind::FieldLoc::ZFACE); + } + amrex::Vector> ref_rho_face( + m_repo.num_active_levels()); + // TODO: remove the or in the if statement for m_has_overset // For now assume variable viscosity for overset // this can be removed once the nsolve overset // masking is implemented in cell based AMReX poisson solvers - if (m_variable_density || m_has_overset || m_mesh_mapping) { + if (m_variable_density || m_has_overset || m_mesh_mapping || + m_is_anelastic) { amrex::Vector> rho_face_const; rho_face_const.reserve(m_repo.num_active_levels()); @@ -153,6 +171,9 @@ void MacProjOp::operator()(const FieldState fstate, const amrex::Real dt) amrex::Vector> rho_face( m_repo.num_active_levels()); + const auto* ref_density = + m_is_anelastic ? &(m_repo.get_field("reference_density")) : nullptr; + for (int lev = 0; lev < m_repo.num_active_levels(); ++lev) { rho_face[lev][0] = &(*rho_xf)(lev); rho_face[lev][1] = &(*rho_yf)(lev); @@ -161,6 +182,19 @@ void MacProjOp::operator()(const FieldState fstate, const amrex::Real dt) amrex::average_cellcenter_to_face( rho_face[lev], density(lev), geom[lev]); + if (m_is_anelastic) { + ref_rho_face[lev][0] = &(*ref_rho_xf)(lev); + ref_rho_face[lev][1] = &(*ref_rho_yf)(lev); + ref_rho_face[lev][2] = &(*ref_rho_zf)(lev); + amrex::average_cellcenter_to_face( + ref_rho_face[lev], (*ref_density)(lev), geom[lev]); + for (int idim = 0; idim < ICNS::ndim; ++idim) { + rho_face[lev][idim]->divide( + *(ref_rho_face[lev][idim]), 0, density.num_comp(), + rho_face[lev][idim]->nGrow()); + } + } + if (m_mesh_mapping) { mac_proj_to_uniform_space( m_repo, u_mac, v_mac, w_mac, rho_face[lev], factor, lev); @@ -187,10 +221,16 @@ void MacProjOp::operator()(const FieldState fstate, const amrex::Real dt) } for (int lev = 0; lev < m_repo.num_active_levels(); ++lev) { - mac_vec[lev][0] = &u_mac(lev); mac_vec[lev][1] = &v_mac(lev); mac_vec[lev][2] = &w_mac(lev); + if (m_is_anelastic) { + for (int idim = 0; idim < ICNS::ndim; ++idim) { + amrex::Multiply( + *(mac_vec[lev][idim]), *(ref_rho_face[lev][idim]), 0, 0, + density.num_comp(), 0); + } + } } m_mac_proj->setUMAC(mac_vec); @@ -209,6 +249,16 @@ void MacProjOp::operator()(const FieldState fstate, const amrex::Real dt) m_mac_proj->project(m_options.rel_tol, m_options.abs_tol); } + if (m_is_anelastic) { + for (int lev = 0; lev < m_repo.num_active_levels(); ++lev) { + for (int idim = 0; idim < ICNS::ndim; ++idim) { + amrex::Divide( + *(mac_vec[lev][idim]), *(ref_rho_face[lev][idim]), 0, 0, + density.num_comp(), 0); + } + } + } + io::print_mlmg_info("MAC_projection", m_mac_proj->getMLMG()); } diff --git a/amr-wind/equation_systems/vof/vof_advection.H b/amr-wind/equation_systems/vof/vof_advection.H index 34068e43a1..fffaee2954 100644 --- a/amr-wind/equation_systems/vof/vof_advection.H +++ b/amr-wind/equation_systems/vof/vof_advection.H @@ -15,7 +15,11 @@ template <> struct AdvectionOp { AdvectionOp( - PDEFields& fields_in, bool /*unused*/, bool /*unused*/, bool /*unused*/) + PDEFields& fields_in, + bool /*unused*/, + bool /*unused*/, + bool /*unused*/, + bool /*unused*/) : fields(fields_in) , u_mac(fields_in.repo.get_field("u_mac")) , v_mac(fields_in.repo.get_field("v_mac")) diff --git a/amr-wind/incflo.H b/amr-wind/incflo.H index 50613d638d..103a5da630 100644 --- a/amr-wind/incflo.H +++ b/amr-wind/incflo.H @@ -173,8 +173,6 @@ private: bool m_do_initial_proj = true; int m_initial_iterations = 3; - bool m_constant_density = true; - bool m_use_godunov = true; // Prescribe advection velocity diff --git a/amr-wind/incflo.cpp b/amr-wind/incflo.cpp index d93ec12467..dfaa553a65 100644 --- a/amr-wind/incflo.cpp +++ b/amr-wind/incflo.cpp @@ -374,7 +374,7 @@ void incflo::init_physics_and_pde() // Register density first so that we can compute its `n+1/2` state before // other scalars attempt to use it in their computations. - if (!m_constant_density) { + if (!pde_mgr.constant_density()) { if (!pde_mgr.scalar_eqns().empty()) { amrex::Abort( "For non-constant density, it must be the first equation " diff --git a/amr-wind/incflo_advance.cpp b/amr-wind/incflo_advance.cpp index b4ef5d95ac..55ec68f25e 100644 --- a/amr-wind/incflo_advance.cpp +++ b/amr-wind/incflo_advance.cpp @@ -288,7 +288,7 @@ void incflo::ApplyPredictor(bool incremental_projection) // ************************************************************************************* // Update density first // ************************************************************************************* - if (m_constant_density) { + if (m_sim.pde_manager().constant_density()) { amr_wind::field_ops::copy(density_nph, density_old, 0, 0, 1, 1); } @@ -544,7 +544,7 @@ void incflo::ApplyCorrector() // ************************************************************************************* // Update density first // ************************************************************************************* - if (m_constant_density) { + if (m_sim.pde_manager().constant_density()) { amr_wind::field_ops::copy(density_nph, density_old, 0, 0, 1, 1); } @@ -694,7 +694,7 @@ void incflo::ApplyPrescribeStep() // ************************************************************************************* // Update density first // ************************************************************************************* - if (m_constant_density) { + if (m_sim.pde_manager().constant_density()) { amr_wind::field_ops::copy(density_nph, density_old, 0, 0, 1, 1); } diff --git a/amr-wind/projection/incflo_apply_nodal_projection.cpp b/amr-wind/projection/incflo_apply_nodal_projection.cpp index fafe4ad5a3..0af66c168e 100644 --- a/amr-wind/projection/incflo_apply_nodal_projection.cpp +++ b/amr-wind/projection/incflo_apply_nodal_projection.cpp @@ -126,7 +126,8 @@ void incflo::ApplyProjection( } } - bool variable_density = + const bool is_anelastic = m_sim.is_anelastic(); + const bool variable_density = (!m_sim.pde_manager().constant_density() || m_sim.physics_manager().contains("MultiPhase")); @@ -143,6 +144,8 @@ void incflo::ApplyProjection( amr_wind::Field const* mesh_detJ = mesh_mapping ? &(m_repo.get_mesh_mapping_detJ(amr_wind::FieldLoc::CELL)) : nullptr; + const auto* ref_density = + is_anelastic ? &(m_repo.get_field("reference_density")) : nullptr; // TODO: Mesh mapping doesn't work with immersed boundaries // Do the pre pressure correction work -- this applies to IB only @@ -266,6 +269,9 @@ void incflo::ApplyProjection( amrex::Array4 detJ = mesh_mapping ? ((*mesh_detJ)(lev).const_array(mfi)) : amrex::Array4(); + const auto& ref_rho = is_anelastic + ? (*ref_density)(lev).const_array(mfi) + : amrex::Array4(); amrex::ParallelFor( bx, ncomp, @@ -276,6 +282,9 @@ void incflo::ApplyProjection( mesh_mapping ? (detJ(i, j, k)) : 1.0; sig(i, j, k, n) = std::pow(fac_cc, -2.) * det_j * scaling_factor / rho(i, j, k); + if (is_anelastic) { + sig(i, j, k, n) *= ref_rho(i, j, k); + } }); } } @@ -296,6 +305,14 @@ void incflo::ApplyProjection( } } + if (is_anelastic) { + for (int lev = 0; lev <= finest_level; ++lev) { + amrex::Multiply( + velocity(lev), (*ref_density)(lev), 0, 0, density[lev]->nComp(), + 0); + } + } + amr_wind::MLMGOptions options("nodal_proj"); if (variable_density || mesh_mapping) { @@ -356,9 +373,18 @@ void incflo::ApplyProjection( } else { nodal_projector->project(options.rel_tol, options.abs_tol); } + amr_wind::io::print_mlmg_info( "Nodal_projection", nodal_projector->getMLMG()); + if (is_anelastic) { + for (int lev = 0; lev <= finest_level; ++lev) { + amrex::Divide( + velocity(lev), (*ref_density)(lev), 0, 0, density[lev]->nComp(), + 0); + } + } + // scale U^* back to -> U = fac/J * U^bar if (mesh_mapping) { velocity.to_stretched_space(); @@ -443,7 +469,8 @@ void incflo::UpdateGradP( // Pressure and sigma are necessary to calculate the pressure gradient - bool variable_density = + const bool is_anelastic = m_sim.is_anelastic(); + const bool variable_density = (!m_sim.pde_manager().constant_density() || m_sim.physics_manager().contains("MultiPhase")); @@ -459,6 +486,11 @@ void incflo::UpdateGradP( amr_wind::Field const* mesh_detJ = mesh_mapping ? &(m_repo.get_mesh_mapping_detJ(amr_wind::FieldLoc::CELL)) : nullptr; + const auto* ref_density = + is_anelastic ? &(m_repo.get_field("reference_density")) : nullptr; + + // ASA does not think we actually need to define sigma here. It + // should not be used by calcGradPhi // Create sigma while accounting for mesh mapping // sigma = 1/(fac^2)*J * dt/rho @@ -482,6 +514,9 @@ void incflo::UpdateGradP( amrex::Array4 detJ = mesh_mapping ? ((*mesh_detJ)(lev).const_array(mfi)) : amrex::Array4(); + const auto& ref_rho = is_anelastic + ? (*ref_density)(lev).const_array(mfi) + : amrex::Array4(); amrex::ParallelFor( bx, ncomp, @@ -492,6 +527,9 @@ void incflo::UpdateGradP( mesh_mapping ? (detJ(i, j, k)) : 1.0; sig(i, j, k, n) = std::pow(fac_cc, -2.) * det_j * scaling_factor / rho(i, j, k); + if (is_anelastic) { + sig(i, j, k, n) *= ref_rho(i, j, k); + } }); } } @@ -510,6 +548,14 @@ void incflo::UpdateGradP( vel.push_back(&(velocity(lev))); } + if (is_anelastic) { + for (int lev = 0; lev <= finest_level; ++lev) { + amrex::Multiply( + velocity(lev), (*ref_density)(lev), 0, 0, density[lev]->nComp(), + 0); + } + } + amr_wind::MLMGOptions options("nodal_proj"); if (variable_density || mesh_mapping) { diff --git a/amr-wind/setup/init.cpp b/amr-wind/setup/init.cpp index 71cf32ad4f..b669a8b2af 100644 --- a/amr-wind/setup/init.cpp +++ b/amr-wind/setup/init.cpp @@ -23,9 +23,6 @@ void incflo::ReadParameters() pp.query("initial_iterations", m_initial_iterations); pp.query("do_initial_proj", m_do_initial_proj); - // Physics - pp.query("constant_density", m_constant_density); - // Godunov-related flags pp.query("use_godunov", m_use_godunov); @@ -84,7 +81,7 @@ void incflo::InitialIterations() vel.copy_state(amr_wind::FieldState::Old, amr_wind::FieldState::New); vel.state(amr_wind::FieldState::Old).fillpatch(m_time.current_time()); - if (m_constant_density) { + if (m_sim.pde_manager().constant_density()) { auto& rho = density(); rho.copy_state( amr_wind::FieldState::Old, amr_wind::FieldState::New); @@ -117,7 +114,7 @@ void incflo::InitialIterations() vel.copy_state( amr_wind::FieldState::New, amr_wind::FieldState::Old); - if (m_constant_density) { + if (m_sim.pde_manager().constant_density()) { auto& rho = density(); rho.copy_state( amr_wind::FieldState::New, amr_wind::FieldState::Old); diff --git a/amr-wind/utilities/MultiLevelVector.cpp b/amr-wind/utilities/MultiLevelVector.cpp index f2354d5eba..3c7e1bbdef 100644 --- a/amr-wind/utilities/MultiLevelVector.cpp +++ b/amr-wind/utilities/MultiLevelVector.cpp @@ -18,7 +18,7 @@ void MultiLevelVector::resize( m_data_d[lev].resize(nc); } copy_host_to_device(); -}; +} void MultiLevelVector::copy_host_to_device() { diff --git a/amr-wind/wind_energy/ABLAnelastic.cpp b/amr-wind/wind_energy/ABLAnelastic.cpp index 52b3509ffe..7c27713e05 100644 --- a/amr-wind/wind_energy/ABLAnelastic.cpp +++ b/amr-wind/wind_energy/ABLAnelastic.cpp @@ -9,10 +9,14 @@ ABLAnelastic::ABLAnelastic(CFDSim& sim) : m_sim(sim) pp.query("anelastic", m_is_anelastic); pp.query("bottom_reference_pressure", m_bottom_reference_pressure); } + std::string godunov_type; + int conserv = 1; { amrex::ParmParse pp("incflo"); pp.queryarr("gravity", m_gravity); pp.query("density", m_rho0_const); + pp.query("godunov_type", godunov_type); + pp.query("icns_conserv", conserv); } if (m_is_anelastic) { { @@ -20,7 +24,16 @@ ABLAnelastic::ABLAnelastic(CFDSim& sim) : m_sim(sim) amrex::ParmParse pp("ICNS"); pp.add("use_perturb_pressure", (bool)true); } - m_sim.repo().declare_field("reference_density", 1, 0, 1); + if (conserv != 1) { + amrex::Abort("ABLAnelastic is not supported for icns_conserv != 1"); + } + if (amrex::toLower(godunov_type) == "bds") { + amrex::Abort("ABLAnelastic is not supported by BDS"); + } + const auto& density = m_sim.repo().get_field("density"); + auto& ref_density = m_sim.repo().declare_field( + "reference_density", 1, density.num_grow()[0], 1); + ref_density.set_default_fillpatch_bc(m_sim.time()); m_sim.repo().declare_nd_field("reference_pressure", 1, 0, 1); } } @@ -47,6 +60,9 @@ void ABLAnelastic::initialize_data() m_density.resize(m_axis, m_sim.mesh().Geom()); m_pressure.resize(m_axis, m_sim.mesh().Geom()); + AMREX_ALWAYS_ASSERT(m_sim.repo().num_active_levels() == m_density.size()); + AMREX_ALWAYS_ASSERT(m_sim.repo().num_active_levels() == m_pressure.size()); + for (int lev = 0; lev < m_density.size(); lev++) { auto& dens = m_density.host_data(lev); auto& pres = m_pressure.host_data(lev); @@ -64,5 +80,7 @@ void ABLAnelastic::initialize_data() auto& p0 = m_sim.repo().get_field("reference_pressure"); m_density.copy_to_field(rho0); m_pressure.copy_to_field(p0); + + rho0.fillpatch(m_sim.time().current_time()); } } // namespace amr_wind diff --git a/test/test_files/abl_anelastic/abl_anelastic.inp b/test/test_files/abl_anelastic/abl_anelastic.inp index 92b08ac78e..d4f33abadb 100644 --- a/test/test_files/abl_anelastic/abl_anelastic.inp +++ b/test/test_files/abl_anelastic/abl_anelastic.inp @@ -27,6 +27,7 @@ incflo.density = 1.0 # Reference density incflo.use_godunov = 1 incflo.godunov_type = "weno_z" +incflo.constant_density = false incflo.diffusion_type = 2 incflo.initial_iterations = 3 incflo.do_initial_proj = 1 diff --git a/unit_tests/equation_systems/test_icns_cstdens.cpp b/unit_tests/equation_systems/test_icns_cstdens.cpp index 1392aceffc..a98c4182ab 100644 --- a/unit_tests/equation_systems/test_icns_cstdens.cpp +++ b/unit_tests/equation_systems/test_icns_cstdens.cpp @@ -23,7 +23,7 @@ class ICNSConstDensTest : public MeshTest // Initialize MAC projection operator const auto& mco = - amr_wind::pde::MacProjOp(sim().repo(), false, false, false); + amr_wind::pde::MacProjOp(sim().repo(), false, false, false, false); // Get background density and check const amrex::Real rho0 = mco.rho0(); EXPECT_EQ(rho0, m_rho_0);