From cd2a82809bc7f30e5ce664015b86dacf0da1dbb7 Mon Sep 17 00:00:00 2001 From: Aaron Ridley Date: Fri, 10 Nov 2023 08:40:25 -0500 Subject: [PATCH 01/23] FEAT: iErr (int) to didWork (bool) --- include/advance.h | 20 ++++---- src/advance.cpp | 119 +++++++++++++++++++++++++--------------------- 2 files changed, 74 insertions(+), 65 deletions(-) diff --git a/include/advance.h b/include/advance.h index b4444072..541dd683 100644 --- a/include/advance.h +++ b/include/advance.h @@ -32,15 +32,15 @@ **/ -int advance(Planets &planet, - Grid &gGrid, - Times &time, - Euv &euv, - Neutrals &neutrals, - Ions &ions, - Chemistry &chemistry, - Electrodynamics &electrodynamics, - Indices &indices, - Logfile &logfile); +bool advance(Planets &planet, + Grid &gGrid, + Times &time, + Euv &euv, + Neutrals &neutrals, + Ions &ions, + Chemistry &chemistry, + Electrodynamics &electrodynamics, + Indices &indices, + Logfile &logfile); #endif // INCLUDE_ADVANCE_H_ diff --git a/src/advance.cpp b/src/advance.cpp index 8b0afa43..355dd294 100644 --- a/src/advance.cpp +++ b/src/advance.cpp @@ -10,19 +10,19 @@ // so many inputs because it alters all of the states in the model. // ----------------------------------------------------------------------------- -int advance(Planets &planet, - Grid &gGrid, - Times &time, - Euv &euv, - Neutrals &neutrals, - Ions &ions, - Chemistry &chemistry, - Electrodynamics &electrodynamics, - Indices &indices, - Logfile &logfile) { - - int iErr = 0; - +bool advance(Planets &planet, + Grid &gGrid, + Times &time, + Euv &euv, + Neutrals &neutrals, + Ions &ions, + Chemistry &chemistry, + Electrodynamics &electrodynamics, + Indices &indices, + Logfile &logfile) { + + bool didWork = true; + std::string function = "advance"; static int iFunction = -1; report.enter(function, iFunction); @@ -38,14 +38,13 @@ int advance(Planets &planet, gGrid.calc_sza(planet, time); neutrals.calc_mass_density(); neutrals.calc_mean_major_mass(); - neutrals.calc_specific_heat(); neutrals.calc_concentration(); neutrals.calc_pressure(); neutrals.calc_bulk_velocity(); neutrals.calc_kappa_eddy(); - neutrals.calc_cMax(); + precision_t dtNeutral = neutrals.calc_dt(gGrid); precision_t dtIon = 100.0; time.calc_dt(dtNeutral, dtIon); @@ -57,7 +56,7 @@ int advance(Planets &planet, // first neutrals.calc_scale_height(gGrid); - neutrals.set_bcs(gGrid, time, indices); + didWork = neutrals.set_bcs(gGrid, time, indices); if (input.get_nAltsGeo() > 1) neutrals.advect_vertical(gGrid, time); @@ -65,60 +64,70 @@ int advance(Planets &planet, // ------------------------------------ // Calculate source terms next: - iErr = calc_euv(planet, - gGrid, - time, - euv, - neutrals, - ions, - indices); + if (didWork) + didWork = calc_euv(planet, + gGrid, + time, + euv, + neutrals, + ions, + indices); - iErr = electrodynamics.update(planet, - gGrid, - time, - ions); - calc_ion_neutral_coll_freq(neutrals, ions); - ions.calc_ion_drift(neutrals, gGrid, time.get_dt()); + if (didWork) + didWork = electrodynamics.update(planet, + gGrid, + time, + ions); + + if (didWork) { + calc_ion_neutral_coll_freq(neutrals, ions); + ions.calc_ion_drift(neutrals, gGrid, time.get_dt()); - calc_aurora(gGrid, neutrals, ions); + calc_aurora(gGrid, neutrals, ions); - // Calculate some neutral source terms: - neutrals.calc_conduction(gGrid, time); - chemistry.calc_chemistry(neutrals, ions, time, gGrid); + // Calculate some neutral source terms: + neutrals.calc_conduction(gGrid, time); + chemistry.calc_chemistry(neutrals, ions, time, gGrid); - if (input.get_O_cooling()) - neutrals.calc_O_cool(); + if (input.get_O_cooling()) + neutrals.calc_O_cool(); - if (input.get_NO_cooling()) - neutrals.calc_NO_cool(); + if (input.get_NO_cooling()) + neutrals.calc_NO_cool(); - neutrals.calc_conduction(gGrid, time); - chemistry.calc_chemistry(neutrals, ions, time, gGrid); + neutrals.calc_conduction(gGrid, time); + chemistry.calc_chemistry(neutrals, ions, time, gGrid); - neutrals.vertical_momentum_eddy(gGrid); - calc_ion_collisions(neutrals, ions); - calc_neutral_friction(neutrals); + neutrals.vertical_momentum_eddy(gGrid); + calc_ion_collisions(neutrals, ions); + calc_neutral_friction(neutrals); - neutrals.add_sources(time); + neutrals.add_sources(time); - ions.calc_ion_temperature(neutrals, gGrid, time); - ions.calc_electron_temperature(neutrals, gGrid); + ions.calc_ion_temperature(neutrals, gGrid, time); + ions.calc_electron_temperature(neutrals, gGrid); - neutrals.exchange(gGrid); + neutrals.exchange(gGrid); - time.increment_time(); + time.increment_time(); - if (time.check_time_gate(input.get_dt_write_restarts())) { - report.print(3, "Writing restart files"); - neutrals.restart_file(input.get_restartout_dir(), DoWrite); - ions.restart_file(input.get_restartout_dir(), DoWrite); - time.restart_file(input.get_restartout_dir(), DoWrite); + if (time.check_time_gate(input.get_dt_write_restarts())) { + report.print(3, "Writing restart files"); + neutrals.restart_file(input.get_restartout_dir(), DoWrite); + ions.restart_file(input.get_restartout_dir(), DoWrite); + time.restart_file(input.get_restartout_dir(), DoWrite); + } } - iErr = output(neutrals, ions, gGrid, time, planet); + if (didWork) + didWork = output(neutrals, ions, gGrid, time, planet); - logfile.write_logfile(indices, neutrals, ions, gGrid, time); + if (didWork) + didWork = logfile.write_logfile(indices, neutrals, ions, gGrid, time); + if (!didWork) + report.error("Error in Advance!"); + report.exit(function); - return iErr; + return didWork; } From edc6d9f72238acc3a40548b113b30625b8ee5376 Mon Sep 17 00:00:00 2001 From: Aaron Ridley Date: Fri, 10 Nov 2023 08:41:26 -0500 Subject: [PATCH 02/23] FEAT: iErr (int) to didWork (bool) --- include/euv.h | 11 +++++------ src/calc_euv.cpp | 41 ++++++++++++++++++++++++----------------- src/euv.cpp | 31 +++++++++++++++---------------- 3 files changed, 44 insertions(+), 39 deletions(-) diff --git a/include/euv.h b/include/euv.h index f91743f0..138e74f3 100644 --- a/include/euv.h +++ b/include/euv.h @@ -107,28 +107,28 @@ class Euv { \param time The times within the model (dt is needed) \param indices Need the F107 and F107a **/ - int euvac(Times time, Indices indices); + bool euvac(Times time, Indices indices); /********************************************************************** \brief Compute the EUV spectrum given F107 and F107a \param time The times within the model (dt is needed) \param indices Need the F107 and F107a **/ - int solomon_hfg(Times time, Indices indices); + bool solomon_hfg(Times time, Indices indices); /********************************************************************** \brief Compute the EUV spectrum given F107 and F107a (new version) \param time The times within the model (dt is needed) \param indices Need the F107 and F107a **/ - int neuvac(Times time, Indices indices); + bool neuvac(Times time, Indices indices); /********************************************************************** \brief Scale the EUV spectrum given the star - planet distance \param planet needed to compute the star - planet distance \param time Needed to compute orbital position around star **/ - int scale_from_1au(Planets planet, Times time); + void scale_from_1au(Planets planet, Times time); /********************************************************************** \brief Pairs rows in the EUV CSV file with neutral and ions @@ -142,8 +142,7 @@ class Euv { \param neutrals Needs names of the neutrals, stores lines in Neutrals \param ions Needs names of the ions **/ - bool pair_euv(Neutrals &neutrals, - Ions ions); + bool pair_euv(Neutrals &neutrals, Ions ions); /********************************************************************** \brief Check to see if internal state of class is ok diff --git a/src/calc_euv.cpp b/src/calc_euv.cpp index bd13ca5e..ed51773e 100644 --- a/src/calc_euv.cpp +++ b/src/calc_euv.cpp @@ -13,15 +13,15 @@ // - calculate ionization and heating // ----------------------------------------------------------------------------- -int calc_euv(Planets planet, - Grid grid, - Times time, - Euv &euv, - Neutrals &neutrals, - Ions &ions, - Indices indices) { +bool calc_euv(Planets planet, + Grid grid, + Times time, + Euv &euv, + Neutrals &neutrals, + Ions &ions, + Indices indices) { - int iErr = 0; + bool didWork; if (time.check_time_gate(input.get_dt_euv())) { std::string function = "Euv::calc_euv"; @@ -37,14 +37,18 @@ int calc_euv(Planets planet, neutrals.calc_chapman(grid); if (euv.doUse) { - if (input.get_euv_model() == "euvac") - iErr = euv.euvac(time, indices); - else if (input.get_euv_model() == "neuvac") - iErr = euv.neuvac(time, indices); - else if (input.get_euv_model() == "hfg") - iErr = euv.solomon_hfg(time, indices); - - iErr = euv.scale_from_1au(planet, time); + // set didWork to false in order to catch bad euv models: + didWork = false; + std::string euvModel = mklower(input.get_euv_model()); + if (euvModel == "euvac") + didWork = euv.euvac(time, indices); + else if (euvModel == "neuvac") + didWork = euv.neuvac(time, indices); + else if (euvModel == "hfg") + didWork = euv.solomon_hfg(time, indices); + + if (didWork) + euv.scale_from_1au(planet, time); calc_ionization_heating(euv, neutrals, ions); } else @@ -53,7 +57,10 @@ int calc_euv(Planets planet, report.exit(function); } - return iErr; + if (!didWork) + report.error("Error in calc_euv! Check euv models."); + + return didWork; } // ----------------------------------------------------------------------------- diff --git a/src/euv.cpp b/src/euv.cpp index 8ddb2c59..0084eaea 100644 --- a/src/euv.cpp +++ b/src/euv.cpp @@ -312,9 +312,8 @@ bool Euv::pair_euv(Neutrals &neutrals, // Scale flux (intensity) at 1 AU to distance from the sun: // -------------------------------------------------------------------------- -int Euv::scale_from_1au(Planets planet, - Times time) { - int iErr = 0; +void Euv::scale_from_1au(Planets planet, + Times time) { precision_t d = planet.get_star_to_planet_dist(time); precision_t scale = 1.0 / (d * d); @@ -324,17 +323,17 @@ int Euv::scale_from_1au(Planets planet, for (int iWave = 0; iWave < nWavelengths; iWave++) wavelengths_intensity_top[iWave] = scale * wavelengths_intensity_1au[iWave]; - return iErr; + return; } // -------------------------------------------------------------------------- // Calculate EUVAC // -------------------------------------------------------------------------- -int Euv::euvac(Times time, - Indices indices) { +bool Euv::euvac(Times time, + Indices indices) { - int iErr = 0; + bool didWork = true; precision_t slope; std::string function = "Euv::euvac"; @@ -368,17 +367,17 @@ int Euv::euvac(Times time, } report.exit(function); - return iErr; + return didWork; } // -------------------------------------------------------------------------- // Calculate EUVAC // -------------------------------------------------------------------------- -int Euv::neuvac(Times time, - Indices indices) { +bool Euv::neuvac(Times time, + Indices indices) { - int iErr = 0; + bool didWork = true; precision_t slope; std::string function = "Euv::neuvac"; @@ -412,21 +411,21 @@ int Euv::neuvac(Times time, } report.exit(function); - return iErr; + return didWork; } // -------------------------------------------------------------------------- // Calculate HFG // -------------------------------------------------------------------------- -int Euv::solomon_hfg(Times time, - Indices indices) { +bool Euv::solomon_hfg(Times time, + Indices indices) { std::string function = "Euv::solomon_hfg"; static int iFunction = -1; report.enter(function, iFunction); - int iErr = 0; + bool didWork = true; precision_t r1; precision_t r2; @@ -455,7 +454,7 @@ int Euv::solomon_hfg(Times time, } report.exit(function); - return iErr; + return didWork; } // -------------------------------------------------------------------------- From c0808cc9e71142026b40dfd4cdc34c7a6dba0a1c Mon Sep 17 00:00:00 2001 From: Aaron Ridley Date: Fri, 10 Nov 2023 08:41:41 -0500 Subject: [PATCH 03/23] FEAT: iErr (int) to didWork (bool) --- include/calc_euv.h | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/include/calc_euv.h b/include/calc_euv.h index e350ade2..561930a2 100644 --- a/include/calc_euv.h +++ b/include/calc_euv.h @@ -20,13 +20,13 @@ // // ------------------------------------------------------------------------- -int calc_euv(Planets planet, - Grid grid, - Times time, - Euv &euv, - Neutrals &neutrals, - Ions &ions, - Indices indices); +bool calc_euv(Planets planet, + Grid grid, + Times time, + Euv &euv, + Neutrals &neutrals, + Ions &ions, + Indices indices); void calc_ionization_heating(Euv euv, Neutrals &neutrals, From b5040913f4c8ed18f9d41d3ffbaf2e74cdc6a94e Mon Sep 17 00:00:00 2001 From: Aaron Ridley Date: Fri, 10 Nov 2023 08:42:34 -0500 Subject: [PATCH 04/23] FEAT: iErr (int) to didWork (bool) --- include/electrodynamics.h | 8 ++++---- src/electrodynamics.cpp | 18 ++++++++---------- 2 files changed, 12 insertions(+), 14 deletions(-) diff --git a/include/electrodynamics.h b/include/electrodynamics.h index aeb1f485..827db32e 100644 --- a/include/electrodynamics.h +++ b/include/electrodynamics.h @@ -45,10 +45,10 @@ class Electrodynamics { \param ions Going to set the potential and aurora **/ - int update(Planets planet, - Grid gGrid, - Times time, - Ions &ions); + bool update(Planets planet, + Grid gGrid, + Times time, + Ions &ions); /************************************************************** \brief used in main.cpp to ensure electrodynamics times and diff --git a/src/electrodynamics.cpp b/src/electrodynamics.cpp index 16776576..14d757fa 100644 --- a/src/electrodynamics.cpp +++ b/src/electrodynamics.cpp @@ -33,10 +33,10 @@ Electrodynamics::Electrodynamics(Times time) { // Update Electrodynamics // ----------------------------------------------------------------------------- -int Electrodynamics::update(Planets planet, - Grid gGrid, - Times time, - Ions &ions) { +bool Electrodynamics::update(Planets planet, + Grid gGrid, + Times time, + Ions &ions) { std::string function = "Electrodynamics::update"; static int iFunction = -1; @@ -60,7 +60,7 @@ int Electrodynamics::update(Planets planet, } report.exit(function); - return 0; + return true; } // ----------------------------------------------------------------------------- @@ -182,12 +182,10 @@ arma_mat Electrodynamics::get_values(arma_mat matToInterpolateOn, // Get potential, electron energy flux, and electron average energy patterns // ----------------------------------------------------------------------------- -//average energy and eflux energy as well call to get_values - std::tuple Electrodynamics::get_electrodynamics(arma_cube magLat, -arma_cube magLocalTime) { + arma_mat, + arma_mat> Electrodynamics::get_electrodynamics(arma_cube magLat, + arma_cube magLocalTime) { arma_cube pot; arma_mat eflux; arma_mat avee; From a259758fa0be8aae0425f44cf88a7efca2778f02 Mon Sep 17 00:00:00 2001 From: Aaron Ridley Date: Fri, 10 Nov 2023 08:43:35 -0500 Subject: [PATCH 05/23] FEAT: iErr (int) to didWork (bool) --- include/neutrals.h | 6 +++--- src/neutrals.cpp | 10 +++++----- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/include/neutrals.h b/include/neutrals.h index 1226f2c8..be61813d 100644 --- a/include/neutrals.h +++ b/include/neutrals.h @@ -257,9 +257,9 @@ class Neutrals { \param time contains information about the current time \param indices used to help set initial conditions **/ - int initial_conditions(Grid grid, - Times time, - Indices indices); + bool initial_conditions(Grid grid, + Times time, + Indices indices); /********************************************************************** \brief temporary function to set neutral densities with in the model diff --git a/src/neutrals.cpp b/src/neutrals.cpp index 17659042..0b5c73c7 100644 --- a/src/neutrals.cpp +++ b/src/neutrals.cpp @@ -72,6 +72,7 @@ Neutrals::Neutrals(Grid grid, Indices indices) { int iErr; + bool didWork = true; species_chars tmp; int64_t nLons = grid.get_nLons(); @@ -139,15 +140,14 @@ Neutrals::Neutrals(Grid grid, // This gets a bunch of the species-dependent characteristics: iErr = read_planet_file(planet); - if (iErr > 0) - std::cout << "Error reading planet file!" << '\n'; + report.error("Error reading planet file!"); // This specifies the initial conditions for the neutrals: - iErr = initial_conditions(grid, time, indices); + didWork = initial_conditions(grid, time, indices); - if (iErr > 0) - std::cout << "Error in setting neutral initial conditions!" << '\n'; + if (!didWork) + report.error("Error in setting neutral initial conditions!"); return; } From 0353a3b5cae88bb51f02d242d1424fe73fe0947e Mon Sep 17 00:00:00 2001 From: Aaron Ridley Date: Fri, 10 Nov 2023 08:44:20 -0500 Subject: [PATCH 06/23] FEAT: use lower case to check for conditions --- src/neutrals_bcs.cpp | 29 ++++++++--- src/neutrals_ics.cpp | 121 +++++++++++++++++++++++-------------------- 2 files changed, 85 insertions(+), 65 deletions(-) diff --git a/src/neutrals_bcs.cpp b/src/neutrals_bcs.cpp index 4939236d..b8d497e0 100644 --- a/src/neutrals_bcs.cpp +++ b/src/neutrals_bcs.cpp @@ -32,10 +32,15 @@ bool Neutrals::set_bcs(Grid grid, if (input.get_nAltsGeo() > 1) { didWork = set_lower_bcs(grid, time, indices); - didWork = set_upper_bcs(grid); - calc_mass_density(); + if (didWork) + didWork = set_upper_bcs(grid); + if (didWork) + calc_mass_density(); } + if (!didWork) + report.error("issue with BCs!"); + report.exit(function); return didWork; } @@ -106,17 +111,19 @@ bool Neutrals::set_lower_bcs(Grid grid, static int iFunction = -1; report.enter(function, iFunction); - bool didWork = true; + bool didWork = false; json bcs = input.get_boundary_condition_types(); int64_t nGCs = grid.get_nGCs(); int64_t iSpecies, iAlt, iDir; + std::string bcsType = mklower(bcs["type"]); + //----------------------------------------------- // MSIS BCs - only works if FORTRAN is enabled! //----------------------------------------------- - if (bcs["type"] == "Msis") { + if (bcsType == "msis") { report.print(2, "Using MSIS for Boundary Conditions"); @@ -124,12 +131,13 @@ bool Neutrals::set_lower_bcs(Grid grid, if (!msis.is_ok()) { didWork = false; - + report.error("MSIS initialization not ok"); if (report.test_verbose(0)) { std::cout << "MSIS Boundary Conditions asked for, "; std::cout << "but MSIS is not compiled! Yikes!\n"; } - } + } else + didWork = true; msis.set_time(time); precision_t f107 = indices.get_f107(time.get_current()); @@ -174,7 +182,7 @@ bool Neutrals::set_lower_bcs(Grid grid, // Planet BCs - set to fixed constant values. //----------------------------------------------- - if (bcs["type"] == "Planet") { + if (bcsType == "planet") { report.print(2, "setting lower bcs to planet"); @@ -183,8 +191,8 @@ bool Neutrals::set_lower_bcs(Grid grid, species[iSpecies].density_scgc.slice(0). fill(species[iSpecies].lower_bc_density); } - temperature_scgc.slice(0).fill(initial_temperatures[0]); + didWork = true; } // fill the second+ grid cells with the bottom temperature: @@ -208,6 +216,11 @@ bool Neutrals::set_lower_bcs(Grid grid, } } + if (!didWork) { + report.error("issue with lower BCs!"); + report.error("maybe check boundaryconditions type : " + bcsType); + } + report.exit(function); return didWork; } diff --git a/src/neutrals_ics.cpp b/src/neutrals_ics.cpp index 9566598e..558a5b07 100644 --- a/src/neutrals_ics.cpp +++ b/src/neutrals_ics.cpp @@ -16,7 +16,7 @@ // file and fill with hydrostatic. // ----------------------------------------------------------------------------- -int Neutrals::initial_conditions(Grid grid, +bool Neutrals::initial_conditions(Grid grid, Times time, Indices indices) { @@ -24,7 +24,10 @@ int Neutrals::initial_conditions(Grid grid, static int iFunction = -1; report.enter(function, iFunction); - int iErr = 0; + // initialize didWork to false, so that we can catch if the initial + // conditions are not actually set + bool didWork = false; + int64_t iLon, iLat, iAlt, iA; precision_t alt, r; int64_t nAlts = grid.get_nZ(); @@ -33,70 +36,71 @@ int Neutrals::initial_conditions(Grid grid, if (input.get_do_restart()) { report.print(1, "Restarting! Reading neutral files!"); - bool DidWork = restart_file(input.get_restartin_dir(), DoRead); - - if (!DidWork) - std::cout << "Reading Restart for Neutrals Failed!!!\n"; + didWork = restart_file(input.get_restartin_dir(), DoRead); + if (!didWork) { + report.error("Reading Restart for Neutrals Failed!!!");\ + } } else { json ics = input.get_initial_condition_types(); + std::string icsType = mklower(ics["type"]); - if (ics["type"] == "Msis") { + if (icsType == "msis") { report.print(2, "Using MSIS for Initial Conditions"); Msis msis; if (!msis.is_ok()) { - iErr = 1; - - if (report.test_verbose(0)) { - std::cout << "MSIS Initial Conditions asked for, "; - std::cout << "but MSIS is not compiled! Yikes!\n"; - } - } - - msis.set_time(time); - precision_t f107 = indices.get_f107(time.get_current()); - precision_t f107a = indices.get_f107a(time.get_current()); - msis.set_f107(f107, f107a); - msis.set_ap(10.0); - msis.set_locations(grid.geoLon_scgc, - grid.geoLat_scgc, - grid.geoAlt_scgc); - - // This is just to check if MSIS is actually working: - if (msis.is_valid_species("Tn")) - // if it is, fill will temperature: - temperature_scgc = msis.get_cube("Tn"); - else - // if it is not, then fill with a value: - temperature_scgc.fill(initial_temperatures[0]); - - for (int iSpecies = 0; iSpecies < nSpecies; iSpecies++) { - if (report.test_verbose(3)) - std::cout << "Setting Species : " << species[iSpecies].cName << "\n"; - - if (msis.is_valid_species(species[iSpecies].cName)) { - if (report.test_verbose(3)) - std::cout << " Found in MSIS!\n"; - - species[iSpecies].density_scgc = - msis.get_cube(species[iSpecies].cName); - } else { - if (report.test_verbose(3)) - std::cout << " NOT Found in MSIS - setting constant\n"; - - species[iSpecies].density_scgc.slice(0). - fill(species[iSpecies].lower_bc_density); - fill_with_hydrostatic(iSpecies, 1, nAlts, grid); - } - - } - - } // type = Msis - - if (ics["type"] == "Planet") { + didWork = false; + report.error("MSIS Initial Conditions asked for, "); + report.error("but MSIS is not compiled/working!"); + } else { + didWork = true; + msis.set_time(time); + precision_t f107 = indices.get_f107(time.get_current()); + precision_t f107a = indices.get_f107a(time.get_current()); + msis.set_f107(f107, f107a); + msis.set_ap(10.0); + msis.set_locations(grid.geoLon_scgc, + grid.geoLat_scgc, + grid.geoAlt_scgc); + + // This is just to check if MSIS is actually working: + if (msis.is_valid_species("Tn")) + // if it is, fill will temperature: + temperature_scgc = msis.get_cube("Tn"); + else + // if it is not, then fill with a value: + temperature_scgc.fill(initial_temperatures[0]); + + for (int iSpecies = 0; iSpecies < nSpecies; iSpecies++) { + if (report.test_verbose(3)) + std::cout << "Setting Species : " + << species[iSpecies].cName << "\n"; + + if (msis.is_valid_species(species[iSpecies].cName)) { + if (report.test_verbose(3)) + std::cout << " Found in MSIS!\n"; + + species[iSpecies].density_scgc = + msis.get_cube(species[iSpecies].cName); + } else { + if (report.test_verbose(3)) + std::cout << " NOT Found in MSIS - setting constant\n"; + + species[iSpecies].density_scgc.slice(0). + fill(species[iSpecies].lower_bc_density); + fill_with_hydrostatic(iSpecies, 1, nAlts, grid); + } + + } // for species + } // msis init worked ok + } // type = msis + + if (icsType == "planet") { + + didWork = true; // --------------------------------------------------------------------- // This section assumes we want a hydrostatic solution given the @@ -165,8 +169,11 @@ int Neutrals::initial_conditions(Grid grid, } // type = planet } + if (!didWork) + report.error("Issue with initial conditions!"); + report.exit(function); - return iErr; + return didWork; } From 6ae7401a8f5e23fed0500f50c8cf022cea3ab348 Mon Sep 17 00:00:00 2001 From: Aaron Ridley Date: Fri, 10 Nov 2023 08:45:03 -0500 Subject: [PATCH 07/23] FEAT: function to make lower case --- include/transform.h | 3 +++ src/transform.cpp | 13 +++++++++++++ 2 files changed, 16 insertions(+) diff --git a/include/transform.h b/include/transform.h index d37f6735..03fa72e3 100644 --- a/include/transform.h +++ b/include/transform.h @@ -10,6 +10,9 @@ #include "aether.h" using namespace arma; +std::string mklower(std::string inString); +std::string mkupper(std::string inString); + void copy_cube_to_array(arma_cube cube_in, float *array_out); diff --git a/src/transform.cpp b/src/transform.cpp index f21caf05..6c5973dd 100644 --- a/src/transform.cpp +++ b/src/transform.cpp @@ -6,6 +6,19 @@ #include "aether.h" +// ----------------------------------------------------------------------- +// transform string to lower case +// ----------------------------------------------------------------------- + +std::string mklower(std::string inString) { + std::string outString = inString; + int64_t nChars = outString.length(); + for (int64_t iChar = 0; iChar < nChars; iChar++) + outString[iChar] = tolower(outString[iChar]); + return outString; +} + + // ----------------------------------------------------------------------- // copy from c++ vector to c-native array // ----------------------------------------------------------------------- From 2a539e7854b65c3b36cb7aa1257b98f6430e61a2 Mon Sep 17 00:00:00 2001 From: Aaron Ridley Date: Fri, 10 Nov 2023 08:45:53 -0500 Subject: [PATCH 08/23] FEAT: add cmember to filename --- include/tools.h | 7 +++++++ src/tools.cpp | 12 ++++++++++++ 2 files changed, 19 insertions(+) diff --git a/include/tools.h b/include/tools.h index b79ababe..91322cf5 100644 --- a/include/tools.h +++ b/include/tools.h @@ -15,6 +15,13 @@ struct mat_2x2{ arma_mat A22; }; +// ----------------------------------------------------------------------------- +// add cMember into a string just before last period +// ----------------------------------------------------------------------------- + +std::string add_cmember(std::string inString); + + // ---------------------------------------------------------------------- // Display an armadillo vector // ---------------------------------------------------------------------- diff --git a/src/tools.cpp b/src/tools.cpp index 841ca678..f8643e17 100644 --- a/src/tools.cpp +++ b/src/tools.cpp @@ -112,6 +112,18 @@ bool compare(precision_t value1, precision_t value2) { return false; } +// ----------------------------------------------------------------------------- +// add cMember into a string just before last period +// ----------------------------------------------------------------------------- + +std::string add_cmember(std::string inString) { + std::string outString = inString; + std::size_t found = outString.rfind("."); + if (found!=std::string::npos) + outString.replace(found, 1, "_" + cMember + "."); + return outString; +} + // ----------------------------------------------------------------------------- // Convert an integer to a zero-padded string // ----------------------------------------------------------------------------- From a971d861f5980d37980f631b48438f382e3aba67 Mon Sep 17 00:00:00 2001 From: Aaron Ridley Date: Fri, 10 Nov 2023 08:47:07 -0500 Subject: [PATCH 09/23] BUG: MSIS has temps below 200 --- src/solver_vertical_rusanov.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solver_vertical_rusanov.cpp b/src/solver_vertical_rusanov.cpp index fdf836a6..0b9258aa 100644 --- a/src/solver_vertical_rusanov.cpp +++ b/src/solver_vertical_rusanov.cpp @@ -306,7 +306,7 @@ void Neutrals::solver_vertical_rusanov(Grid grid, species[iSpecies].newVelocity_vcgc[2](iX, iY, iZ) = -100.0; } - if (newTemperature_scgc(iX, iY, iZ) < 200.0 || + if (newTemperature_scgc(iX, iY, iZ) < 100.0 || std::isnan(newTemperature_scgc(iX, iY, iZ))) { std::cout << "low temp found : " << iX << " " From 89e8e13696e854962c756542fc568ce647b607e6 Mon Sep 17 00:00:00 2001 From: Aaron Ridley Date: Fri, 10 Nov 2023 08:48:02 -0500 Subject: [PATCH 10/23] FEAT: add get integer settings from 2 keys --- include/inputs.h | 1 + src/inputs.cpp | 16 +++++++++++++++- 2 files changed, 16 insertions(+), 1 deletion(-) diff --git a/include/inputs.h b/include/inputs.h index 0fb03d75..0e3d990b 100644 --- a/include/inputs.h +++ b/include/inputs.h @@ -118,6 +118,7 @@ class Inputs { std::string get_settings_str(std::string key1); std::string get_settings_str(std::string key1, std::string key2); + int get_settings(std::string key1, std::string key2); bool check_settings(std::string key1, std::string key2); bool check_settings(std::string key1); std::string check_settings_str(std::string key1, std::string key2); diff --git a/src/inputs.cpp b/src/inputs.cpp index 46295d39..ce84b4f7 100644 --- a/src/inputs.cpp +++ b/src/inputs.cpp @@ -90,7 +90,10 @@ bool Inputs::write_restart() { // ----------------------------------------------------------------------- std::string Inputs::get_logfile() { - return check_settings_str("Logfile", "name"); + std::string logfile = check_settings_str("Logfile", "name"); + if (nMembers > 1) + logfile = add_cmember(logfile); + return logfile; } // ----------------------------------------------------------------------- @@ -226,6 +229,17 @@ std::string Inputs::get_settings_str(std::string key1, return value; } +int Inputs::get_settings(std::string key1, + std::string key2) { + int value = -1; + + if (settings.find(key1) != settings.end()) + if (settings.at(key1).find(key2) != settings.at(key1).end()) + value = settings.at(key1).at(key2); + + return value; +} + // ----------------------------------------------------------------------- // Check for missing settings // ----------------------------------------------------------------------- From 600ca02c045694b8934902d6b31dc8dcfa3268c6 Mon Sep 17 00:00:00 2001 From: Aaron Ridley Date: Fri, 10 Nov 2023 08:49:14 -0500 Subject: [PATCH 11/23] FEAT: iErr (int) to didWork (bool) --- include/output.h | 18 +++++------ src/main/main.cpp | 69 +++++++++++++++++++++++-------------------- src/output.cpp | 58 ++++++++++++++++++------------------ src/output_netcdf.cpp | 19 ++++++------ 4 files changed, 85 insertions(+), 79 deletions(-) diff --git a/include/output.h b/include/output.h index 8ae9d0a3..b328825f 100644 --- a/include/output.h +++ b/include/output.h @@ -122,22 +122,22 @@ class OutputContainer { /********************************************************************** \brief write a file with the information in the container **/ - void write(); + bool write(); /********************************************************************** \brief write a json header file with the information in the container **/ - int write_container_header(); + bool write_container_header(); /********************************************************************** \brief write a binary file with the information in the container **/ - int write_container_binary(); + bool write_container_binary(); /********************************************************************** \brief write a netcdf file with the information in the container **/ - int write_container_netcdf(); + bool write_container_netcdf(); /********************************************************************** \brief read from a file an load into the container @@ -222,11 +222,11 @@ class OutputContainer { \param planet information about the planet **/ -int output(const Neutrals &neutrals, - const Ions &ions, - const Grid &grid, - Times time, - const Planets &planet); +bool output(const Neutrals &neutrals, + const Ions &ions, + const Grid &grid, + Times time, + const Planets &planet); void output_binary_3d(std::ofstream &binary, arma_cube value); diff --git a/src/main/main.cpp b/src/main/main.cpp index 7d0f0c56..f10d5025 100644 --- a/src/main/main.cpp +++ b/src/main/main.cpp @@ -12,7 +12,7 @@ int main() { int iErr = 0; std::string sError; - bool DidWork = true; + bool didWork = true; Times time; @@ -36,13 +36,13 @@ int main() { throw std::string("quadtree initialization failed!"); // Initialize MPI and parallel aspects of the code: - DidWork = init_parallel(quadtree); - if (!DidWork) + didWork = init_parallel(quadtree); + if (!didWork) throw std::string("init_parallel failed!"); // Everything should be set for the inputs now, so write a restart file: - DidWork = input.write_restart(); - if (!DidWork) + didWork = input.write_restart(); + if (!didWork) throw std::string("input.write_restart failed!"); // Initialize the EUV system: @@ -58,9 +58,9 @@ int main() { // Initialize the indices, read the files, and perturb: Indices indices; - DidWork = read_and_store_indices(indices); + didWork = read_and_store_indices(indices); MPI_Barrier(aether_comm); - if (!DidWork) + if (!didWork) throw std::string("read_and_store_indices failed!"); // Perturb the inputs if user has asked for this @@ -72,9 +72,9 @@ int main() { input.get_nLatsGeo(), input.get_nAltsGeo(), nGeoGhosts); - DidWork = gGrid.init_geo_grid(quadtree, planet); + didWork = gGrid.init_geo_grid(quadtree, planet); MPI_Barrier(aether_comm); - if (!DidWork) + if (!didWork) throw std::string("init_geo_grid failed!"); // Calculate centripetal acceleration, since this is a constant @@ -103,8 +103,8 @@ int main() { } if (input.get_check_for_nans()) { - DidWork = neutrals.check_for_nonfinites(); - DidWork = ions.check_for_nonfinites(); + didWork = neutrals.check_for_nonfinites(); + didWork = ions.check_for_nonfinites(); } // ----------------------------------------------------------------- @@ -130,14 +130,16 @@ int main() { // If the user wants to restart, then get the time of the restart if (input.get_do_restart()) { report.print(1, "Restarting! Reading time file!"); - DidWork = time.restart_file(input.get_restartin_dir(), DoRead); - if (!DidWork) + didWork = time.restart_file(input.get_restartin_dir(), DoRead); + if (!didWork) throw std::string("Reading Restart for time Failed!!!\n"); } // This is for the initial output. If it is not a restart, this will go: if (time.check_time_gate(input.get_dt_output(0))) - iErr = output(neutrals, ions, gGrid, time, planet); + didWork = output(neutrals, ions, gGrid, time, planet); + if (!didWork) + throw std::string("output failed!"); // This is advancing now... We are not coupling, so set dt_couple to the // end of the simulation @@ -157,18 +159,21 @@ int main() { time.increment_intermediate(dt_couple); // Increment until the intermediate time: - while (time.get_current() < time.get_intermediate()) - iErr = advance(planet, - gGrid, - time, - euv, - neutrals, - ions, - chemistry, - electrodynamics, - indices, - logfile); - + while (time.get_current() < time.get_intermediate()) { + didWork = advance(planet, + gGrid, + time, + euv, + neutrals, + ions, + chemistry, + electrodynamics, + indices, + logfile); + if (!didWork) + throw std::string("Error in advance!"); + } + // Should write out some restart files every time we are done with // intermediate times. Just so when we restart, we know that we can // couple first thing and everything should be good. (Not sure if @@ -182,16 +187,16 @@ int main() { if (!time.check_time_gate(input.get_dt_write_restarts())) { report.print(3, "Writing restart files"); - DidWork = neutrals.restart_file(input.get_restartout_dir(), DoWrite); - if (!DidWork) + didWork = neutrals.restart_file(input.get_restartout_dir(), DoWrite); + if (!didWork) throw std::string("Writing Restart for Neutrals Failed!!!\n"); - DidWork = ions.restart_file(input.get_restartout_dir(), DoWrite); - if (!DidWork) + didWork = ions.restart_file(input.get_restartout_dir(), DoWrite); + if (!didWork) throw std::string("Writing Restart for Ions Failed!!!\n"); - DidWork = time.restart_file(input.get_restartout_dir(), DoWrite); - if (!DidWork) + didWork = time.restart_file(input.get_restartout_dir(), DoWrite); + if (!didWork) throw std::string("Writing Restart for time Failed!!!\n"); } diff --git a/src/output.cpp b/src/output.cpp index f23e0cfe..f35b7ee5 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -4,7 +4,6 @@ // #include #include "aether.h" -//#include "output.cpp" /* --------------------------------------------------------------------- @@ -20,11 +19,11 @@ // Fills output containers and outputs them for common output types // ----------------------------------------------------------------------------- -int output(const Neutrals &neutrals, - const Ions &ions, - const Grid &grid, - Times time, - const Planets &planet) { +bool output(const Neutrals &neutrals, + const Ions &ions, + const Grid &grid, + Times time, + const Planets &planet) { std::string function = "output"; static int iFunction = -1; @@ -32,7 +31,7 @@ int output(const Neutrals &neutrals, static bool IsFirstTime = true; - int iErr = 0; + bool didWork = true; int nOutputs = input.get_n_outputs(); static std::vector AllOutputContainers; @@ -311,7 +310,10 @@ int output(const Neutrals &neutrals, // ------------------------------------------------------------ // write output container - AllOutputContainers[iOutput].write(); + if (!AllOutputContainers[iOutput].write()) { + report.error("Error in writing output container!"); + didWork = false; + } // ------------------------------------------------------------ // Clear variables for next time @@ -321,7 +323,7 @@ int output(const Neutrals &neutrals, } report.exit(function); - return iErr; + return didWork; } /* --------------------------------------------------------------------- @@ -509,19 +511,21 @@ OutputContainer::OutputContainer() { // the user wants, and then output that particular type. // ----------------------------------------------------------------------------- -void OutputContainer::write() { +bool OutputContainer::write() { - int iErr = 0; + bool didWork = true; if (output_type == binary_type) { - iErr = write_container_header(); + didWork = write_container_header(); - if (iErr == 0) - iErr = write_container_binary(); + if (didWork) + didWork = write_container_binary(); } if (output_type == netcdf_type) - iErr = write_container_netcdf(); + didWork = write_container_netcdf(); + + return didWork; } // ----------------------------------------------------------------------------- @@ -558,9 +562,9 @@ void OutputContainer::display() { // formatted file. // ----------------------------------------------------------------------------- -int OutputContainer::write_container_header() { +bool OutputContainer::write_container_header() { - int iErr = 0; + bool didWork = true; int64_t nVars = elements.size(); int64_t nX = elements[0].value.n_rows; int64_t nY = elements[0].value.n_cols; @@ -592,12 +596,11 @@ int OutputContainer::write_container_header() { file << header.dump(4) << "\n"; file.close(); } catch (...) { - std::cout << "Error writing header file : " - << whole_filename << "\n"; - iErr = 1; + report.error("Error writing header file : " + whole_filename); + didWork = false; } - return iErr; + return didWork; } //---------------------------------------------------------------------- @@ -643,9 +646,9 @@ void output_binary_3d(std::ofstream &binary, // dump the contents of the container out into a binary file // ----------------------------------------------------------------------------- -int OutputContainer::write_container_binary() { +bool OutputContainer::write_container_binary() { - int iErr = 0; + bool didWork = true; std::ofstream binary; std::string whole_filename = directory + "/" + filename + ".bin"; @@ -657,14 +660,13 @@ int OutputContainer::write_container_binary() { for (int64_t iVar = 0; iVar < nVars; iVar++) output_binary_3d(binary, elements[iVar].value); - return iErr; + return didWork; } catch (...) { - std::cout << "Error writing header file : " - << whole_filename << "\n"; - iErr = 1; + report.error("Error writing binary file : " + whole_filename); + didWork = false; } - return iErr; + return didWork; } // ----------------------------------------------------------------------------- diff --git a/src/output_netcdf.cpp b/src/output_netcdf.cpp index 5227c8a9..0544ecc1 100644 --- a/src/output_netcdf.cpp +++ b/src/output_netcdf.cpp @@ -158,9 +158,9 @@ int OutputContainer::read_container_netcdf() { // dump the contents of the container out into a binary file // ----------------------------------------------------------------------------- -int OutputContainer::write_container_netcdf() { +bool OutputContainer::write_container_netcdf() { - int iErr = 0; + bool didWork = true; std::string whole_filename = directory + "/" + filename + ".nc"; std::string UNITS = "units"; std::string LONG_NAME = "long_name"; @@ -204,12 +204,11 @@ int OutputContainer::write_container_netcdf() { ncdf_file.close(); } catch (...) { - std::cout << "Error writing netcdf container file : " - << whole_filename << "\n"; - iErr = 1; + report.error("Error writing netcdf container file : " + whole_filename); + didWork = false; } - return iErr; + return didWork; } #else @@ -226,10 +225,10 @@ int OutputContainer::read_container_netcdf() { return iErr; } -int OutputContainer::write_container_netcdf() { - int iErr = 1; - std::cout << "write_container_netcdf is not working!\n"; - return iErr; +bool OutputContainer::write_container_netcdf() { + bool didWork = false; + report.error("write_container_netcdf is not working!"); + return didWork; } #endif From 6cb1f151e15fd18f3ce005aae15b93b750a12e7a Mon Sep 17 00:00:00 2001 From: Aaron Ridley Date: Fri, 10 Nov 2023 08:50:05 -0500 Subject: [PATCH 12/23] FEAT: use get_settings to check for bad keys --- src/read_input_file.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/read_input_file.cpp b/src/read_input_file.cpp index 262400a1..e7b9f51d 100644 --- a/src/read_input_file.cpp +++ b/src/read_input_file.cpp @@ -68,12 +68,12 @@ bool Inputs::read_inputs_json(Times &time) { report.print(1, "Using planet file : " + planet_filename); // Debug Stuff: - report.set_verbose(settings["Debug"]["iVerbose"]); - report.set_DefaultVerbose(settings["Debug"]["iVerbose"]); + report.set_verbose(get_settings("Debug", "iVerbose")); + report.set_DefaultVerbose(get_settings("Debug", "iVerbose")); report.set_doInheritVerbose(settings["Debug"]["doInheritVerbose"]); - report.set_timing_depth(settings["Debug"]["iTimingDepth"]); + report.set_timing_depth(get_settings("Debug", "iTimingDepth")); report.set_timing_percent(settings["Debug"]["TimingPercent"]); - report.set_iProc(settings["Debug"]["iProc"]); + report.set_iProc(get_settings("Debug", "iProc")); for (auto &item : settings["Debug"]["iFunctionVerbose"].items()) report.set_FunctionVerbose(item.key(), item.value()); From 8fe24bbc7c805d5282cdfcc46582bbbc48118c3c Mon Sep 17 00:00:00 2001 From: Aaron Ridley Date: Fri, 10 Nov 2023 08:51:47 -0500 Subject: [PATCH 13/23] BUG: only 0th grid proc output and add cMember --- src/logfile.cpp | 86 ++++++++++++++++--------------------------------- 1 file changed, 27 insertions(+), 59 deletions(-) diff --git a/src/logfile.cpp b/src/logfile.cpp index ff657a6f..3864d1ed 100644 --- a/src/logfile.cpp +++ b/src/logfile.cpp @@ -252,20 +252,22 @@ Logfile::Logfile(Indices &indices) { std::vector sat_names = input.get_satellite_names(); std::vector sat_dts = input.get_satellite_dts(); - // Open the general log file stream - if (doAppend) - logfilestream.open(logfileName, std::ofstream::app); + // Only the 0th processor in each member needs to open the logfile - else { - logfilestream.open(logfileName, std::ofstream::trunc); - logfilestream.precision(4); - } - - // Report error if can not open the log file stream - if (!logfilestream.is_open()) { - // TRY TO EXIT GRACEFULLY HERE. ALL THE FOLLOWING CODE SHOULD - // NOT BE EXECUTED - throw std::string("Can not open log file"); + if (iGrid == 0) { + // Open the general log file stream + if (doAppend) + logfilestream.open(logfileName, std::ofstream::app); + else { + logfilestream.open(logfileName, std::ofstream::trunc); + logfilestream.precision(4); + } + // Report error if can not open the log file stream + if (!logfilestream.is_open()) { + // TRY TO EXIT GRACEFULLY HERE. ALL THE FOLLOWING CODE SHOULD + // NOT BE EXECUTED + throw std::string("Can not open log file"); + } } // The header of time. They are placed at the beginning of all values @@ -288,12 +290,14 @@ Logfile::Logfile(Indices &indices) { for (int i = 0; i < indices.all_indices_array_size(); ++i) header_log += ' ' + indices.get_name(i); - // Write the header to log - logfilestream << header_time << header_log << '\n'; - - // Close the file stream if append - if (doAppend) - logfilestream.close(); + // only the 0th processor writes to the logfile + if (iGrid == 0) { + // Write the header to log + logfilestream << header_time << header_log << '\n'; + // Close the file stream if append + if (doAppend) + logfilestream.close(); + } // Check whether the input settings contain Satellites if (!sat_files.empty()) { @@ -329,7 +333,8 @@ Logfile::Logfile(Indices &indices) { //------------------------------------------------------------- Logfile::~Logfile() { - if (!doAppend) + // only the 0th processor in each ensemble opened the file + if (!doAppend & iGrid == 0) logfilestream.close(); } @@ -410,7 +415,7 @@ bool Logfile::write_logfile(Indices &indices, } // Check if the time gate for general log file is open - if (iProc == 0 && time.check_time_gate(dt)) { + if (iGrid == 0 && time.check_time_gate(dt)) { // Open the file stream if append if (doAppend) { logfilestream.open(logfileName, std::ofstream::app); @@ -419,7 +424,7 @@ bool Logfile::write_logfile(Indices &indices, // Report error if the file stream is not open if (!logfilestream.is_open()) { - std::cout << "Logfile: Can not open output file stream!\n"; + report.error("Logfile: Can not open output file stream!"); report.exit(function); return false; } @@ -456,43 +461,6 @@ bool Logfile::write_logfile(Indices &indices, logfilestream.close(); } - // Report error if the file stream is not open - if (!logfilestream.is_open()) { - std::cout << "Logfile: Can not open output file stream!\n"; - return false; - } - - std::vector values_log; - std::vector min_mean_max; - // Temperature - min_mean_max = get_min_mean_max(neutrals.temperature_scgc); - values_log.insert(values_log.end(), - min_mean_max.begin(), - min_mean_max.end()); - // Temperature at the chosen point - values_log.push_back(neutrals.temperature_scgc(lla[0], - lla[1], - lla[2])); - - // Specified variables - for (auto it : species) { - min_mean_max = get_min_mean_max_density(it, neutrals, ions); - values_log.insert(values_log.end(), - min_mean_max.begin(), - min_mean_max.end()); - } - - // Indices - for (int i = 0; i < indices.all_indices_array_size(); ++i) - values_log.push_back(indices.get_index(current, i)); - - // Output - write_logfile_line(logfilestream, iCurrent, values_log); - - // Close the file stream if append - if (doAppend) - logfilestream.close(); - report.exit(function); return true; } From ff9ba215c5fc9b080bee185292d920f3a6a42fd3 Mon Sep 17 00:00:00 2001 From: Aaron James Ridley Date: Fri, 10 Nov 2023 14:47:41 -0500 Subject: [PATCH 14/23] FEAT: better restart test with ensembles --- tests/restart_ensembles/aether.json.restart1 | 67 ++++++++++++++++++++ tests/restart_ensembles/aether.json.restart2 | 67 ++++++++++++++++++++ tests/restart_ensembles/aether.json.start | 67 ++++++++++++++++++++ tests/restart_ensembles/aether.json.whole | 67 ++++++++++++++++++++ tests/restart_ensembles/run_all.sh | 38 +++++++++++ 5 files changed, 306 insertions(+) create mode 100644 tests/restart_ensembles/aether.json.restart1 create mode 100644 tests/restart_ensembles/aether.json.restart2 create mode 100644 tests/restart_ensembles/aether.json.start create mode 100644 tests/restart_ensembles/aether.json.whole create mode 100755 tests/restart_ensembles/run_all.sh diff --git a/tests/restart_ensembles/aether.json.restart1 b/tests/restart_ensembles/aether.json.restart1 new file mode 100644 index 00000000..ecf11375 --- /dev/null +++ b/tests/restart_ensembles/aether.json.restart1 @@ -0,0 +1,67 @@ + +{ + "Ensembles" : { + "nMembers" : 5}, + + "Perturb": { + "f107" : { "Mean" : 1.0, + "Std" : 0.1, + "Add" : false, + "Constant" : true}, + "f107a" : { "Mean" : 1.0, + "Std" : 0.02, + "Add" : false, + "Constant" : true}}, + + "Restart" : { + "do" : true, + "OutDir" : "UA/restartOut", + "InDir" : "UA/restartIn", + "dt" : 3600.0}, + + "Debug" : { + "iVerbose" : 0, + "iFunctionVerbose" : { + "Grid::create_altitudes": 0}, + "dt" : 10.0, + "check_for_nans" : false + }, + + "EndTime" : [2011, 3, 20, 0, 10, 0], + + "GeoBlockSize" : { + "nLons" : 18, + "nLats" : 18, + "nAlts" : 50}, + + "GeoGrid" : { + "dAlt" : 0.25, + "IsUniformAlt" : false}, + + "InitialConditions" : { + "type" : "msis"}, + + "BoundaryConditions" : { + "type" : "msis"}, + + "Advection" : { + "Neutrals" : { + "Vertical" : "hydro", + "Horizontal" : "default"}, + "Ions" : { + "Along" : "rusanov", + "Across" : "default"} }, + + "OmniwebFiles" : ["UA/inputs/omni_20110319.txt"], + + "ElectrodynamicsFile" : "UA/inputs/b20110320n_omni.bin", + + "Outputs" : { + "type" : ["states"], + "dt" : [900] }, + + "DoCalcBulkIonTemp" : false, + + "PlanetFile" : "UA/inputs/earth.in" + +} diff --git a/tests/restart_ensembles/aether.json.restart2 b/tests/restart_ensembles/aether.json.restart2 new file mode 100644 index 00000000..2d83d738 --- /dev/null +++ b/tests/restart_ensembles/aether.json.restart2 @@ -0,0 +1,67 @@ + +{ + "Ensembles" : { + "nMembers" : 5}, + + "Perturb": { + "f107" : { "Mean" : 0.9, + "Std" : 0.1, + "Add" : false, + "Constant" : true}, + "f107a" : { "Mean" : 1.0, + "Std" : 0.02, + "Add" : false, + "Constant" : true}}, + + "Restart" : { + "do" : true, + "OutDir" : "UA/restartOut", + "InDir" : "UA/restartIn", + "dt" : 3600.0}, + + "Debug" : { + "iVerbose" : 0, + "iFunctionVerbose" : { + "Grid::create_altitudes": 0}, + "dt" : 10.0, + "check_for_nans" : false + }, + + "EndTime" : [2011, 3, 20, 0, 10, 0], + + "GeoBlockSize" : { + "nLons" : 18, + "nLats" : 18, + "nAlts" : 50}, + + "GeoGrid" : { + "dAlt" : 0.25, + "IsUniformAlt" : false}, + + "InitialConditions" : { + "type" : "msis"}, + + "BoundaryConditions" : { + "type" : "msis"}, + + "Advection" : { + "Neutrals" : { + "Vertical" : "hydro", + "Horizontal" : "default"}, + "Ions" : { + "Along" : "rusanov", + "Across" : "default"} }, + + "OmniwebFiles" : ["UA/inputs/omni_20110319.txt"], + + "ElectrodynamicsFile" : "UA/inputs/b20110320n_omni.bin", + + "Outputs" : { + "type" : ["states"], + "dt" : [900] }, + + "DoCalcBulkIonTemp" : false, + + "PlanetFile" : "UA/inputs/earth.in" + +} diff --git a/tests/restart_ensembles/aether.json.start b/tests/restart_ensembles/aether.json.start new file mode 100644 index 00000000..da15c5ee --- /dev/null +++ b/tests/restart_ensembles/aether.json.start @@ -0,0 +1,67 @@ + +{ + "Ensembles" : { + "nMembers" : 5}, + + "Perturb": { + "f107" : { "Mean" : 1.0, + "Std" : 0.1, + "Add" : false, + "Constant" : true}, + "f107a" : { "Mean" : 1.0, + "Std" : 0.02, + "Add" : false, + "Constant" : true}}, + + "Restart" : { + "do" : false, + "OutDir" : "UA/restartOut", + "InDir" : "UA/restartIn", + "dt" : 3600.0}, + + "Debug" : { + "iVerbose" : 0, + "iFunctionVerbose" : { + "Grid::create_altitudes": 0}, + "dt" : 10.0, + "check_for_nans" : false + }, + + "EndTime" : [2011, 3, 20, 0, 5, 0], + + "GeoBlockSize" : { + "nLons" : 18, + "nLats" : 18, + "nAlts" : 50}, + + "GeoGrid" : { + "dAlt" : 0.25, + "IsUniformAlt" : false}, + + "InitialConditions" : { + "type" : "msis"}, + + "BoundaryConditions" : { + "type" : "msis"}, + + "Advection" : { + "Neutrals" : { + "Vertical" : "hydro", + "Horizontal" : "default"}, + "Ions" : { + "Along" : "rusanov", + "Across" : "default"} }, + + "OmniwebFiles" : ["UA/inputs/omni_20110319.txt"], + + "ElectrodynamicsFile" : "UA/inputs/b20110320n_omni.bin", + + "Outputs" : { + "type" : ["states"], + "dt" : [900] }, + + "DoCalcBulkIonTemp" : false, + + "PlanetFile" : "UA/inputs/earth.in" + +} diff --git a/tests/restart_ensembles/aether.json.whole b/tests/restart_ensembles/aether.json.whole new file mode 100644 index 00000000..5805cf9b --- /dev/null +++ b/tests/restart_ensembles/aether.json.whole @@ -0,0 +1,67 @@ + +{ + "Ensembles" : { + "nMembers" : 5}, + + "Perturb": { + "f107" : { "Mean" : 1.0, + "Std" : 0.1, + "Add" : false, + "Constant" : true}, + "f107a" : { "Mean" : 1.0, + "Std" : 0.02, + "Add" : false, + "Constant" : true}}, + + "Restart" : { + "do" : false, + "OutDir" : "UA/restartOut", + "InDir" : "UA/restartIn", + "dt" : 3600.0}, + + "Debug" : { + "iVerbose" : 0, + "iFunctionVerbose" : { + "Grid::create_altitudes": 0}, + "dt" : 10.0, + "check_for_nans" : false + }, + + "EndTime" : [2011, 3, 20, 0, 10, 0], + + "GeoBlockSize" : { + "nLons" : 18, + "nLats" : 18, + "nAlts" : 50}, + + "GeoGrid" : { + "dAlt" : 0.25, + "IsUniformAlt" : false}, + + "InitialConditions" : { + "type" : "msis"}, + + "BoundaryConditions" : { + "type" : "msis"}, + + "Advection" : { + "Neutrals" : { + "Vertical" : "hydro", + "Horizontal" : "default"}, + "Ions" : { + "Along" : "rusanov", + "Across" : "default"} }, + + "OmniwebFiles" : ["UA/inputs/omni_20110319.txt"], + + "ElectrodynamicsFile" : "UA/inputs/b20110320n_omni.bin", + + "Outputs" : { + "type" : ["states"], + "dt" : [900] }, + + "DoCalcBulkIonTemp" : false, + + "PlanetFile" : "UA/inputs/earth.in" + +} diff --git a/tests/restart_ensembles/run_all.sh b/tests/restart_ensembles/run_all.sh new file mode 100755 index 00000000..64de4caa --- /dev/null +++ b/tests/restart_ensembles/run_all.sh @@ -0,0 +1,38 @@ +#!/bin/sh + +rm -rf run.test +cp -R ../../share/run ./run.test +cd run.test + +rm -f aether.json UA/output/* +rm -rf UA/restart* ; mkdir UA/restartOut +rm -rf UA/output_whole +rm -rf UA/output_nonperturbed +rm -rf UA/output_perturbed +rm -rf UA/output_save + +cp ../aether.json.whole aether.json +mpirun -np 20 ./aether +cd UA/output ; ../../../../../srcPython/plot_logfiles.py -vars=10 log_m00*.txt ; cd - +cd UA ; mv output output_whole ; mkdir output ; cd - +rm -f aether.json + +rm -f aether.json UA/output/* UA/restartOut/* +cp ../aether.json.start aether.json +mpirun -np 20 ./aether +cd UA ; mv restartOut restartIn ; mkdir restartOut ; cd - +cd UA ; cp -R output output_save ; cd .. + +rm -f aether.json +cp ../aether.json.restart1 aether.json +mpirun -np 20 ./aether +cd UA/output ; ../../../../../srcPython/plot_logfiles.py -vars=10 log_m00*.txt ; cd - +cd UA ; mv output output_nonperturbed ; cp -R output_save output ; cd - + +rm -f aether.json +cp ../aether.json.restart2 aether.json +mpirun -np 20 ./aether +cd UA/output ; ../../../../../srcPython/plot_logfiles.py -vars=10 log_m00*.txt ; cd - +cd UA ; mv output output_perturbed ; mkdir output ; cd - + + From 7f73374cd311cb59842cd75ef5c9e42b9f81dedb Mon Sep 17 00:00:00 2001 From: Aaron James Ridley Date: Fri, 10 Nov 2023 14:51:17 -0500 Subject: [PATCH 15/23] FEAT: read and plot log files --- srcPython/plot_logfiles.py | 342 +++++++++++++++++++++++++++++++++++++ 1 file changed, 342 insertions(+) create mode 100755 srcPython/plot_logfiles.py diff --git a/srcPython/plot_logfiles.py b/srcPython/plot_logfiles.py new file mode 100755 index 00000000..5483b90a --- /dev/null +++ b/srcPython/plot_logfiles.py @@ -0,0 +1,342 @@ +#!/usr/bin/env python3 + +import matplotlib as mpl +mpl.use('Agg') +import matplotlib.pyplot as plt +import numpy as np +import re +import argparse +import datetime as dt + +# ---------------------------------------------------------------------- +# Function to parse input arguments +# ---------------------------------------------------------------------- + +def get_args_timeline(): + + parser = argparse.ArgumentParser(description = + 'Post process and move model results') + parser.add_argument('-plotfile', + help = 'output file for plot', + default = 'timeline.png') + + parser.add_argument('-vars', + help = 'var(s) to plot (e.g. -vars 0 or -vars 0 1 2', + default = [0], nargs = '+', type = int) + + parser.add_argument('filelist', nargs='+', \ + help = 'list files to use for generating plots') + + args = parser.parse_args() + + return args + +# ---------------------------------------------------------------------- +# Read space seperated file with header containing variable names +# ---------------------------------------------------------------------- + +def read_ssv_file(file): + + data = {"times" : [], + 'vars': [], + 'alt' : 0.0, + "integral" : 'values', + "file": file} + + fpin = open(file, 'r') + lines = fpin.readlines() + fpin.close() + + vars = lines[0].split(' ') + if (vars[6] == 'milli'): + iOff = 7 + else: + iOff = 6 + + data['vars'] = vars[iOff:] + nVars = len(data['vars']) + + iStart = 1 + iEnd = len(lines) + iLine = iStart + allValues = np.zeros((nVars, iEnd - iStart)) + + while (iLine < iEnd): + aline = lines[iLine].split(' ') + + year = int(aline[0]) + month = int(aline[1]) + day = int(aline[2]) + hour = int(aline[3]) + minute = int(aline[4]) + second = int(aline[5]) + if (iOff == 6): + milli = 0 + else: + milli = int(aline[6]) + t = dt.datetime(year, month, day, hour, minute, second, milli) + + data["times"].append(t) + for iVar in range(nVars): + allValues[iVar, iLine - iStart] = float(aline[iOff + iVar]) + iLine += 1 + + data['values'] = allValues + + return data + + +# ---------------------------------------------------------------------- +# Read file +# ---------------------------------------------------------------------- + +def read_timeline_file(file): + + data = {"times" : [], + 'vars': [], + 'alt' : 0.0, + "integral" : 'values', + "file": file} + + fpin = open(file, 'r') + + lines = fpin.readlines() + + iStart = 0 + iEnd = len(lines) + iLine = iStart + + while (iLine < iEnd): + line = lines[iLine] + + m = re.match(r'#VAR',line) + if m: + # skip year, month, day, hour, minute, second: + if (lines[iLine+1].strip().lower() == 'year'): + iLine += 7 + else: + iLine += 1 + while (len(lines[iLine].strip()) > 1): + data["vars"].append(lines[iLine].strip()) + iLine += 1 + data['var'] = data["vars"][0] + + m = re.match(r'#INTEGRAL',line) + if m: + iLine += 1 + data["integral"] = lines[iLine].strip() + + m = re.match(r'#ALTITUDE',line) + if m: + iLine += 1 + data["alt"] = float(lines[iLine].strip()) + + m = re.match(r'#DIRECTORY',line) + if m: + iLine += 1 + data["dir"] = lines[iLine].strip() + + m = re.match(r'#START',line) + if m: + iStart = iLine + 1 + break + + iLine += 1 + + nVars = len(data['vars']) + nTimes = iEnd - iStart + allValues = np.zeros([nVars, nTimes]) + for i in range(iStart, iEnd): + aline = lines[i].split() + year = int(aline[0]) + month = int(aline[1]) + day = int(aline[2]) + hour = int(aline[3]) + minute = int(aline[4]) + second = int(aline[5]) + t = dt.datetime(year, month, day, hour, minute, second) + + data["times"].append(t) + for iVar in range(nVars): + allValues[iVar, i - iStart] = float(aline[6 + iVar]) + + data['values'] = allValues + + return data + +# ---------------------------------------------------------------------- +# match filename with colors and linestyles +# ---------------------------------------------------------------------- + +def assign_var_to_color(var): + + color = 'black' + line = 'solid' + label = file + + m = re.match('.*GOCE.*', var) + if m: + color = 'black' + line = 'solid' + label = 'GOCE' + + m = re.match('.*GITM.*', var) + if m: + color = 'blue' + line = 'solid' + label = 'GITM' + + m = re.match('.*Smooth.*', var) + if m: + line = 'dashed' + label = label + ' (smoothed)' + + return color, line, label + +# ---------------------------------------------------------------------- +# match filename with colors and linestyles +# ---------------------------------------------------------------------- + +def assign_file_to_color(file): + + color = 'black' + line = 'solid' + label = file + + m = re.match('gdcA.*', file) + if m: + color = 'black' + + m = re.match('gdcB.*', file) + if m: + color = 'blue' + + m = re.match('gdcC.*', file) + if m: + color = 'cyan' + + m = re.match('gdcD.*', file) + if m: + color = 'purple' + + m = re.match('gdcE.*', file) + if m: + color = 'orange' + + m = re.match('gdcF.*', file) + if m: + color = 'red' + + m = re.match(r'base', file) + if m: + color = 'black' + line = 'solid' + label = 'Idealized Model' + + m = re.match(r'fre', file) + if m: + color = 'grey' + line = 'solid' + label = 'FRE Model' + + m = re.match(r'ovation', file) + if m: + color = 'blue' + line = 'dashed' + label = 'Ovation' + + m = re.match(r'avee_050', file) + if m: + color = 'darkblue' + line = 'dashed' + m = re.match(r'avee_075', file) + if m: + color = 'blue' + line = 'dashed' + + m = re.match(r'avee_150', file) + if m: + color = 'firebrick' + line = 'dashdot' + + m = re.match(r'avee_200', file) + if m: + color = 'red' + line = 'dashdot' + + m = re.match(r'fta', file) + if m: + color = 'forestgreen' + line = 'dashed' + label = 'FTA Model' + + m = re.match(r'ae', file) + if m: + color = 'plum' + line = 'dashed' + label = 'FRE, AE Power' + + m = re.match(r'ions', file) + if m: + color = 'plum' + line = 'dashdot' + label = 'Ions' + + m = re.match(r'.*no', file) + if m: + color = 'cyan' + label = 'NO Cooling' + m = re.match(r'.*ch', file) + if m: + color = 'red' + label = 'Chemical Heating' + m = re.match(r'.*ht', file) + if m: + color = 'blue' + label = 'Heat Transfer' + m = re.match(r'.*jh', file) + if m: + label = 'Joule Heating' + + return color, line, label + +# Needed to run main script as the default executable from the command line +if __name__ == '__main__': + + args = get_args_timeline() + + iVar = args.vars[0] + + nFiles = len(args.filelist) + nVars = len(args.vars) + + fig = plt.figure(figsize = (10,10)) + ax = fig.add_subplot(111) + + for file in args.filelist: + print("Reading file : ",file) + #data = read_timeline_file(file) + data = read_ssv_file(file) + for iVar in args.vars: + label = data["vars"][iVar] + color = 'black' + line = '--' + if (nFiles > 1): + color,line,label = assign_file_to_color(file) + if (nVars > 1): + color,line,label = assign_var_to_color(data['vars'][iVar]) + + ax.plot(data["times"], data["values"][iVar], label = label, + color = color, linestyle = line, linewidth = 2.0) + + iVar = args.vars[0] + ytitle = data["integral"] + " of " + data["vars"][iVar] + if (data["alt"] > 0): + ytitle += " at %5.1f" % data["alt"] + " km Alt" + ax.set_ylabel(ytitle) + + ax.legend() + plotfile = args.plotfile + print('writing : ',plotfile) + fig.savefig(plotfile) + plt.close() From 09024c715f1cca92652206e19f494e9ad03885bb Mon Sep 17 00:00:00 2001 From: Aaron James Ridley Date: Fri, 10 Nov 2023 14:53:44 -0500 Subject: [PATCH 16/23] BUG: didWork needs to be initialized --- src/calc_euv.cpp | 73 ++++++++++++++++++++++++------------------------ 1 file changed, 36 insertions(+), 37 deletions(-) diff --git a/src/calc_euv.cpp b/src/calc_euv.cpp index ed51773e..c4e400ed 100644 --- a/src/calc_euv.cpp +++ b/src/calc_euv.cpp @@ -21,45 +21,44 @@ bool calc_euv(Planets planet, Ions &ions, Indices indices) { - bool didWork; - - if (time.check_time_gate(input.get_dt_euv())) { - std::string function = "Euv::calc_euv"; - static int iFunction = -1; - report.enter(function, iFunction); - - if (input.get_is_student()) - report.print(-1, "(2) What function is this " + - input.get_student_name() + "?"); - - // Chapman integrals for EUV energy deposition: - // Need chapman integrals for aurora too! - neutrals.calc_chapman(grid); - - if (euv.doUse) { - // set didWork to false in order to catch bad euv models: - didWork = false; - std::string euvModel = mklower(input.get_euv_model()); - if (euvModel == "euvac") - didWork = euv.euvac(time, indices); - else if (euvModel == "neuvac") - didWork = euv.neuvac(time, indices); - else if (euvModel == "hfg") - didWork = euv.solomon_hfg(time, indices); - - if (didWork) - euv.scale_from_1au(planet, time); - - calc_ionization_heating(euv, neutrals, ions); - } else - neutrals.heating_euv_scgc.zeros(); - - report.exit(function); - } + bool didWork = true; + + if (!time.check_time_gate(input.get_dt_euv())) + return true; + + std::string function = "Euv::calc_euv"; + static int iFunction = -1; + report.enter(function, iFunction); + + if (input.get_is_student()) + report.print(-1, "(2) What function is this " + + input.get_student_name() + "?"); + + // Chapman integrals for EUV energy deposition: + // Need chapman integrals for aurora too! + neutrals.calc_chapman(grid); + + if (euv.doUse) { + // set didWork to false in order to catch bad euv models: + didWork = false; + std::string euvModel = mklower(input.get_euv_model()); + if (euvModel == "euvac") + didWork = euv.euvac(time, indices); + else if (euvModel == "neuvac") + didWork = euv.neuvac(time, indices); + else if (euvModel == "hfg") + didWork = euv.solomon_hfg(time, indices); + + if (didWork) + euv.scale_from_1au(planet, time); + + calc_ionization_heating(euv, neutrals, ions); + } else + neutrals.heating_euv_scgc.zeros(); if (!didWork) - report.error("Error in calc_euv! Check euv models."); - + report.error("Error in calc_euv! Check euv models."); + report.exit(function); return didWork; } From dcd15ffbab64e6129ffa8f6ddd9712ae5e99b6fd Mon Sep 17 00:00:00 2001 From: Aaron James Ridley Date: Fri, 10 Nov 2023 14:54:39 -0500 Subject: [PATCH 17/23] BUG: catch bad dts --- src/calc_neutral_derived.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/calc_neutral_derived.cpp b/src/calc_neutral_derived.cpp index 926886e1..cdd60a7d 100644 --- a/src/calc_neutral_derived.cpp +++ b/src/calc_neutral_derived.cpp @@ -359,6 +359,8 @@ precision_t Neutrals::calc_dt(Grid grid) { dt = dta.min(); + if (dt < 0.0) dt = 2.0; + if (report.test_verbose(3)) std::cout << "dt for neutrals : " << dt << "\n"; From 63e682b6c25f12e360d22d46228310a675e9c8aa Mon Sep 17 00:00:00 2001 From: Aaron James Ridley Date: Fri, 10 Nov 2023 14:55:18 -0500 Subject: [PATCH 18/23] BUG: return 0, since this is catch on other side --- src/inputs.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/inputs.cpp b/src/inputs.cpp index ce84b4f7..7385764d 100644 --- a/src/inputs.cpp +++ b/src/inputs.cpp @@ -514,9 +514,8 @@ int Inputs::get_original_seed() { if (settings.find("Seed") == settings.end()) { IsOk = false; std::cout << "Error in getting seed!\n"; - return dummy_int; + return 0; } - return settings.at("Seed"); } From b162cfc5d079fe79f2b1fa37f6bf06502ced1099 Mon Sep 17 00:00:00 2001 From: Aaron James Ridley Date: Fri, 10 Nov 2023 14:55:47 -0500 Subject: [PATCH 19/23] BUG: don't write header on restart --- src/logfile.cpp | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/logfile.cpp b/src/logfile.cpp index 3864d1ed..1d43d326 100644 --- a/src/logfile.cpp +++ b/src/logfile.cpp @@ -252,6 +252,8 @@ Logfile::Logfile(Indices &indices) { std::vector sat_names = input.get_satellite_names(); std::vector sat_dts = input.get_satellite_dts(); + bool isRestart = input.get_do_restart(); + // Only the 0th processor in each member needs to open the logfile if (iGrid == 0) { @@ -263,7 +265,7 @@ Logfile::Logfile(Indices &indices) { logfilestream.precision(4); } // Report error if can not open the log file stream - if (!logfilestream.is_open()) { + if (!logfilestream.is_open() & !isRestart) { // TRY TO EXIT GRACEFULLY HERE. ALL THE FOLLOWING CODE SHOULD // NOT BE EXECUTED throw std::string("Can not open log file"); @@ -291,9 +293,10 @@ Logfile::Logfile(Indices &indices) { header_log += ' ' + indices.get_name(i); // only the 0th processor writes to the logfile - if (iGrid == 0) { + if (logfilestream.is_open()) { // Write the header to log - logfilestream << header_time << header_log << '\n'; + if (!isRestart) + logfilestream << header_time << header_log << '\n'; // Close the file stream if append if (doAppend) logfilestream.close(); From dd21954496fa2d74e222dda21c70b493a6311c48 Mon Sep 17 00:00:00 2001 From: Aaron James Ridley Date: Fri, 10 Nov 2023 14:57:02 -0500 Subject: [PATCH 20/23] BUG: on restart append to the log file --- src/read_input_file.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/read_input_file.cpp b/src/read_input_file.cpp index e7b9f51d..6acf4990 100644 --- a/src/read_input_file.cpp +++ b/src/read_input_file.cpp @@ -44,14 +44,15 @@ bool Inputs::read_inputs_json(Times &time) { // - This is BEFORE the user inputs are merged!!! if (user_inputs.contains("Restart")) { - cout << "Contains restart" << endl; - if (user_inputs["Restart"].contains("do")) { if (user_inputs["Restart"]["do"]) { std::string restart_file = settings["Restart"]["InDir"]; restart_file = restart_file + "/settings.json"; json restart_inputs; restart_inputs = read_json(restart_file); + // This forces the logfile to append. User can override + // if they really want: + restart_inputs["Logfile"]["append"] = true; settings.merge_patch(restart_inputs); } } From 772726cc6bbe3dbe043f7fee044476f854619760 Mon Sep 17 00:00:00 2001 From: Aaron James Ridley Date: Fri, 10 Nov 2023 15:01:26 -0500 Subject: [PATCH 21/23] BUG: new exchange works with cubesphere --- src/advance.cpp | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/advance.cpp b/src/advance.cpp index 355dd294..5d59c03f 100644 --- a/src/advance.cpp +++ b/src/advance.cpp @@ -107,8 +107,11 @@ bool advance(Planets &planet, ions.calc_ion_temperature(neutrals, gGrid, time); ions.calc_electron_temperature(neutrals, gGrid); - neutrals.exchange(gGrid); - + if (input.get_is_cubesphere()) + neutrals.exchange(gGrid); + else + neutrals.exchange_old(gGrid); + time.increment_time(); if (time.check_time_gate(input.get_dt_write_restarts())) { From df231a02e729721e356429a3740fb09d2d2d40b6 Mon Sep 17 00:00:00 2001 From: Aaron James Ridley Date: Fri, 10 Nov 2023 15:02:19 -0500 Subject: [PATCH 22/23] BUG: crashes if trying to write std --- srcPython/postAether.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/srcPython/postAether.py b/srcPython/postAether.py index deeeb228..0282138e 100755 --- a/srcPython/postAether.py +++ b/srcPython/postAether.py @@ -984,11 +984,11 @@ def write_and_plot_data(dataToWrite, write_and_plot_data(ensembleData, fileInfo['ensembleFile'], '_mean', iVar, iAlt, output_netcdf) - stdData = calc_std_of_ensembles(filesInfo, - ensembleIndexList, - ensembleData) - write_and_plot_data(stdData, fileInfo['ensembleFile'], - '_std', iVar, iAlt, output_netcdf) + #stdData = calc_std_of_ensembles(filesInfo, + # ensembleIndexList, + # ensembleData) + #write_and_plot_data(stdData, fileInfo['ensembleFile'], + # '_std', iVar, iAlt, output_netcdf) if (args.rm): print(' --> Removing files ...') From 196e9382882446cad518d804f384e813eb0c9418 Mon Sep 17 00:00:00 2001 From: Aaron Ridley Date: Tue, 14 Nov 2023 10:51:51 -0500 Subject: [PATCH 23/23] STY: astyle --- src/advance.cpp | 44 +- src/calc_euv.cpp | 22 +- src/calc_neutral_derived.cpp | 17 +- src/electrodynamics.cpp | 12 +- src/euv.cpp | 8 +- src/exchange_messages_v2.cpp | 869 +++++++++++++++-------------- src/init_geo_grid.cpp | 8 +- src/inputs.cpp | 19 +- src/logfile.cpp | 4 +- src/neutrals.cpp | 1 + src/neutrals_bcs.cpp | 6 +- src/neutrals_ics.cpp | 92 +-- src/neutrals_momentum_friction.cpp | 2 +- src/output.cpp | 12 +- src/read_input_file.cpp | 8 +- src/solver_gradients.cpp | 37 +- src/tools.cpp | 4 +- src/transform.cpp | 2 + 18 files changed, 611 insertions(+), 556 deletions(-) diff --git a/src/advance.cpp b/src/advance.cpp index 5d59c03f..ff78c827 100644 --- a/src/advance.cpp +++ b/src/advance.cpp @@ -11,18 +11,18 @@ // ----------------------------------------------------------------------------- bool advance(Planets &planet, - Grid &gGrid, - Times &time, - Euv &euv, - Neutrals &neutrals, - Ions &ions, - Chemistry &chemistry, - Electrodynamics &electrodynamics, - Indices &indices, - Logfile &logfile) { + Grid &gGrid, + Times &time, + Euv &euv, + Neutrals &neutrals, + Ions &ions, + Chemistry &chemistry, + Electrodynamics &electrodynamics, + Indices &indices, + Logfile &logfile) { bool didWork = true; - + std::string function = "advance"; static int iFunction = -1; report.enter(function, iFunction); @@ -66,19 +66,19 @@ bool advance(Planets &planet, if (didWork) didWork = calc_euv(planet, - gGrid, - time, - euv, - neutrals, - ions, - indices); + gGrid, + time, + euv, + neutrals, + ions, + indices); if (didWork) didWork = electrodynamics.update(planet, - gGrid, - time, - ions); - + gGrid, + time, + ions); + if (didWork) { calc_ion_neutral_coll_freq(neutrals, ions); ions.calc_ion_drift(neutrals, gGrid, time.get_dt()); @@ -111,7 +111,7 @@ bool advance(Planets &planet, neutrals.exchange(gGrid); else neutrals.exchange_old(gGrid); - + time.increment_time(); if (time.check_time_gate(input.get_dt_write_restarts())) { @@ -130,7 +130,7 @@ bool advance(Planets &planet, if (!didWork) report.error("Error in Advance!"); - + report.exit(function); return didWork; } diff --git a/src/calc_euv.cpp b/src/calc_euv.cpp index c4e400ed..09be872e 100644 --- a/src/calc_euv.cpp +++ b/src/calc_euv.cpp @@ -14,25 +14,25 @@ // ----------------------------------------------------------------------------- bool calc_euv(Planets planet, - Grid grid, - Times time, - Euv &euv, - Neutrals &neutrals, - Ions &ions, - Indices indices) { + Grid grid, + Times time, + Euv &euv, + Neutrals &neutrals, + Ions &ions, + Indices indices) { bool didWork = true; if (!time.check_time_gate(input.get_dt_euv())) return true; - + std::string function = "Euv::calc_euv"; static int iFunction = -1; report.enter(function, iFunction); if (input.get_is_student()) report.print(-1, "(2) What function is this " + - input.get_student_name() + "?"); + input.get_student_name() + "?"); // Chapman integrals for EUV energy deposition: // Need chapman integrals for aurora too! @@ -42,6 +42,7 @@ bool calc_euv(Planets planet, // set didWork to false in order to catch bad euv models: didWork = false; std::string euvModel = mklower(input.get_euv_model()); + if (euvModel == "euvac") didWork = euv.euvac(time, indices); else if (euvModel == "neuvac") @@ -51,13 +52,14 @@ bool calc_euv(Planets planet, if (didWork) euv.scale_from_1au(planet, time); - + calc_ionization_heating(euv, neutrals, ions); } else neutrals.heating_euv_scgc.zeros(); if (!didWork) - report.error("Error in calc_euv! Check euv models."); + report.error("Error in calc_euv! Check euv models."); + report.exit(function); return didWork; } diff --git a/src/calc_neutral_derived.cpp b/src/calc_neutral_derived.cpp index b8d28423..9e1c64a4 100644 --- a/src/calc_neutral_derived.cpp +++ b/src/calc_neutral_derived.cpp @@ -349,7 +349,7 @@ precision_t Neutrals::calc_dt(Grid grid) { int64_t nAlts = grid.get_nAlts(); int64_t nXs = grid.get_nLons(); int64_t nYs = grid.get_nLats(); - + // dtx dty for reference coordinate system arma_cube dtx(size(cMax_vcgc[0])); arma_cube dty(size(cMax_vcgc[0])); @@ -360,15 +360,16 @@ precision_t Neutrals::calc_dt(Grid grid) { // Loop through altitudes for (int iAlt = 0; iAlt < nAlts; iAlt++) { // Conver cMax to contravariant velocity first - arma_mat u1 = cMax_vcgc[0].slice(iAlt) % grid.A11_inv_scgc.slice(iAlt) - + cMax_vcgc[1].slice(iAlt) % grid.A12_inv_scgc.slice(iAlt); - arma_mat u2 = cMax_vcgc[0].slice(iAlt) % grid.A21_inv_scgc.slice(iAlt) - + cMax_vcgc[1].slice(iAlt) % grid.A22_inv_scgc.slice(iAlt); + arma_mat u1 = cMax_vcgc[0].slice(iAlt) % grid.A11_inv_scgc.slice(iAlt) + + cMax_vcgc[1].slice(iAlt) % grid.A12_inv_scgc.slice(iAlt); + arma_mat u2 = cMax_vcgc[0].slice(iAlt) % grid.A21_inv_scgc.slice(iAlt) + + cMax_vcgc[1].slice(iAlt) % grid.A22_inv_scgc.slice(iAlt); // Extract a scalar dx and divide by contravariant velocity - dtx.slice(iAlt) = grid.drefx(iAlt)*dummy_1 / u1; - dty.slice(iAlt) = grid.drefy(iAlt)*dummy_1 / u2; + dtx.slice(iAlt) = grid.drefx(iAlt) * dummy_1 / u1; + dty.slice(iAlt) = grid.drefy(iAlt) * dummy_1 / u2; } + dta(0) = dtx.min(); dta(1) = dty.min(); } else { @@ -378,7 +379,7 @@ precision_t Neutrals::calc_dt(Grid grid) { arma_cube dty = grid.dlat_center_dist_scgc / cMax_vcgc[1]; dta(1) = dty.min(); } - + if (input.get_nAltsGeo() > 1) { arma_cube dtz = grid.dalt_center_scgc / cMax_vcgc[2]; dta(2) = dtz.min(); diff --git a/src/electrodynamics.cpp b/src/electrodynamics.cpp index 14d757fa..476a34f0 100644 --- a/src/electrodynamics.cpp +++ b/src/electrodynamics.cpp @@ -34,9 +34,9 @@ Electrodynamics::Electrodynamics(Times time) { // ----------------------------------------------------------------------------- bool Electrodynamics::update(Planets planet, - Grid gGrid, - Times time, - Ions &ions) { + Grid gGrid, + Times time, + Ions &ions) { std::string function = "Electrodynamics::update"; static int iFunction = -1; @@ -183,9 +183,9 @@ arma_mat Electrodynamics::get_values(arma_mat matToInterpolateOn, // ----------------------------------------------------------------------------- std::tuple Electrodynamics::get_electrodynamics(arma_cube magLat, - arma_cube magLocalTime) { + arma_mat, + arma_mat> Electrodynamics::get_electrodynamics(arma_cube magLat, +arma_cube magLocalTime) { arma_cube pot; arma_mat eflux; arma_mat avee; diff --git a/src/euv.cpp b/src/euv.cpp index 0084eaea..c4a4069a 100644 --- a/src/euv.cpp +++ b/src/euv.cpp @@ -313,7 +313,7 @@ bool Euv::pair_euv(Neutrals &neutrals, // -------------------------------------------------------------------------- void Euv::scale_from_1au(Planets planet, - Times time) { + Times time) { precision_t d = planet.get_star_to_planet_dist(time); precision_t scale = 1.0 / (d * d); @@ -331,7 +331,7 @@ void Euv::scale_from_1au(Planets planet, // -------------------------------------------------------------------------- bool Euv::euvac(Times time, - Indices indices) { + Indices indices) { bool didWork = true; precision_t slope; @@ -375,7 +375,7 @@ bool Euv::euvac(Times time, // -------------------------------------------------------------------------- bool Euv::neuvac(Times time, - Indices indices) { + Indices indices) { bool didWork = true; precision_t slope; @@ -419,7 +419,7 @@ bool Euv::neuvac(Times time, // -------------------------------------------------------------------------- bool Euv::solomon_hfg(Times time, - Indices indices) { + Indices indices) { std::string function = "Euv::solomon_hfg"; static int iFunction = -1; diff --git a/src/exchange_messages_v2.cpp b/src/exchange_messages_v2.cpp index d64e9c73..220b4bf5 100644 --- a/src/exchange_messages_v2.cpp +++ b/src/exchange_messages_v2.cpp @@ -20,35 +20,35 @@ int pos_iProc_surface(precision_t oLR, precision_t LR, precision_t DU, int numProcs) { - // The iProc of the lower left corner - int retProc = 0; - - // Repeat until there is only one processor in the described area - while (numProcs != 1) { - // Divide into 4 parts - dLR /= 2; - dDU /= 2; - - if (LR >= oLR + dLR) { - // The point is on the right, and the iProc of the lower left - // corner where the point is located is [1/4] or [3/4] - // * remainProc, depending on DU - retProc += numProcs / 4; - // Update the origin - oLR += dLR; - } - - if (DU >= oDU + dDU) { - // The point is above - retProc += numProcs / 2; - // Update the origin - oDU += dDU; - } + // The iProc of the lower left corner + int retProc = 0; + + // Repeat until there is only one processor in the described area + while (numProcs != 1) { + // Divide into 4 parts + dLR /= 2; + dDU /= 2; + + if (LR >= oLR + dLR) { + // The point is on the right, and the iProc of the lower left + // corner where the point is located is [1/4] or [3/4] + // * remainProc, depending on DU + retProc += numProcs / 4; + // Update the origin + oLR += dLR; + } - numProcs /= 4; + if (DU >= oDU + dDU) { + // The point is above + retProc += numProcs / 2; + // Update the origin + oDU += dDU; } - return retProc; + numProcs /= 4; + } + + return retProc; } // -------------------------------------------------------------------------- @@ -56,13 +56,13 @@ int pos_iProc_surface(precision_t oLR, // -------------------------------------------------------------------------- int pos_iProc_sphere(precision_t lon_in, precision_t lat_in) { - return pos_iProc_surface(Sphere::ORIGINS(0) * cPI, - Sphere::ORIGINS(1) * cPI, - Sphere::RIGHTS(0) * cPI, - Sphere::UPS(1) * cPI, - lon_in, - lat_in, - nGrids); + return pos_iProc_surface(Sphere::ORIGINS(0) * cPI, + Sphere::ORIGINS(1) * cPI, + Sphere::RIGHTS(0) * cPI, + Sphere::UPS(1) * cPI, + lon_in, + lat_in, + nGrids); } // -------------------------------------------------------------------------- @@ -70,11 +70,11 @@ int pos_iProc_sphere(precision_t lon_in, precision_t lat_in) { // -------------------------------------------------------------------------- bool inRange(precision_t a, precision_t b, precision_t c) { - if (b > 0) { - return c >= a && c < a + b; - } else { - return c <= a && c > a + b; - } + if (b > 0) + return c >= a && c < a + b; + + else + return c <= a && c > a + b; } // -------------------------------------------------------------------------- @@ -82,60 +82,61 @@ bool inRange(precision_t a, precision_t b, precision_t c) { // -------------------------------------------------------------------------- int pos_iProc_cube(precision_t lon_in, precision_t lat_in) { - // Transfer polar coordinate to cartesian coordinate - arma_vec point = sphere_to_cube(lon_in, lat_in); - - // Store the index of surface, LR and DU - // INITIAL VALUES ARE HARD-CODED BASED ON SURFACE 5(6) TO - // SOLVE BOUNDARY PROBLEMS - int iFace = 0, iLR = 1, iDU = 0; - bool on_surface = false; - // Find the surface where the point is located - for (;iFace < 6 && !on_surface; ++iFace) { - // Assume that the point is on current surface - on_surface = true; - int sLR = -1, sDU = -1; - - // Check three dimensions - for (int j = 0; j < 3 && on_surface; ++j) { - if (CubeSphere::RIGHTS(iFace, j) != 0) { - sLR = j; - on_surface = inRange(CubeSphere::ORIGINS(iFace, j), - CubeSphere::RIGHTS(iFace, j), - point(j)); - } else if (CubeSphere::UPS(iFace, j) != 0) { - sDU = j; - on_surface = inRange(CubeSphere::ORIGINS(iFace, j), - CubeSphere::UPS(iFace, j), - point(j)); - } else { - on_surface = (CubeSphere::ORIGINS(iFace, j) == point(j)); - } - } + // Transfer polar coordinate to cartesian coordinate + arma_vec point = sphere_to_cube(lon_in, lat_in); + + // Store the index of surface, LR and DU + // INITIAL VALUES ARE HARD-CODED BASED ON SURFACE 5(6) TO + // SOLVE BOUNDARY PROBLEMS + int iFace = 0, iLR = 1, iDU = 0; + bool on_surface = false; + + // Find the surface where the point is located + for (; iFace < 6 && !on_surface; ++iFace) { + // Assume that the point is on current surface + on_surface = true; + int sLR = -1, sDU = -1; + + // Check three dimensions + for (int j = 0; j < 3 && on_surface; ++j) { + if (CubeSphere::RIGHTS(iFace, j) != 0) { + sLR = j; + on_surface = inRange(CubeSphere::ORIGINS(iFace, j), + CubeSphere::RIGHTS(iFace, j), + point(j)); + } else if (CubeSphere::UPS(iFace, j) != 0) { + sDU = j; + on_surface = inRange(CubeSphere::ORIGINS(iFace, j), + CubeSphere::UPS(iFace, j), + point(j)); + } else + on_surface = (CubeSphere::ORIGINS(iFace, j) == point(j)); + } - // Assign iLR and iDU if the point is indeed on current surface - if (on_surface) { - iLR = sLR; - iDU = sDU; - } + // Assign iLR and iDU if the point is indeed on current surface + if (on_surface) { + iLR = sLR; + iDU = sDU; } - // Minus 1 due to for loop mechanism - --iFace; - - // Change the sign if delta is negative - precision_t coefLR, coefDU; - coefLR = CubeSphere::RIGHTS(iFace, iLR) > 0 ? 1.0 : -1.0; - coefDU = CubeSphere::UPS(iFace, iDU) > 0 ? 1.0 : -1.0; - // Compute iProc in the surface - int ip = pos_iProc_surface(CubeSphere::ORIGINS(iFace, iLR) * coefLR, - CubeSphere::ORIGINS(iFace, iDU) * coefDU, - CubeSphere::RIGHTS(iFace, iLR) * coefLR, - CubeSphere::UPS(iFace, iDU) * coefDU, - point(iLR) * coefLR, - point(iDU) * coefDU, - nGrids / 6); - // Consider surface number - return ip + iFace * nGrids / 6; + } + + // Minus 1 due to for loop mechanism + --iFace; + + // Change the sign if delta is negative + precision_t coefLR, coefDU; + coefLR = CubeSphere::RIGHTS(iFace, iLR) > 0 ? 1.0 : -1.0; + coefDU = CubeSphere::UPS(iFace, iDU) > 0 ? 1.0 : -1.0; + // Compute iProc in the surface + int ip = pos_iProc_surface(CubeSphere::ORIGINS(iFace, iLR) * coefLR, + CubeSphere::ORIGINS(iFace, iDU) * coefDU, + CubeSphere::RIGHTS(iFace, iLR) * coefLR, + CubeSphere::UPS(iFace, iDU) * coefDU, + point(iLR) * coefLR, + point(iDU) * coefDU, + nGrids / 6); + // Consider surface number + return ip + iFace * nGrids / 6; } /******** Functions to find iproc ends ********/ @@ -149,16 +150,16 @@ int pos_iProc_cube(precision_t lon_in, precision_t lat_in) { template struct mpi_msg_t { - // The message to send - std::vector buf; - // The processor to send - int rank; - // Default constructor - mpi_msg_t() - : buf(), rank() {} - // Move constructor - mpi_msg_t(std::vector &&_buf, const int _rank) - : buf(std::move(_buf)), rank(_rank) {} + // The message to send + std::vector buf; + // The processor to send + int rank; + // Default constructor + mpi_msg_t() + : buf(), rank() {} + // Move constructor + mpi_msg_t(std::vector &&_buf, const int _rank) + : buf(std::move(_buf)), rank(_rank) {} }; // -------------------------------------------------------------------------- @@ -168,15 +169,17 @@ struct mpi_msg_t { template std::ostream& operator<<(std::ostream &out, const std::vector> &msg) { - for (auto& it : msg) { - out << "\tTo " << it.rank << ":\n"; - for (auto& it2 : it.buf) { - out << "\t\t" << it2 << '\n'; - } - out << '\n'; - } + for (auto& it : msg) { + out << "\tTo " << it.rank << ":\n"; + + for (auto& it2 : it.buf) + out << "\t\t" << it2 << '\n'; + out << '\n'; - return out; + } + + out << '\n'; + return out; } /** @@ -194,146 +197,158 @@ template bool CUSTOM_MPI_SENDRECV(const std::vector> &msg_in, std::vector> &msg_out, const MPI_Comm &comm) { - // Promise previous send and receive have all ended - MPI_Barrier(comm); - - // Clear the output variables - msg_out.clear(); - - // Get the rank the size of communicator - // Processor with rank 0 checks whether all ends - const int num = msg_in.size(); - int comm_rank, comm_size; - MPI_Comm_rank(comm, &comm_rank); - MPI_Comm_size(comm, &comm_size); - - // Define constants - const int TAG_RESERVED = 3; - const char ALL_FINISH_TAG = 0; - const char PROC_FINISH_TAG = 1; - const char ACK_RECV_TAG = 2; - - // Only used to check correctness - std::vector request; - - // Send - for (int i = 0; i < num; ++i) { - request.emplace_back(); - MPI_Isend(msg_in[i].buf.data(), - msg_in[i].buf.size() * sizeof(T), - MPI_BYTE, - msg_in[i].rank, - i + TAG_RESERVED, - comm, - &request.back()); + // Promise previous send and receive have all ended + MPI_Barrier(comm); + + // Clear the output variables + msg_out.clear(); + + // Get the rank the size of communicator + // Processor with rank 0 checks whether all ends + const int num = msg_in.size(); + int comm_rank, comm_size; + MPI_Comm_rank(comm, &comm_rank); + MPI_Comm_size(comm, &comm_size); + + // Define constants + const int TAG_RESERVED = 3; + const char ALL_FINISH_TAG = 0; + const char PROC_FINISH_TAG = 1; + const char ACK_RECV_TAG = 2; + + // Only used to check correctness + std::vector request; + + // Send + for (int i = 0; i < num; ++i) { + request.emplace_back(); + MPI_Isend(msg_in[i].buf.data(), + msg_in[i].buf.size() * sizeof(T), + MPI_BYTE, + msg_in[i].rank, + i + TAG_RESERVED, + comm, + &request.back()); + } + + // Receive + // Record how many messages are sent successfully and + // how many processors finish sending messages + // finished_proc is only used by iProc 0 + int success_sent = 0, finished_proc = 0; + bool done = false; + + // If nothing needs to send, then the processor already finish sending + if (!num) { + if (!comm_rank) { + // Increase finished_proc if iProc = 0 + ++finished_proc; + } else { + // Send signal to iProc 0 otherwise + request.emplace_back(); + MPI_Isend(NULL, 0, MPI_BYTE, 0, PROC_FINISH_TAG, + comm, &request.back()); } - - // Receive - // Record how many messages are sent successfully and - // how many processors finish sending messages - // finished_proc is only used by iProc 0 - int success_sent = 0, finished_proc = 0; - bool done = false; - // If nothing needs to send, then the processor already finish sending - if (!num) { - if (!comm_rank) { - // Increase finished_proc if iProc = 0 + } + + while (!done) { + MPI_Status status; + int count; + // Probe the size of next message + MPI_Probe(MPI_ANY_SOURCE, MPI_ANY_TAG, comm, &status); + MPI_Get_count(&status, MPI_BYTE, &count); + + if (count) { + // Receiving actual messages + mpi_msg_t msg; + // Set rank and count, create space to receive message + msg.buf.resize(count / sizeof(T)); + msg.rank = status.MPI_SOURCE; + MPI_Recv(msg.buf.data(), + count, + MPI_BYTE, + status.MPI_SOURCE, + status.MPI_TAG, + comm, + MPI_STATUS_IGNORE); + msg_out.push_back(msg); + // Send ack message to the origin of this message + request.emplace_back(); + MPI_Isend(NULL, 0, MPI_BYTE, status.MPI_SOURCE, + ACK_RECV_TAG, comm, &request.back()); + } else { + // The message is a control signal with size = 0 + MPI_Recv(NULL, 0, MPI_BYTE, status.MPI_SOURCE, + status.MPI_TAG, comm, &status); + + switch (status.MPI_TAG) { + case ACK_RECV_TAG: { + // Increase the success_sent + ++success_sent; + + if (success_sent == num) { + // The processor finish sending + if (!comm_rank) { + // Increase finished_proc for rank 0 ++finished_proc; - } else { - // Send signal to iProc 0 otherwise + } else { + // Send proc finish to rank 0 for other processors request.emplace_back(); MPI_Isend(NULL, 0, MPI_BYTE, 0, PROC_FINISH_TAG, comm, &request.back()); + } } - } - while (!done) { - MPI_Status status; - int count; - // Probe the size of next message - MPI_Probe(MPI_ANY_SOURCE, MPI_ANY_TAG, comm, &status); - MPI_Get_count(&status, MPI_BYTE, &count); - if (count) { - // Receiving actual messages - mpi_msg_t msg; - // Set rank and count, create space to receive message - msg.buf.resize(count / sizeof(T)); - msg.rank = status.MPI_SOURCE; - MPI_Recv(msg.buf.data(), - count, - MPI_BYTE, - status.MPI_SOURCE, - status.MPI_TAG, - comm, - MPI_STATUS_IGNORE); - msg_out.push_back(msg); - // Send ack message to the origin of this message - request.emplace_back(); - MPI_Isend(NULL, 0, MPI_BYTE, status.MPI_SOURCE, - ACK_RECV_TAG, comm, &request.back()); - } else { - // The message is a control signal with size = 0 - MPI_Recv(NULL, 0, MPI_BYTE, status.MPI_SOURCE, - status.MPI_TAG, comm, &status); - switch (status.MPI_TAG) { - case ACK_RECV_TAG: { - // Increase the success_sent - ++success_sent; - if (success_sent == num) { - // The processor finish sending - if (!comm_rank) { - // Increase finished_proc for rank 0 - ++finished_proc; - } else { - // Send proc finish to rank 0 for other processors - request.emplace_back(); - MPI_Isend(NULL, 0, MPI_BYTE, 0, PROC_FINISH_TAG, - comm, &request.back()); - } - } - break; - } - case PROC_FINISH_TAG: { - // Only processor with rank 0 will get this tag - ++finished_proc; - break; - } - case ALL_FINISH_TAG: { - // Exit the while loop - done = true; - break; - } - default: { - // This should never be executed - return false; - } - } - // Check whether everything is done by processor with rank 0 - if (finished_proc == comm_size) { - // Send all finish to all other processors - for (int i = 1; i < comm_size; ++i) { - request.emplace_back(); - MPI_Isend(NULL, 0, MPI_BYTE, i, ALL_FINISH_TAG, - comm, &request.back()); - } - done = true; - } + + break; + } + + case PROC_FINISH_TAG: { + // Only processor with rank 0 will get this tag + ++finished_proc; + break; + } + + case ALL_FINISH_TAG: { + // Exit the while loop + done = true; + break; + } + + default: { + // This should never be executed + return false; + } + } + + // Check whether everything is done by processor with rank 0 + if (finished_proc == comm_size) { + // Send all finish to all other processors + for (int i = 1; i < comm_size; ++i) { + request.emplace_back(); + MPI_Isend(NULL, 0, MPI_BYTE, i, ALL_FINISH_TAG, + comm, &request.back()); } + + done = true; + } } + } - // Correctness check using request - for (auto& it : request) { - int flag; - MPI_Test(&it, &flag, MPI_STATUS_IGNORE); - if (!flag) { - // This should never be executed - return false; - } + // Correctness check using request + for (auto& it : request) { + int flag; + MPI_Test(&it, &flag, MPI_STATUS_IGNORE); + + if (!flag) { + // This should never be executed + return false; } + } - // Promise all send and receive end - MPI_Barrier(comm); + // Promise all send and receive end + MPI_Barrier(comm); - return true; + return true; } /** @@ -350,49 +365,50 @@ template bool CUSTOM_MPI_SENDRECV2(const std::vector> &msg_in, std::vector> &msg_out, const MPI_Comm &comm) { - // Promise previous send and receive have all ended - MPI_Barrier(comm); - - // Only used to check correctness - std::vector request; - - // Send - for (auto& it : msg_in) { - request.emplace_back(); - MPI_Isend(it.buf.data(), - it.buf.size() * sizeof(T), - MPI_BYTE, - it.rank, - 0, - comm, - &request.back()); + // Promise previous send and receive have all ended + MPI_Barrier(comm); + + // Only used to check correctness + std::vector request; + + // Send + for (auto& it : msg_in) { + request.emplace_back(); + MPI_Isend(it.buf.data(), + it.buf.size() * sizeof(T), + MPI_BYTE, + it.rank, + 0, + comm, + &request.back()); + } + + // Receive + for (auto& it : msg_out) { + MPI_Recv(it.buf.data(), + it.buf.size() * sizeof(T), + MPI_BYTE, + it.rank, + 0, + comm, + MPI_STATUS_IGNORE); + } + + // Promise all send and receive end + MPI_Barrier(comm); + + // Correctness check using request + for (auto& it : request) { + int flag; + MPI_Test(&it, &flag, MPI_STATUS_IGNORE); + + if (!flag) { + // This should never be executed + return false; } + } - // Receive - for (auto& it : msg_out) { - MPI_Recv(it.buf.data(), - it.buf.size() * sizeof(T), - MPI_BYTE, - it.rank, - 0, - comm, - MPI_STATUS_IGNORE); - } - - // Promise all send and receive end - MPI_Barrier(comm); - - // Correctness check using request - for (auto& it : request) { - int flag; - MPI_Test(&it, &flag, MPI_STATUS_IGNORE); - if (!flag) { - // This should never be executed - return false; - } - } - - return true; + return true; } /******** Functions related to MPI ends ********/ @@ -402,8 +418,8 @@ bool CUSTOM_MPI_SENDRECV2(const std::vector> &msg_in, // -------------------------------------------------------------------------- struct lon_lat_t { - precision_t lon; - precision_t lat; + precision_t lon; + precision_t lat; }; // -------------------------------------------------------------------------- @@ -411,8 +427,8 @@ struct lon_lat_t { // -------------------------------------------------------------------------- std::ostream& operator<<(std::ostream &out, const lon_lat_t &lonlat) { - out << "lon = " << lonlat.lon << "\tlat = " << lonlat.lat << '\n'; - return out; + out << "lon = " << lonlat.lon << "\tlat = " << lonlat.lat << '\n'; + return out; } // -------------------------------------------------------------------------- @@ -421,23 +437,25 @@ std::ostream& operator<<(std::ostream &out, const lon_lat_t &lonlat) { // -------------------------------------------------------------------------- int wrap_point_sphere(precision_t &lon, precision_t &lat) { - int inverse = 1; - if (lat < -0.5 * cPI) { - lat = -cPI - lat; - inverse = -1; - lon += cPI; - } else if (lat > 0.5 * cPI) { - lat = cPI - lat; - inverse = -1; - lon += cPI; - } - while (lon < 0) { - lon += cTWOPI; - } - while (lon > cTWOPI) { - lon -= cTWOPI; - } - return inverse; + int inverse = 1; + + if (lat < -0.5 * cPI) { + lat = -cPI - lat; + inverse = -1; + lon += cPI; + } else if (lat > 0.5 * cPI) { + lat = cPI - lat; + inverse = -1; + lon += cPI; + } + + while (lon < 0) + lon += cTWOPI; + + while (lon > cTWOPI) + lon -= cTWOPI; + + return inverse; } // -------------------------------------------------------------------------- @@ -445,85 +463,92 @@ int wrap_point_sphere(precision_t &lon, precision_t &lat) { // -------------------------------------------------------------------------- void Grid::init_connection() { - // Initialize communicator. Message exchange is limited in each - // ensemble and rank 0 in each ensemble has special purpose - MPI_Comm_split(aether_comm, iMember, iGrid, &grid_comm); - - // Classify the {lon, lat} for all ghost cells based on destination - std::unordered_map> lonlatmap; - for (int64_t iLon = 0; iLon < nLons; ++iLon) { - for (int64_t iLat = 0; iLat < nLats; ++iLat) { - // Read the lon, lat of ghost cell - precision_t lon = geoLon_scgc(iLon, iLat, 0); - precision_t lat = geoLat_scgc(iLon, iLat, 0); - // Find which processor can handle it and whether - // the message crosses the north/south pole - int dest, inverse = 1; - if (IsCubeSphereGrid) { - dest = pos_iProc_cube(lon, lat); - } else { - // Spherical grid needs wrap - inverse = wrap_point_sphere(lon, lat); - dest = pos_iProc_sphere(lon, lat); - } - - // Prepare to send the lon, lat to dest processor - lonlatmap[dest].push_back({lon, lat}); - // Record where to place data when receiving from that processor - exch_recv[dest].push_back({iLon, iLat, inverse}); - // Skip non ghost cells - if (iLat == nGCs - 1 && iLon >= nGCs && iLon < nLons - nGCs) { - iLat = nLats - nGCs - 1; - } - } + // Initialize communicator. Message exchange is limited in each + // ensemble and rank 0 in each ensemble has special purpose + MPI_Comm_split(aether_comm, iMember, iGrid, &grid_comm); + + // Classify the {lon, lat} for all ghost cells based on destination + std::unordered_map> lonlatmap; + + for (int64_t iLon = 0; iLon < nLons; ++iLon) { + for (int64_t iLat = 0; iLat < nLats; ++iLat) { + // Read the lon, lat of ghost cell + precision_t lon = geoLon_scgc(iLon, iLat, 0); + precision_t lat = geoLat_scgc(iLon, iLat, 0); + // Find which processor can handle it and whether + // the message crosses the north/south pole + int dest, inverse = 1; + + if (IsCubeSphereGrid) + dest = pos_iProc_cube(lon, lat); + + else { + // Spherical grid needs wrap + inverse = wrap_point_sphere(lon, lat); + dest = pos_iProc_sphere(lon, lat); + } + + // Prepare to send the lon, lat to dest processor + lonlatmap[dest].push_back({lon, lat}); + // Record where to place data when receiving from that processor + exch_recv[dest].push_back({iLon, iLat, inverse}); + + // Skip non ghost cells + if (iLat == nGCs - 1 && iLon >= nGCs && iLon < nLons - nGCs) + iLat = nLats - nGCs - 1; } - - // Construct send messages - std::vector> msg_in; - for (auto& it : lonlatmap) { - msg_in.emplace_back(std::move(it.second), it.first); + } + + // Construct send messages + std::vector> msg_in; + + for (auto& it : lonlatmap) + msg_in.emplace_back(std::move(it.second), it.first); + + // Create space to receive messages + std::vector> msg_out; + + // Send all longitude and latitude to where it can be handled + // Receive all longitude and latitude this processor needs to handle + CUSTOM_MPI_SENDRECV(msg_in, msg_out, grid_comm); + + // Build interpolation coefficients based on lon and lat + for (auto& it : msg_out) { + // Initialize the inputs to set_interp_coefs + std::vector Lons; + std::vector Lats; + std::vector Alts; + + // Combine each altitude with lat_lon into points + for (int64_t iAlt = nGCs; iAlt < nAlts - nGCs; ++iAlt) { + for (auto& it2 : it.buf) { + Lons.push_back(it2.lon); + Lats.push_back(it2.lat); + Alts.push_back(geoAlt_scgc(0, 0, iAlt)); + } } - // Create space to receive messages - std::vector> msg_out; - - // Send all longitude and latitude to where it can be handled - // Receive all longitude and latitude this processor needs to handle - CUSTOM_MPI_SENDRECV(msg_in, msg_out, grid_comm); - - // Build interpolation coefficients based on lon and lat - for (auto& it : msg_out) { - // Initialize the inputs to set_interp_coefs - std::vector Lons; - std::vector Lats; - std::vector Alts; - - // Combine each altitude with lat_lon into points - for (int64_t iAlt = nGCs; iAlt < nAlts - nGCs; ++iAlt) { - for (auto& it2 : it.buf) { - Lons.push_back(it2.lon); - Lats.push_back(it2.lat); - Alts.push_back(geoAlt_scgc(0, 0, iAlt)); - } - } - // Compute interpolation coefficients - set_interpolation_coefs(Lons, Lats, Alts); - - // Fix precision issues - for (auto& it : interp_coefs) { - if (it.rRow < 0.0001) { - it.rRow = 0; - } else if (it.rRow > 0.9999) { - it.rRow = 1; - } - if (it.rCol < 0.0001) { - it.rCol = 0; - } else if (it.rCol > 0.9999) { - it.rCol = 1; - } - } - // Store the coefficients - exch_send[it.rank] = std::move(interp_coefs); + + // Compute interpolation coefficients + set_interpolation_coefs(Lons, Lats, Alts); + + // Fix precision issues + for (auto& it : interp_coefs) { + if (it.rRow < 0.0001) + it.rRow = 0; + + else if (it.rRow > 0.9999) + it.rRow = 1; + + if (it.rCol < 0.0001) + it.rCol = 0; + + else if (it.rCol > 0.9999) + it.rCol = 1; } + + // Store the coefficients + exch_send[it.rank] = std::move(interp_coefs); + } } // -------------------------------------------------------------------------- @@ -531,41 +556,43 @@ void Grid::init_connection() { // -------------------------------------------------------------------------- void Grid::exchange(arma_cube &data, const bool pole_inverse) { - // Construct send and receive buffer - std::vector> msg_send; - std::vector> msg_recv; - - // Do the interpolation for other processors - for (auto& it : exch_send) { - std::swap(interp_coefs, it.second); - msg_send.emplace_back(get_interpolation_values(data), - it.first); - std::swap(interp_coefs, it.second); - } - - // Initialize the size of receive buffer - const int64_t numalt = nAlts - 2 * nGCs; - for (auto& it : exch_recv) { - msg_recv.emplace_back(std::vector(it.second.size()*numalt), - it.first); - } - - // Exchange message - CUSTOM_MPI_SENDRECV2(msg_send, msg_recv, grid_comm); - - // Fill the ghost cell with received message - for (auto& it : msg_recv) { - const std::vector &idx = exch_recv[it.rank]; - for (int64_t iAlt = nGCs; iAlt < nAlts - nGCs; ++iAlt) { - for (size_t i = 0; i < idx.size(); ++i) { - if (pole_inverse) { - it.buf[(iAlt - nGCs) * idx.size() + i] *= idx[i].inverse; - } - data(idx[i].ilon, idx[i].ilat, iAlt) = - it.buf[(iAlt - nGCs) * idx.size() + i]; - } - } + // Construct send and receive buffer + std::vector> msg_send; + std::vector> msg_recv; + + // Do the interpolation for other processors + for (auto& it : exch_send) { + std::swap(interp_coefs, it.second); + msg_send.emplace_back(get_interpolation_values(data), + it.first); + std::swap(interp_coefs, it.second); + } + + // Initialize the size of receive buffer + const int64_t numalt = nAlts - 2 * nGCs; + + for (auto& it : exch_recv) { + msg_recv.emplace_back(std::vector(it.second.size()*numalt), + it.first); + } + + // Exchange message + CUSTOM_MPI_SENDRECV2(msg_send, msg_recv, grid_comm); + + // Fill the ghost cell with received message + for (auto& it : msg_recv) { + const std::vector &idx = exch_recv[it.rank]; + + for (int64_t iAlt = nGCs; iAlt < nAlts - nGCs; ++iAlt) { + for (size_t i = 0; i < idx.size(); ++i) { + if (pole_inverse) + it.buf[(iAlt - nGCs) * idx.size() + i] *= idx[i].inverse; + + data(idx[i].ilon, idx[i].ilat, iAlt) = + it.buf[(iAlt - nGCs) * idx.size() + i]; + } } + } } // ----------------------------------------------------------------------------- @@ -573,17 +600,17 @@ void Grid::exchange(arma_cube &data, const bool pole_inverse) { // ----------------------------------------------------------------------------- bool Neutrals::exchange(Grid &grid) { - // For each species, exchange if its DoAdvect is true - for (int i = 0; i < nSpecies; ++i) { - if (species[i].DoAdvect) { - grid.exchange(species[i].density_scgc, false); - } - } - // Exchange temperature - grid.exchange(temperature_scgc, false); - // Exchange velocity - grid.exchange(velocity_vcgc[0], true); - grid.exchange(velocity_vcgc[1], true); - grid.exchange(velocity_vcgc[2], false); - return true; + // For each species, exchange if its DoAdvect is true + for (int i = 0; i < nSpecies; ++i) { + if (species[i].DoAdvect) + grid.exchange(species[i].density_scgc, false); + } + + // Exchange temperature + grid.exchange(temperature_scgc, false); + // Exchange velocity + grid.exchange(velocity_vcgc[0], true); + grid.exchange(velocity_vcgc[1], true); + grid.exchange(velocity_vcgc[2], false); + return true; } diff --git a/src/init_geo_grid.cpp b/src/init_geo_grid.cpp index c7e317a2..e50c1e7f 100644 --- a/src/init_geo_grid.cpp +++ b/src/init_geo_grid.cpp @@ -279,7 +279,8 @@ void transformation_metrics(Quadtree quadtree, A12_inv(i, j) = 0; A21_inv(i, j) = p2 * tan(latp) * tan(lonp - cPI / 4. - 3 * cPI / 2.); A22_inv(i, j) = p2 / cos(latp); - } else if (quadtree.iSide == 6 - 1) { // Face 5 and 6 are flipped than Nair's formulation + } else if (quadtree.iSide == 6 - + 1) { // Face 5 and 6 are flipped than Nair's formulation double p1 = R * sin(latp) / a; A11(i, j) = p1 * cos(lonp - 3 * cPI / 4.); A12(i, j) = p1 * sin(lonp - 3 * cPI / 4.); @@ -291,7 +292,8 @@ void transformation_metrics(Quadtree quadtree, A12_inv(i, j) = -p2 * sin(lonp - 3 * cPI / 4.); A21_inv(i, j) = p2 * sin(latp) * sin(lonp - 3 * cPI / 4.); A22_inv(i, j) = p2 * cos(lonp - 3 * cPI / 4.); - } else if (quadtree.iSide == 5 - 1) { // Face 5 and 6 are flipped than Nair's formulation + } else if (quadtree.iSide == 5 - + 1) { // Face 5 and 6 are flipped than Nair's formulation double p1 = R * sin(latp) / a; A11(i, j) = -p1 * cos(lonp - 3 * cPI / 4.); A12(i, j) = p1 * sin(lonp - 3 * cPI / 4.); @@ -862,7 +864,7 @@ void Grid::correct_xy_grid(Planets planet) { int64_t iAlt; // initialize grid drefx drefy - drefx = arma_vec(nAlts); + drefx = arma_vec(nAlts); drefy = arma_vec(nAlts); // Planet.get_radius() takes in latitude diff --git a/src/inputs.cpp b/src/inputs.cpp index 7385764d..ac912957 100644 --- a/src/inputs.cpp +++ b/src/inputs.cpp @@ -91,8 +91,10 @@ bool Inputs::write_restart() { std::string Inputs::get_logfile() { std::string logfile = check_settings_str("Logfile", "name"); + if (nMembers > 1) logfile = add_cmember(logfile); + return logfile; } @@ -230,7 +232,7 @@ std::string Inputs::get_settings_str(std::string key1, } int Inputs::get_settings(std::string key1, - std::string key2) { + std::string key2) { int value = -1; if (settings.find(key1) != settings.end()) @@ -255,8 +257,8 @@ bool Inputs::check_settings(std::string key1, if (report.test_verbose(1)) std::cout << "checking setting : " - << key1 << " and " - << key2 << "\n"; + << key1 << " and " + << key2 << "\n"; //try to find the keys first if (settings.find(key1) != settings.end()) { @@ -516,6 +518,7 @@ int Inputs::get_original_seed() { std::cout << "Error in getting seed!\n"; return 0; } + return settings.at("Seed"); } @@ -582,20 +585,24 @@ bool Inputs::set_verbose(json in) { bool DidWork = true; int iVerbose = -1; + // Want to set verbose level ASAP: if (in.contains("Debug")) { if (in.at("Debug").contains("iVerbose")) { iVerbose = in.at("Debug").at("iVerbose"); + if (in.at("Debug").contains("iProc")) { - if (iProc != in.at("Debug").at("iProc")) - iVerbose = -1; + if (iProc != in.at("Debug").at("iProc")) + iVerbose = -1; } - } + } } + if (iVerbose > 0) { std::cout << "Setting iVerbose : " << iVerbose << "\n"; report.set_verbose(iVerbose); } + return DidWork; } diff --git a/src/logfile.cpp b/src/logfile.cpp index 1d43d326..7d5db759 100644 --- a/src/logfile.cpp +++ b/src/logfile.cpp @@ -253,7 +253,7 @@ Logfile::Logfile(Indices &indices) { std::vector sat_dts = input.get_satellite_dts(); bool isRestart = input.get_do_restart(); - + // Only the 0th processor in each member needs to open the logfile if (iGrid == 0) { @@ -264,6 +264,7 @@ Logfile::Logfile(Indices &indices) { logfilestream.open(logfileName, std::ofstream::trunc); logfilestream.precision(4); } + // Report error if can not open the log file stream if (!logfilestream.is_open() & !isRestart) { // TRY TO EXIT GRACEFULLY HERE. ALL THE FOLLOWING CODE SHOULD @@ -297,6 +298,7 @@ Logfile::Logfile(Indices &indices) { // Write the header to log if (!isRestart) logfilestream << header_time << header_log << '\n'; + // Close the file stream if append if (doAppend) logfilestream.close(); diff --git a/src/neutrals.cpp b/src/neutrals.cpp index 0b5c73c7..519a9711 100644 --- a/src/neutrals.cpp +++ b/src/neutrals.cpp @@ -140,6 +140,7 @@ Neutrals::Neutrals(Grid grid, // This gets a bunch of the species-dependent characteristics: iErr = read_planet_file(planet); + if (iErr > 0) report.error("Error reading planet file!"); diff --git a/src/neutrals_bcs.cpp b/src/neutrals_bcs.cpp index b8d497e0..0c7c4017 100644 --- a/src/neutrals_bcs.cpp +++ b/src/neutrals_bcs.cpp @@ -32,8 +32,10 @@ bool Neutrals::set_bcs(Grid grid, if (input.get_nAltsGeo() > 1) { didWork = set_lower_bcs(grid, time, indices); + if (didWork) didWork = set_upper_bcs(grid); + if (didWork) calc_mass_density(); } @@ -132,6 +134,7 @@ bool Neutrals::set_lower_bcs(Grid grid, if (!msis.is_ok()) { didWork = false; report.error("MSIS initialization not ok"); + if (report.test_verbose(0)) { std::cout << "MSIS Boundary Conditions asked for, "; std::cout << "but MSIS is not compiled! Yikes!\n"; @@ -191,6 +194,7 @@ bool Neutrals::set_lower_bcs(Grid grid, species[iSpecies].density_scgc.slice(0). fill(species[iSpecies].lower_bc_density); } + temperature_scgc.slice(0).fill(initial_temperatures[0]); didWork = true; } @@ -220,7 +224,7 @@ bool Neutrals::set_lower_bcs(Grid grid, report.error("issue with lower BCs!"); report.error("maybe check boundaryconditions type : " + bcsType); } - + report.exit(function); return didWork; } diff --git a/src/neutrals_ics.cpp b/src/neutrals_ics.cpp index 558a5b07..9f67e005 100644 --- a/src/neutrals_ics.cpp +++ b/src/neutrals_ics.cpp @@ -17,8 +17,8 @@ // ----------------------------------------------------------------------------- bool Neutrals::initial_conditions(Grid grid, - Times time, - Indices indices) { + Times time, + Indices indices) { std::string function = "Neutrals::initial_conditions"; static int iFunction = -1; @@ -27,7 +27,7 @@ bool Neutrals::initial_conditions(Grid grid, // initialize didWork to false, so that we can catch if the initial // conditions are not actually set bool didWork = false; - + int64_t iLon, iLat, iAlt, iA; precision_t alt, r; int64_t nAlts = grid.get_nZ(); @@ -37,8 +37,10 @@ bool Neutrals::initial_conditions(Grid grid, if (input.get_do_restart()) { report.print(1, "Restarting! Reading neutral files!"); didWork = restart_file(input.get_restartin_dir(), DoRead); + if (!didWork) { - report.error("Reading Restart for Neutrals Failed!!!");\ + report.error("Reading Restart for Neutrals Failed!!!"); + \ } } else { @@ -53,48 +55,48 @@ bool Neutrals::initial_conditions(Grid grid, if (!msis.is_ok()) { didWork = false; - report.error("MSIS Initial Conditions asked for, "); + report.error("MSIS Initial Conditions asked for, "); report.error("but MSIS is not compiled/working!"); } else { - didWork = true; - msis.set_time(time); - precision_t f107 = indices.get_f107(time.get_current()); - precision_t f107a = indices.get_f107a(time.get_current()); - msis.set_f107(f107, f107a); - msis.set_ap(10.0); - msis.set_locations(grid.geoLon_scgc, - grid.geoLat_scgc, - grid.geoAlt_scgc); - - // This is just to check if MSIS is actually working: - if (msis.is_valid_species("Tn")) - // if it is, fill will temperature: - temperature_scgc = msis.get_cube("Tn"); - else - // if it is not, then fill with a value: - temperature_scgc.fill(initial_temperatures[0]); - - for (int iSpecies = 0; iSpecies < nSpecies; iSpecies++) { - if (report.test_verbose(3)) - std::cout << "Setting Species : " - << species[iSpecies].cName << "\n"; - - if (msis.is_valid_species(species[iSpecies].cName)) { - if (report.test_verbose(3)) - std::cout << " Found in MSIS!\n"; - - species[iSpecies].density_scgc = - msis.get_cube(species[iSpecies].cName); - } else { - if (report.test_verbose(3)) - std::cout << " NOT Found in MSIS - setting constant\n"; - - species[iSpecies].density_scgc.slice(0). - fill(species[iSpecies].lower_bc_density); - fill_with_hydrostatic(iSpecies, 1, nAlts, grid); - } - - } // for species + didWork = true; + msis.set_time(time); + precision_t f107 = indices.get_f107(time.get_current()); + precision_t f107a = indices.get_f107a(time.get_current()); + msis.set_f107(f107, f107a); + msis.set_ap(10.0); + msis.set_locations(grid.geoLon_scgc, + grid.geoLat_scgc, + grid.geoAlt_scgc); + + // This is just to check if MSIS is actually working: + if (msis.is_valid_species("Tn")) + // if it is, fill will temperature: + temperature_scgc = msis.get_cube("Tn"); + else + // if it is not, then fill with a value: + temperature_scgc.fill(initial_temperatures[0]); + + for (int iSpecies = 0; iSpecies < nSpecies; iSpecies++) { + if (report.test_verbose(3)) + std::cout << "Setting Species : " + << species[iSpecies].cName << "\n"; + + if (msis.is_valid_species(species[iSpecies].cName)) { + if (report.test_verbose(3)) + std::cout << " Found in MSIS!\n"; + + species[iSpecies].density_scgc = + msis.get_cube(species[iSpecies].cName); + } else { + if (report.test_verbose(3)) + std::cout << " NOT Found in MSIS - setting constant\n"; + + species[iSpecies].density_scgc.slice(0). + fill(species[iSpecies].lower_bc_density); + fill_with_hydrostatic(iSpecies, 1, nAlts, grid); + } + + } // for species } // msis init worked ok } // type = msis @@ -171,7 +173,7 @@ bool Neutrals::initial_conditions(Grid grid, if (!didWork) report.error("Issue with initial conditions!"); - + report.exit(function); return didWork; } diff --git a/src/neutrals_momentum_friction.cpp b/src/neutrals_momentum_friction.cpp index 78a738e7..30493882 100644 --- a/src/neutrals_momentum_friction.cpp +++ b/src/neutrals_momentum_friction.cpp @@ -94,7 +94,7 @@ void calc_neutral_friction(Neutrals &neutrals) { for (iSpecies = 0; iSpecies < neutrals.nSpecies; iSpecies++) for (iDir = 0; iDir < 3; iDir++) neutrals.species[iSpecies].acc_neutral_friction[iDir].zeros(); - + if (input.get_advection_neutrals_vertical() != "hydro") { arma_vec vels(neutrals.nSpeciesAdvect, fill::zeros); diff --git a/src/output.cpp b/src/output.cpp index f35b7ee5..01ea5f86 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -20,10 +20,10 @@ // ----------------------------------------------------------------------------- bool output(const Neutrals &neutrals, - const Ions &ions, - const Grid &grid, - Times time, - const Planets &planet) { + const Ions &ions, + const Grid &grid, + Times time, + const Planets &planet) { std::string function = "output"; static int iFunction = -1; @@ -311,8 +311,8 @@ bool output(const Neutrals &neutrals, // write output container if (!AllOutputContainers[iOutput].write()) { - report.error("Error in writing output container!"); - didWork = false; + report.error("Error in writing output container!"); + didWork = false; } // ------------------------------------------------------------ diff --git a/src/read_input_file.cpp b/src/read_input_file.cpp index 6acf4990..fbc3786d 100644 --- a/src/read_input_file.cpp +++ b/src/read_input_file.cpp @@ -29,7 +29,7 @@ bool Inputs::read_inputs_json(Times &time) { // Set the default values first: settings = read_json("UA/inputs/defaults.json"); DidWork = set_verbose(settings); - + // Set the planet-specific file (user can change this in aether.in file!): settings["PlanetSpeciesFile"] = settings["Planet"]["file"]; @@ -50,9 +50,9 @@ bool Inputs::read_inputs_json(Times &time) { restart_file = restart_file + "/settings.json"; json restart_inputs; restart_inputs = read_json(restart_file); - // This forces the logfile to append. User can override - // if they really want: - restart_inputs["Logfile"]["append"] = true; + // This forces the logfile to append. User can override + // if they really want: + restart_inputs["Logfile"]["append"] = true; settings.merge_patch(restart_inputs); } } diff --git a/src/solver_gradients.cpp b/src/solver_gradients.cpp index 0b587a98..f03f4228 100644 --- a/src/solver_gradients.cpp +++ b/src/solver_gradients.cpp @@ -186,35 +186,38 @@ std::vector calc_gradient_cubesphere(arma_cube value, Grid grid) { arma_mat grad_x_curr(nXs, nYs); arma_mat grad_y_curr(nXs, nYs); - + // Calc gradient on x and y direction (since reference grid) // Only update interior cells // May vectorize for future improvements - if (nGCs >=2) { // if more than 1 nGCs, we do fourth order, some foolproofing in case we go into debug hell - for (int j = nGCs; j < nYs - nGCs; j++) - { - for (int i = nGCs; i < nXs - nGCs; i++) - { - grad_x_curr(i, j) = (-curr_value(i+2,j) + 8*curr_value(i+1,j) - 8*curr_value(i-1,j) + curr_value(i-2,j))*(1./12./dx); - grad_y_curr(i, j) = (-curr_value(i,j+2) + 8*curr_value(i,j+1) - 8*curr_value(i,j-1) + curr_value(i,j-2))*(1./12./dy); + if (nGCs >= + 2) { // if more than 1 nGCs, we do fourth order, some foolproofing in case we go into debug hell + for (int j = nGCs; j < nYs - nGCs; j++) { + for (int i = nGCs; i < nXs - nGCs; i++) { + grad_x_curr(i, j) = (-curr_value(i + 2, j) + 8 * curr_value(i + 1, + j) - 8 * curr_value(i - 1, j) + curr_value(i - 2, j)) * (1. / 12. / dx); + grad_y_curr(i, j) = (-curr_value(i, j + 2) + 8 * curr_value(i, + j + 1) - 8 * curr_value(i, j - 1) + curr_value(i, j - 2)) * (1. / 12. / dy); } } } else { // otherwise we do second order - for (int j = nGCs; j < nYs - nGCs; j++) - { - for (int i = nGCs; i < nXs - nGCs; i++) - { - grad_x_curr(i, j) = (curr_value(i+1,j)-curr_value(i-1,j))*(1./2./dx); - grad_y_curr(i, j) = (curr_value(i,j+1)-curr_value(i,j-1))*(1./2./dy); + for (int j = nGCs; j < nYs - nGCs; j++) { + for (int i = nGCs; i < nXs - nGCs; i++) { + grad_x_curr(i, j) = (curr_value(i + 1, j) - curr_value(i - 1, + j)) * (1. / 2. / dx); + grad_y_curr(i, j) = (curr_value(i, j + 1) - curr_value(i, + j - 1)) * (1. / 2. / dy); } } } - + // We then use A transformation matrices to convert grad_xy to grad_latlon // Ref -> Physical, we use A matrix - grad_lon.slice(iAlt) = grad_x_curr%grid.A11_inv_scgc.slice(iAlt)+grad_y_curr%grid.A21_inv_scgc.slice(iAlt); - grad_lat.slice(iAlt) = grad_x_curr%grid.A12_inv_scgc.slice(iAlt)+grad_y_curr%grid.A22_inv_scgc.slice(iAlt); + grad_lon.slice(iAlt) = grad_x_curr % grid.A11_inv_scgc.slice( + iAlt) + grad_y_curr % grid.A21_inv_scgc.slice(iAlt); + grad_lat.slice(iAlt) = grad_x_curr % grid.A12_inv_scgc.slice( + iAlt) + grad_y_curr % grid.A22_inv_scgc.slice(iAlt); } // Not initializing with array like procedure in case I get bugs diff --git a/src/tools.cpp b/src/tools.cpp index f8643e17..92b1fb38 100644 --- a/src/tools.cpp +++ b/src/tools.cpp @@ -119,8 +119,10 @@ bool compare(precision_t value1, precision_t value2) { std::string add_cmember(std::string inString) { std::string outString = inString; std::size_t found = outString.rfind("."); - if (found!=std::string::npos) + + if (found != std::string::npos) outString.replace(found, 1, "_" + cMember + "."); + return outString; } diff --git a/src/transform.cpp b/src/transform.cpp index 6c5973dd..33a23bbe 100644 --- a/src/transform.cpp +++ b/src/transform.cpp @@ -13,8 +13,10 @@ std::string mklower(std::string inString) { std::string outString = inString; int64_t nChars = outString.length(); + for (int64_t iChar = 0; iChar < nChars; iChar++) outString[iChar] = tolower(outString[iChar]); + return outString; }