Skip to content

Commit

Permalink
fix compilation of the exact Riemann solver (#2745)
Browse files Browse the repository at this point in the history
  • Loading branch information
zingale authored Feb 10, 2024
1 parent f1a879b commit c5c466b
Show file tree
Hide file tree
Showing 5 changed files with 127 additions and 111 deletions.
18 changes: 11 additions & 7 deletions Util/exact_riemann/exact_riemann.H
Original file line number Diff line number Diff line change
@@ -1,12 +1,16 @@
#ifndef EXACT_RIEMANN_H
#define EXACT_RIEMANN_H

#include <AMReX_REAL.H>

#include <fstream>

#include <riemann_sample.H>
#include <riemann_star_state.H>
#include <riemann_support.H>

using namespace amrex::literals;

AMREX_INLINE
void
exact_riemann() {
Expand All @@ -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;


Expand Down Expand Up @@ -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,
Expand All @@ -62,7 +66,7 @@ exact_riemann() {

// This follows from the discussion around C&G Eq. 15

Real dx = (problem::xmax - problem::xmin) / static_cast<Real>(problem::npts);
amrex::Real dx = (problem::xmax - problem::xmin) / static_cast<amrex::Real>(problem::npts);

// loop over xi space

Expand All @@ -72,9 +76,9 @@ exact_riemann() {

// compute xi = x/t -- this is the similarity variable for the
// solution
Real x = problem::xmin + (static_cast<Real>(i) - 0.5_rt) * dx;
amrex::Real x = problem::xmin + (static_cast<amrex::Real>(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,
Expand All @@ -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;
}
Expand Down
24 changes: 12 additions & 12 deletions Util/exact_riemann/main.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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();

}
36 changes: 20 additions & 16 deletions Util/exact_riemann/riemann_sample.H
Original file line number Diff line number Diff line change
@@ -1,16 +1,20 @@
#ifndef RIEMANN_SAMPLE_MODULE_H
#define RIEMANN_SAMPLE_MODULE_H

#include <AMReX_REAL.H>

#include <riemann_support.H>

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

Expand All @@ -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;
Expand All @@ -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
Expand All @@ -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;
Expand Down Expand Up @@ -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
Expand All @@ -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 {
Expand All @@ -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;
Expand Down
46 changes: 25 additions & 21 deletions Util/exact_riemann/riemann_star_state.H
Original file line number Diff line number Diff line change
@@ -1,17 +1,21 @@
#ifndef RIEMANN_STAR_STATE_H
#define RIEMANN_STAR_STATE_H

#include <AMReX_REAL.H>

#include <riemann_support.H>

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
Expand All @@ -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;
Expand All @@ -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
Expand Down Expand Up @@ -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;
Expand All @@ -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

Expand All @@ -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);
}

Expand All @@ -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) {
Expand Down
Loading

0 comments on commit c5c466b

Please sign in to comment.