Skip to content

Commit

Permalink
Finish riemann cleaning (#2897)
Browse files Browse the repository at this point in the history 
we now rely on Source/hydro/riemann_constants.H for all of the necessary
tolerances.
  • Loading branch information
@zingale
zingale authored Jul 10, 2024
1 parent dd89b8e commit 5590794
Show file tree
Hide file tree
Showing 10 changed files with 89 additions and 99 deletions.
3 changes: 3 additions & 0 deletions Source/hydro/riemann_constants.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,9 @@ namespace riemann_constants {
constexpr amrex::Real smlp1 = 1.e-10_rt;
constexpr amrex::Real small = 1.e-8_rt;
constexpr amrex::Real smallu = 1.e-12_rt;
constexpr amrex::Real riemann_integral_tol = 1.e-8_rt;
constexpr amrex::Real riemann_u_tol = 1.e-6_rt;
constexpr amrex::Real riemann_p_tol = 1.e-8_rt;
constexpr int HISTORY_SIZE=40;
constexpr int PSTAR_BISECT_FACTOR = 5;
}
Expand Down
7 changes: 4 additions & 3 deletions Util/exact_riemann/Make.package
View file
Original file line number Diff line number Diff line change
@@ -1,9 +1,10 @@
CEXE_sources += main.cpp

CEXE_headers += exact_riemann.H
CEXE_headers += riemann_star_state.H
CEXE_headers += riemann_sample.H
CEXE_headers += riemann_support.H
CEXE_headers += exact_riemann_star_state.H
CEXE_headers += exact_riemann_sample.H
CEXE_headers += exact_riemann_shock.H
CEXE_headers += exact_riemann_rarefaction.H
CEXE_sources += extern_parameters.cpp
CEXE_headers += extern_parameters.H

Expand Down
90 changes: 45 additions & 45 deletions Util/exact_riemann/ci-benchmarks/test1-riemann.out
View file

Large diffs are not rendered by default.

13 changes: 8 additions & 5 deletions Util/exact_riemann/exact_riemann.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,8 +7,8 @@

#include <castro_params.H>

#include <riemann_sample.H>
#include <riemann_star_state.H>
#include <exact_riemann_sample.H>
#include <exact_riemann_star_state.H>

using namespace amrex::literals;

Expand Down Expand Up @@ -50,10 +50,13 @@ exact_riemann() {

amrex::Real ustar, pstar, W_l, W_r;

riemann_star_state(problem::rho_l, problem::u_l, problem::p_l, xn,
problem::rho_r, problem::u_r, problem::p_r, xn,
ustar, pstar, W_l, W_r);
std::cout << "solving for star state: " << std::endl;

auto niter = riemann_star_state(problem::rho_l, problem::u_l, problem::p_l, xn,
problem::rho_r, problem::u_r, problem::p_r, xn,
ustar, pstar, W_l, W_r);

std::cout << "found star state, niter = " << niter << "; pstar, ustar = " << pstar << " " << ustar << std::endl;

// find the solution as a function of xi = x/t
std::ofstream ofile;
Expand Down
14 changes: 5 additions & 9 deletions Util/exact_riemann/riemann_rarefaction.H → ...exact_riemann/exact_riemann_rarefaction.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -174,8 +174,6 @@ rarefaction(const amrex::Real pstar,
const double S1{0.9};
const double S2{4.0};

constexpr amrex::Real tol = 1.e-5;

// initial conditions

amrex::Real tau = 1.0_rt / rho_s;
Expand Down Expand Up @@ -206,7 +204,7 @@ rarefaction(const amrex::Real pstar,

// take a step

while (rel_error > tol) {
while (rel_error > riemann_constants::riemann_integral_tol) {
dp = dp_new;
if (p + dp < pstar) {
dp = pstar - p;
Expand All @@ -232,7 +230,7 @@ rarefaction(const amrex::Real pstar,

// adaptive step estimate

double dp_est = S1 * dp * std::pow(std::abs(tol/rel_error), 0.2);
double dp_est = S1 * dp * std::pow(std::abs(riemann_constants::riemann_integral_tol/rel_error), 0.2);
dp_new = dp_est; //std::clamp(S1*dp_est, dp/S2, S2*dp);

}
Expand Down Expand Up @@ -295,8 +293,6 @@ rarefaction_to_u(const amrex::Real rho_s, const amrex::Real u_s, const amrex::Re

const double S1{0.9};

constexpr amrex::Real tol = 1.e-8;

// initial conditions

amrex::Real tau = 1.0_rt / rho_s;
Expand Down Expand Up @@ -343,7 +339,7 @@ rarefaction_to_u(const amrex::Real rho_s, const amrex::Real u_s, const amrex::Re

// take a step

while (rel_error > tol) {
while (rel_error > riemann_constants::riemann_integral_tol) {
du = du_new;

// take 2 half steps
Expand All @@ -368,7 +364,7 @@ rarefaction_to_u(const amrex::Real rho_s, const amrex::Real u_s, const amrex::Re

amrex::Real du_est{};
if (rel_error > 0) {
du_est = S1 * du * std::pow(std::abs(tol / rel_error), 0.2);
du_est = S1 * du * std::pow(std::abs(riemann_constants::riemann_integral_tol / rel_error), 0.2);
} else {
du_est = 10.0 * du;
}
Expand Down Expand Up @@ -417,7 +413,7 @@ rarefaction_to_u(const amrex::Real rho_s, const amrex::Real u_s, const amrex::Re
}
}

if (std::abs(du_new) < riemann_solver::tol * std::abs(u) + riemann_solver::tol) {
if (std::abs(du_new) < riemann_constants::riemann_u_tol * std::abs(u) + riemann_constants::riemann_u_tol) {
finished = true;
}

Expand Down
4 changes: 2 additions & 2 deletions Util/exact_riemann/riemann_sample.H → Util/exact_riemann/exact_riemann_sample.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,8 +3,8 @@

#include <AMReX_REAL.H>

#include <riemann_support.H>
#include <riemann_rarefaction.H>
#include <riemann_constants.H>
#include <exact_riemann_rarefaction.H>

using namespace amrex::literals;

Expand Down
23 changes: 12 additions & 11 deletions Util/exact_riemann/riemann_shock.H → Util/exact_riemann/exact_riemann_shock.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,18 +40,18 @@ W_s_shock(const amrex::Real W_s, const amrex::Real pstar,
amrex::Real dedrho_p = eos_state.dedr - eos_state.dedT * eos_state.dpdr / eos_state.dpdT;

fprime = 2.0_rt * W_s * (eos_state.e - e_s) -
2.0_rt * dedrho_p * (pstar - p_s) * std::pow(rhostar_s, 2) / W_s;
2.0_rt * dedrho_p * (pstar - p_s) * amrex::Math::powi<2>(rhostar_s) / W_s;

}


AMREX_INLINE
void
bool
newton_shock(amrex::Real& W_s, const amrex::Real pstar,
const amrex::Real rho_s, const amrex::Real p_s, const amrex::Real e_s,
const amrex::Real* xn,
amrex::Real* rhostar_hist, amrex::Real* Ws_hist,
eos_rep_t& eos_state, bool& converged) {
eos_rep_t& eos_state) {

// Newton iterations -- we are zeroing the energy R-H jump condition
// W^2 [e] = 1/2 [p^2]
Expand All @@ -60,7 +60,7 @@ newton_shock(amrex::Real& W_s, const amrex::Real pstar,
//
// and then compute f'

converged = false;
bool converged = false;

int iter = 1;
while (! converged && iter < castro::riemann_shock_maxiter) {
Expand All @@ -76,7 +76,7 @@ newton_shock(amrex::Real& W_s, const amrex::Real pstar,
converged = true;
}

W_s = amrex::min(2.0_rt * W_s, amrex::max(0.5_rt * W_s, W_s + dW));
W_s = std::clamp(W_s + dW, 0.5_rt * W_s, 2.0_rt * W_s);

// store some history

Expand All @@ -85,6 +85,8 @@ newton_shock(amrex::Real& W_s, const amrex::Real pstar,

iter++;
}

return converged;
}

AMREX_INLINE
Expand Down Expand Up @@ -134,7 +136,7 @@ shock(const amrex::Real pstar,
// just doesn't work. In this case, we use the ideas from CG, Eq. 35, and
// take W = sqrt(gamma p rho)

if (pstar - p_s < riemann_solver::tol_p * p_s) {
if (pstar - p_s < riemann_constants::riemann_p_tol * p_s) {
W_s = std::sqrt(eos_state.gam1 * p_s * rho_s);
} else {
W_s = std::sqrt((pstar - p_s) *
Expand All @@ -154,10 +156,9 @@ shock(const amrex::Real pstar,

// newton

bool converged;
newton_shock(W_s, pstar, rho_s, p_s, e_s, xn,
rhostar_hist, Ws_hist,
eos_state, converged);
auto converged = newton_shock(W_s, pstar, rho_s, p_s, e_s, xn,
rhostar_hist, Ws_hist,
eos_state);


// now did we converge?
Expand Down Expand Up @@ -185,7 +186,7 @@ shock(const amrex::Real pstar,
amrex::Real p_e = eos_state.dpdT / eos_state.dedT;
amrex::Real p_rho = eos_state.dpdr - eos_state.dpdT * eos_state.dedr / eos_state.dedT;

amrex::Real p_tau = -std::pow(rhostar_s, 2) * p_rho;
amrex::Real p_tau = -amrex::Math::powi<2>(rhostar_s) * p_rho;

amrex::Real dW2dpstar = (C*C - W_s*W_s) * W_s * W_s /
((0.5_rt * (pstar + p_s) * p_e - p_tau) * (pstar - p_s));
Expand Down
19 changes: 9 additions & 10 deletions Util/exact_riemann/riemann_star_state.H → .../exact_riemann/exact_riemann_star_state.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,14 +3,14 @@

#include <AMReX_REAL.H>

#include <riemann_support.H>
#include <riemann_shock.H>
#include <riemann_rarefaction.H>
#include <riemann_constants.H>
#include <exact_riemann_shock.H>
#include <exact_riemann_rarefaction.H>

using namespace amrex::literals;

AMREX_INLINE
void
int
riemann_star_state(const amrex::Real rho_l, const amrex::Real u_l, const amrex::Real p_l, const amrex::Real* xn_l,
const amrex::Real rho_r, const amrex::Real u_r, const amrex::Real p_r, const amrex::Real* xn_r,
amrex::Real& ustar, amrex::Real& pstar, amrex::Real& W_l, amrex::Real& W_r) {
Expand Down Expand Up @@ -74,13 +74,11 @@ riemann_star_state(const amrex::Real rho_l, const amrex::Real u_l, const amrex::
pstar = ((W_r * p_l + W_l * p_r) + W_l * W_r * (u_l - u_r)) / (W_l + W_r);
}

pstar = amrex::max(pstar, smallp);
pstar = std::max(pstar, smallp);


// find the exact pstar and ustar

std::cout << "solving for star state: " << std::endl;

// this procedure follows directly from Colella & Glaz 1985, section 1
// the basic idea is that we want to find the pstar that satisfies:
//
Expand Down Expand Up @@ -108,7 +106,7 @@ riemann_star_state(const amrex::Real rho_l, const amrex::Real u_l, const amrex::
bool converged = false;

int iter = 1;
while (! converged && iter < riemann_solver::max_iters) {
while (! converged && iter < castro::riemann_shock_maxiter) {

// compute Z_l, W_l and Z_r, W_r -- the form of these depend
// on whether the wave is a shock or a rarefaction
Expand Down Expand Up @@ -151,7 +149,7 @@ riemann_star_state(const amrex::Real rho_l, const amrex::Real u_l, const amrex::
}

// get ready for the next iteration
pstar = pstar_new;
pstar = std::clamp(pstar_new, 0.5 * pstar, 2.0 * pstar);

iter++;
}
Expand All @@ -162,6 +160,7 @@ riemann_star_state(const amrex::Real rho_l, const amrex::Real u_l, const amrex::

ustar = 0.5_rt * (ustar_l + ustar_r);

std::cout << "found pstar, ustar: " << pstar << " " << ustar << std::endl;
return iter;

}
#endif
1 change: 1 addition & 0 deletions Util/exact_riemann/inputs.bug8
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,3 +16,4 @@ problem.t = 0.0008

problem.npts = 512

castro.riemann_shock_maxiter = 100
14 changes: 0 additions & 14 deletions Util/exact_riemann/riemann_support.H
View file

This file was deleted.

0 comments on commit 5590794

Please sign in to comment.