Skip to content

Commit

Permalink
Merge branch 'development' into fix_geom_source_flag
Browse files Browse the repository at this point in the history
  • Loading branch information
@zingale
zingale authored Nov 21, 2024
2 parents abf1a66 + f02444d commit 1ef1731
Show file tree
Hide file tree
Showing 4 changed files with 104 additions and 41 deletions.
8 changes: 4 additions & 4 deletions Exec/science/xrb_spherical/analysis/slice.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -66,12 +66,12 @@ def slice(fname:str, field:str,
if field in ["x_velocity", "y_velocity", "z_velocity"]:
sp.set_cmap(field, "coolwarm")
elif field == "Temp":
sp.set_zlim(f, 5.e7, 2.5e9)
sp.set_cmap(f, "magma_r")
sp.set_zlim(field. 5.e7, 2.5e9)
sp.set_cmap(field, "magma_r")
elif field == "enuc":
sp.set_zlim(f, 1.e18, 1.e20)
sp.set_zlim(field, 1.e18, 1.e20)
elif field == "density":
sp.set_zlim(f, 1.e-3, 5.e8)
sp.set_zlim(field, 1.e-3, 5.e8)

# sp.annotate_text((0.05, 0.05), f"{currentTime.in_cgs():8.5f} s")
sp.save(f"{ds}_{loc}")
Expand Down
13 changes: 7 additions & 6 deletions Source/driver/Castro.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -851,12 +851,13 @@ Castro::buildMetrics ()
dLogArea[dir].clear();
geom.GetDLogA(dLogArea[dir],grids,dmap,dir,NUM_GROW);
}
for (int dir = AMREX_SPACEDIM; dir < 2; dir++)
{
dLogArea[dir].clear();
dLogArea[dir].define(grids, dmap, 1, 0);
dLogArea[dir].setVal(0.0);
}

#if (AMREX_SPACEDIM == 1)
dLogArea[1].clear();
dLogArea[1].define(grids, dmap, 1, 0);
dLogArea[1].setVal(0.0);
#endif

#endif

wall_time_start = 0.0;
Expand Down
69 changes: 54 additions & 15 deletions Source/hydro/reconstruction.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -103,11 +103,24 @@ add_geometric_rho_source(amrex::Array4<amrex::Real const> const& q_arr,

// 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);
if (ncomp == QU) {

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);

} else if (ncomp == QV) {

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

}

}

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
Expand All @@ -131,11 +144,24 @@ add_geometric_rhoe_source(amrex::Array4<amrex::Real const> const& q_arr,

// 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,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);
if (ncomp == 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);

} else if (ncomp == QV) {

s[im2] += -dloga(i,j-2,k) * (q_arr(i,j-2,k,QREINT) + q_arr(i,j-2,k,QPRES)) * q_arr(i,j-2,k,ncomp);
s[im1] += -dloga(i,j-1,k) * (q_arr(i,j-1,k,QREINT) + q_arr(i,j-1,k,QPRES)) * q_arr(i,j-1,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,j+1,k) * (q_arr(i,j+1,k,QREINT) + q_arr(i,j+1,k,QPRES)) * q_arr(i,j+1,k,ncomp);
s[ip2] += -dloga(i,j+2,k) * (q_arr(i,j+2,k,QREINT) + q_arr(i,j+2,k,QPRES)) * q_arr(i,j+2,k,ncomp);

}

}

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
Expand All @@ -160,11 +186,24 @@ add_geometric_p_source(amrex::Array4<amrex::Real const> const& q_arr,

// 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,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);

if (ncomp == 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);

} else if (ncomp == QV) {

s[im2] += -dloga(i,j-2,k) * q_arr(i,j-2,k,QPRES) * qaux_arr(i,j-2,k,QGAMC) * q_arr(i,j-2,k,ncomp);
s[im1] += -dloga(i,j-1,k) * q_arr(i,j-1,k,QPRES) * qaux_arr(i,j-1,k,QGAMC) * q_arr(i,j-1,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,j+1,k) * q_arr(i,j+1,k,QPRES) * qaux_arr(i,j+1,k,QGAMC) * q_arr(i,j+1,k,ncomp);
s[ip2] += -dloga(i,j+2,k) * q_arr(i,j+2,k,QPRES) * qaux_arr(i,j+2,k,QGAMC) * q_arr(i,j+2,k,ncomp);
}

}


Expand Down
55 changes: 39 additions & 16 deletions Source/hydro/trace_plm.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,6 +33,8 @@ Castro::trace_plm(const Box& bx, const int idir,
// vbx is the valid box (no ghost cells)

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

const int* lo_bc = phys_bc.lo();
const int* hi_bc = phys_bc.hi();
Expand All @@ -43,8 +45,6 @@ Castro::trace_plm(const Box& bx, const int idir,
const auto domlo = geom.Domain().loVect3d();
const auto domhi = geom.Domain().hiVect3d();

const Real dtdx = dt/dx[idir];

auto vlo = vbx.loVect3d();
auto vhi = vbx.hiVect3d();

Expand Down Expand Up @@ -100,6 +100,16 @@ Castro::trace_plm(const Box& bx, const int idir,
[=] 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;
}

bool lo_bc_test = lo_symm && ((idir == 0 && i == domlo[0]) ||
(idir == 1 && j == domlo[1]) ||
(idir == 2 && k == domlo[2]));
Expand Down Expand Up @@ -166,7 +176,7 @@ Castro::trace_plm(const Box& bx, const int idir,
load_stencil(srcQ, idir, i, j, k, QUN, src);

Real dp = dq[IEIGN_P];
pslope(trho, s, src, flat, lo_bc_test, hi_bc_test, dx[idir], dp);
pslope(trho, s, src, flat, lo_bc_test, hi_bc_test, dL, dp);
dq[IEIGN_P] = dp;

}
Expand All @@ -185,7 +195,7 @@ Castro::trace_plm(const Box& bx, const int idir,

// construct the right state on the i interface

Real ref_fac = 0.5_rt*(1.0_rt + dtdx*amrex::min(e[0], 0.0_rt));
Real ref_fac = 0.5_rt*(1.0_rt + dtdL*amrex::min(e[0], 0.0_rt));
Real rho_ref = rho - ref_fac*dq[IEIGN_RHO];
Real un_ref = un - ref_fac*dq[IEIGN_UN];
Real ut_ref = ut - ref_fac*dq[IEIGN_UT];
Expand All @@ -195,8 +205,8 @@ Castro::trace_plm(const Box& bx, const int idir,

// this is -(1/2) ( 1 + dt/dx lambda) (l . dq) r
Real trace_fac0 = 0.0_rt; // FOURTH*dtdx*(e(1) - e(1))*(1.0_rt - sign(1.0_rt, e(1)))
Real trace_fac1 = 0.25_rt*dtdx*(e[0] - e[1])*(1.0_rt - std::copysign(1.0_rt, e[1]));
Real trace_fac2 = 0.25_rt*dtdx*(e[0] - e[2])*(1.0_rt - std::copysign(1.0_rt, e[2]));
Real trace_fac1 = 0.25_rt*dtdL*(e[0] - e[1])*(1.0_rt - std::copysign(1.0_rt, e[1]));
Real trace_fac2 = 0.25_rt*dtdL*(e[0] - e[2])*(1.0_rt - std::copysign(1.0_rt, e[2]));

Real apright = trace_fac2*alphap;
Real amright = trace_fac0*alpham;
Expand Down Expand Up @@ -230,16 +240,16 @@ Castro::trace_plm(const Box& bx, const int idir,

// now construct the left state on the i+1 interface

ref_fac = 0.5_rt*(1.0_rt - dtdx*amrex::max(e[2], 0.0_rt));
ref_fac = 0.5_rt*(1.0_rt - dtdL*amrex::max(e[2], 0.0_rt));
rho_ref = rho + ref_fac*dq[IEIGN_RHO];
un_ref = un + ref_fac*dq[IEIGN_UN];
ut_ref = ut + ref_fac*dq[IEIGN_UT];
utt_ref = utt + ref_fac*dq[IEIGN_UTT];
p_ref = p + ref_fac*dq[IEIGN_P];
rhoe_ref = rhoe + ref_fac*dq[IEIGN_RE];

trace_fac0 = 0.25_rt*dtdx*(e[2] - e[0])*(1.0_rt + std::copysign(1.0_rt, e[0]));
trace_fac1 = 0.25_rt*dtdx*(e[2] - e[1])*(1.0_rt + std::copysign(1.0_rt, e[1]));
trace_fac0 = 0.25_rt*dtdL*(e[2] - e[0])*(1.0_rt + std::copysign(1.0_rt, e[0]));
trace_fac1 = 0.25_rt*dtdL*(e[2] - e[1])*(1.0_rt + std::copysign(1.0_rt, e[1]));
trace_fac2 = 0.0_rt; // FOURTH*dtdx*(e(3) - e(3))*(1.0_rt + sign(1.0_rt, e(3)))

Real apleft = trace_fac2*alphap;
Expand Down Expand Up @@ -301,29 +311,42 @@ Castro::trace_plm(const Box& bx, const int idir,
}

#if (AMREX_SPACEDIM < 3)
// geometry source terms -- these only apply to the x-states
if (idir == 0 && dloga(i,j,k) != 0.0_rt) {
Real courn = dtdx*(cc + std::abs(un));
// geometry source terms
// these only apply to the x-states for cylindrical and spherical
// or y-states for spherical

if (dloga(i,j,k) != 0.0_rt) {
Real courn = dtdL*(cc + std::abs(un));
Real eta = (1.0_rt-courn)/(cc*dt*std::abs(dloga(i,j,k)));
Real dlogatmp = amrex::min(eta, 1.0_rt)*dloga(i,j,k);
Real sourcr = -0.5_rt*dt*rho*dlogatmp*un;
Real sourcp = sourcr*csq;
Real source = sourcp*enth;

if (i <= vhi[0]) {
if (idir == 0 && i <= vhi[0]) {
qm(i+1,j,k,QRHO) += sourcr;
qm(i+1,j,k,QRHO) = amrex::max(qm(i+1,j,k,QRHO), lsmall_dens);
qm(i+1,j,k,QPRES) += sourcp;
qm(i+1,j,k,QREINT) += source;
}

if (i >= vlo[0]) {
if (idir == 1 && j <= vhi[1]) {
qm(i,j+1,k,QRHO) += sourcr;
qm(i,j+1,k,QRHO) = amrex::max(qm(i,j+1,k,QRHO), lsmall_dens);
qm(i,j+1,k,QPRES) += sourcp;
qm(i,j+1,k,QREINT) += source;
}

if ((idir == 0 && i >= vlo[0]) ||
(idir == 1 && j >= vlo[1])) {
qp(i,j,k,QRHO) += sourcr;
qp(i,j,k,QRHO) = amrex::max(qp(i,j,k,QRHO), lsmall_dens);
qp(i,j,k,QPRES) += sourcp;
qp(i,j,k,QREINT) += source;
}

}

#endif

for (int ipassive = 0; ipassive < npassive; ipassive++) {
Expand All @@ -340,12 +363,12 @@ Castro::trace_plm(const Box& bx, const int idir,
(idir == 1 && j >= vlo[1]) ||
(idir == 2 && k >= vlo[2])) {

Real spzero = un >= 0.0_rt ? -1.0_rt : un*dtdx;
Real spzero = un >= 0.0_rt ? -1.0_rt : un*dtdL;
qp(i,j,k,n) = s[i0] + 0.5_rt*(-1.0_rt - spzero)*dX;
}

// Left state
Real spzero = un >= 0.0_rt ? un*dtdx : 1.0_rt;
Real spzero = un >= 0.0_rt ? un*dtdL : 1.0_rt;
Real acmpleft = 0.5_rt*(1.0_rt - spzero )*dX;

if (idir == 0 && i <= vhi[0]) {
Expand Down

0 comments on commit 1ef1731

Please sign in to comment.