Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

update primitive variable reconstruction with spherical 2d geometry. #2964

Merged
merged 25 commits into from
Sep 26, 2024
Merged

update primitive variable reconstruction with spherical 2d geometry. #2964

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
25 commits
Select commit Hold shift + click to select a range
4eabeab
add geometric source terms to the PPM tracing
zingale May 30, 2023
2a44c10
finish implementation
zingale May 30, 2023
f0647a6
fix 3d compilation
zingale May 30, 2023
8516265
Merge branch 'development' into ppm_source_tracing
zingale Jun 7, 2023
d41376e
Merge branch 'development' into ppm_source_tracing
zingale Jun 17, 2023
dd08912
Merge branch 'development' into ppm_source_tracing
zingale Jun 17, 2023
12d552d
Merge branch 'development' into ppm_source_tracing
zingale Jun 18, 2023
7b0343c
Merge branch 'development' into ppm_source_tracing
zingale Jul 10, 2023
01f092c
Merge branch 'development' into ppm_source_tracing
zingale Nov 8, 2023
4c0b75d
Merge branch 'development' into ppm_source_tracing
zingale Jul 26, 2024
02e3a9e
Merge branch 'development' into ppm_source_tracing
zingale Jul 27, 2024
8828c4b
fix pressure source
zingale Jul 27, 2024
ec64d4f
Merge branch 'development' into ppm_source_tracing
zingale Aug 21, 2024
2c57bc8
remove unused var
zingale Aug 21, 2024
83e735b
Merge branch 'ppm_source_tracing' of github.com:zingale/Castro into p…
zingale Aug 21, 2024
8e914c3
update CI
zingale Aug 21, 2024
3a4662b
Merge branch 'development' into ppm_source_tracing
zingale Sep 16, 2024
267d57c
Merge branch 'development' into ppm_source_tracing
zhichen3 Sep 23, 2024
be16896
update trace ppm to be compatible with spherical 2d
zhichen3 Sep 23, 2024
a2794fb
missing comma
zhichen3 Sep 23, 2024
592fc9c
misspell
zhichen3 Sep 23, 2024
12803d5
fix compilation
zhichen3 Sep 23, 2024
0ea7f1b
delete comment
zhichen3 Sep 23, 2024
0517453
Merge branch 'development' into 2d_spherical_trace_geom
zhichen3 Sep 26, 2024
4d10d82
Merge branch 'development' into 2d_spherical_trace_geom
zhichen3 Sep 26, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
65 changes: 43 additions & 22 deletions Source/hydro/reconstruction.H
View file
Open in desktop
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
Open in desktop
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