Skip to content

Commit

Permalink
manifold-specific transport, derives, divu
Browse files Browse the repository at this point in the history
  • Loading branch information
@baperry2
baperry2 committed Sep 18, 2024
1 parent 2670304 commit 94cf2e4
Show file tree
Hide file tree
Showing 7 changed files with 239 additions and 47 deletions.
16 changes: 16 additions & 0 deletions Source/PeleLMeX_DeriveFunc.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -290,6 +290,22 @@ void pelelmex_deruserdef(
const amrex::Vector<amrex::BCRec>& bcrec,
int level);

#ifdef USE_MANIFOLD_EOS
void pelelmex_dermaniout(
PeleLM* a_pelelm,
const amrex::Box& bx,
amrex::FArrayBox& derfab,
int dcomp,
int ncomp,
const amrex::FArrayBox& statefab,
const amrex::FArrayBox& reactfab,
const amrex::FArrayBox& pressfab,
const amrex::Geometry& geom,
amrex::Real time,
const amrex::Vector<amrex::BCRec>& bcrec,
int level);
#endif

#ifdef PELE_USE_EFIELD
#include <PeleLMeX_EFDeriveFunc.H>
#endif
Expand Down
55 changes: 52 additions & 3 deletions Source/PeleLMeX_DeriveFunc.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@
#include "PeleLMeX.H"
#include "PeleLMeX_K.H"
#include "PeleLMeX_DeriveFunc.H"
#include "PeleLMeX_Utils.H"

#include <PelePhysics.H>
#include <mechanism.H>
Expand Down Expand Up @@ -1382,7 +1383,7 @@ pelelmex_derdiffc(
bx, [do_fixed_Le, do_fixed_Pr, do_soret, LeInv, PrInv, rhoY, T, rhoD,
rhotheta, lambda, mu, ltransparm,
leosparm] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
getTransportCoeff(
getTransportCoeff<pele::physics::PhysicsType::eos_type>(
i, j, k, do_fixed_Le, do_fixed_Pr, do_soret, LeInv, PrInv, rhoY, T,
rhoD, rhotheta, lambda, mu, ltransparm, leosparm);
});
Expand Down Expand Up @@ -1428,7 +1429,7 @@ pelelmex_derlambda(
bx, [do_fixed_Le, do_fixed_Pr, do_soret, LeInv, PrInv, rhoY, T, rhoD,
rhotheta, lambda, mu, ltransparm,
leosparm] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
getTransportCoeff(
getTransportCoeff<pele::physics::PhysicsType::eos_type>(
i, j, k, do_fixed_Le, do_fixed_Pr, do_soret, LeInv, PrInv, rhoY, T,
rhoD, rhotheta, lambda, mu, ltransparm, leosparm);
});
Expand All @@ -1451,7 +1452,6 @@ pelelmex_derdmap(
Real /*time*/,
const Vector<BCRec>& /*bcrec*/,
int /*level*/)

{
AMREX_ASSERT(derfab.box().contains(bx));
auto der = derfab.array(dcomp);
Expand All @@ -1460,3 +1460,52 @@ pelelmex_derdmap(
der(i, j, k) = myrank;
});
}

//
// Derive manifold output quantities
//
#ifdef USE_MANIFOLD_EOS
void
pelelmex_dermaniout(
PeleLM* a_pelelm,
const Box& bx,
FArrayBox& derfab,
int dcomp,
int ncomp,
const FArrayBox& statefab,
const FArrayBox& /*reactfab*/,
const FArrayBox& /*pressfab*/,
const Geometry& /*geom*/,
Real /*time*/,
const Vector<BCRec>& /*bcrec*/,
int /*level*/)
{
auto& h_manf_data =
a_pelelm->eos_parms.host_only_parm().manfunc_par->host_parm();
auto* d_manf_data =
a_pelelm->eos_parms.host_only_parm().manfunc_par->device_parm();
int nmanivar = h_manf_data.Nvar;

amrex::ignore_unused(ncomp);
AMREX_ASSERT(derfab.box().contains(bx));
AMREX_ASSERT(statefab.box().contains(bx));
AMREX_ASSERT(derfab.nComp() >= dcomp + ncomp);
AMREX_ASSERT(statefab.nComp() >= NUM_SPECIES + 1);
AMREX_ASSERT(ncomp == nmanivar);
AMREX_ASSERT(!a_pelelm->m_incompressible);

auto const in_dat = statefab.array();
auto der = derfab.array(dcomp);
amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
amrex::Real rho, rhoinv, maniparm[NUM_SPECIES];
// pele::physics::PhysicsType::eos_type::RY2RRinvY(array4_to_array(i, j, k,
// in_dat, FIRSTSPEC).data(), rho, rhoinv, maniparm);
pele::physics::PhysicsType::eos_type::RY2RRinvY(
in_dat.ptr(i, j, k, FIRSTSPEC), rho, rhoinv, maniparm);
pele::physics::BlackBoxFunctionFactory<
pele::physics::eos::ManifoldFunctionType>
manfunc{d_manf_data};
manfunc.get_func()->get_all_values(maniparm, der.ptr(i, j, k));
});
}
#endif
4 changes: 2 additions & 2 deletions Source/PeleLMeX_Eos.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -139,7 +139,7 @@ PeleLM::calcDivU(
if (flag(i, j, k).isCovered()) {
divu(i, j, k) = 0.0;
} else {
compute_divu(
compute_divu<pele::physics::PhysicsType::eos_type>(
i, j, k, rhoY, T, SpecD, Fourier, DiffDiff, r, extRhoY, extRhoH,
divu, use_react, leosparm);
}
Expand All @@ -151,7 +151,7 @@ PeleLM::calcDivU(
bx,
[rhoY, T, SpecD, Fourier, DiffDiff, r, extRhoY, extRhoH, divu,
use_react, leosparm] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
compute_divu(
compute_divu<pele::physics::PhysicsType::eos_type>(
i, j, k, rhoY, T, SpecD, Fourier, DiffDiff, r, extRhoY, extRhoH,
divu, use_react, leosparm);
});
Expand Down
136 changes: 124 additions & 12 deletions Source/PeleLMeX_K.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,9 +4,8 @@
#include <mechanism.H>
#include <PelePhysics.H>

AMREX_GPU_DEVICE
AMREX_FORCE_INLINE
void
template <typename EOSType>
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void
getTransportCoeff(
int i,
int j,
Expand Down Expand Up @@ -113,6 +112,70 @@ getTransportCoeff(
}
}

template <>
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void
getTransportCoeff<pele::physics::eos::Manifold>(
int i,
int j,
int k,
const bool /*do_fixed_Le*/,
const bool /*do_fixed_Pr*/,
const bool do_soret,
amrex::Real /*LeInv*/,
amrex::Real /*PrInv*/,
amrex::Array4<const amrex::Real> const& rhoY,
amrex::Array4<const amrex::Real> const& T,
amrex::Array4<amrex::Real> const& rhoDi,
amrex::Array4<amrex::Real> const& /*rhotheta*/,
amrex::Array4<amrex::Real> const& lambda,
amrex::Array4<amrex::Real> const& mu,
pele::physics::transport::TransParm<
pele::physics::PhysicsType::eos_type,
pele::physics::PhysicsType::transport_type> const* trans_parm,
const pele::physics::eos::EosParm<pele::physics::PhysicsType::eos_type>*
eosparm) noexcept
{
using namespace amrex::literals;

// No Soret with manifold moels
AMREX_ASSERT(!do_soret);

// Get rho & Y from rhoY
amrex::Real rho, rhoinv;
amrex::Real y[NUM_SPECIES] = {0.0};
amrex::Real massdens[NUM_SPECIES];
for (int n = 0; n < NUM_SPECIES; n++) {
massdens[n] = rhoY(i, j, k, n);
}
auto eos = pele::physics::PhysicsType::eos(eosparm);
eos.RY2RRinvY(massdens, rho, rhoinv, y);

rho *= 1.0e-3_rt; // MKS -> CGS conversion
amrex::Real rhoDi_cgs[NUM_SPECIES] = {0.0};
amrex::Real lambda_cgs = 0.0_rt;
amrex::Real mu_cgs = 0.0_rt;
amrex::Real dummy_xi = 0.0_rt;
amrex::Real* dummy_chi = nullptr;
amrex::Real Tloc = T(i, j, k);

bool get_xi = false;
bool get_mu = true;
bool get_lam = false;
bool get_Ddiag = true;
bool get_chi = do_soret;
auto trans = pele::physics::PhysicsType::transport();
trans.transport(
get_xi, get_mu, get_lam, get_Ddiag, get_chi, Tloc, rho, y, rhoDi_cgs,
dummy_chi, mu_cgs, dummy_xi, lambda_cgs, trans_parm);

// Do CGS -> MKS conversions
for (int n = 0; n < NUM_SPECIES; n++) {
rhoDi(i, j, k, n) = rhoDi_cgs[n] * 1.0e-1_rt;
}
lambda(i, j, k) = 0.0; // No need to carry lambda for manifold
mu(i, j, k) = mu_cgs * 1.0e-1_rt;
}

AMREX_GPU_DEVICE
AMREX_FORCE_INLINE
void
Expand All @@ -130,15 +193,14 @@ getVelViscosity(
using namespace amrex::literals;

// Get rho & Y from rhoY
amrex::Real rho = 0.0_rt;
amrex::Real massdens[NUM_SPECIES];
for (int n = 0; n < NUM_SPECIES; n++) {
rho += rhoY(i, j, k, n);
massdens[n] = rhoY(i, j, k, n);
}
amrex::Real rho = 0.0_rt;
amrex::Real rhoinv = 1.0_rt / rho;
amrex::Real y[NUM_SPECIES] = {0.0};
for (int n = 0; n < NUM_SPECIES; n++) {
y[n] = rhoY(i, j, k, n) * rhoinv;
}
pele::physics::PhysicsType::eos_type::RY2RRinvY(massdens, rho, rhoinv, y);

rho *= 1.0e-3_rt; // MKS -> CGS conversion
amrex::Real temp = T(i, j, k);
Expand Down Expand Up @@ -191,9 +253,8 @@ getPGivenRTY(
P(i, j, k) = P(i, j, k) * 0.1_rt; // CGS -> MKS conversion
}

AMREX_GPU_DEVICE
AMREX_FORCE_INLINE
void
template <typename EOSType>
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void
compute_divu(
int i,
int j,
Expand Down Expand Up @@ -251,6 +312,56 @@ compute_divu(
}
}

template <>
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void
compute_divu<pele::physics::eos::Manifold>(
int i,
int j,
int k,
amrex::Array4<const amrex::Real> const& rhoY,
amrex::Array4<const amrex::Real> const& /*T*/,
amrex::Array4<const amrex::Real> const& specDiff,
amrex::Array4<const amrex::Real> const& /*tempDiff*/,
amrex::Array4<const amrex::Real> const& /*specEnthDiff*/,
amrex::Array4<const amrex::Real> const& rhoYdot,
amrex::Array4<const amrex::Real> const& extRhoY,
amrex::Array4<const amrex::Real> const& /*extRhoH*/,
amrex::Array4<amrex::Real> const& divu,
int do_react,
pele::physics::eos::EosParm<pele::physics::PhysicsType::eos_type> const*
eosparm) noexcept
{
using namespace amrex::literals;

// Get rho & Ys from rhoYs.
amrex::Real rho, rhoinv;
amrex::Real y[NUM_SPECIES] = {0.0_rt};
amrex::Real massdens[NUM_SPECIES] = {0.0_rt};
for (int n = 0; n < NUM_SPECIES; n++) {
massdens[n] = rhoY(i, j, k, n);
}

auto eos = pele::physics::PhysicsType::eos(eosparm);
eos.RY2RRinvY(massdens, rho, rhoinv, y);
amrex::Real derivs[NUM_SPECIES - 1] = {0.0};
rho /= 1.0e3; // Unit conversion MKS -> CGS
eos.RY2dRdY(rho, y, derivs);

divu(i, j, k) = 0.0;

for (int n = 0; n < NUM_SPECIES - 1; n++) {
derivs[n] *= 1.0e3; // Unit conversion CGS -> MKS
divu(i, j, k) += derivs[n] * specDiff(i, j, k, n);
if (do_react != 0)
divu(i, j, k) += derivs[n] * (rhoYdot(i, j, k, n));
divu(i, j, k) += derivs[n] * (extRhoY(i, j, k, n));
}
// First -rhoinv: from divu = -1/rho Drho/Dt
// Second rhoinv: convert from rho*Dxi/Dt (which RHS terms are evaluated for)
// to Dxi/Dt
divu(i, j, k) *= -rhoinv * rhoinv;
}

AMREX_GPU_DEVICE
AMREX_FORCE_INLINE
void
Expand Down Expand Up @@ -624,7 +735,8 @@ makeVelForce(
const amrex::Real& rho_incomp,
int pseudo_gravity,
int pseudo_gravity_dir,
const amrex::Real& /*time*/,
const amrex::Real& /*time*/
,
amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const gravity,
amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const gp0,
const amrex::Real& dV_control,
Expand Down
Loading

0 comments on commit 94cf2e4

Please sign in to comment.