Skip to content

Commit

Permalink
initialupdate
Browse files Browse the repository at this point in the history
  • Loading branch information
@zhichen3
zhichen3 committed Nov 13, 2024
1 parent b38a018 commit 84ffcff
Show file tree
Hide file tree
Showing 6 changed files with 54 additions and 82 deletions.
14 changes: 0 additions & 14 deletions Source/hydro/Castro_ctu.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -64,20 +64,6 @@ Castro::consup_hydro(const Box& bx,
pdu = 0.5 * pdu * volinv;

U_new(i,j,k,n) = U_new(i,j,k,n) - dt * pdu;

} else if (n == UMX && !mom_flux_has_p(0, 0, geomdata.Coord())) {
// Add gradp term to momentum equation -- only for axisymmetric
// coords (and only for the radial flux).

U_new(i,j,k,UMX) += - dt * (qx(i+1,j,k,GDPRES) - qx(i,j,k,GDPRES)) / geomdata.CellSize(0);

#if AMREX_SPACEDIM >= 2
} else if (n == UMY && !mom_flux_has_p(1, 1, geomdata.Coord())) {
// Add gradp term to polar(theta) momentum equation for Spherical 2D geometry

Real r = geomdata.ProbLo(0) + (static_cast<Real>(i) + 0.5_rt) * geomdata.CellSize(0);
U_new(i,j,k,UMY) += - dt * (qy(i,j+1,k,GDPRES) - qy(i,j,k,GDPRES)) / (r * geomdata.CellSize(1));
#endif
}
});
}
Expand Down
25 changes: 1 addition & 24 deletions Source/hydro/Castro_mol.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -259,23 +259,6 @@ Castro::mol_consup(const Box& bx, // NOLINT(readability-convert-member-function
flux2(i,j,k,n) * area2(i,j,k) - flux2(i,j,k+1,n) * area2(i,j,k+1) ) / vol(i,j,k);
#endif

#if AMREX_SPACEDIM <= 2
if (do_hydro == 1) {
if (n == UMX && !mom_flux_has_p(0, 0, coord)) {
// Add gradp term to radial momentum equation -- only for axisymmetric
// coords.

update(i,j,k,UMX) -= (q0(i+1,j,k,GDPRES) - q0(i,j,k,GDPRES)) / dx[0];

} else if (n == UMY && !mom_flux_has_p(1, 1, coord)) {
// Add gradp term to polar(theta) momentum equation for Spherical 2D geometry

Real r = prob_lo[0] + (static_cast<Real>(i) + 0.5_rt) * dx[0];
update(i,j,k,UMY) -= (q1(i,j+1,k,GDPRES) - q1(i,j,k,GDPRES)) / (r * dx[1]);
}
}
#endif

// this assumes that the species are at the end of the conserved state
if (n < NSRC) {
update(i,j,k,n) += srcU(i,j,k,n);
Expand Down Expand Up @@ -347,9 +330,6 @@ Castro::compute_flux_from_q(const Box& bx,
int iu, iv1, iv2;
int im1, im2, im3;

auto coord = geom.Coord();
auto mom_check = mom_flux_has_p(idir, idir, coord);

if (idir == 0) {
iu = QU;
iv1 = QV;
Expand Down Expand Up @@ -389,10 +369,7 @@ Castro::compute_flux_from_q(const Box& bx,
// Compute fluxes, order as conserved state (not q)
F(i,j,k,URHO) = qint(i,j,k,QRHO)*u_adv;

F(i,j,k,im1) = F(i,j,k,URHO)*qint(i,j,k,iu);
if (mom_check) {
F(i,j,k,im1) += qint(i,j,k,QPRES);
}
F(i,j,k,im1) = F(i,j,k,URHO)*qint(i,j,k,iu) + qint(i,j,k,QPRES);
F(i,j,k,im2) = F(i,j,k,URHO)*qint(i,j,k,iv1);
F(i,j,k,im3) = F(i,j,k,URHO)*qint(i,j,k,iv2);

Expand Down
22 changes: 4 additions & 18 deletions Source/hydro/HLL_solvers.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -79,18 +79,10 @@ namespace HLL {

F[URHO] = U[URHO]*u_flx;

F[UMX] = U[UMX]*u_flx;
F[UMX] = U[UMX]*u_flx + p;
F[UMY] = U[UMY]*u_flx;
F[UMZ] = U[UMZ]*u_flx;

auto mom_check = mom_flux_has_p(idir, idir, coord);

if (mom_check) {
// we do not include the pressure term in any non-Cartesian
// coordinate directions
F[UMX+idir] = F[UMX+idir] + p;
}

F[UEINT] = U[UEINT]*u_flx;
F[UEDEN] = (U[UEDEN] + p)*u_flx;

Expand Down Expand Up @@ -257,15 +249,9 @@ namespace HLL {

flux_hll[URHO] = (bp*fl_tmp - bm*fr_tmp)*bd + bp*bm*bd*(qr[QRHO] - ql[QRHO]);

// normal momentum flux. Note for 1-d and 2-d non cartesian
// r-coordinate, we leave off the pressure term and handle that
// separately in the update, to accommodate different geometries
fl_tmp = ql[QRHO]*ql[ivel]*ql[ivel];
fr_tmp = qr[QRHO]*qr[ivel]*qr[ivel];
if (mom_flux_has_p(idir, idir, coord)) {
fl_tmp = fl_tmp + ql[QPRES];
fr_tmp = fr_tmp + qr[QPRES];
}
// normal momentum flux.
fl_tmp = ql[QRHO]*ql[ivel]*ql[ivel] + ql[QPRES];
fr_tmp = qr[QRHO]*qr[ivel]*qr[ivel] + qr[QPRES];

flux_hll[imom] = (bp*fl_tmp - bm*fr_tmp)*bd + bp*bm*bd*(qr[QRHO]*qr[ivel] - ql[QRHO]*ql[ivel]);

Expand Down
11 changes: 2 additions & 9 deletions Source/hydro/riemann_solvers.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,9 +28,6 @@ compute_flux_q(const int i, const int j, const int k, const int idir,

int im1, im2, im3;

auto coord = geomdata.Coord();
auto mom_check = mom_flux_has_p(idir, idir, coord);

if (idir == 0) {
im1 = UMX;
im2 = UMY;
Expand All @@ -55,9 +52,7 @@ compute_flux_q(const int i, const int j, const int k, const int idir,
F(i,j,k,im3) = F(i,j,k,URHO) * qint.utt;

#ifndef RADIATION
if (mom_check) {
F(i,j,k,im1) += qint.p;
}
F(i,j,k,im1) += qint.p;

Real rhoetot = qint.rhoe + 0.5_rt * qint.rho *
(qint.un * qint.un +
Expand All @@ -67,9 +62,7 @@ compute_flux_q(const int i, const int j, const int k, const int idir,
F(i,j,k,UEDEN) = qint.un * (rhoetot + qint.p);
F(i,j,k,UEINT) = qint.un * qint.rhoe;
#else
if (mom_check) {
F(i,j,k,im1) += qint.p_g;
}
F(i,j,k,im1) += qint.p_g;

Real rhoetot = qint.rhoe_g + 0.5_rt * qint.rho *
(qint.un * qint.un +
Expand Down
7 changes: 0 additions & 7 deletions Source/hydro/trans.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -307,15 +307,8 @@ Castro::actual_trans_single(const Box& bx, // NOLINT(readability-convert-member

Real runewn = run - hdt * (area_t(ir,jr,kr) * flux_t(ir,jr,kr,UMX) -
area_t(il,jl,kl) * flux_t(il,jl,kl,UMX)) * volinv;
if (idir_t == 0 && !mom_flux_has_p(0, idir_t, coord)) {
runewn = runewn - cdtdx * (pgp - pgm);
}
Real rvnewn = rvn - hdt * (area_t(ir,jr,kr) * flux_t(ir,jr,kr,UMY) -
area_t(il,jl,kl) * flux_t(il,jl,kl,UMY)) * volinv;
if (idir_t == 1 && !mom_flux_has_p(1, idir_t, coord)) {
Real r = problo[0] + static_cast<Real>(il + 0.5_rt) * dx[0];
rvnewn = rvnewn - cdtdx / r * (pgp - pgm);
}
Real rwnewn = rwn - hdt * (area_t(ir,jr,kr) * flux_t(ir,jr,kr,UMZ) -
area_t(il,jl,kl) * flux_t(il,jl,kl,UMZ)) * volinv;
Real renewn = ren - hdt * (area_t(ir,jr,kr) * flux_t(ir,jr,kr,UEDEN) -
Expand Down
57 changes: 47 additions & 10 deletions Source/sources/Castro_geom.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -167,14 +167,32 @@ Castro::fill_RZ_geom_source (Real time, Real dt, MultiFab& cons_state, MultiFab&
[=] 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 r = prob_lo[0] + (static_cast<Real>(i) + 0.5_rt)*dx[0];
Real rinv = 1.0_rt / (prob_lo[0] + (static_cast<Real>(i) + 0.5_rt)*dx[0]);

// 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) / (U_arr(i,j,k,URHO) * r);
src(i,j,k,UMX) = U_arr(i,j,k,UMZ) * U_arr(i,j,k,UMZ) * rhoinv * rinv;

// 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) / (U_arr(i,j,k,URHO) * r);
src(i,j,k,UMZ) = - U_arr(i,j,k,UMX) * U_arr(i,j,k,UMZ) * rhoinv * rinv;

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

});
}
Expand Down Expand Up @@ -209,23 +227,42 @@ Castro::fill_RTheta_geom_source (Real time, Real dt, MultiFab& cons_state, Multi
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 r = prob_lo[0] + (static_cast<Real>(i) + 0.5_rt)*dx[0];
Real theta = prob_lo[1] + (static_cast<Real>(j) + 0.5_rt)*dx[1];
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]);

// 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)) / (U_arr(i,j,k,URHO) * r);
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) * cot(theta) -
U_arr(i,j,k,UMX) * U_arr(i,j,k,UMY)) / (U_arr(i,j,k,URHO) * 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;

// 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) * cot(theta) -
U_arr(i,j,k,UMX) * U_arr(i,j,k,UMZ)) / (U_arr(i,j,k,URHO) * 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;

// 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 = 2 p / r
src(i,j,k,UMX) = 2.0_rt * eos_state.p * rinv;

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

0 comments on commit 84ffcff

Please sign in to comment.