Skip to content

Commit

Permalink
some more diagnostics for the subch_planar problem
Browse files Browse the repository at this point in the history
  • Loading branch information
@zingale
zingale committed Jun 24, 2024
1 parent 830d937 commit be9fc30
Show file tree
Hide file tree
Showing 3 changed files with 328 additions and 0 deletions.
15 changes: 15 additions & 0 deletions Exec/science/subch_planar/Problem_Derive.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,3 +8,18 @@ void ca_dergradpoverp1
(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*/, int /*level*/);

void ca_dergradpx
(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*/, int /*level*/);

void ca_dergradpy
(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*/, int /*level*/);

void ca_dergradrhooverrho
(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*/, int /*level*/);
304 changes: 304 additions & 0 deletions Exec/science/subch_planar/Problem_Derive.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -210,6 +210,7 @@ void ca_dergradpoverp(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*nco

}


void ca_dergradpoverp1(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/,
const FArrayBox& datfab, const Geometry& geomdata,
Real /*time*/, const int* /*bcrec*/, int /*level*/)
Expand Down Expand Up @@ -411,3 +412,306 @@ void ca_dergradpoverp1(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*nc
});

}


void ca_dergradpx(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/,
const FArrayBox& datfab, const Geometry& geomdata,
Real /*time*/, const int* /*bcrec*/, int /*level*/)
{

// Compute grad p / p . xhat

const auto dx = geomdata.CellSizeArray();
const int coord_type = geomdata.Coord();

auto const dat = datfab.array();
auto const der = derfab.array();

#if AMREX_SPACEDIM == 3
amrex::Error("3D not supported");
#endif
#if AMREX_SPACEDIM == 1
amrex::Error("1D not supported");
#endif

Real dxinv = 1.0_rt / dx[0];
Real dyinv = 1.0_rt / dx[1];

amrex::ParallelFor(bx,
[=] AMREX_GPU_HOST_DEVICE (int i, int j, int k)
{

// get the pressure
eos_t eos_state;

Real rhoInv = 1.0 / dat(i,j,k,URHO);

eos_state.rho = dat(i,j,k,URHO);
eos_state.T = dat(i,j,k,UTEMP);
eos_state.e = dat(i,j,k,UEINT) * rhoInv;
for (int n = 0; n < NumSpec; n++) {
eos_state.xn[n] = dat(i,j,k,UFS+n) * rhoInv;
}
#if NAUX_NET > 0
for (int n = 0; n < NumAux; n++) {
eos_state.aux[n] = dat(i,j,k,UFX+n) * rhoInv;
}
#endif

if (eos_state.e <= 0.0_rt) {
eos(eos_input_rt, eos_state);
} else {
eos(eos_input_re, eos_state);
}
Real p_zone = eos_state.p;


// we need to compute p in the full stencil

Real p_ip1{};
{
Real lrhoInv = 1.0 / dat(i+1,j,k,URHO);

eos_state.rho = dat(i+1,j,k,URHO);
eos_state.T = dat(i+1,j,k,UTEMP);
eos_state.e = dat(i+1,j,k,UEINT) * lrhoInv;
for (int n = 0; n < NumSpec; n++) {
eos_state.xn[n] = dat(i+1,j,k,UFS+n) * lrhoInv;
}
#if NAUX_NET > 0
for (int n = 0; n < NumAux; n++) {
eos_state.aux[n] = dat(i+1,j,k,UFX+n) * lrhoInv;
}
#endif

if (eos_state.e <= 0.0_rt) {
eos(eos_input_rt, eos_state);
} else {
eos(eos_input_re, eos_state);
}
p_ip1 = eos_state.p;
}

Real p_im1{};
{
Real lrhoInv = 1.0 / dat(i-1,j,k,URHO);

eos_state.rho = dat(i-1,j,k,URHO);
eos_state.T = dat(i-1,j,k,UTEMP);
eos_state.e = dat(i-1,j,k,UEINT) * lrhoInv;
for (int n = 0; n < NumSpec; n++) {
eos_state.xn[n] = dat(i-1,j,k,UFS+n) * lrhoInv;
}
#if NAUX_NET > 0
for (int n = 0; n < NumAux; n++) {
eos_state.aux[n] = dat(i-1,j,k,UFX+n) * lrhoInv;
}
#endif

if (eos_state.e <= 0.0_rt) {
eos(eos_input_rt, eos_state);
} else {
eos(eos_input_re, eos_state);
}
p_im1 = eos_state.p;
}

Real dp_x = 0.5_rt * (p_ip1 - p_im1);

der(i,j,k,0) = std::abs(dp_x) / p_zone;

});

}


void ca_dergradpy(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/,
const FArrayBox& datfab, const Geometry& geomdata,
Real /*time*/, const int* /*bcrec*/, int /*level*/)
{

// compute grad p / p . yhat

const auto dx = geomdata.CellSizeArray();
const int coord_type = geomdata.Coord();

auto const dat = datfab.array();
auto const der = derfab.array();

#if AMREX_SPACEDIM == 3
amrex::Error("3D not supported");
#endif
#if AMREX_SPACEDIM == 1
amrex::Error("1D not supported");
#endif

Real dxinv = 1.0_rt / dx[0];
Real dyinv = 1.0_rt / dx[1];


amrex::ParallelFor(bx,
[=] AMREX_GPU_HOST_DEVICE (int i, int j, int k)
{

// get the pressure
eos_t eos_state;

Real rhoInv = 1.0 / dat(i,j,k,URHO);

eos_state.rho = dat(i,j,k,URHO);
eos_state.T = dat(i,j,k,UTEMP);
eos_state.e = dat(i,j,k,UEINT) * rhoInv;
for (int n = 0; n < NumSpec; n++) {
eos_state.xn[n] = dat(i,j,k,UFS+n) * rhoInv;
}
#if NAUX_NET > 0
for (int n = 0; n < NumAux; n++) {
eos_state.aux[n] = dat(i,j,k,UFX+n) * rhoInv;
}
#endif

if (eos_state.e <= 0.0_rt) {
eos(eos_input_rt, eos_state);
} else {
eos(eos_input_re, eos_state);
}
Real p_zone = eos_state.p;



// we need to compute p in the full stencil

Real p_jp1{};
{
Real lrhoInv = 1.0 / dat(i,j+1,k,URHO);

eos_state.rho = dat(i,j+1,k,URHO);
eos_state.T = dat(i,j+1,k,UTEMP);
eos_state.e = dat(i,j+1,k,UEINT) * lrhoInv;
for (int n = 0; n < NumSpec; n++) {
eos_state.xn[n] = dat(i,j+1,k,UFS+n) * lrhoInv;
}
#if NAUX_NET > 0
for (int n = 0; n < NumAux; n++) {
eos_state.aux[n] = dat(i,j+1,k,UFX+n) * lrhoInv;
}
#endif

if (eos_state.e <= 0.0_rt) {
eos(eos_input_rt, eos_state);
} else {
eos(eos_input_re, eos_state);
}
p_jp1 = eos_state.p;
}

Real p_jm1{};
{
Real lrhoInv = 1.0 / dat(i,j-1,k,URHO);

eos_state.rho = dat(i,j-1,k,URHO);
eos_state.T = dat(i,j-1,k,UTEMP);
eos_state.e = dat(i,j-1,k,UEINT) * lrhoInv;
for (int n = 0; n < NumSpec; n++) {
eos_state.xn[n] = dat(i,j-1,k,UFS+n) * lrhoInv;
}
#if NAUX_NET > 0
for (int n = 0; n < NumAux; n++) {
eos_state.aux[n] = dat(i,j-1,k,UFX+n) * lrhoInv;
}
#endif

if (eos_state.e <= 0.0_rt) {
eos(eos_input_rt, eos_state);
} else {
eos(eos_input_re, eos_state);
}
p_jm1 = eos_state.p;
}

Real dp_y = 0.5_rt * (p_jp1 - p_jm1);

der(i,j,k,0) = std::abs(dp_y) / p_zone;

});

}


void ca_dergradrhooverrho(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/,
const FArrayBox& datfab, const Geometry& geomdata,
Real /*time*/, const int* /*bcrec*/, int /*level*/)
{

// compute grad rho / rho

const auto dx = geomdata.CellSizeArray();
const int coord_type = geomdata.Coord();

auto const dat = datfab.array();
auto const der = derfab.array();

#if AMREX_SPACEDIM == 3
amrex::Error("3D not supported");
#endif
#if AMREX_SPACEDIM == 1
amrex::Error("1D not supported");
#endif

Real dxinv = 1.0_rt / dx[0];
Real dyinv = 1.0_rt / dx[1];

amrex::ParallelFor(bx,
[=] AMREX_GPU_HOST_DEVICE (int i, int j, int k)
{

Real div_u = 0.0_rt;

Real up = dat(i+1,j,k,UMX) / dat(i+1,j,k,URHO);
Real um = dat(i-1,j,k,UMX) / dat(i-1,j,k,URHO);
Real u0 = dat(i,j,k,UMX) / dat(i,j,k,URHO);

Real vp = dat(i,j+1,k,UMY) / dat(i,j+1,k,URHO);
Real vm = dat(i,j-1,k,UMY) / dat(i,j-1,k,URHO);
Real v0 = dat(i,j,k,UMY) / dat(i,j,k,URHO);

// construct div{U}
if (coord_type == 0) {

// Cartesian
div_u += 0.5_rt * (up - um) * dxinv;
div_u += 0.5_rt * (vp - vm) * dyinv;

} else if (coord_type == 1) {

// r-z
Real rc = (i + 0.5_rt) * dx[0];
Real rm = (i - 1 + 0.5_rt) * dx[0];
Real rp = (i + 1 + 0.5_rt) * dx[0];

div_u += 0.5_rt * (rp * up - rm * um) / (rc * dx[0]) +
0.5_rt * (vp - vm) * dyinv;

#ifndef AMREX_USE_GPU
} else {
amrex::Error("ERROR: invalid coord_type in shock");
#endif
}


Real drho_x = 0.5_rt * (dat(i+1,j,k,URHO) - dat(i-1,j,k,URHO));
Real drho_y = 0.5_rt * (dat(i,j+1,k,URHO) - dat(i,j-1,k,URHO));

Real vel = std::sqrt(u0 * u0 + v0 * v0);

Real gradrhodx_over_rho{0.0_rt};
if (vel != 0.0) {
gradrhodx_over_rho = std::abs(drho_x * u0 + drho_y * v0) / vel;
}
gradrhodx_over_rho /= dat(i,j,k,URHO);

der(i,j,k,0) = gradrhodx_over_rho;

});

}
9 changes: 9 additions & 0 deletions Exec/science/subch_planar/Problem_Derives.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,3 +6,12 @@

derive_lst.add("gradp_over_p1", IndexType::TheCellType(), 1, ca_dergradpoverp1, grow_box_by_one);
derive_lst.addComponent("gradp_over_p1", desc_lst, State_Type, URHO, NUM_STATE);

derive_lst.add("gradp_x", IndexType::TheCellType(), 1, ca_dergradpx, grow_box_by_one);
derive_lst.addComponent("gradp_x", desc_lst, State_Type, URHO, NUM_STATE);

derive_lst.add("gradp_y", IndexType::TheCellType(), 1, ca_dergradpy, grow_box_by_one);
derive_lst.addComponent("gradp_y", desc_lst, State_Type, URHO, NUM_STATE);

derive_lst.add("gradrho_over_rho", IndexType::TheCellType(), 1, ca_dergradrhooverrho, grow_box_by_one);
derive_lst.addComponent("gradrho_over_rho", desc_lst, State_Type, URHO, NUM_STATE);

0 comments on commit be9fc30

Please sign in to comment.