Skip to content

Commit

Permalink
move thin body source terms into ERF_make_mom_sources.cpp (#1572)
Browse files Browse the repository at this point in the history
* move thin body source terms into ERF_make_mom_sources.cpp

* thin body forces need to be constructed after the rest of the RHS is done

* fix oops in last commit
  • Loading branch information
asalmgren authored Apr 9, 2024
1 parent e6bb94f commit 9ac4f3a
Show file tree
Hide file tree
Showing 15 changed files with 333 additions and 228 deletions.
1 change: 1 addition & 0 deletions CMake/BuildERFExe.cmake
Original file line number Diff line number Diff line change
Expand Up @@ -153,6 +153,7 @@ function(build_erf_lib erf_lib_name)
${SRC_DIR}/IO/console_io.cpp
${SRC_DIR}/SourceTerms/ERF_ApplySpongeZoneBCs.cpp
${SRC_DIR}/SourceTerms/ERF_make_buoyancy.cpp
${SRC_DIR}/SourceTerms/ERF_add_thin_body_sources.cpp
${SRC_DIR}/SourceTerms/ERF_make_mom_sources.cpp
${SRC_DIR}/SourceTerms/ERF_make_sources.cpp
${SRC_DIR}/SourceTerms/ERF_moist_set_rhs.cpp
Expand Down
4 changes: 3 additions & 1 deletion Source/BoundaryConditions/ERF_FillPatch.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -561,7 +561,9 @@ ERF::FillBdyCCVels (Vector<MultiFab>& mf_cc_vel)

for (MFIter mfi(mf_cc_vel[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi)
{
const Box& bx = mfi.growntilebox(1);
// Note that we don't fill corners here -- only the cells that share a face
// with interior cells -- this is all that is needed to calculate vorticity
const Box& bx = mfi.tilebox();
const Array4<Real>& vel_arr = mf_cc_vel[lev].array(mfi);

if (!Geom(lev).isPeriodic(0)) {
Expand Down
10 changes: 10 additions & 0 deletions Source/Derive.H
Original file line number Diff line number Diff line change
Expand Up @@ -122,5 +122,15 @@ void erf_dervortz (
const int* bcrec,
const int level);

void erf_dermagvel (
const amrex::Box& bx,
amrex::FArrayBox& derfab,
int dcomp,
int ncomp,
const amrex::FArrayBox& datfab,
const amrex::Geometry& geomdata,
amrex::Real time,
const int* bcrec,
const int level);
}
#endif
26 changes: 26 additions & 0 deletions Source/Derive.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -278,5 +278,31 @@ erf_dervortz (
});
}

void
erf_dermagvel (
const amrex::Box& bx,
amrex::FArrayBox& derfab,
int dcomp,
int ncomp,
const amrex::FArrayBox& datfab,
const amrex::Geometry& /*geomdata*/,
amrex::Real /*time*/,
const int* /*bcrec*/,
const int /*level*/)
{
AMREX_ALWAYS_ASSERT(dcomp == 0);
AMREX_ALWAYS_ASSERT(ncomp == 1);

auto const dat = datfab.array(); // cell-centered velocity
auto tfab = derfab.array(); // cell-centered magvel

ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
{
Real u = dat(i,j,k,0);
Real v = dat(i,j,k,1);
Real w = dat(i,j,k,2);
tfab(i,j,k,dcomp) = std::sqrt(u*u + v*v + w*w);
});
}

} // namespace
2 changes: 1 addition & 1 deletion Source/ERF.H
Original file line number Diff line number Diff line change
Expand Up @@ -767,7 +767,7 @@ private:
// Note that the order of variable names here must match the order in Derive.cpp
const amrex::Vector<std::string> derived_names {"soundspeed", "temp", "theta", "KE", "QKE", "scalar",
"vorticity_x","vorticity_y","vorticity_z",
"divU",
"magvel", "divU",
"pres_hse", "dens_hse", "pressure", "pert_pres", "pert_dens", "eq_pot_temp", "num_turb",
"dpdx", "dpdy", "pres_hse_x", "pres_hse_y",
"z_phys", "detJ" , "mapfac",
Expand Down
138 changes: 70 additions & 68 deletions Source/IO/Plotfile.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -198,6 +198,7 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
if (containerHasElement(plot_var_names, "x_velocity" ) ||
containerHasElement(plot_var_names, "y_velocity" ) ||
containerHasElement(plot_var_names, "z_velocity" ) ||
containerHasElement(plot_var_names, "magvel" ) ||
containerHasElement(plot_var_names, "vorticity_x") ||
containerHasElement(plot_var_names, "vorticity_y") ||
containerHasElement(plot_var_names, "vorticity_z") ) {
Expand Down Expand Up @@ -294,6 +295,7 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
calculate_derived("vorticity_x", mf_cc_vel[lev] , derived::erf_dervortx);
calculate_derived("vorticity_y", mf_cc_vel[lev] , derived::erf_dervorty);
calculate_derived("vorticity_z", mf_cc_vel[lev] , derived::erf_dervortz);
calculate_derived("magvel" , mf_cc_vel[lev] , derived::erf_dermagvel);

if (containerHasElement(plot_var_names, "divU"))
{
Expand Down Expand Up @@ -1256,10 +1258,10 @@ ERF::WriteMultiLevelPlotfileWithTerrain (const std::string& plotfilename, int nl
{
BL_PROFILE("WriteMultiLevelPlotfileWithTerrain()");

BL_ASSERT(nlevels <= mf.size());
BL_ASSERT(nlevels <= ref_ratio.size()+1);
BL_ASSERT(nlevels <= level_steps.size());
BL_ASSERT(mf[0]->nComp() == varnames.size());
AMREX_ALWAYS_ASSERT(nlevels <= mf.size());
AMREX_ALWAYS_ASSERT(nlevels <= ref_ratio.size()+1);
AMREX_ALWAYS_ASSERT(nlevels <= level_steps.size());
AMREX_ALWAYS_ASSERT(mf[0]->nComp() == varnames.size());

bool callBarrier(false);
PreBuildDirectorHierarchy(plotfilename, levelPrefix, nlevels, callBarrier);
Expand Down Expand Up @@ -1338,78 +1340,78 @@ ERF::WriteGenericPlotfileHeaderWithTerrain (std::ostream &HeaderFile,
const std::string &levelPrefix,
const std::string &mfPrefix) const
{
BL_ASSERT(nlevels <= bArray.size());
BL_ASSERT(nlevels <= ref_ratio.size()+1);
BL_ASSERT(nlevels <= level_steps.size());
AMREX_ALWAYS_ASSERT(nlevels <= bArray.size());
AMREX_ALWAYS_ASSERT(nlevels <= ref_ratio.size()+1);
AMREX_ALWAYS_ASSERT(nlevels <= level_steps.size());

HeaderFile.precision(17);
HeaderFile.precision(17);

// ---- this is the generic plot file type name
HeaderFile << versionName << '\n';
// ---- this is the generic plot file type name
HeaderFile << versionName << '\n';

HeaderFile << varnames.size() << '\n';
HeaderFile << varnames.size() << '\n';

for (int ivar = 0; ivar < varnames.size(); ++ivar) {
HeaderFile << varnames[ivar] << "\n";
}
HeaderFile << AMREX_SPACEDIM << '\n';
HeaderFile << time << '\n';
HeaderFile << finest_level << '\n';
for (int i = 0; i < AMREX_SPACEDIM; ++i) {
HeaderFile << geom[0].ProbLo(i) << ' ';
}
HeaderFile << '\n';
for (int i = 0; i < AMREX_SPACEDIM; ++i) {
HeaderFile << geom[0].ProbHi(i) << ' ';
}
HeaderFile << '\n';
for (int i = 0; i < finest_level; ++i) {
HeaderFile << ref_ratio[i][0] << ' ';
}
HeaderFile << '\n';
for (int i = 0; i <= finest_level; ++i) {
HeaderFile << geom[i].Domain() << ' ';
}
HeaderFile << '\n';
for (int i = 0; i <= finest_level; ++i) {
HeaderFile << level_steps[i] << ' ';
for (int ivar = 0; ivar < varnames.size(); ++ivar) {
HeaderFile << varnames[ivar] << "\n";
}
HeaderFile << AMREX_SPACEDIM << '\n';
HeaderFile << time << '\n';
HeaderFile << finest_level << '\n';
for (int i = 0; i < AMREX_SPACEDIM; ++i) {
HeaderFile << geom[0].ProbLo(i) << ' ';
}
HeaderFile << '\n';
for (int i = 0; i < AMREX_SPACEDIM; ++i) {
HeaderFile << geom[0].ProbHi(i) << ' ';
}
HeaderFile << '\n';
for (int i = 0; i < finest_level; ++i) {
HeaderFile << ref_ratio[i][0] << ' ';
}
HeaderFile << '\n';
for (int i = 0; i <= finest_level; ++i) {
HeaderFile << geom[i].Domain() << ' ';
}
HeaderFile << '\n';
for (int i = 0; i <= finest_level; ++i) {
HeaderFile << level_steps[i] << ' ';
}
HeaderFile << '\n';
for (int i = 0; i <= finest_level; ++i) {
for (int k = 0; k < AMREX_SPACEDIM; ++k) {
HeaderFile << geom[i].CellSize()[k] << ' ';
}
HeaderFile << '\n';
for (int i = 0; i <= finest_level; ++i) {
for (int k = 0; k < AMREX_SPACEDIM; ++k) {
HeaderFile << geom[i].CellSize()[k] << ' ';
}
HeaderFile << '\n';
}
HeaderFile << (int) geom[0].Coord() << '\n';
HeaderFile << "0\n";
}
HeaderFile << (int) geom[0].Coord() << '\n';
HeaderFile << "0\n";

for (int level = 0; level <= finest_level; ++level) {
HeaderFile << level << ' ' << bArray[level].size() << ' ' << time << '\n';
HeaderFile << level_steps[level] << '\n';
for (int level = 0; level <= finest_level; ++level) {
HeaderFile << level << ' ' << bArray[level].size() << ' ' << time << '\n';
HeaderFile << level_steps[level] << '\n';

const IntVect& domain_lo = geom[level].Domain().smallEnd();
for (int i = 0; i < bArray[level].size(); ++i)
{
// Need to shift because the RealBox ctor we call takes the
// physical location of index (0,0,0). This does not affect
// the usual cases where the domain index starts with 0.
const Box& b = shift(bArray[level][i], -domain_lo);
RealBox loc = RealBox(b, geom[level].CellSize(), geom[level].ProbLo());
for (int n = 0; n < AMREX_SPACEDIM; ++n) {
HeaderFile << loc.lo(n) << ' ' << loc.hi(n) << '\n';
}
const IntVect& domain_lo = geom[level].Domain().smallEnd();
for (int i = 0; i < bArray[level].size(); ++i)
{
// Need to shift because the RealBox ctor we call takes the
// physical location of index (0,0,0). This does not affect
// the usual cases where the domain index starts with 0.
const Box& b = shift(bArray[level][i], -domain_lo);
RealBox loc = RealBox(b, geom[level].CellSize(), geom[level].ProbLo());
for (int n = 0; n < AMREX_SPACEDIM; ++n) {
HeaderFile << loc.lo(n) << ' ' << loc.hi(n) << '\n';
}

HeaderFile << MultiFabHeaderPath(level, levelPrefix, mfPrefix) << '\n';
}
HeaderFile << "1" << "\n";
HeaderFile << "3" << "\n";
HeaderFile << "amrexvec_nu_x" << "\n";
HeaderFile << "amrexvec_nu_y" << "\n";
HeaderFile << "amrexvec_nu_z" << "\n";
std::string mf_nodal_prefix = "Nu_nd";
for (int level = 0; level <= finest_level; ++level) {
HeaderFile << MultiFabHeaderPath(level, levelPrefix, mf_nodal_prefix) << '\n';
}

HeaderFile << MultiFabHeaderPath(level, levelPrefix, mfPrefix) << '\n';
}
HeaderFile << "1" << "\n";
HeaderFile << "3" << "\n";
HeaderFile << "amrexvec_nu_x" << "\n";
HeaderFile << "amrexvec_nu_y" << "\n";
HeaderFile << "amrexvec_nu_z" << "\n";
std::string mf_nodal_prefix = "Nu_nd";
for (int level = 0; level <= finest_level; ++level) {
HeaderFile << MultiFabHeaderPath(level, levelPrefix, mf_nodal_prefix) << '\n';
}
}
131 changes: 131 additions & 0 deletions Source/SourceTerms/ERF_add_thin_body_sources.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,131 @@
#include <AMReX_MultiFab.H>
#include <AMReX_ArrayLim.H>
#include <AMReX_BCRec.H>
#include <AMReX_TableData.H>
#include <AMReX_GpuContainers.H>

#include <NumericalDiffusion.H>
#include <PlaneAverage.H>
#include <TI_slow_headers.H>
#include <Src_headers.H>
#include <Utils.H>

using namespace amrex;

/**
* Function for computing the slow RHS for the evolution equations for the density, potential temperature and momentum.
*
* @param[in] xmom_src source terms for x-momentum
* @param[in] ymom_src source terms for y-momentum
* @param[in] zmom_src source terms for z-momentum
* @param[in] xflux_imask_lev thin-body mask on x-faces
* @param[in] yflux_imask_lev thin-body mask on y-faces
* @param[in] zflux_imask_lev thin-body mask on z-faces
* @param[in] thin_xforce_lev x-component of forces on thin immersed bodies
* @param[in] thin_yforce_lev y-component of forces on thin immersed bodies
* @param[in] thin_zforce_lev z-component of forces on thin immersed bodies
*/

void add_thin_body_sources ( MultiFab & xmom_src,
MultiFab & ymom_src,
MultiFab & zmom_src,
std::unique_ptr<iMultiFab>& xflux_imask_lev,
std::unique_ptr<iMultiFab>& yflux_imask_lev,
std::unique_ptr<iMultiFab>& zflux_imask_lev,
std::unique_ptr<MultiFab>& thin_xforce_lev,
std::unique_ptr<MultiFab>& thin_yforce_lev,
std::unique_ptr<MultiFab>& thin_zforce_lev)
{
BL_PROFILE_REGION("erf_add_thin_body_sources()");

const bool l_have_thin_xforce = (thin_xforce_lev != nullptr);
const bool l_have_thin_yforce = (thin_yforce_lev != nullptr);
const bool l_have_thin_zforce = (thin_zforce_lev != nullptr);

// *****************************************************************************
// If a thin immersed body is present, add forcing terms
// *****************************************************************************
if (l_have_thin_xforce) {
MultiFab::Copy(*thin_xforce_lev, xmom_src, 0, 0, 1, 0);
thin_xforce_lev->mult(-1., 0, 1, 0);
ApplyInvertedMask(*thin_xforce_lev, *xflux_imask_lev, 0);
MultiFab::Add(xmom_src, *thin_xforce_lev, 0, 0, 1, 0);
}

if (l_have_thin_yforce) {
MultiFab::Copy(*thin_yforce_lev, ymom_src, 0, 0, 1, 0);
thin_yforce_lev->mult(-1., 0, 1, 0);
ApplyInvertedMask(*thin_yforce_lev, *yflux_imask_lev, 0);
MultiFab::Add(ymom_src, *thin_yforce_lev, 0, 0, 1, 0);
}

if (l_have_thin_zforce) {
MultiFab::Copy(*thin_zforce_lev, zmom_src, 0, 0, 1, 0);
thin_zforce_lev->mult(-1., 0, 1, 0);
ApplyInvertedMask(*thin_zforce_lev, *zflux_imask_lev, 0);
MultiFab::Add(zmom_src, *thin_zforce_lev, 0, 0, 1, 0);
}

#if 0
#ifndef AMREX_USE_GPU
if (l_have_thin_xforce) {
// TODO: Implement particles to better track and output these data
if (nrk==2) {
for ( MFIter mfi(S_data[IntVars::cons],TileNoZ()); mfi.isValid(); ++mfi)
{
const Box& tbx = mfi.nodaltilebox(0);
const Array4<const Real> & fx = thin_xforce_lev->const_array(mfi);
const Array4<const int> & mask = xflux_imask_lev->const_array(mfi);
ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
if (mask(i,j,k)==0) {
amrex::AllPrint() << "thin body fx"<<IntVect(i,j,k)<<" = " << fx(i,j,k) << std::endl;
}
});
}
}
}
#endif
#endif

#if 0
#ifndef AMREX_USE_GPU
if (l_have_thin_yforce) {
// TODO: Implement particles to better track and output these data
if (nrk==2) {
for ( MFIter mfi(S_data[IntVars::cons],TileNoZ()); mfi.isValid(); ++mfi)
{
const Box& tby = mfi.nodaltilebox(1);
const Array4<const Real> & fy = thin_yforce_lev->const_array(mfi);
const Array4<const int> & mask = yflux_imask_lev->const_array(mfi);
ParallelFor(tby, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
if (mask(i,j,k)==0) {
amrex::AllPrint() << "thin body fy"<<IntVect(i,j,k)<<" = " << fy(i,j,k) << std::endl;
}
});
}
}
}
#endif
#endif

#if 0
#ifndef AMREX_USE_GPU
if (l_have_thin_zforce) {
// TODO: Implement particles to better track and output these data
if (nrk==2) {
for ( MFIter mfi(S_data[IntVars::cons],TileNoZ()); mfi.isValid(); ++mfi)
{
const Box& tbz = mfi.nodaltilebox(2);
const Array4<const Real> & fz = thin_zforce_lev->const_array(mfi);
const Array4<const int> & mask = zflux_imask_lev->const_array(mfi);
ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
if (mask(i,j,k)==0) {
amrex::AllPrint() << "thin body fz"<<IntVect(i,j,k)<<" = " << fz(i,j,k) << std::endl;
}
});
}
}
}
#endif
#endif
}
Loading

0 comments on commit 9ac4f3a

Please sign in to comment.