diff --git a/Util/exact_riemann/exact_riemann.H b/Util/exact_riemann/exact_riemann.H index 021120a322..d7e7c4831b 100644 --- a/Util/exact_riemann/exact_riemann.H +++ b/Util/exact_riemann/exact_riemann.H @@ -1,12 +1,16 @@ #ifndef EXACT_RIEMANN_H #define EXACT_RIEMANN_H +#include + #include #include #include #include +using namespace amrex::literals; + AMREX_INLINE void exact_riemann() { @@ -15,7 +19,7 @@ exact_riemann() { // exploring composition jumps here. We'll take a constant // composition. - Real xn[NumSpec] = {0.0}; + amrex::Real xn[NumSpec] = {0.0}; xn[0] = 1.0_rt; @@ -44,12 +48,12 @@ exact_riemann() { } if (problem::co_moving_frame) { - Real W_avg = 0.5_rt * (problem::u_l + problem::u_r); + amrex::Real W_avg = 0.5_rt * (problem::u_l + problem::u_r); problem::u_l -= W_avg; problem::u_r -= W_avg; } - Real ustar, pstar, W_l, W_r; + 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, @@ -62,7 +66,7 @@ exact_riemann() { // This follows from the discussion around C&G Eq. 15 - Real dx = (problem::xmax - problem::xmin) / static_cast(problem::npts); + amrex::Real dx = (problem::xmax - problem::xmin) / static_cast(problem::npts); // loop over xi space @@ -72,9 +76,9 @@ exact_riemann() { // compute xi = x/t -- this is the similarity variable for the // solution - Real x = problem::xmin + (static_cast(i) - 0.5_rt) * dx; + amrex::Real x = problem::xmin + (static_cast(i) - 0.5_rt) * dx; - Real rho, u, p, xn_s[NumSpec]; + amrex::Real rho, u, p, xn_s[NumSpec]; riemann_sample(problem::rho_l, problem::u_l, problem::p_l, xn, problem::rho_r, problem::u_r, problem::p_r, xn, @@ -83,7 +87,7 @@ exact_riemann() { rho, u, p, xn_s); if (problem::co_moving_frame) { - Real W_avg = 0.5_rt * (problem::u_l + problem::u_r); + amrex::Real W_avg = 0.5_rt * (problem::u_l + problem::u_r); u += W_avg; x += problem::t * W_avg; } diff --git a/Util/exact_riemann/main.cpp b/Util/exact_riemann/main.cpp index 20d8447738..e91b337308 100644 --- a/Util/exact_riemann/main.cpp +++ b/Util/exact_riemann/main.cpp @@ -17,25 +17,25 @@ int main(int argc, char *argv[]) { amrex::Initialize(argc, argv); - std::cout << "starting the exact Riemann solver..." << std::endl; - std::cout << argv[1] << std::endl; - std::cout << strlen(argv[1]) << std::endl; + std::cout << "starting the exact Riemann solver..." << std::endl; + std::cout << argv[1] << std::endl; + std::cout << strlen(argv[1]) << std::endl; - // initialize the external runtime parameters in C++ + // initialize the external runtime parameters in C++ - init_prob_parameters(); + init_prob_parameters(); - init_extern_parameters(); + init_extern_parameters(); - // now initialize the C++ Microphysics + // now initialize the C++ Microphysics #ifdef REACTIONS - network_init(); + network_init(); #endif - Real small_temp = 1.e-200; - Real small_dens = 1.e-200; - eos_init(small_temp, small_dens); + amrex::Real small_temp = 1.e-200; + amrex::Real small_dens = 1.e-200; + eos_init(small_temp, small_dens); - exact_riemann(); + exact_riemann(); } diff --git a/Util/exact_riemann/riemann_sample.H b/Util/exact_riemann/riemann_sample.H index dcadceb8a9..c2580dd258 100644 --- a/Util/exact_riemann/riemann_sample.H +++ b/Util/exact_riemann/riemann_sample.H @@ -1,16 +1,20 @@ #ifndef RIEMANN_SAMPLE_MODULE_H #define RIEMANN_SAMPLE_MODULE_H +#include + #include +using namespace amrex::literals; + AMREX_INLINE void -riemann_sample(const Real rho_l, const Real u_l, const Real p_l, const Real* xn_l, - const Real rho_r, const Real u_r, const Real p_r, const Real* xn_r, - const Real ustar, const Real pstar, - const Real W_l, const Real W_r, - const Real x, const Real xjump, const Real time, - Real& rho, Real& u, Real& p, Real* xn) { +riemann_sample(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, + const amrex::Real ustar, const amrex::Real pstar, + const amrex::Real W_l, const amrex::Real W_r, + const amrex::Real x, const amrex::Real xjump, const amrex::Real time, + amrex::Real& rho, amrex::Real& u, amrex::Real& p, amrex::Real* xn) { // get the initial sound speeds @@ -24,7 +28,7 @@ riemann_sample(const Real rho_l, const Real u_l, const Real p_l, const Real* xn_ eos(eos_input_rp, eos_state); - Real cs_l = std::sqrt(eos_state.gam1 * p_l / rho_l); + amrex::Real cs_l = std::sqrt(eos_state.gam1 * p_l / rho_l); eos_state.rho = rho_r; eos_state.p = p_r; @@ -35,7 +39,7 @@ riemann_sample(const Real rho_l, const Real u_l, const Real p_l, const Real* xn_ eos(eos_input_rp, eos_state); - Real cs_r = std::sqrt(eos_state.gam1 * p_r / rho_r); + amrex::Real cs_r = std::sqrt(eos_state.gam1 * p_r / rho_r); // find the solution as a function of xi = x/t @@ -45,14 +49,14 @@ riemann_sample(const Real rho_l, const Real u_l, const Real p_l, const Real* xn_ // compute xi = x/t -- this is the similarity variable for the // solution - Real xi = (x - xjump) / time; + amrex::Real xi = (x - xjump) / time; // check which side of the contact we need to worry about - Real chi = std::copysign(1.0_rt, xi - ustar); + amrex::Real chi = std::copysign(1.0_rt, xi - ustar); - Real rho_s, u_s, p_s, W_s, cs_s; - Real uhat_s, xihat, uhat_star; + amrex::Real rho_s, u_s, p_s, W_s, cs_s; + amrex::Real uhat_s, xihat, uhat_star; if (chi == -1.0_rt) { rho_s = rho_l; @@ -97,7 +101,7 @@ riemann_sample(const Real rho_l, const Real u_l, const Real p_l, const Real* xn_ // are we a shock or rarefaction? - Real rhostar, Z_temp; + amrex::Real rhostar, Z_temp; if (pstar > p_s) { // shock @@ -107,7 +111,7 @@ riemann_sample(const Real rho_l, const Real u_l, const Real p_l, const Real* xn_ } else { // rarefaction. Here we need to integrate the Riemann // invariant curves to our pstar to get rhostar and cs_star - Real Z_temp, W_temp; + amrex::Real Z_temp, W_temp; if (chi == -1.0_rt) { rarefaction(pstar, rho_s, u_s, p_s, xn, 1, Z_temp, W_temp, rhostar); } else { @@ -128,11 +132,11 @@ riemann_sample(const Real rho_l, const Real u_l, const Real p_l, const Real* xn_ eos(eos_input_rp, eos_state); - Real cs_star = std::sqrt(eos_state.gam1 * pstar / rhostar); + amrex::Real cs_star = std::sqrt(eos_state.gam1 * pstar / rhostar); // now deal with the cases where we are not spanning a rarefaction - Real lambdahat_s, lambdahat_star; + amrex::Real lambdahat_s, lambdahat_star; if (pstar <= p_s) { lambdahat_s = uhat_s + cs_s; diff --git a/Util/exact_riemann/riemann_star_state.H b/Util/exact_riemann/riemann_star_state.H index 9cc23a8158..b29c81dd50 100644 --- a/Util/exact_riemann/riemann_star_state.H +++ b/Util/exact_riemann/riemann_star_state.H @@ -1,17 +1,21 @@ #ifndef RIEMANN_STAR_STATE_H #define RIEMANN_STAR_STATE_H +#include + #include +using namespace amrex::literals; + AMREX_INLINE void -riemann_star_state(const Real rho_l, const Real u_l, const Real p_l, const Real* xn_l, - const Real rho_r, const Real u_r, const Real p_r, const Real* xn_r, - Real& ustar, Real& pstar, Real& W_l, Real& W_r) { +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) { - const Real tol = 1.e-10_rt; - const Real smallp = 1.e-8_rt; - const Real SMALL = 1.e-13_rt; + const amrex::Real tol = 1.e-10_rt; + const amrex::Real smallp = 1.e-8_rt; + const amrex::Real SMALL = 1.e-13_rt; // get the initial sound speeds @@ -26,10 +30,10 @@ riemann_star_state(const Real rho_l, const Real u_l, const Real p_l, const Real* eos(eos_input_rp, eos_state); - Real cs_l = std::sqrt(eos_state.gam1 * p_l / rho_l); + amrex::Real cs_l = std::sqrt(eos_state.gam1 * p_l / rho_l); - Real gammaE_l = p_l / (rho_l * eos_state.e) + 1.0_rt; - Real gammaC_l = eos_state.gam1; + amrex::Real gammaE_l = p_l / (rho_l * eos_state.e) + 1.0_rt; + amrex::Real gammaC_l = eos_state.gam1; eos_state.rho = rho_r; eos_state.p = p_r; @@ -40,13 +44,13 @@ riemann_star_state(const Real rho_l, const Real u_l, const Real p_l, const Real* eos(eos_input_rp, eos_state); - Real cs_r = std::sqrt(eos_state.gam1 * p_r / rho_r); + amrex::Real cs_r = std::sqrt(eos_state.gam1 * p_r / rho_r); - Real gammaE_r = p_r / (rho_r * eos_state.e) + 1.0_rt; - Real gammaC_r = eos_state.gam1; + amrex::Real gammaE_r = p_r / (rho_r * eos_state.e) + 1.0_rt; + amrex::Real gammaC_r = eos_state.gam1; - Real gammaE_bar = 0.5_rt * (gammaE_l + gammaE_r); - Real gammaC_bar = 0.5_rt * (gammaC_l + gammaC_r); + amrex::Real gammaE_bar = 0.5_rt * (gammaE_l + gammaE_r); + amrex::Real gammaC_bar = 0.5_rt * (gammaC_l + gammaC_r); // create an initial guess for pstar using a primitive variable @@ -78,7 +82,7 @@ riemann_star_state(const Real rho_l, const Real u_l, const Real p_l, const Real* // this procedure follows directly from Colella & Glaz 1985, section 1 - Real ustar_l, ustar_r; + amrex::Real ustar_l, ustar_r; bool converged = false; int iter = 1; @@ -87,7 +91,7 @@ riemann_star_state(const Real rho_l, const Real u_l, const Real p_l, const Real* // compute Z_l and Z_r -- the form of these depend on whether the // wave is a shock or a rarefaction - Real Z_l, Z_r; + amrex::Real Z_l, Z_r; // left wave @@ -96,7 +100,7 @@ riemann_star_state(const Real rho_l, const Real u_l, const Real p_l, const Real* shock(pstar, rho_l, u_l, p_l, xn_l, gammaE_bar, gammaC_bar, Z_l, W_l); } else { // left rarefaction - Real rhostar_dummy; + amrex::Real rhostar_dummy; rarefaction(pstar, rho_l, u_l, p_l, xn_l, 1, Z_l, W_l, rhostar_dummy); } @@ -107,21 +111,21 @@ riemann_star_state(const Real rho_l, const Real u_l, const Real p_l, const Real* shock(pstar, rho_r, u_r, p_r, xn_r, gammaE_bar, gammaC_bar, Z_r, W_r); } else { // right rarefaction - Real rhostar_dummy; + amrex::Real rhostar_dummy; rarefaction(pstar, rho_r, u_r, p_r, xn_r, 3, Z_r, W_r, rhostar_dummy); } ustar_l = u_l - (pstar - p_l) / W_l; ustar_r = u_r + (pstar - p_r) / W_r; - Real pstar_new = pstar - Z_l * Z_r * (ustar_r - ustar_l) / (Z_l + Z_r); + amrex::Real pstar_new = pstar - Z_l * Z_r * (ustar_r - ustar_l) / (Z_l + Z_r); std::cout << "done with iteration " << iter << std::endl; std::cout << "ustar_l/r, pstar: " << ustar_l << " " << ustar_r << " " << pstar_new << std::endl; // estimate the error in the current star solution - Real err1 = std::abs(ustar_r - ustar_l); - Real err2 = pstar_new - pstar; + amrex::Real err1 = std::abs(ustar_r - ustar_l); + amrex::Real err2 = pstar_new - pstar; if (err1 < tol * amrex::max(std::abs(ustar_l), std::abs(ustar_r)) && err2 < tol * pstar) { diff --git a/Util/exact_riemann/riemann_support.H b/Util/exact_riemann/riemann_support.H index 002e183815..80f3391492 100644 --- a/Util/exact_riemann/riemann_support.H +++ b/Util/exact_riemann/riemann_support.H @@ -1,22 +1,26 @@ #ifndef RIEMANN_SUPPORT_H #define RIEMANN_SUPPORT_H +#include + #include #include -const Real smallrho = 1.e-5_rt; +using namespace amrex::literals; + +const amrex::Real smallrho = 1.e-5_rt; const int max_iters = 100; -const Real tol = 1.e-6_rt; +const amrex::Real tol = 1.e-6_rt; AMREX_INLINE void -W_s_shock(const Real W_s, const Real pstar, const Real rho_s, const Real p_s, const Real e_s, const Real* xn, - Real& rhostar_s, eos_t& eos_state, Real& f, Real& fprime) { +W_s_shock(const 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_s, eos_t& eos_state, amrex::Real& f, amrex::Real& fprime) { // we need rhostar -- get it from the R-H conditions - Real taustar_s = (1.0_rt / rho_s) - (pstar - p_s) / (W_s * W_s); + amrex::Real taustar_s = (1.0_rt / rho_s) - (pstar - p_s) / (W_s * W_s); rhostar_s = 1.0_rt / taustar_s; // get the thermodynamics @@ -36,7 +40,7 @@ W_s_shock(const Real W_s, const Real pstar, const Real rho_s, const Real p_s, co // we need de/drho at constant p -- this is not returned by the EOS - Real dedrho_p = eos_state.dedr - eos_state.dedT * eos_state.dpdr / eos_state.dpdT; + 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; @@ -45,8 +49,8 @@ W_s_shock(const Real W_s, const Real pstar, const Real rho_s, const Real p_s, co AMREX_INLINE void -newton_shock(Real& W_s, const Real pstar, const Real rho_s, const Real p_s, const Real e_s, const Real* xn, - Real* rhostar_hist, Real* Ws_hist, +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_t& eos_state, bool& converged) { // Newton iterations -- we are zeroing the energy R-H jump condition @@ -61,12 +65,12 @@ newton_shock(Real& W_s, const Real pstar, const Real rho_s, const Real p_s, cons int iter = 1; while (! converged && iter < max_iters) { - Real rhostar_s, f, fprime; + amrex::Real rhostar_s, f, fprime; W_s_shock(W_s, pstar, rho_s, p_s, e_s, xn, rhostar_s, eos_state, f, fprime); - Real dW = -f / fprime; + amrex::Real dW = -f / fprime; if (std::abs(dW) < tol * W_s) { converged = true; @@ -85,14 +89,14 @@ newton_shock(Real& W_s, const Real pstar, const Real rho_s, const Real p_s, cons AMREX_INLINE void -shock(const Real pstar, const Real rho_s, const Real u_s, const Real p_s, const Real* xn, - const Real gammaE_bar, const Real gammaC_bar, - Real& Z_s, Real& W_s) { +shock(const amrex::Real pstar, const amrex::Real rho_s, const amrex::Real u_s, const amrex::Real p_s, const amrex::Real* xn, + const amrex::Real gammaE_bar, const amrex::Real gammaC_bar, + amrex::Real& Z_s, amrex::Real& W_s) { - const Real tol_p = 1.e-6_rt; + const amrex::Real tol_p = 1.e-6_rt; - Real rhostar_hist[max_iters], Ws_hist[max_iters]; + amrex::Real rhostar_hist[max_iters], Ws_hist[max_iters]; // compute the Z_s function for a shock following C&G Eq. 20 and // 23. Here the "_s" variables are the state (L or R) that we are @@ -113,15 +117,15 @@ shock(const Real pstar, const Real rho_s, const Real u_s, const Real p_s, const eos(eos_input_rp, eos_state); - Real e_s = eos_state.e; + amrex::Real e_s = eos_state.e; // to kick things off, we need a guess for W_s. We'll use the // approximation from Colella & Glaz (Eq. 34), which in turn // makes an approximation for gammaE_star, using equation 31. - Real gammaE_s = p_s / (rho_s * e_s) + 1.0_rt; + amrex::Real gammaE_s = p_s / (rho_s * e_s) + 1.0_rt; - Real gammaE_star = gammaE_s + + amrex::Real gammaE_star = gammaE_s + 2.0_rt * (1.0_rt - gammaE_bar / gammaC_bar) * (gammaE_bar - 1.0_rt) * (pstar - p_s) / (pstar + p_s); @@ -143,15 +147,15 @@ shock(const Real pstar, const Real rho_s, const Real u_s, const Real p_s, const // we need rhostar -- get it from the R-H conditions - Real taustar_s = (1.0_rt / rho_s) - (pstar - p_s) / (W_s * W_s); - Real rhostar_s; + amrex::Real taustar_s = (1.0_rt / rho_s) - (pstar - p_s) / (W_s * W_s); + amrex::Real rhostar_s; if (taustar_s < 0.0_rt) { rhostar_s = smallrho; W_s = std::sqrt((pstar - p_s) / (1.0_rt / rho_s - 1.0_rt / rhostar_s)); } - Real W_s_guess = W_s; + amrex::Real W_s_guess = W_s; // newton @@ -181,17 +185,17 @@ shock(const Real pstar, const Real rho_s, const Real u_s, const Real p_s, const // next we compute the derivative dW_s/dpstar -- the paper gives // dW**2/dpstar (Eq. 23), so we take 1/2W of that - Real C = std::sqrt(eos_state.gam1 * pstar * rhostar_s); + amrex::Real C = std::sqrt(eos_state.gam1 * pstar * rhostar_s); - Real p_e = eos_state.dpdT / eos_state.dedT; - Real p_rho = eos_state.dpdr - eos_state.dpdT * eos_state.dedr / eos_state.dedT; + 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; - Real p_tau = -std::pow(rhostar_s, 2) * p_rho; + amrex::Real p_tau = -std::pow(rhostar_s, 2) * p_rho; - Real dW2dpstar = (C*C - W_s*W_s) * W_s * W_s / + 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)); - Real dWdpstar = 0.5_rt * dW2dpstar / W_s; + amrex::Real dWdpstar = 0.5_rt * dW2dpstar / W_s; // finally compute Z_s @@ -202,8 +206,8 @@ shock(const Real pstar, const Real rho_s, const Real u_s, const Real p_s, const AMREX_INLINE void -riemann_invariant_rhs(const Real p, const Real tau, const Real u, const Real* xn, const int iwave, - Real& dtaudp, Real& dudp) { +riemann_invariant_rhs(const amrex::Real p, const amrex::Real tau, const amrex::Real u, const amrex::Real* xn, const int iwave, + amrex::Real& dtaudp, amrex::Real& dudp) { // here, p is out independent variable, and tau, u are the // dependent variables. We return the derivatives of these @@ -221,7 +225,7 @@ riemann_invariant_rhs(const Real p, const Real tau, const Real u, const Real* xn eos(eos_input_rp, eos_state); - Real C = std::sqrt(eos_state.gam1 * p / tau); + amrex::Real C = std::sqrt(eos_state.gam1 * p / tau); dtaudp = -1.0_rt / (C * C); @@ -236,8 +240,8 @@ riemann_invariant_rhs(const Real p, const Real tau, const Real u, const Real* xn AMREX_INLINE void -riemann_invariant_rhs2(const Real u, const Real tau, const Real p, const Real* xn, const int iwave, - Real& dtaudu, Real& dpdu) { +riemann_invariant_rhs2(const amrex::Real u, const amrex::Real tau, const amrex::Real p, const amrex::Real* xn, const int iwave, + amrex::Real& dtaudu, amrex::Real& dpdu) { // here, u is out independent variable, and tau, p are the // dependent variables. We return the derivatives of these @@ -255,7 +259,7 @@ riemann_invariant_rhs2(const Real u, const Real tau, const Real p, const Real* x eos(eos_input_rp, eos_state); - Real C = std::sqrt(eos_state.gam1 * p / tau); + amrex::Real C = std::sqrt(eos_state.gam1 * p / tau); if (iwave == 3) { dpdu = C; @@ -270,8 +274,8 @@ riemann_invariant_rhs2(const Real u, const Real tau, const Real p, const Real* x AMREX_INLINE void -rarefaction(const Real pstar, const Real rho_s, const Real u_s, const Real p_s, - const Real* xn, const int iwave, Real& Z_s, Real& W_s, Real& rhostar) { +rarefaction(const amrex::Real pstar, const amrex::Real rho_s, const amrex::Real u_s, const amrex::Real p_s, + const amrex::Real* xn, const int iwave, amrex::Real& Z_s, amrex::Real& W_s, amrex::Real& rhostar) { const int npts = 1000; @@ -279,27 +283,27 @@ rarefaction(const Real pstar, const Real rho_s, const Real u_s, const Real p_s, // region by integrating the Riemann invariant from p_s to pstar. // This means solving a system of ODEs. We use 4th-order R-K. - Real tau = 1.0_rt / rho_s; - Real u = u_s; - Real p = p_s; + amrex::Real tau = 1.0_rt / rho_s; + amrex::Real u = u_s; + amrex::Real p = p_s; - Real dp = (pstar - p_s) / static_cast(npts); - Real dp2 = 0.5_rt * dp; + amrex::Real dp = (pstar - p_s) / static_cast(npts); + amrex::Real dp2 = 0.5_rt * dp; for (int i = 1; i <= npts; ++i) { // do 4th-order RT - Real dtaudp1, dudp1; + amrex::Real dtaudp1, dudp1; riemann_invariant_rhs(p, tau, u, xn, iwave, dtaudp1, dudp1); - Real dtaudp2, dudp2; + amrex::Real dtaudp2, dudp2; riemann_invariant_rhs(p+dp2, tau+dp2*dtaudp1, u+dp2*dudp1, xn, iwave, dtaudp2, dudp2); - Real dtaudp3, dudp3; + amrex::Real dtaudp3, dudp3; riemann_invariant_rhs(p+dp2, tau+dp2*dtaudp2, u+dp2*dudp2, xn, iwave, dtaudp3, dudp3); - Real dtaudp4, dudp4; + amrex::Real dtaudp4, dudp4; riemann_invariant_rhs(p+dp, tau+dp*dtaudp3, u+dp*dudp3, xn, iwave, dtaudp4, dudp4); p += dp; @@ -336,8 +340,8 @@ rarefaction(const Real pstar, const Real rho_s, const Real u_s, const Real p_s, AMREX_INLINE void -rarefaction_to_u(const Real rho_s, const Real u_s, const Real p_s, const Real* xn, const int iwave, const Real xi, - Real& rho, Real& p, Real& u) { +rarefaction_to_u(const amrex::Real rho_s, const amrex::Real u_s, const amrex::Real p_s, const amrex::Real* xn, const int iwave, const amrex::Real xi, + amrex::Real& rho, amrex::Real& p, amrex::Real& u) { const int npts = 1000; @@ -355,7 +359,7 @@ rarefaction_to_u(const Real rho_s, const Real u_s, const Real p_s, const Real* x // dp/du = C; dtau/du = -1/C for the 1-wave // dp/du = -C; dtau/du = 1/C for the 3-wave - Real tau = 1.0_rt / rho_s; + amrex::Real tau = 1.0_rt / rho_s; u = u_s; p = p_s; @@ -372,37 +376,37 @@ rarefaction_to_u(const Real rho_s, const Real u_s, const Real p_s, const Real* x eos(eos_input_rp, eos_state); - Real c = std::sqrt(eos_state.gam1 * p * tau); + amrex::Real c = std::sqrt(eos_state.gam1 * p * tau); - Real ustop; + amrex::Real ustop; if (iwave == 1) { ustop = xi + c; } else if (iwave == 3) { ustop = xi - c; } - Real du = (ustop - u_s) / static_cast(npts); + amrex::Real du = (ustop - u_s) / static_cast(npts); bool finished = false; std::cout << "integrating from u: " << u << " " << ustop << " " << xi << " " << c << std::endl; - Real du2 = 0.5_rt * du; + amrex::Real du2 = 0.5_rt * du; while (! finished) { // do 4th-order RT - Real dtaudu1, dpdu1; + amrex::Real dtaudu1, dpdu1; riemann_invariant_rhs2(u, tau, p, xn, iwave, dtaudu1, dpdu1); - Real dtaudu2, dpdu2; + amrex::Real dtaudu2, dpdu2; riemann_invariant_rhs2(u+du2, tau+du2*dtaudu1, p+du2*dpdu1, xn, iwave, dtaudu2, dpdu2); - Real dtaudu3, dpdu3; + amrex::Real dtaudu3, dpdu3; riemann_invariant_rhs2(u+du2, tau+du2*dtaudu2, p+du2*dpdu2, xn, iwave, dtaudu3, dpdu3); - Real dtaudu4, dpdu4; + amrex::Real dtaudu4, dpdu4; riemann_invariant_rhs2(u+du, tau+du*dtaudu3, p+du*dpdu3, xn, iwave, dtaudu4, dpdu4); u += du;