Skip to content

Commit

Permalink
Merge branch 'development' into ProfileOutFix
Browse files Browse the repository at this point in the history
  • Loading branch information
@AMLattanzi
AMLattanzi authored Dec 9, 2024
2 parents 03d3d70 + 5a88eda commit 01108d8
Show file tree
Hide file tree
Showing 4 changed files with 40 additions and 49 deletions.
4 changes: 2 additions & 2 deletions Source/ERF.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -164,8 +164,8 @@ public:
void WriteMyEBSurface ();
#endif

// Compute the divergence of the momenta
void compute_divergence (int lev, amrex::MultiFab& rhs, amrex::Vector<amrex::MultiFab>& mom_mf,
// Compute the divergence -- whether EB, no-terrain, flat terrain or general terrain
void compute_divergence (int lev, amrex::MultiFab& rhs, amrex::Array<amrex::MultiFab const*,AMREX_SPACEDIM> rho0_u_const,
amrex::Geometry const& geom_at_lev);

// Project the velocities to be divergence-free -- this is only relevant if anelastic == 1
Expand Down
3 changes: 1 addition & 2 deletions Source/IO/ERF_Plotfile.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,6 @@
#include <ERF_EOS.H>
#include <ERF.H>
#include "AMReX_Interp_3D_C.H"
#include "AMReX_PlotFileUtil.H"
#include "ERF_TerrainMetrics.H"
#include "ERF_Constants.H"

Expand Down Expand Up @@ -368,7 +367,7 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector<std::string> p
u[0] = &(vars_new[lev][Vars::xvel]);
u[1] = &(vars_new[lev][Vars::yvel]);
u[2] = &(vars_new[lev][Vars::zvel]);
computeDivergence(dmf, u, geom[lev]);
compute_divergence (lev, dmf, u, geom[lev]);
mf_comp += 1;
}

Expand Down
73 changes: 30 additions & 43 deletions Source/LinearSolvers/ERF_ComputeDivergence.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,18 +7,13 @@ using namespace amrex;
* Project the single-level velocity field to enforce incompressibility
* Note that the level may or may not be level 0.
*/
void ERF::compute_divergence (int lev, MultiFab& rhs, Vector<MultiFab>& mom_mf, Geometry const& geom_at_lev)
void ERF::compute_divergence (int lev, MultiFab& rhs, Array<MultiFab const*,AMREX_SPACEDIM> rho0_u_const, Geometry const& geom_at_lev)
{
BL_PROFILE("ERF::compute_divergence()");

bool l_use_terrain = (solverChoice.terrain_type != TerrainType::None);

auto dxInv = geom[lev].InvCellSizeArray();

Array<MultiFab const*, AMREX_SPACEDIM> rho0_u_const;
rho0_u_const[0] = &mom_mf[IntVars::xmom];
rho0_u_const[1] = &mom_mf[IntVars::ymom];
rho0_u_const[2] = &mom_mf[IntVars::zmom];
auto dxInv = geom_at_lev.InvCellSizeArray();

// ****************************************************************************
// Compute divergence which will form RHS
Expand All @@ -33,45 +28,37 @@ void ERF::compute_divergence (int lev, MultiFab& rhs, Vector<MultiFab>& mom_mf,
for ( MFIter mfi(rhs,TilingIfNotGPU()); mfi.isValid(); ++mfi)
{
Box bx = mfi.tilebox();
const Array4<Real const>& rho0u_arr = mom_mf[IntVars::xmom].const_array(mfi);
const Array4<Real const>& rho0v_arr = mom_mf[IntVars::ymom].const_array(mfi);
const Array4<Real const>& rho0w_arr = mom_mf[IntVars::zmom].const_array(mfi);
const Array4<Real const>& rho0u_arr = rho0_u_const[0]->const_array(mfi);
const Array4<Real const>& rho0v_arr = rho0_u_const[1]->const_array(mfi);
const Array4<Real const>& rho0w_arr = rho0_u_const[2]->const_array(mfi);
const Array4<Real >& rhs_arr = rhs.array(mfi);

Real* stretched_dz_d_ptr = stretched_dz_d[lev].data();
ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
Real dz = stretched_dz_d_ptr[k];
rhs_arr(i,j,k) = (rho0u_arr(i+1,j,k) - rho0u_arr(i,j,k)) * dxInv[0]
+(rho0v_arr(i,j+1,k) - rho0v_arr(i,j,k)) * dxInv[1]
+(rho0w_arr(i,j,k+1) - rho0w_arr(i,j,k)) / dz;
});
} // mfi
}
else if (l_use_terrain) // terrain is not flat
{
//
// Note we compute the divergence using "rho0w" == Omega
//
for ( MFIter mfi(rhs,TilingIfNotGPU()); mfi.isValid(); ++mfi)
{
Box bx = mfi.tilebox();
const Array4<Real >& rho0u_arr = mom_mf[IntVars::xmom].array(mfi);
const Array4<Real >& rho0v_arr = mom_mf[IntVars::ymom].array(mfi);
const Array4<Real >& rho0w_arr = mom_mf[IntVars::zmom].array(mfi);
const Array4<Real >& rhs_arr = rhs.array(mfi);
if (SolverChoice::terrain_is_flat) { // flat terrain
Real* stretched_dz_d_ptr = stretched_dz_d[lev].data();
ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
Real dz = stretched_dz_d_ptr[k];
rhs_arr(i,j,k) = (rho0u_arr(i+1,j,k) - rho0u_arr(i,j,k)) * dxInv[0]
+(rho0v_arr(i,j+1,k) - rho0v_arr(i,j,k)) * dxInv[1]
+(rho0w_arr(i,j,k+1) - rho0w_arr(i,j,k)) / dz;
});
} else { // terrain is not flat

const Array4<Real const>& ax_arr = ax[lev]->const_array(mfi);
const Array4<Real const>& ay_arr = ay[lev]->const_array(mfi);
const Array4<Real const>& dJ_arr = detJ_cc[lev]->const_array(mfi);
//
// az == 1 for terrain-fitted coordinates
//
ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
rhs_arr(i,j,k) = ((ax_arr(i+1,j,k)*rho0u_arr(i+1,j,k) - ax_arr(i,j,k)*rho0u_arr(i,j,k)) * dxInv[0]
+(ay_arr(i,j+1,k)*rho0v_arr(i,j+1,k) - ay_arr(i,j,k)*rho0v_arr(i,j,k)) * dxInv[1]
+( rho0w_arr(i,j,k+1) - rho0w_arr(i,j,k)) * dxInv[2]) / dJ_arr(i,j,k);
});
//
// Note we compute the divergence using "rho0w" == Omega
//
const Array4<Real const>& ax_arr = ax[lev]->const_array(mfi);
const Array4<Real const>& ay_arr = ay[lev]->const_array(mfi);
const Array4<Real const>& dJ_arr = detJ_cc[lev]->const_array(mfi);
//
// az == 1 for terrain-fitted coordinates
//
ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
rhs_arr(i,j,k) = ((ax_arr(i+1,j,k)*rho0u_arr(i+1,j,k) - ax_arr(i,j,k)*rho0u_arr(i,j,k)) * dxInv[0]
+(ay_arr(i,j+1,k)*rho0v_arr(i,j+1,k) - ay_arr(i,j,k)*rho0v_arr(i,j,k)) * dxInv[1]
+( rho0w_arr(i,j,k+1) - rho0w_arr(i,j,k)) * dxInv[2]) / dJ_arr(i,j,k);
});
}
} // mfi

}
Expand Down
9 changes: 7 additions & 2 deletions Source/LinearSolvers/ERF_PoissonSolve.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -72,7 +72,12 @@ void ERF::project_velocities (int lev, Real l_dt, Vector<MultiFab>& mom_mf, Mult
// Compute divergence which will form RHS
// Note that we replace "rho0w" with the contravariant momentum, Omega
// ****************************************************************************
compute_divergence(lev, rhs[0], mom_mf, geom_tmp[0]);
Array<MultiFab const*, AMREX_SPACEDIM> rho0_u_const;
rho0_u_const[0] = &mom_mf[IntVars::xmom];
rho0_u_const[1] = &mom_mf[IntVars::ymom];
rho0_u_const[2] = &mom_mf[IntVars::zmom];

compute_divergence(lev, rhs[0], rho0_u_const, geom_tmp[0]);

Real rhsnorm = rhs[0].norm0();

Expand Down Expand Up @@ -212,7 +217,7 @@ void ERF::project_velocities (int lev, Real l_dt, Vector<MultiFab>& mom_mf, Mult
//
if (mg_verbose > 0)
{
compute_divergence(lev, rhs[0], mom_mf, geom_tmp[0]);
compute_divergence(lev, rhs[0], rho0_u_const, geom_tmp[0]);

Print() << "Max/L2 norm of divergence after solve at level " << lev << " : " << rhs[0].norm0() << " " << rhs[0].norm2() << std::endl;

Expand Down

0 comments on commit 01108d8

Please sign in to comment.