Skip to content

Commit

Permalink
update primitive vairable geometric pres source and conserved variabl…
Browse files Browse the repository at this point in the history 
…e geometric pres source
  • Loading branch information
@zhichen3
zhichen3 committed Nov 18, 2024
1 parent 265ea33 commit 7bf9636
Show file tree
Hide file tree
Showing 3 changed files with 94 additions and 40 deletions.
37 changes: 37 additions & 0 deletions Source/hydro/reconstruction.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -110,6 +110,43 @@ add_geometric_rho_source(amrex::Array4<amrex::Real const> const& q_arr,
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
delete_UN_geometric_pres_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) {

using namespace reconstruction;

// Depending on the direction of reconstruction, delete corresponding geometric
// pressure source in the velocity equations.

// If idir == 0, i.e. r-direction, we should delete the geometric pressure source
// that were obtained by grouping pressure into x (r) momentum flux.

// If idir == 1, i.e. theta-direction, delete the geometric pressure source in V
// that were obtained by grouping pressure into y (theta) momentum flux.

// For idir == 0, i.e. r-direction:
// this takes the form: -alpha p / 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: - p 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,QPRES);
s[im1] += -dloga(i-1,j,k) * q_arr(i-1,j,k,QPRES);
s[i0] += -dloga(i,j,k) * q_arr(i,j,k,QPRES);
s[ip1] += -dloga(i+1,j,k) * q_arr(i+1,j,k,QPRES);
s[ip2] += -dloga(i+2,j,k) * q_arr(i+2,j,k,QPRES);
}


AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
void
add_geometric_rhoe_source(amrex::Array4<amrex::Real const> const& q_arr,
Expand Down
20 changes: 18 additions & 2 deletions Source/hydro/trace_ppm.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -86,12 +86,19 @@ Castro::trace_ppm(const Box& bx,
// geometric source terms in r-direction for nonCartesian coordinate
// or theta-direction in spherical coordinate need tracing

if ((coord == 2 || (idir == 0 && coord == 1))
if ((idir == 0 && coord == 1)
&& (n == QRHO || n == QPRES || n == QREINT)) {
do_source_trace[n] = 1;
continue;
}

if (coord == 2 &&
(n == QRHO || n == QPRES || n == QREINT
n == QU || n == QV)) {
do_source_trace[n] = 1;
continue;
}

for (int k = lo[2]-2*dg2; k <= hi[2]+2*dg2; k++) {
for (int j = lo[1]-2*dg1; j <= hi[1]+2*dg1; j++) {
for (int i = lo[0]-2; i <= hi[0]+2; i++) {
Expand Down Expand Up @@ -337,11 +344,20 @@ Castro::trace_ppm(const Box& bx,
#ifndef AMREX_USE_GPU
do_trace = do_source_trace[QUN];
#else
do_trace = check_trace_source(srcQ, idir, i, j, k, QUN);
if (coord == 2 || (idir == 0 && coord = 1)) {
do_trace = 1;
} else {
do_trace = check_trace_source(srcQ, idir, i, j, k, QUN);
}
#endif

if (do_trace) {
load_stencil(srcQ, idir, i, j, k, QUN, s);
#if AMREX_SPACEDIM <= 2
if (coord == 2 || (idir == 0 && coord == 1)) {
delete_UN_geometric_pres_source(q_arr, dloga, i, j, k, s);
}
#endif
ppm_reconstruct(s, flat, sm, sp);
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);
Expand Down
77 changes: 39 additions & 38 deletions Source/sources/Castro_geom.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -159,41 +159,40 @@ Castro::fill_RZ_geom_source (Real time, Real dt, MultiFab& cons_state, MultiFab&

const Box& bx = mfi.tilebox();

Array4<Real const> const U_arr = cons_state.array(mfi);
// dlogaX = 1 / r for Cylindrical

Array4<Real const> const dlogaX = (dLogArea[0]).array(mfi);

Array4<Real const> const U_arr = cons_state.array(mfi);
Array4<Real> const src = geom_src.array(mfi);

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

Real rhoinv = 1.0_rt / U_arr(i,j,k,QRHO);

// radius for non-Cartesian
Real rinv = 1.0_rt / (prob_lo[0] + (static_cast<Real>(i) + 0.5_rt)*dx[0]);
Real rhoinv = 1.0_rt / U_arr(i,j,k,URHO);

// radial momentum: F = rho v_phi**2 / r
src(i,j,k,UMX) = U_arr(i,j,k,UMZ) * U_arr(i,j,k,UMZ) * rhoinv * rinv;
src(i,j,k,UMX) = U_arr(i,j,k,UMZ) * U_arr(i,j,k,UMZ) * dlogaX(i,j,k) * rhoinv ;

// azimuthal momentum: F = - rho v_r v_phi / r
src(i,j,k,UMZ) = - U_arr(i,j,k,UMX) * U_arr(i,j,k,UMZ) * rhoinv * rinv;
src(i,j,k,UMZ) = - U_arr(i,j,k,UMX) * U_arr(i,j,k,UMZ) * dlogaX(i,j,k) * rhoinv;

// We absorb pressure gradient into div Flux term,
// which causes an extra geometric pressure term. Address that here.

// Find pressure first:

// eos_rep_t eos_state;
// eos_state.rho = U_arr(i,j,k,URHO);
// eos_state.e = U_arr(i,j,k,UEINT) * rhoinv;
// for (int n = 0; n < NumSpec; n++) {
// eos_state.xn[n] = U_arr(i,j,k,UFS+n) * rhoinv;
// }
// eos(eos_input_re, eos_state);

// // radial geometric pres term: F = p/r
// src(i,j,k,UMX) = eos_state.p * rinv;
eos_rep_t eos_state;
eos_state.rho = U_arr(i,j,k,URHO);
eos_state.e = U_arr(i,j,k,UEINT) * rhoinv;
for (int n = 0; n < NumSpec; n++) {
eos_state.xn[n] = U_arr(i,j,k,UFS+n) * rhoinv;
}
eos(eos_input_re, eos_state);

// radial geometric pres term: F = p/r
src(i,j,k,UMX) = eos_state.p * dlogaX(i,j,k);
});
}
}
Expand All @@ -220,49 +219,51 @@ Castro::fill_RTheta_geom_source (Real time, Real dt, MultiFab& cons_state, Multi

const Box& bx = mfi.tilebox();

Array4<Real const> const U_arr = cons_state.array(mfi);
// dlogaX = 2 / r for Spherical
// dlogaY = cot(theta) / r for Spherical

Array4<Real const> const dlogaX = (dLogArea[0]).array(mfi);
Array4<Real const> const dlogaY = (dLogArea[1]).array(mfi);

Array4<Real const> const U_arr = cons_state.array(mfi);
Array4<Real> const src = geom_src.array(mfi);

amrex::ParallelFor(bx,
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
Real rhoinv = 1.0_rt / U_arr(i,j,k,QRHO);

// Cell-centered Spherical Radius and Theta
Real rinv = 1.0_rt / (prob_lo[0] + (static_cast<Real>(i) + 0.5_rt)*dx[0]);
Real cotTheta = cot(prob_lo[1] + (static_cast<Real>(j) + 0.5_rt)*dx[1]);
Real rhoinv = 1.0_rt / U_arr(i,j,k,URHO);
Real rinv = 0.5_rt * dlogaX(i,j,k);

// radial momentum: F = rho (v_theta**2 + v_phi**2) / r
src(i,j,k,UMX) = (U_arr(i,j,k,UMY) * U_arr(i,j,k,UMY) +
U_arr(i,j,k,UMZ) * U_arr(i,j,k,UMZ)) * rhoinv * rinv;

// Theta momentum F = rho v_phi**2 cot(theta) / r - rho v_r v_theta / r
src(i,j,k,UMY) = (U_arr(i,j,k,UMZ) * U_arr(i,j,k,UMZ) * cotTheta -
U_arr(i,j,k,UMX) * U_arr(i,j,k,UMY)) * rhoinv * rinv;
src(i,j,k,UMY) = (U_arr(i,j,k,UMZ) * U_arr(i,j,k,UMZ) * dlogaY(i,j,k) +
U_arr(i,j,k,UMX) * U_arr(i,j,k,UMY) * rinv) * rhoinv;

// Phi momentum: F = - rho v_r v_phi / r - rho v_theta v_phi cot(theta) / r
src(i,j,k,UMZ) = (- U_arr(i,j,k,UMY) * U_arr(i,j,k,UMZ) * cotTheta -
U_arr(i,j,k,UMX) * U_arr(i,j,k,UMZ)) * rhoinv * rinv;
src(i,j,k,UMZ) = (-U_arr(i,j,k,UMY) * U_arr(i,j,k,UMZ) * dlogaY(i,j,k) -
U_arr(i,j,k,UMX) * U_arr(i,j,k,UMZ) * rinv) * rhoinv;

// We absorb pressure gradient into div Flux term,
// which causes an extra geometric pressure term. Address that here.

// Find pressure first:

// eos_rep_t eos_state;
// eos_state.rho = U_arr(i,j,k,URHO);
// eos_state.e = U_arr(i,j,k,UEINT) * rhoinv;
// for (int n = 0; n < NumSpec; n++) {
// eos_state.xn[n] = U_arr(i,j,k,UFS+n) * rhoinv;
// }
// eos(eos_input_re, eos_state);
eos_rep_t eos_state;
eos_state.rho = U_arr(i,j,k,URHO);
eos_state.e = U_arr(i,j,k,UEINT) * rhoinv;
for (int n = 0; n < NumSpec; n++) {
eos_state.xn[n] = U_arr(i,j,k,UFS+n) * rhoinv;
}
eos(eos_input_re, eos_state);

// // radial geometric pres term: F = 2 p / r
// src(i,j,k,UMX) = 2.0_rt * eos_state.p * rinv;
// radial geometric pres term: F = 2 p / r
src(i,j,k,UMX) += eos_state.p * dlogaX(i,j,k);

// // Theta geometric pres term: F = cot(theta) p / r
// src(i,j,k,UMY) = cotTheta * eos_state.p * rinv;
// Theta geometric pres term: F = cot(theta) p / r
src(i,j,k,UMY) += eos_state.p * dlogaY(i,j,k);
});
}
}

0 comments on commit 7bf9636

Please sign in to comment.