Update anelastic (#961)

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;
 

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<ABL>();
+ return abl.anelastic().is_anelastic();
+ }
+ return false;
+}
 } // namespace amr_wind

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"))

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"))

amr-wind/equation_systems/PDE.H

@@ -63,13 +63,13 @@ public:
  new TurbulenceOp<PDE>(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<PDE, Scheme>(
  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<PDE, Scheme>(
  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());

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();
 }

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<ICNS, fvm::Godunov>
  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<ICNS, fvm::MOL>
  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(

amr-wind/equation_systems/icns/icns_advection.cpp

@@ -41,12 +41,17 @@ amrex::Array<amrex::LinOpBCType, AMREX_SPACEDIM> 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<ScratchField> 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<amrex::Array<amrex::MultiFab*, ICNS::ndim>> 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<amrex::Array<amrex::MultiFab const*, ICNS::ndim>>
  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<amrex::Array<amrex::MultiFab*, ICNS::ndim>> 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());
 }
 

amr-wind/equation_systems/vof/vof_advection.H

@@ -15,7 +15,11 @@ template <>
 struct AdvectionOp<VOF, fvm::Godunov>
 {
  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"))

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

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 "

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);
  }