#include "nuc_eos.h"
#include "../con2prim.h"
// constants (required by EOS_Omni)
const double kBerg = 1.38064852e-16;
const double amu = 1.6605402e-24;
// conversion factors between cgs and M_Sun = c = G = 1
// (temperature between K and MeV)
// see EOS_Omni/doc/units.py
const double rho_gf = 1.61887093132742e-18;
const double press_gf = 1.80123683248503e-39;
const double eps_gf = 1.11265005605362e-21;
const double temp_gf = 8.617330337217213e-11;
// Inverses of the numbers above, calculated manually instead of by
// the compiler
const double inv_rho_gf = 6.17714470405638e17;
const double inv_press_gf = 5.55174079257738e38;
const double inv_eps_gf = 8.98755178736818e20;
const double inv_temp_gf = 1.1604522060401007e10;
// EOS quantities for ideal gas, required by EOS_Omni
double EOS_const[5];
// con2prim parameters
struct c2p_steer c2p;
extern void eos_omni_helmholtzreadtable_(char * );
void collect_metric(struct metric *h, const double hxx, const double hxy, const double hxz,
const double hyy, const double hyz, const double hzz){
// get the metric and inverse metric components
for (int i=0;i<4;i++) {
h->lo[0][i] = 0.0;
h->lo[i][0] = 0.0;
h->up[0][i] = 0.0;
h->up[i][0] = 0.0;
}
h->lo[1][1] = hxx;
h->lo[1][2] = hxy;
h->lo[1][3] = hxz;
h->lo[2][2] = hyy;
h->lo[2][3] = hyz;
h->lo[3][3] = hzz;
h->lo[2][1] = h->lo[1][2];
h->lo[3][1] = h->lo[1][3];
h->lo[3][2] = h->lo[2][3];
h->lo_det = hxx*hyy*hzz + 2.0*hxy*hyz*hxz - hxz*hxz*hyy - hxx*hyz*hyz - hxy*hxy*hzz;
h->lo_sqrt_det = sqrt(h->lo_det);
h->up_det = 1.0 / h->lo_det;
h->up[1][1] = (-hyz*hyz + hyy*hzz) * h->up_det;
h->up[1][2] = (hxz*hyz - hxy*hzz) * h->up_det;
h->up[2][2] = (-hxz*hxz + hxx*hzz) * h->up_det;
h->up[1][3] = (-hxz*hyy + hxy*hyz) * h->up_det;
h->up[2][3] = (hxy*hxz - hxx*hyz) * h->up_det;
h->up[3][3] = (-hxy*hxy + hxx*hyy) * h->up_det;
h->up[2][1] = h->up[1][2];
h->up[3][1] = h->up[1][3];
h->up[3][2] = h->up[2][3];
}
int main(int argc, char *argv[]){
if (argc != 7) {
printf("size of argv= %i\n", argc );
printf("Usage: name_executable Recovery_scheme EOS param1 param2 Ye rhoT_key");
printf("Incorrect number of arguments. Terminating...\n");
exit(1);
}
// general EOS settings
c2p.eoskey_polytrope = 1;
c2p.eos_prec = 1e-10;
// tolerance for convergence in recovery scheme
c2p.tol_x = 5e-9;
c2p.tol_x_retry = 1e2*c2p.tol_x; // reduced tolerance for retry if first attempt failed
// maximum number of iterations for c2p scheme
c2p.max_iterations = 300;
c2p.extra_iterations = 0;
// atmosphere settings
c2p.rho_atmo = 6e-16; // in M_Sun = c = G = 1
c2p.rho_atmo_tol = 0.01; // relative tolerance
c2p.T_atmo = 9e-07; // in MeV, 8.62e-06 MeV = 1e5K
c2p.Ye_atmo = 0.1; // electron fraction
c2p.retain_B_atmo = false; // keep the magnetic field in the atmo?
// number of primitive and conserved quantities
c2p.numprims = 14; // rho, vel123, eps, B123, ye, temp, press, W, ent, abar
c2p.numcons = 10; // D, S123, tau, B123, ye_con, temp
// stricly enforce v^2 < 1
c2p.enforce_v2 = 1;
// define relative amplitude for random perturbations
// applied to primitives before starting recovery process
const double perturbation=0.05;
const double perturbation_W=0.05; // only for Lorentz factor
// tolerance defining successful recovery of primitive variables
// abs((recovered-original)/original)<tol_prim
double tol_prim = 10.*c2p.tol_x;
// Choose angle d between B and v
// See Section 4.1 of Cerda-Duran et al. 2008
// ALIGNED = 0: d = pi/2
// ALIGNED = 1: d = 0
const int ALIGNED = 1;
// EOS parameters for ideal gas and Hybrid EOS
for(int i=0;i<5;i++)
EOS_const[i] = 0.0;
EOS_const[2] = 2.5; //gamma2
EOS_const[3] = 1.5; //gammath
EOS_const[4] = ((2.0e14)*rho_gf); //rho_nuc
// set recovery method
c2p.c2p_method = atoi(argv[1]);
c2p.c2p_method_backup = 5;
c2p.use_c2p_method_backup = false;
if ((c2p.c2p_method >= 8) || (c2p.c2p_method_backup >= 8)){
printf("Invalid c2p method(s). Terminating...\n");
exit(1);
}
// set EOS
c2p.eoskey = atoi(argv[2]);
if (c2p.eoskey<4) {
// ideal gas and Hybrid EOS
c2p.evolve_T = 0;
c2p.evolve_Ye = 0;
EOS_const[0] = 4./3.;
} else {
c2p.evolve_T = 1;
c2p.evolve_Ye = 1;
}
// set electron fraction
double ye = strtod(argv[5], NULL);
// set recovery test type
// plot rho vs. T, or p_mag/p vs. W-1
double logWminus1 = 0.0;
double W = 0.0;
double logPmagP = 0.0;
double rho_cgs = 0.0;
double temp_mev = 0.0;
bool rhoT_key = atoi(argv[6]);
if (rhoT_key){
logWminus1 = strtod(argv[3], NULL);
W = 1.0 + pow(10,logWminus1);
logPmagP = strtod(argv[4], NULL);
printf("RS= %i, EOS= %i, logWminus1= %e, logPmagP= %e, ye= %e\n",c2p.c2p_method, c2p.eoskey, logWminus1, logPmagP, ye);
} else {
rho_cgs = strtod(argv[3], NULL);
temp_mev = strtod(argv[4], NULL);
printf("RS= %i, EOS= %i, rho= %e, temp= %e, ye= %e\n",c2p.c2p_method, c2p.eoskey, rho_cgs, temp_mev, ye);
}
// choose size of code test
// sqrt(npoints) points in each variable (rho,T or W-1,p_mag/p)
int npoints;
#if DEBUG
npoints = pow(2,4);
#else
npoints = pow(2,12);
#endif
// choose range for variables
double rho_min, rho_max, temp_min, temp_max;
// set rho in g/cm^3
// vary: 1e3 g/cm^3 < rho < 1e15 g/cm^3
rho_min = 1.0e3;
rho_min *= 1.1; // add 10% to not avoid problems at table boundary
if (c2p.eoskey==8)
rho_max = 9.0e11;
else
rho_max = 1.0e15;
// set temp in K
// vary 1e5*Kelvin < temp < 50 MeV
if (c2p.eoskey==4) {
temp_min = 1.16e8;
temp_min *= 1.1; // add 10% to not avoid problems at table boundary
} else {
temp_min = 1e5;
}
temp_max = 50.0*inv_temp_gf;
// set Wminus1, PmagP
double logWminus1_min, logWminus1_max, logPmagP_min, logPmagP_max;
// vary W-1 from 1e-4 to 1e4
logWminus1_min = -4.0;
logWminus1_max = 4.0;
// vary P_mag/P from 1e-8 to 1e10
logPmagP_min = -8.0;
logPmagP_max = 10.0;
printf("rho_min = %g\n", rho_min);
printf("rho_max = %g\n", rho_max);
printf("temp_min = %g\n", temp_min);
printf("temp_max = %g\n", temp_max);
printf("logWminus1_min = %g\n", logWminus1_min);
printf("logWminus1_max = %g\n", logWminus1_max);
printf("logPmagP_min = %g\n", logPmagP_min);
printf("logPmagP_max = %g\n", logPmagP_max);
// READ EOS TABLES
// tabulated nuclear EOS
if(c2p.eoskey == 4){
nuc_eos_C_ReadTable("eostable.h5");
printf("nrho: %d\n",nrho);
printf("ntemp: %d\n",ntemp);
printf("nye: %d\n\n",nye);
}
// Helmholtz EOS
if(c2p.eoskey == 8){
char helmeos_table_name[50];
strcpy(helmeos_table_name, "helm_table_FLASH.dat");
eos_omni_helmholtzreadtable_(helmeos_table_name);
}
// further EOS settings
double rhomi,rhoma,tmi,tma,yemi,yema,emi,ema,pmi,pma;
EOS_Omni_get_eos_limits(c2p.eoskey,&rhomi,&rhoma,&tmi,&tma,&yemi,&yema,&emi,&ema,&pmi,&pma);
c2p.eos_rho_min = rhomi;
c2p.eos_rho_max = rhoma;
c2p.eos_temp_min = tmi;
c2p.eos_temp_max = tma;
c2p.eos_eps_min = emi; //-1.518793e-03; //emi;
c2p.eos_eps_max = ema;
c2p.eos_press_min = pmi;
c2p.eos_press_max = pma;
printf("rma = %g,rmi = %g,temp_ma = %g,temp_mi = %g,eps_ma = %g,eps_mi = %g,pmi=%g,pma=%g\n",rhoma,rhomi,tma,tmi,ema,emi,pmi,pma);
// Open output file
FILE *fp;
char fname[256];
if (rhoT_key)
sprintf(fname, "c2p_output_A%i_RS-%i_EOS-%i_logWminus1-%.2g_logPmagP-%.2g_YE-%.1f.dat", ALIGNED, c2p.c2p_method, c2p.eoskey, logWminus1, logPmagP, ye);
else
sprintf(fname, "c2p_output_A%i_RS-%i_EOS-%i_RHO-%.2g_T-%.2g_YE-%.1f.dat", ALIGNED, c2p.c2p_method, c2p.eoskey, rho_cgs, temp_mev, ye);
fp = fopen(fname, "w");
if (fp == NULL) {
printf("Couldn't open output file for writing.\n");
exit(1);
}
// Write setup parameters to output file
if(c2p.eoskey==1)
fprintf(fp,"# Polytropic EOS\n");
else if(c2p.eoskey==2)
fprintf(fp,"# Ideal gas EOS\n");
else if(c2p.eoskey==3)
fprintf(fp,"# Hybrid EOS\n");
else if(c2p.eoskey==4)
fprintf(fp,"# Tabulated EOS\n");
else if(c2p.eoskey==8)
fprintf(fp,"# Helmholtz EOS\n");
else{
printf("EOS not implemented. Terminating...\n");
exit(1);
}
if(ALIGNED==0)
fprintf(fp,"# NOT ALIGNED: cos^2(delta) = 0\n");
else if(ALIGNED==1)
fprintf(fp,"# ALIGNED: cos^2(delta) = 1\n");
else {
printf("ALIGNED = %i parameter not known. Aborting.", ALIGNED);
exit(1);
}
if (rhoT_key) {
fprintf(fp,"# rhoT_key = 1");
fprintf(fp,"# logWminus1 = %e\n",logWminus1);
fprintf(fp,"# logPmagP = %e\n",logPmagP);
} else {
fprintf(fp,"# rhoT_key = 0");
fprintf(fp,"# rho_cgs = %e\n",rho_cgs);
fprintf(fp,"# temp_mev = %e\n",temp_mev);
}
// define the metric
struct metric g;
collect_metric(&g, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0);
if(g.lo_det <= 0.0){
printf("det(g) <= 0. Con2prim does not make sense at this point.");
}
// allocate primitive and conservative arrays
double* prim = malloc(sizeof(double) * c2p.numprims);
double* prim_all = malloc(sizeof(double) * npoints * c2p.numprims);
double* original_prim = malloc(sizeof(double) * npoints * c2p.numprims);
if (prim_all == NULL) {
printf("Could not allocate enough much memory for prims...\n");
exit(1);
}
double* con = malloc(sizeof(double) * c2p.numcons);
double* con_all= malloc(sizeof(double) * npoints * c2p.numcons);
if (con_all == NULL) {
printf("Could not allocate enough memory for cons...\n");
exit(1);
}
// initialize prims and cons
for(int i=0;i<(npoints*c2p.numprims);i++){
prim_all[i]=0.0;
}
for(int i=0;i<(c2p.numprims);i++){
prim[i]=0.0;
}
for(int i=0;i<(npoints*c2p.numcons);i++){
con_all[i]=0.0;
}
for(int i=0;i<(c2p.numcons);i++){
con[i]=0.0;
}
//------------------------------------------------------------
// Generate grid of primitive quantities and compute conserved
// quantities for subsequent inversion
//------------------------------------------------------------
double n1,n2;
double v,vx,vy,vz;
double rhoexp;
double delrho;
double exp1;
double del1;
double tempexp;
double deltemp;
double tempexp_i;
double rhoexp_i;
// define range of rho or W-1 depending on rhoT_key
if (rhoT_key) {
rhoexp_i = log10(rho_min);
rhoexp = rhoexp_i;
delrho = (log10(rho_max)-rhoexp_i)/(sqrt(1.0*npoints)-1.0);
} else {
exp1 = logWminus1_min;
del1 = (logWminus1_max-logWminus1_min)/(sqrt(1.0*npoints)-1.0);
}
printf("\n---------------------------------------------------\n");
printf("Generating conservatives: prim2con\n");
printf("---------------------------------------------------\n\n");
for(n1=0;n1<sqrt(1.0*npoints);n1++){
double exp2;
double del2;
if (rhoT_key) {
tempexp_i = log10(temp_min*temp_gf);
tempexp = tempexp_i;
deltemp = (log10(temp_max*temp_gf)-tempexp_i)/(sqrt(1.0*npoints)-1.0);
} else {
exp2 = logPmagP_min;
del2 = (logPmagP_max-logPmagP_min)/(sqrt(1.0*npoints)-1.0);
}
for(n2=0;n2<sqrt(1.0*npoints);n2++){
// Translate 2D indexing to 1D index
int i = n1*sqrt(1.0*npoints)+n2;
printf("**************************************************\n");
printf("Creating data set No.: %d\n", i+1);
printf("n1,n2 = %g,%g\n", n1,n2);
if (rhoT_key) {
rho_cgs = pow(10, rhoexp);
temp_mev = pow(10,tempexp);
} else {
// (W-1) = 10^(exp1)
W = 1.0+pow(10.0,exp1);
}
// Velocity magnitude
v = sqrt(1.0-1.0/(W*W));
// Randomly choose vx, vy, vz
vx = v*((double)rand())/((double)RAND_MAX);
vy = sqrt(v*v - vx*vx)*((double)rand())/((double)RAND_MAX);
vz = sqrt(v*v - vx*vx - vy*vy);
// Verify the magnitude
if(fabs(v*v - (vx*vx+vy*vy+vz*vz))>1.0e-10){
printf("v problem!!!\n");
exit(1);
}
// Set primitive variables
prim_all[i*c2p.numprims+RHO] = rho_cgs*rho_gf;
prim_all[i*c2p.numprims+v1_cov] = vx;
prim_all[i*c2p.numprims+v2_cov] = vy;
prim_all[i*c2p.numprims+v3_cov] = vz;
prim_all[i*c2p.numprims+TEMP] = temp_mev;
prim_all[i*c2p.numprims+YE] = ye;
prim_all[i*c2p.numprims+WLORENTZ] = W;
for(int k=0;k<c2p.numprims;k++){
prim[k] = prim_all[i*c2p.numprims+k];
}
// compute eps, press
double epsilon, pressure;
epsilon = eps_from_prim(prim);
prim_all[i*c2p.numprims+EPS] = epsilon;
prim[EPS] = epsilon;
pressure = P_from_prim(prim);
prim_all[i*c2p.numprims+PRESS] = pressure;
prim[PRESS] = pressure;
// check P > 0
if(pressure<0.0){
printf("negative initial pressure: press = %e\n",pressure);
printf("Terminating...\n");
exit(1);
}
// set up magnetic field
double B;
if (rhoT_key) {
B = sqrt(2.0*pow(10.0,logPmagP)*pressure);
} else {
B = sqrt(2.0*pow(10.0,exp2)*pressure);
}
double Bxhat;
double Byhat;
double Bzhat;
if(ALIGNED==1){
double rand_temp = ((double)rand())/((double)RAND_MAX);
if(rand_temp<0.5){
Bxhat = vx/v;
Byhat = vy/v;
Bzhat = vz/v;
}
else{
Bxhat = -1.0*vx/v;
Byhat = -1.0*vy/v;
Bzhat = -1.0*vz/v;
}
}
else{
Bxhat = ((double)rand())/((double)RAND_MAX);
Byhat = sqrt(1.0 - Bxhat*Bxhat)*((double)rand())/((double)RAND_MAX);
Bzhat = (-Bxhat*vx-Byhat*vy)/vz;
double normB = sqrt(Bxhat*Bxhat + Byhat*Byhat + Bzhat*Bzhat);
Bxhat = Bxhat/normB;
Byhat = Byhat/normB;
Bzhat = Bzhat/normB;
// Verify orthogonality
if(fabs((Bxhat*vx+Byhat*vy+Bzhat*vz))>1.0e-10){
printf("Orthogonality problem!\n");
printf("Terminating...\n");
exit(1);
}
}
prim_all[i*c2p.numprims+B1_con] = Bxhat*B;
prim_all[i*c2p.numprims+B2_con] = Byhat*B;
prim_all[i*c2p.numprims+B3_con] = Bzhat*B;
// Verify the magnitude
if(fabs(1.0 -(Bxhat*Bxhat+Byhat*Byhat+Bzhat*Bzhat))>1.0e-10){
printf("B problem!!!\n");
printf("%e\n",Bxhat*Bxhat+Byhat*Byhat+Bzhat*Bzhat);
printf("Terminating...\n");
exit(1);
}
// Save original primitive variables
// and set primitive variables for present prim2con conversion
for(int k=0;k<c2p.numprims;k++){
original_prim[i*c2p.numprims+k] = prim_all[i*c2p.numprims+k];
prim[k] = prim_all[i*c2p.numprims+k];
printf("Save prims:\n");
printf("%e\t%e\n",prim_all[i*c2p.numprims+k],original_prim[i*c2p.numprims+k]);
}
// Set up conservatives by converting prim to con
prim2con_MHD_(prim, con, g.up, g.lo);
// update con_all
for(int k=0;k<c2p.numcons;k++){
con_all[i*c2p.numcons+k] = con[k];
}
#if DEBUG
// print primitive and conservative data sets
printf("\nprim2con Results (i=%d):\nprims\t\tcons\n",i);
for(int k=0;k<c2p.numprims;k++){
if (k>=c2p.numcons) {
printf("%e\t%e\n",prim[k],0.0);
} else{
printf("%e\t%e\n",prim[k],con[k]);
}
}
for(int k=0;k<c2p.numcons;k++){
printf("final2 p2c:\n");
printf("con[%i] \t con_all[%i]\n", k, i*c2p.numcons+k);
printf("%g \t %g\n", con[k], con_all[i*c2p.numcons+k]);
}
#endif
if (rhoT_key) {
tempexp = tempexp + deltemp;
} else {
exp2 = exp2 + del2;
}
}
if (rhoT_key) {
rhoexp = rhoexp + delrho;
} else {
exp1 = exp1 + del1;
}
}
//------------------------------------------------------------
// Conservative to Primitive Recovery
//------------------------------------------------------------
// Array to store number of iterations required for con2prim
// at each point
int count[npoints]; // number of iterations
int nEOScalls[npoints]; // number of EOS calls
double error[npoints]; // relative accuracy
bool retry[npoints]; // record whether c2p has been carried out with reduced tolerance
int failed[npoints]; // record whether c2p has failed
int c2p_keyerr[npoints]; // record c2p error code
// Reset primitive variables/setting up initial guesses
for(int i=0;i<npoints;i++){
for(int k=0;k<c2p.numprims;k++){
int randnum = ((int)rand() % 2); // 0 or 1
randnum = randnum*2; // 0 or 2
randnum = randnum - 1; // -1 or 1
if ((k == RHO) && (c2p.c2p_method == 3)){
prim_all[i*c2p.numprims+k] *= (1.0 + ((double)randnum)*perturbation_W);
} else if (k == WLORENTZ){
// perturb W-1
prim_all[i*c2p.numprims+k] = 1.0 + (prim_all[i*c2p.numprims+k] - 1.0) * (1.0 + ((double)randnum)*perturbation_W);
} else {
// perturb all other primitives
prim_all[i*c2p.numprims+k] *= (1.0 + ((double)randnum)*perturbation);
}
}
// checks on initial guesses
// ensure they are physical
// density
if (prim_all[i*c2p.numprims+RHO] < c2p.eos_rho_min){
prim_all[i*c2p.numprims+RHO] = (1.0+perturbation)*c2p.eos_rho_min;
}
if (prim_all[i*c2p.numprims+RHO] > c2p.eos_rho_max){
prim_all[i*c2p.numprims+RHO] = (1.0-perturbation)*c2p.eos_rho_max;
}
// velocities
// make sure that initial guesses satisfy v<c
vx = prim_all[i*c2p.numprims+v1_cov];
vy = prim_all[i*c2p.numprims+v2_cov];
vz = prim_all[i*c2p.numprims+v3_cov];
if(sqrt(vx*vx+vy*vy+vz*vz) >= 1.0){
printf("v problem with initial guesses...!\n");
printf("...reset to original values!\n");
prim_all[i*c2p.numprims+v1_cov] = original_prim[i*c2p.numprims+v1_cov];
prim_all[i*c2p.numprims+v2_cov] = original_prim[i*c2p.numprims+v2_cov];
prim_all[i*c2p.numprims+v3_cov] = original_prim[i*c2p.numprims+v3_cov];
vx = prim_all[i*c2p.numprims+v1_cov];
vy = prim_all[i*c2p.numprims+v2_cov];
vz = prim_all[i*c2p.numprims+v3_cov];
if(sqrt(vx*vx+vy*vy+vz*vz) >= 1.0){
printf("Reset of v did not work... something went wrong. Terminating...\n");
exit(1);
}
}
// specific internal energy
if (prim_all[i*c2p.numprims+EPS] < c2p.eos_eps_min){
prim_all[i*c2p.numprims+EPS] = fabs((1.0+perturbation)*c2p.eos_eps_min - c2p.eos_eps_min) + c2p.eos_eps_min;
}
if (prim_all[i*c2p.numprims+EPS] > c2p.eos_eps_max){
prim_all[i*c2p.numprims+EPS] = (1.0-perturbation)*c2p.eos_eps_max;
}
}
// Start the con2prim process
printf("\n---------------------------------------------------\n");
printf("Starting con2prim\n");
printf("---------------------------------------------------\n\n");
for(int i=0;i<npoints;i++){
printf("**************************************************\n");
printf("Starting inversion No. %d\n", i+1);
// set primitives (initial guesses) and conservatives
for(int k=0;k<c2p.numprims;k++){
prim[k] = prim_all[i*c2p.numprims+k];
}
for(int k=0;k<c2p.numcons;k++){
con[k] = con_all[i*c2p.numcons+k];
}
// recover Ye from conservatives
prim[YE]=con[YE]/con[D];
#if DEBUG
// check initial guesses
printf("primitives and conservatives:\n");
printf("prim guesses\tprims\t\tcons\n");
for(int k=0;k<c2p.numprims;k++)
{
if (k< c2p.numcons)
{
printf("%e\t%e\t%e\n",prim[k],original_prim[i*c2p.numprims+k],con_all[i*c2p.numcons+k]);
} else {
printf("%e\t%e\n",prim[k],original_prim[i*c2p.numprims+k]);
}
}
#endif
// con2prim process
struct c2p_report report;
bool excise = false;
bool c2p_grace = false;
failed[i] = 0;
c2p_keyerr[i] = 0;
retry[i] = false;
// call con2prim process
con2prim_MHD_(prim,con,g.up,g.lo,excise,c2p_grace,&report);
count[i] = report.count;
retry[i] = report.retry;
nEOScalls[i] = report.nEOScalls;
if (report.failed) {
failed[i] = 1;
c2p_keyerr[i] = report.c2p_keyerr;
printf("FAILED: con2prim reconstruction failed for the following reason:\n");
printf("Report: %s\n", report.err_msg);
} else {
printf("SUCCESS: con2prim converged in %i iterations!\n",count[i]);
}
// update primitives
for(int k=0;k<c2p.numprims;k++){
prim_all[i*c2p.numprims+k] = prim[k];
}
}
// free storage of EOS tables
if(c2p.eoskey==4){
free(alltables);
free(epstable);
free(presstable);
free(logrho);
free(logtemp);
free(yes);
}
//------------------------------------------------------------
// Analysis
//------------------------------------------------------------
// check for convergence and successful recovery of primitives
// count number of iterations for successful and failed recoveries
double nofailaverage = 0.0;
double average = 0.0;
int nofailtotal = 0;
int nofailconv = 0;
int converged;
int recovery;
double tol_prim_loc;
printf("\n---------------------------------------------------\n");
printf("Analysis of con2prim\n");
printf("---------------------------------------------------\n\n");
for(int i=0;i<npoints;i++){
// For each point, check for discrepancy in all primitive variables
// check that appropriate tolerance is used
if (retry[i]){
tol_prim_loc = tol_prim*c2p.tol_x_retry/c2p.tol_x;
}else{
tol_prim_loc = tol_prim;
}
// check for discrepancy
recovery = 1;
error[i] = 0.0;
for(int k=0;k<c2p.numprims;k++){
if((!(c2p.eoskey==1 && k == EPS))&& ( (k==RHO) || (k==v1_cov) || (k==v2_cov) || (k==v3_cov) || (k==EPS) ))
{
//tolerance enforced on rho, vel, eps
if((fabs((original_prim[i*c2p.numprims+k] -
prim_all[i*c2p.numprims+k])/original_prim[i*c2p.numprims+k])>tol_prim_loc) ||
(prim_all[i*c2p.numprims+k]!=prim_all[i*c2p.numprims+k])){
// If there is a discrepancy
#if DEBUG
printf("i,k, discrepancy: %i, %i, %e\n", i,k,fabs((original_prim[i*c2p.numprims+k] -
prim_all[i*c2p.numprims+k])/original_prim[i*c2p.numprims+k]));
printf("%g\t%g\n",prim_all[i*c2p.numprims+k],original_prim[i*c2p.numprims+k]);
#endif
recovery = 0;
}
error[i] += fabs((original_prim[i*c2p.numprims+k]-prim_all[i*c2p.numprims+k]) / original_prim[i*c2p.numprims+k]);
}
}
if (!(c2p.eoskey==1)) {
error[i] = 0.2 * error[i];
} else {
error[i] = 0.25 * error[i];
}
// count the total number of iterations
average += count[i];
// check if scheme converged and count the corresponding number of iterations
if (failed[i] != 1) {
nofailconv++;
}
// count the successful recoveries and the corresponding number of iterations
if(recovery==1){
nofailtotal++;
nofailaverage += count[i];
}
// check if eps<0 and v>c
int epsneg = 0;
if(prim_all[i*c2p.numprims+EPS]<0.0)
epsneg = 1;
int gtc = 0;
if((sqrt(prim_all[i*c2p.numprims+v1_cov]*prim_all[i*c2p.numprims+v1_cov]+prim_all[i*c2p.numprims+v2_cov]*prim_all[i*c2p.numprims+v2_cov]+prim_all[i*c2p.numprims+v3_cov]*prim_all[i*c2p.numprims+v3_cov]) -1.0)>=0.0)
gtc = 1;
// for each point, store its parameters (log(rho), log(temp)) or (log(W-1), log(Pmag/P))
// an recovery characteristics in output file
int n2 = (i%((int)sqrt(1.0*npoints)));
int n1 = (i-n2)/sqrt(1.0*npoints);
double del1;
double del2;
double exp1;
double exp2;
if (rhoT_key)
{
tempexp = tempexp_i+n2*deltemp;
rhoexp = rhoexp_i + n1*delrho;
} else {
del1 = (logWminus1_max-logWminus1_min)/(sqrt(1.0*npoints)-1.0);
del2 = (logPmagP_max-logPmagP_min)/(sqrt(1.0*npoints)-1.0);
exp1 = logWminus1_min+n1*del1;
exp2 = logPmagP_min+n2*del2;
}
if (rhoT_key) {
fprintf(fp,"%e %e %d %i %i %d %d %i %e %e %e %i\n",logWminus1,logPmagP,recovery,failed[i],count[i],epsneg,gtc, c2p_keyerr[i], rhoexp, tempexp, error[i], nEOScalls[i]);
} else {
fprintf(fp,"%e %e %d %i %i %d %d %i %e %e %e %i\n",exp1,exp2,recovery,failed[i],count[i],epsneg,gtc, c2p_keyerr[i], rhoexp, tempexp, error[i], nEOScalls[i]);
}
#if DEBUG
// print recovery results
printf("**************************************\n");
printf("Con2prim results, inversion No.: %d\n", i+1);
if(recovery==0)
{
if (rhoT_key) {
printf("\n*****FAILED recovery at\nrho\ttemp\n%e\t%e\nafter %d iterations\n\n",pow(10,rhoexp), pow(10,tempexp),count[i]);
} else {
printf("\n*****FAILED recovery at\nlog(W-1)\tlog(Pmag/P)\n%e\t%e\nafter %d iterations\n\n",exp1,exp2,count[i]);
}
}
else
{
if (rhoT_key) {
printf("\nSUCCESSFUL recovery at\nrho\ttemp\n%e\t%e\nafter %d iterations\n\n",pow(10,rhoexp), pow(10,tempexp),count[i]);
} else {
printf("\nSUCCESSFUL recovery at\nlog(W-1)\tlog(Pmag/P)\n%e\t%e\nafter %d iterations\n\n",exp1,exp2,count[i]);
}
}
printf("recovered prims\toriginal prims\tcons\n");
for(int k=0;k<c2p.numprims;k++)
{
if (k< c2p.numcons) {
printf("%e\t%e\t%e\n",prim_all[i*c2p.numprims+k],original_prim[i*c2p.numprims+k],con_all[i*c2p.numcons+k]);
} else {
printf("%e\t%e\n",prim_all[i*c2p.numprims+k],original_prim[i*c2p.numprims+k]);
}
}
#endif
}
// print final report
printf("\nFINAL REPORT:\n");
printf("%d/%d recovered (%f percent)\n",nofailtotal,npoints,nofailtotal/(1.0*npoints)*100.0);
printf("%d/%d converged (%f percent)\n",nofailconv,npoints,nofailconv/(1.0*npoints)*100.0);
printf("Average number of iterations (all): %f\n",average/(1.0*npoints));
printf("Average number of iterations (recovered only): %f\n\n",nofailaverage/(1.0*nofailtotal));
if (rhoT_key==1){
printf("W=%e, logPmag/P= %e\n",W, logPmagP );
} else {
printf("rho=%e, temp_mev= %e\n",rho_cgs, temp_mev );
}
// free arrays
free(prim);
free(con);
free(prim_all);
free(con_all);
free(original_prim);
return 0;
}