Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

fix compilation of the exact Riemann solver #2745

Merged
merged 1 commit into from
Feb 10, 2024
Merged

fix compilation of the exact Riemann solver #2745

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
18 changes: 11 additions & 7 deletions Util/exact_riemann/exact_riemann.H
View file
Open in desktop
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
View file
Open in desktop
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
View file
Open in desktop
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
View file
Open in desktop
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