Skip to content

Commit

Permalink
update primitive variable reconstruction with spherical 2d geometry. (#…
Browse files Browse the repository at this point in the history 
…2964)

this adds both the geometry source terms to the interface tracing and makes sure that the d/dtheta terms have the r weighting
  • Loading branch information
@zhichen3
zhichen3 authored Sep 26, 2024
1 parent ae6b17c commit 12666ac
Show file tree
Hide file tree
Showing 2 changed files with 85 additions and 52 deletions.
65 changes: 43 additions & 22 deletions Source/hydro/reconstruction.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -87,41 +87,55 @@ void
add_geometric_rho_source(amrex::Array4<amrex::Real const> const& q_arr,
amrex::Array4<amrex::Real const> const& dloga,
const int i, const int j, const int k,
amrex::Real* s) {
const int ncomp, amrex::Real* s) {

using namespace reconstruction;

// For idir == 0, i.e. r-direction:
// this takes the form: -alpha rho u / r
// where alpha = 1 for cylindrical and 2 for spherical
// note dloga(idir==0) = alpha/r

// note: this is assumed to be working only in the x-direction
// For idir == 1,
// i.e. theta-direction for spherical and z-direction for cylindrical
// this takes the form: -rho v cot(theta) / r for spherical and 0 for cylindrical
// note: dloga(idir==1) = cot(theta)/r for spherical and 0 for cylindrical.

s[im2] += -dloga(i-2,j,k) * q_arr(i-2,j,k,QRHO) * q_arr(i-2,j,k,QU);
s[im1] += -dloga(i-1,j,k) * q_arr(i-1,j,k,QRHO) * q_arr(i-1,j,k,QU);
s[i0] += -dloga(i,j,k) * q_arr(i,j,k,QRHO) * q_arr(i,j,k,QU);
s[ip1] += -dloga(i+1,j,k) * q_arr(i+1,j,k,QRHO) * q_arr(i+1,j,k,QU);
s[ip2] += -dloga(i+2,j,k) * q_arr(i+2,j,k,QRHO) * q_arr(i+2,j,k,QU);
// ncomp should be QU for idir == 0 and QV for idir == 1

s[im2] += -dloga(i-2,j,k) * q_arr(i-2,j,k,QRHO) * q_arr(i-2,j,k,ncomp);
s[im1] += -dloga(i-1,j,k) * q_arr(i-1,j,k,QRHO) * q_arr(i-1,j,k,ncomp);
s[i0] += -dloga(i,j,k) * q_arr(i,j,k,QRHO) * q_arr(i,j,k,ncomp);
s[ip1] += -dloga(i+1,j,k) * q_arr(i+1,j,k,QRHO) * q_arr(i+1,j,k,ncomp);
s[ip2] += -dloga(i+2,j,k) * q_arr(i+2,j,k,QRHO) * q_arr(i+2,j,k,ncomp);
}

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
void
add_geometric_rhoe_source(amrex::Array4<amrex::Real const> const& q_arr,
amrex::Array4<amrex::Real const> const& dloga,
const int i, const int j, const int k,
amrex::Real* s) {
const int ncomp, amrex::Real* s) {

using namespace reconstruction;

// this takes the form: -alpha (rho e + p) u / r
// For idir == 0, i.e. r-direction:
// this takes the form: -alpha (rho e + p ) u / r
// where alpha = 1 for cylindrical and 2 for spherical
// note dloga(idir==0) = alpha/r

// For idir == 1,
// i.e. theta-direction for spherical and z-direction for cylindrical
// this takes the form: - (rho e + p) v cot(theta) / r for spherical and 0 for cylindrical
// note: dloga(idir==1) = cot(theta)/r for spherical and 0 for cylindrical.

// note: this is assumed to be working only in the x-direction
// ncomp should be QU for idir == 0 and QV for idir == 1

s[im2] += -dloga(i-2,j,k) * (q_arr(i-2,j,k,QREINT) + q_arr(i-2,j,k,QPRES)) * q_arr(i-2,j,k,QU);
s[im1] += -dloga(i-1,j,k) * (q_arr(i-1,j,k,QREINT) + q_arr(i-1,j,k,QPRES)) * q_arr(i-1,j,k,QU);
s[i0] += -dloga(i,j,k) * (q_arr(i,j,k,QREINT) + q_arr(i,j,k,QPRES)) * q_arr(i,j,k,QU);
s[ip1] += -dloga(i+1,j,k) * (q_arr(i+1,j,k,QREINT) + q_arr(i+1,j,k,QPRES)) * q_arr(i+1,j,k,QU);
s[ip2] += -dloga(i+2,j,k) * (q_arr(i+2,j,k,QREINT) + q_arr(i+2,j,k,QPRES)) * q_arr(i+2,j,k,QU);
s[im2] += -dloga(i-2,j,k) * (q_arr(i-2,j,k,QREINT) + q_arr(i-2,j,k,QPRES)) * q_arr(i-2,j,k,ncomp);
s[im1] += -dloga(i-1,j,k) * (q_arr(i-1,j,k,QREINT) + q_arr(i-1,j,k,QPRES)) * q_arr(i-1,j,k,ncomp);
s[i0] += -dloga(i,j,k) * (q_arr(i,j,k,QREINT) + q_arr(i,j,k,QPRES)) * q_arr(i,j,k,ncomp);
s[ip1] += -dloga(i+1,j,k) * (q_arr(i+1,j,k,QREINT) + q_arr(i+1,j,k,QPRES)) * q_arr(i+1,j,k,ncomp);
s[ip2] += -dloga(i+2,j,k) * (q_arr(i+2,j,k,QREINT) + q_arr(i+2,j,k,QPRES)) * q_arr(i+2,j,k,ncomp);
}

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
Expand All @@ -130,20 +144,27 @@ add_geometric_p_source(amrex::Array4<amrex::Real const> const& q_arr,
amrex::Array4<amrex::Real const> const& qaux_arr,
amrex::Array4<amrex::Real const> const& dloga,
const int i, const int j, const int k,
amrex::Real* s) {
const int ncomp, amrex::Real* s) {

using namespace reconstruction;

// For idir == 0, i.e. r-direction:
// this takes the form: -alpha Gamma1 p u / r
// where alpha = 1 for cylindrical and 2 for spherical
// note: dloga(idir==0) = alpha/r

// For idir == 1,
// i.e. theta-direction for spherical and z-direction for cylindrical
// this takes the form: - Gamma1 p v cot(theta) / r for spherical and 0 for cylindrical
// note: dloga(idir==1) = cot(theta)/r for spherical and 0 for cylindrical.

// note: this is assumed to be working only in the x-direction
// ncomp should be QU for idir == 0 and QV for idir == 1

s[im2] += -dloga(i-2,j,k) * q_arr(i-2,j,k,QPRES) * qaux_arr(i-2,j,k,QGAMC) * q_arr(i-2,j,k,QU);
s[im1] += -dloga(i-1,j,k) * q_arr(i-1,j,k,QPRES) * qaux_arr(i-1,j,k,QGAMC) * q_arr(i-1,j,k,QU);
s[i0] += -dloga(i,j,k) * q_arr(i,j,k,QPRES) * qaux_arr(i,j,k,QGAMC) * q_arr(i,j,k,QU);
s[ip1] += -dloga(i+1,j,k) * q_arr(i+1,j,k,QPRES) * qaux_arr(i+1,j,k,QGAMC) * q_arr(i+1,j,k,QU);
s[ip2] += -dloga(i+2,j,k) * q_arr(i+2,j,k,QPRES) * qaux_arr(i+2,j,k,QGAMC) * q_arr(i+2,j,k,QU);
s[im2] += -dloga(i-2,j,k) * q_arr(i-2,j,k,QPRES) * qaux_arr(i-2,j,k,QGAMC) * q_arr(i-2,j,k,ncomp);
s[im1] += -dloga(i-1,j,k) * q_arr(i-1,j,k,QPRES) * qaux_arr(i-1,j,k,QGAMC) * q_arr(i-1,j,k,ncomp);
s[i0] += -dloga(i,j,k) * q_arr(i,j,k,QPRES) * qaux_arr(i,j,k,QGAMC) * q_arr(i,j,k,ncomp);
s[ip1] += -dloga(i+1,j,k) * q_arr(i+1,j,k,QPRES) * qaux_arr(i+1,j,k,QGAMC) * q_arr(i+1,j,k,ncomp);
s[ip2] += -dloga(i+2,j,k) * q_arr(i+2,j,k,QPRES) * qaux_arr(i+2,j,k,QGAMC) * q_arr(i+2,j,k,ncomp);
}


Expand Down
72 changes: 42 additions & 30 deletions Source/hydro/trace_ppm.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,10 +55,10 @@ Castro::trace_ppm(const Box& bx,


const auto dx = geom.CellSizeArray();
const auto problo = geom.ProbLoArray();
const int coord = geom.Coord();

Real hdt = 0.5_rt * dt;
Real dtdx = dt / dx[idir];

#ifndef AMREX_USE_GPU
auto lo = bx.loVect3d();
Expand All @@ -83,8 +83,11 @@ Castro::trace_ppm(const Box& bx,
for (int n = 0; n < NQSRC; n++) {
do_source_trace[n] = 0;

// geometric source terms in r-direction need tracing
if (coord > 0 && idir == 0 && (n == QRHO || n == QPRES || n == QREINT)) {
// geometric source terms in r-direction for nonCartesian coordinate
// or theta-direction in spherical coordinate need tracing

if ((coord == 2 || (idir == 0 && coord == 1))
&& (n == QRHO || n == QPRES || n == QREINT)) {
do_source_trace[n] = 1;
continue;
}
Expand Down Expand Up @@ -154,6 +157,15 @@ Castro::trace_ppm(const Box& bx,
amrex::ParallelFor(bx,
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
Real dtdL = dt / dx[idir];
Real dL = dx[idir];

// Want dt/(rdtheta) instead of dt/dtheta for 2d Spherical
if (coord == 2 && idir == 1) {
Real r = problo[0] + static_cast<Real>(i + 0.5_rt) * dx[0];
dL = r * dx[1];
dtdL = dt / dL;
}

Real cc = qaux_arr(i,j,k,QC);
Real un = q_arr(i,j,k,QUN);
Expand Down Expand Up @@ -196,7 +208,7 @@ Castro::trace_ppm(const Box& bx,

load_stencil(q_arr, idir, i, j, k, QRHO, s);
ppm_reconstruct(s, flat, sm, sp);
ppm_int_profile(sm, sp, s[i0], un, cc, dtdx, Ip_rho, Im_rho);
ppm_int_profile(sm, sp, s[i0], un, cc, dtdL, Ip_rho, Im_rho);

// reconstruct normal velocity

Expand All @@ -207,8 +219,8 @@ Castro::trace_ppm(const Box& bx,

load_stencil(q_arr, idir, i, j, k, QUN, s);
ppm_reconstruct(s, flat, sm, sp);
ppm_int_profile_single(sm, sp, s[i0], un-cc, dtdx, Ip_un_0, Im_un_0);
ppm_int_profile_single(sm, sp, s[i0], un+cc, dtdx, Ip_un_2, Im_un_2);
ppm_int_profile_single(sm, sp, s[i0], un-cc, dtdL, Ip_un_0, Im_un_0);
ppm_int_profile_single(sm, sp, s[i0], un+cc, dtdL, Ip_un_2, Im_un_2);

// reconstruct pressure

Expand All @@ -224,12 +236,12 @@ Castro::trace_ppm(const Box& bx,
load_stencil(q_arr, idir, i, j, k, QRHO, trho);
load_stencil(srcQ, idir, i, j, k, QUN, src);

ppm_reconstruct_pslope(trho, s, src, flat, dx[idir], sm, sp);
ppm_reconstruct_pslope(trho, s, src, flat, dL, sm, sp);

} else {
ppm_reconstruct(s, flat, sm, sp);
}
ppm_int_profile(sm, sp, s[i0], un, cc, dtdx, Ip_p, Im_p);
ppm_int_profile(sm, sp, s[i0], un, cc, dtdL, Ip_p, Im_p);

// reconstruct rho e

Expand All @@ -238,7 +250,7 @@ Castro::trace_ppm(const Box& bx,

load_stencil(q_arr, idir, i, j, k, QREINT, s);
ppm_reconstruct(s, flat, sm, sp);
ppm_int_profile(sm, sp, s[i0], un, cc, dtdx, Ip_rhoe, Im_rhoe);
ppm_int_profile(sm, sp, s[i0], un, cc, dtdL, Ip_rhoe, Im_rhoe);

// reconstruct transverse velocities

Expand All @@ -249,11 +261,11 @@ Castro::trace_ppm(const Box& bx,

load_stencil(q_arr, idir, i, j, k, QUT, s);
ppm_reconstruct(s, flat, sm, sp);
ppm_int_profile_single(sm, sp, s[i0], un, dtdx, Ip_ut_1, Im_ut_1);
ppm_int_profile_single(sm, sp, s[i0], un, dtdL, Ip_ut_1, Im_ut_1);

load_stencil(q_arr, idir, i, j, k, QUTT, s);
ppm_reconstruct(s, flat, sm, sp);
ppm_int_profile_single(sm, sp, s[i0], un, dtdx, Ip_utt_1, Im_utt_1);
ppm_int_profile_single(sm, sp, s[i0], un, dtdL, Ip_utt_1, Im_utt_1);

// gamma_c

Expand All @@ -264,8 +276,8 @@ Castro::trace_ppm(const Box& bx,

load_stencil(qaux_arr, idir, i, j, k, QGAMC, s);
ppm_reconstruct(s, flat, sm, sp);
ppm_int_profile_single(sm, sp, s[i0], un-cc, dtdx, Ip_gc_0, Im_gc_0);
ppm_int_profile_single(sm, sp, s[i0], un+cc, dtdx, Ip_gc_2, Im_gc_2);
ppm_int_profile_single(sm, sp, s[i0], un-cc, dtdL, Ip_gc_0, Im_gc_0);
ppm_int_profile_single(sm, sp, s[i0], un+cc, dtdL, Ip_gc_2, Im_gc_2);


// source terms -- we only trace if the terms in the stencil are non-zero
Expand All @@ -279,7 +291,7 @@ Castro::trace_ppm(const Box& bx,
#ifndef AMREX_USE_GPU
do_trace = do_source_trace[QRHO];
#else
if (idir == 0 && coord > 0) {
if (coord == 2 || (idir == 0 && coord == 1)) {
do_trace = 1;
} else {
do_trace = check_trace_source(srcQ, idir, i, j, k, QRHO);
Expand All @@ -289,12 +301,12 @@ Castro::trace_ppm(const Box& bx,
if (do_trace) {
load_stencil(srcQ, idir, i, j, k, QRHO, s);
#if AMREX_SPACEDIM <= 2
if (idir == 0 && coord > 0) {
add_geometric_rho_source(q_arr, dloga, i, j, k, s);
if (coord == 2 || (idir == 0 && coord == 1)) {
add_geometric_rho_source(q_arr, dloga, i, j, k, QUN, s);
}
#endif
ppm_reconstruct(s, flat, sm, sp);
ppm_int_profile(sm, sp, s[i0], un, cc, dtdx, Ip_src_rho, Im_src_rho);
ppm_int_profile(sm, sp, s[i0], un, cc, dtdL, Ip_src_rho, Im_src_rho);
}

// normal velocity
Expand All @@ -313,8 +325,8 @@ Castro::trace_ppm(const Box& bx,
if (do_trace) {
load_stencil(srcQ, idir, i, j, k, QUN, s);
ppm_reconstruct(s, flat, sm, sp);
ppm_int_profile_single(sm, sp, s[i0], un-cc, dtdx, Ip_src_un_0, Im_src_un_0);
ppm_int_profile_single(sm, sp, s[i0], un+cc, dtdx, Ip_src_un_2, Im_src_un_2);
ppm_int_profile_single(sm, sp, s[i0], un-cc, dtdL, Ip_src_un_0, Im_src_un_0);
ppm_int_profile_single(sm, sp, s[i0], un+cc, dtdL, Ip_src_un_2, Im_src_un_2);
}

// pressure
Expand All @@ -325,7 +337,7 @@ Castro::trace_ppm(const Box& bx,
#ifndef AMREX_USE_GPU
do_trace = do_source_trace[QPRES];
#else
if (idir == 0 && coord > 0) {
if (coord == 2 || (idir == 0 && coord == 1)) {
do_trace = 1;
} else {
do_trace = check_trace_source(srcQ, idir, i, j, k, QPRES);
Expand All @@ -335,12 +347,12 @@ Castro::trace_ppm(const Box& bx,
if (do_trace) {
load_stencil(srcQ, idir, i, j, k, QPRES, s);
#if AMREX_SPACEDIM <= 2
if (idir == 0 && coord > 0) {
add_geometric_p_source(q_arr, qaux_arr, dloga, i, j, k, s);
if (coord == 2 || (idir == 0 && coord == 1)) {
add_geometric_p_source(q_arr, qaux_arr, dloga, i, j, k, QUN, s);
}
#endif
ppm_reconstruct(s, flat, sm, sp);
ppm_int_profile(sm, sp, s[i0], un, cc, dtdx, Ip_src_p, Im_src_p);
ppm_int_profile(sm, sp, s[i0], un, cc, dtdL, Ip_src_p, Im_src_p);
}

// rho e
Expand All @@ -351,7 +363,7 @@ Castro::trace_ppm(const Box& bx,
#ifndef AMREX_USE_GPU
do_trace = do_source_trace[QREINT];
#else
if (idir == 0 && coord > 0) {
if (coord == 2 || (idir == 0 && coord == 1)) {
do_trace = 1;
} else {
do_trace = check_trace_source(srcQ, idir, i, j, k, QREINT);
Expand All @@ -361,12 +373,12 @@ Castro::trace_ppm(const Box& bx,
if (do_trace) {
load_stencil(srcQ, idir, i, j, k, QREINT, s);
#if AMREX_SPACEDIM <= 2
if (idir == 0 && coord > 0) {
add_geometric_rhoe_source(q_arr, dloga, i, j, k, s);
if (coord == 2 || (idir == 0 && coord == 1)) {
add_geometric_rhoe_source(q_arr, dloga, i, j, k, QUN, s);
}
#endif
ppm_reconstruct(s, flat, sm, sp);
ppm_int_profile(sm, sp, s[i0], un, cc, dtdx, Ip_src_rhoe, Im_src_rhoe);
ppm_int_profile(sm, sp, s[i0], un, cc, dtdL, Ip_src_rhoe, Im_src_rhoe);
}

// transverse velocities
Expand All @@ -383,7 +395,7 @@ Castro::trace_ppm(const Box& bx,
if (do_trace) {
load_stencil(srcQ, idir, i, j, k, QUT, s);
ppm_reconstruct(s, flat, sm, sp);
ppm_int_profile_single(sm, sp, s[i0], un, dtdx, Ip_src_ut_1, Im_src_ut_1);
ppm_int_profile_single(sm, sp, s[i0], un, dtdL, Ip_src_ut_1, Im_src_ut_1);
}

Real Ip_src_utt_1 = 0.0_rt;
Expand All @@ -398,7 +410,7 @@ Castro::trace_ppm(const Box& bx,
if (do_trace) {
load_stencil(srcQ, idir, i, j, k, QUTT, s);
ppm_reconstruct(s, flat, sm, sp);
ppm_int_profile_single(sm, sp, s[i0], un, dtdx, Ip_src_utt_1, Im_src_utt_1);
ppm_int_profile_single(sm, sp, s[i0], un, dtdL, Ip_src_utt_1, Im_src_utt_1);
}


Expand All @@ -416,7 +428,7 @@ Castro::trace_ppm(const Box& bx,

load_passive_stencil(U_arr, rho_inv_arr, idir, i, j, k, nc, s);
ppm_reconstruct(s, flat, sm, sp);
ppm_int_profile_single(sm, sp, s[i0], un, dtdx, Ip_passive, Im_passive);
ppm_int_profile_single(sm, sp, s[i0], un, dtdL, Ip_passive, Im_passive);

// Plus state on face i

Expand Down

0 comments on commit 12666ac

Please sign in to comment.