diff --git a/CMakeLists.txt b/CMakeLists.txt index 0da83cfd..a479a1fe 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -31,6 +31,8 @@ elseif(TEST_COORD) elseif(TEST_EXCHANGE) add_executable(aether ${SRC_FILES} ${MSIS_FILES} ${IE_FILES} ${MAIN_DIR}/main_test_exchange.cpp) set(USE_DOUBLE_PRECISION True) +elseif(TEST_GRADIENT) + add_executable(aether ${SRC_FILES} ${MSIS_FILES} ${MAIN_DIR}/main_test_gradient.cpp) else() add_executable(aether ${SRC_FILES} ${MSIS_FILES} ${IE_FILES} ${MAIN_DIR}/main.cpp) endif() 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/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, diff --git a/include/electrodynamics.h b/include/electrodynamics.h index 14504d41..13a631ab 100644 --- a/include/electrodynamics.h +++ b/include/electrodynamics.h @@ -67,11 +67,12 @@ class Electrodynamics { \param ions Going to set the potential and aurora **/ - int update(Planets planet, - Grid gGrid, - Times time, - Indices &indices, - Ions &ions); + bool update(Planets planet, + Grid gGrid, + Times time, + Indices &indices, + Ions &ions); + /************************************************************** \brief used in main.cpp to ensure electrodynamics times and 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/include/inputs.h b/include/inputs.h index 5dbf7354..4fad713c 100644 --- a/include/inputs.h +++ b/include/inputs.h @@ -123,6 +123,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/include/neutrals.h b/include/neutrals.h index 1226f2c8..47436e1b 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 @@ -339,15 +339,8 @@ class Neutrals { \param grid The grid to define the neutrals on **/ precision_t calc_dt(Grid grid); - - /********************************************************************** - \brief Calculate dt (cell size / cMax) in each direction, and take min - \param dt returns the neutral time-step - \param grid The grid to define the neutrals on - This function is for cubesphere dt calculation only - **/ precision_t calc_dt_cubesphere(Grid grid); - + /********************************************************************** \brief Calculate the chapman integrals for the individual species \param grid The grid to define the neutrals on 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/include/solvers.h b/include/solvers.h index ed4f36ef..3a20e3d4 100644 --- a/include/solvers.h +++ b/include/solvers.h @@ -56,6 +56,7 @@ arma_cube calc_gradient_lon(arma_cube value, Grid grid); arma_cube calc_gradient_lat(arma_cube value, Grid grid); arma_cube calc_gradient_alt(arma_cube value, Grid grid); std::vector calc_gradient_vector(arma_cube value_scgc, Grid grid); +std::vector calc_gradient_cubesphere(arma_cube value, Grid grid); arma_cube calc_gradient_alt_4th(arma_cube value, Grid grid); // interpolation in 1D 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/include/transform.h b/include/transform.h index 8db3ccab..a56313b7 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); void copy_mat_to_array(arma_mat mat_in, diff --git a/src/advance.cpp b/src/advance.cpp index d27f9cde..ce478c4f 100644 --- a/src/advance.cpp +++ b/src/advance.cpp @@ -10,18 +10,18 @@ // 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; @@ -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,61 +64,75 @@ 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, - indices, - 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, + indices, + ions); - calc_aurora(gGrid, neutrals, ions); - // Calculate some neutral source terms: - neutrals.calc_conduction(gGrid, time); - chemistry.calc_chemistry(neutrals, ions, time, gGrid); + if (didWork) { + calc_ion_neutral_coll_freq(neutrals, ions); + ions.calc_ion_drift(neutrals, gGrid, time.get_dt()); - if (input.get_O_cooling()) - neutrals.calc_O_cool(); + calc_aurora(gGrid, neutrals, ions); - if (input.get_NO_cooling()) - neutrals.calc_NO_cool(); + // Calculate some neutral source terms: + neutrals.calc_conduction(gGrid, time); + chemistry.calc_chemistry(neutrals, ions, time, gGrid); - neutrals.calc_conduction(gGrid, time); - chemistry.calc_chemistry(neutrals, ions, time, gGrid); + if (input.get_O_cooling()) + neutrals.calc_O_cool(); - neutrals.vertical_momentum_eddy(gGrid); - calc_ion_collisions(neutrals, ions); - calc_neutral_friction(neutrals); + if (input.get_NO_cooling()) + neutrals.calc_NO_cool(); - neutrals.add_sources(time); + neutrals.calc_conduction(gGrid, time); + chemistry.calc_chemistry(neutrals, ions, time, gGrid); - ions.calc_ion_temperature(neutrals, gGrid, time); - ions.calc_electron_temperature(neutrals, gGrid); + neutrals.vertical_momentum_eddy(gGrid); + calc_ion_collisions(neutrals, ions); + calc_neutral_friction(neutrals); - neutrals.exchange(gGrid); + neutrals.add_sources(time); - time.increment_time(); + ions.calc_ion_temperature(neutrals, gGrid, time); + ions.calc_electron_temperature(neutrals, gGrid); - 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 (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())) { + 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); + + if (didWork) + didWork = logfile.write_logfile(indices, neutrals, ions, gGrid, time); - logfile.write_logfile(indices, neutrals, ions, gGrid, time); + if (!didWork) + report.error("Error in Advance!"); report.exit(function); - return iErr; + return didWork; } diff --git a/src/calc_euv.cpp b/src/calc_euv.cpp index bd13ca5e..09be872e 100644 --- a/src/calc_euv.cpp +++ b/src/calc_euv.cpp @@ -13,47 +13,55 @@ // - calculate ionization and heating // ----------------------------------------------------------------------------- -int calc_euv(Planets planet, - Grid grid, - Times time, - Euv &euv, - Neutrals &neutrals, - Ions &ions, - Indices indices) { - - int iErr = 0; - - 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) { - 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); - - calc_ionization_heating(euv, neutrals, ions); - } else - neutrals.heating_euv_scgc.zeros(); - - report.exit(function); - } - - return iErr; +bool calc_euv(Planets planet, + 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() + "?"); + + // 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.exit(function); + return didWork; } // ----------------------------------------------------------------------------- diff --git a/src/calc_neutral_derived.cpp b/src/calc_neutral_derived.cpp index 6d28076c..d101f567 100644 --- a/src/calc_neutral_derived.cpp +++ b/src/calc_neutral_derived.cpp @@ -328,8 +328,7 @@ void Neutrals::calc_cMax() { } // ---------------------------------------------------------------------- -// Calculate cMax, which is the sound speed + velocity in each -// direction +// Calculate dt primarily for the spherical grid // ---------------------------------------------------------------------- precision_t Neutrals::calc_dt(Grid grid) { @@ -373,6 +372,10 @@ precision_t Neutrals::calc_dt(Grid grid) { return dt; } +// ---------------------------------------------------------------------- +// Calculate dt for the cubesphere grid. +// ---------------------------------------------------------------------- + precision_t Neutrals::calc_dt_cubesphere(Grid grid) { std::string function = "Neutrals::calc_dt_cubesphere"; diff --git a/src/electrodynamics.cpp b/src/electrodynamics.cpp index b9d8d408..c52eee44 100644 --- a/src/electrodynamics.cpp +++ b/src/electrodynamics.cpp @@ -132,11 +132,12 @@ void Electrodynamics::set_all_indices_for_ie(Times time, // Update Electrodynamics // ----------------------------------------------------------------------------- -int Electrodynamics::update(Planets planet, - Grid gGrid, - Times time, - Indices &indices, - Ions &ions) { +bool Electrodynamics::update(Planets planet, + Grid gGrid, + Times time, + Indices &indices, + Ions &ions) { + std::string function = "Electrodynamics::update"; static int iFunction = -1; @@ -209,7 +210,7 @@ int Electrodynamics::update(Planets planet, } report.exit(function); - return 0; + return true; } // ----------------------------------------------------------------------------- @@ -331,8 +332,6 @@ 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, diff --git a/src/euv.cpp b/src/euv.cpp index 8ddb2c59..c4a4069a 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; } // -------------------------------------------------------------------------- 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 2f29467c..e50c1e7f 100644 --- a/src/init_geo_grid.cpp +++ b/src/init_geo_grid.cpp @@ -194,7 +194,7 @@ void transformation_metrics(Quadtree quadtree, int64_t nY = lat2d.n_cols; // Assume R = 1 (since lat-lon/ xy generation assumes unit vect) double R = 1; - double a = 1 / sqrt(3); + double a = R / sqrt(3); double xref, yref, rref; double latp, lonp; double g; @@ -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 == 5 - 1) { + } 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 == 6 - 1) { + } 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 @@ -879,19 +881,19 @@ void Grid::correct_xy_grid(Planets planet) { precision_t R = R_Alts(iAlt); refx_scgc.slice(iAlt) *= R; refy_scgc.slice(iAlt) *= R; - A11_scgc.slice(iAlt) *= R; - A12_scgc.slice(iAlt) *= R; - A21_scgc.slice(iAlt) *= R; - A22_scgc.slice(iAlt) *= R; - A11_inv_scgc.slice(iAlt) /= R; - A12_inv_scgc.slice(iAlt) /= R; - A21_inv_scgc.slice(iAlt) /= R; - A22_inv_scgc.slice(iAlt) /= R; - g11_upper_scgc.slice(iAlt) /= R * R; - g12_upper_scgc.slice(iAlt) /= R * R; - g21_upper_scgc.slice(iAlt) /= R * R; - g22_upper_scgc.slice(iAlt) /= R * R; - sqrt_g_scgc.slice(iAlt) *= R * R; + //A11_scgc.slice(iAlt) *= R; + //A12_scgc.slice(iAlt) *= R; + //A21_scgc.slice(iAlt) *= R; + //A22_scgc.slice(iAlt) *= R; + //A11_inv_scgc.slice(iAlt) /= R; + //A12_inv_scgc.slice(iAlt) /= R; + //A21_inv_scgc.slice(iAlt) /= R; + //A22_inv_scgc.slice(iAlt) /= R; + //g11_upper_scgc.slice(iAlt) /= R * R; + //g12_upper_scgc.slice(iAlt) /= R * R; + //g21_upper_scgc.slice(iAlt) /= R * R; + //g22_upper_scgc.slice(iAlt) /= R * R; + //sqrt_g_scgc.slice(iAlt) *= R * R; // Addition: Get a copy of dx dy arma_mat curr_refx = refx_scgc.slice(iAlt); @@ -902,35 +904,35 @@ void Grid::correct_xy_grid(Planets planet) { refx_Left.slice(iAlt) *= R; refy_Left.slice(iAlt) *= R; - A11_Left.slice(iAlt) *= R; - A12_Left.slice(iAlt) *= R; - A21_Left.slice(iAlt) *= R; - A22_Left.slice(iAlt) *= R; - A11_inv_Left.slice(iAlt) /= R; - A12_inv_Left.slice(iAlt) /= R; - A21_inv_Left.slice(iAlt) /= R; - A22_inv_Left.slice(iAlt) /= R; - g11_upper_Left.slice(iAlt) /= R * R; - g12_upper_Left.slice(iAlt) /= R * R; - g21_upper_Left.slice(iAlt) /= R * R; - g22_upper_Left.slice(iAlt) /= R * R; - sqrt_g_Left.slice(iAlt) *= R * R; + //A11_Left.slice(iAlt) *= R; + //A12_Left.slice(iAlt) *= R; + //A21_Left.slice(iAlt) *= R; + //A22_Left.slice(iAlt) *= R; + //A11_inv_Left.slice(iAlt) /= R; + //A12_inv_Left.slice(iAlt) /= R; + //A21_inv_Left.slice(iAlt) /= R; + //A22_inv_Left.slice(iAlt) /= R; + //g11_upper_Left.slice(iAlt) /= R * R; + //g12_upper_Left.slice(iAlt) /= R * R; + //g21_upper_Left.slice(iAlt) /= R * R; + //g22_upper_Left.slice(iAlt) /= R * R; + //sqrt_g_Left.slice(iAlt) *= R * R; refx_Down.slice(iAlt) *= R; refy_Down.slice(iAlt) *= R; - A11_Down.slice(iAlt) *= R; - A12_Down.slice(iAlt) *= R; - A21_Down.slice(iAlt) *= R; - A22_Down.slice(iAlt) *= R; - A11_inv_Down.slice(iAlt) /= R; - A12_inv_Down.slice(iAlt) /= R; - A21_inv_Down.slice(iAlt) /= R; - A22_inv_Down.slice(iAlt) /= R; - g11_upper_Down.slice(iAlt) /= R * R; - g12_upper_Down.slice(iAlt) /= R * R; - g21_upper_Down.slice(iAlt) /= R * R; - g22_upper_Down.slice(iAlt) /= R * R; - sqrt_g_Down.slice(iAlt) *= R * R; + //A11_Down.slice(iAlt) *= R; + //A12_Down.slice(iAlt) *= R; + //A21_Down.slice(iAlt) *= R; + //A22_Down.slice(iAlt) *= R; + //A11_inv_Down.slice(iAlt) /= R; + //A12_inv_Down.slice(iAlt) /= R; + //A21_inv_Down.slice(iAlt) /= R; + //A22_inv_Down.slice(iAlt) /= R; + //g11_upper_Down.slice(iAlt) /= R * R; + //g12_upper_Down.slice(iAlt) /= R * R; + //g21_upper_Down.slice(iAlt) /= R * R; + //g22_upper_Down.slice(iAlt) /= R * R; + //sqrt_g_Down.slice(iAlt) *= R * R; refx_Corner.slice(iAlt) *= R; refy_Corner.slice(iAlt) *= R; diff --git a/src/inputs.cpp b/src/inputs.cpp index 09ad677e..186bed3a 100644 --- a/src/inputs.cpp +++ b/src/inputs.cpp @@ -90,7 +90,12 @@ 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 +231,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 // ----------------------------------------------------------------------- @@ -241,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()) { @@ -552,7 +568,7 @@ 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"); @@ -621,20 +637,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 ff657a6f..7d5db759 100644 --- a/src/logfile.cpp +++ b/src/logfile.cpp @@ -252,20 +252,25 @@ 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); + bool isRestart = input.get_do_restart(); - else { - logfilestream.open(logfileName, std::ofstream::trunc); - logfilestream.precision(4); - } + // Only the 0th processor in each member needs to open the logfile - // 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() & !isRestart) { + // 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 +293,16 @@ 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'; + // only the 0th processor writes to the logfile + if (logfilestream.is_open()) { + // Write the header to log + if (!isRestart) + logfilestream << header_time << header_log << '\n'; - // Close the file stream if append - if (doAppend) - logfilestream.close(); + // Close the file stream if append + if (doAppend) + logfilestream.close(); + } // Check whether the input settings contain Satellites if (!sat_files.empty()) { @@ -329,7 +338,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 +420,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 +429,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 +466,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; } 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/main/main_test_gradient.cpp b/src/main/main_test_gradient.cpp new file mode 100644 index 00000000..467bd53d --- /dev/null +++ b/src/main/main_test_gradient.cpp @@ -0,0 +1,188 @@ +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) +// Full license can be found in License.md + +#include "../include/aether.h" + +/** + * Output function for debugging + * + * @param values Values + * @param filename FileName + * @param DoAppend + */ +void output(arma_mat &values, + std::string filename, + bool DoAppend) +{ + + std::ofstream outfile; + if (DoAppend) + outfile.open(filename, std::ios_base::app); + else + { + outfile.open(filename); + int64_t nX = values.n_rows; + int64_t nY = values.n_cols; + outfile << nX << " " << nY << "\n"; + } + outfile << values; + outfile << "----"; + outfile.close(); +} + +int main() { + + int iErr = 0; + std::string sError; + bool DidWork = true; + + Times time; + + // Define the function and report: + std::string function = "main"; + static int iFunction = -1; + report.enter(function, iFunction); + + try { + + // Create inputs (reading the input file): + input = Inputs(time); + if (!input.is_ok()) + throw std::string("input initialization failed!"); + + Quadtree quadtree; + if (!quadtree.is_ok()) + throw std::string("quadtree initialization failed!"); + + // Initialize MPI and parallel aspects of the code: + DidWork = init_parallel(quadtree); + if (!DidWork) + throw std::string("init_parallel failed!"); + + // Initialize the planet: + Planets planet; + MPI_Barrier(aether_comm); + if (!planet.is_ok()) + throw std::string("planet initialization failed!"); + + // Initialize Geographic grid: + Grid gGrid(input.get_nLonsGeo(), + input.get_nLatsGeo(), + input.get_nAltsGeo(), + nGeoGhosts); + DidWork = gGrid.init_geo_grid(quadtree, planet); + + // First check whether the initialization uses exactly 6 processes. + // The exactly 6 requirements is due to the checking of the range of reference coordinate system + int world_size; + MPI_Comm_size(aether_comm, &world_size); + if (world_size != 6) { + throw std::string("Comm size must be 6!!!"); + } + + // Gradient Test + { + // Set tolerance limit + precision_t tol = 1e-5; + + // Print current side number + std::string side_num = std::to_string(quadtree.iSide + 1); + std::cout << "Initiating Test 1 for Side Number (1-based index): " << side_num << std::endl; + + /** + * Extract some test data generated by Aether Model + */ + + // Cell center coordinates + arma_mat aether_lon_cc = gGrid.geoLon_scgc.slice(0); + arma_mat aether_lat_cc = gGrid.geoLat_scgc.slice(0); + + int64_t nXs = gGrid.get_nY(); + int64_t nYs = gGrid.get_nX(); + int64_t nGCs = gGrid.get_nGCs(); + int64_t nAlts = gGrid.get_nAlts(); + + // Test scalar field and gradients + arma_cube scgc(nXs, nYs, nAlts); + arma_cube grad_lon_analytical(nXs, nYs, nAlts); + arma_cube grad_lat_analytical(nXs, nYs, nAlts); + + // Radius Information + precision_t planet_R = planet.get_radius(0); + // radius of planet + altitude + // just pick alt at (0,0) loction + arma_vec R_Alts = gGrid.geoAlt_scgc.tube(0, 0) + planet_R; + + for (int iAlt = 0; iAlt < nAlts; iAlt++) { + arma_mat curr_scalar(nXs, nYs, arma::fill::zeros); // setup zero mat + arma_mat curr_grad_lon(nXs, nYs); + arma_mat curr_grad_lat(nXs, nYs); + precision_t A = 1; + precision_t B = 1; + + for (int j = 0; j < nYs; j++) + { + for (int i = 0; i < nXs; i++) + { + precision_t curr_lat = aether_lat_cc(i, j); + precision_t curr_lon = aether_lon_cc(i, j); + + curr_scalar(i,j) = std::sin(curr_lat); + curr_grad_lon(i,j) = 0.; + curr_grad_lat(i,j) = std::cos(curr_lat); // Assume R=1, we will scale the numerical result + } + } + scgc.slice(iAlt) = curr_scalar; + grad_lon_analytical.slice(iAlt) = curr_grad_lon; + grad_lat_analytical.slice(iAlt) = curr_grad_lat; + } + + std::vector test_res = calc_gradient_cubesphere(scgc, gGrid); + + // Perform Tests + for (int iAlt = 0; iAlt < nAlts; iAlt++) { + arma_mat curr_grad_lon = grad_lon_analytical.slice(iAlt); + arma_mat curr_grad_lat = grad_lat_analytical.slice(iAlt); + arma_mat curr_numgrad_lon = test_res[0].slice(iAlt); + arma_mat curr_numgrad_lat = test_res[1].slice(iAlt); + + + // Evaluate actual cells only + for (int j = nGCs; j < nYs - nGCs; j++) + { + for (int i = nGCs; i < nXs - nGCs; i++) + { + if (std::abs(curr_grad_lat(i,j) - curr_numgrad_lat(i,j) * R_Alts(iAlt)) > 1e-4) { // For float precision + std::cout << "Found Incorrect latitudinal gradient for face " + side_num + ", test f = sin(lat)" << std::endl; + std::cout << std::abs(curr_grad_lat(i,j) - curr_numgrad_lat(i,j)* R_Alts(iAlt)) << std::endl; + std::cout << iAlt << std::endl; + goto endloop1; + } + + if (std::abs(curr_grad_lon(i,j) - curr_numgrad_lon(i,j) * R_Alts(iAlt)) > 1e-4) { // For float precision + std::cout << "Found Incorrect longitudinal gradient for face " + side_num + ", test f = sin(lat)" << std::endl; + goto endloop1; + } + } + } + } + } + endloop1: + + report.exit(function); + report.times(); + + } catch (std::string error) { + //if (iProc == 0) { + std::cout << error << "\n"; + std::cout << "---- Must Exit! ----\n"; + //} + } + MPI_Barrier(aether_comm); + + + // End parallel tasks: + iErr = MPI_Finalize(); + + return iErr; +} \ No newline at end of file diff --git a/src/neutrals.cpp b/src/neutrals.cpp index 17659042..519a9711 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(); @@ -141,13 +142,13 @@ Neutrals::Neutrals(Grid grid, 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; } diff --git a/src/neutrals_bcs.cpp b/src/neutrals_bcs.cpp index 4939236d..0c7c4017 100644 --- a/src/neutrals_bcs.cpp +++ b/src/neutrals_bcs.cpp @@ -32,10 +32,17 @@ 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 +113,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 +133,14 @@ 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 +185,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"); @@ -185,6 +196,7 @@ bool Neutrals::set_lower_bcs(Grid grid, } temperature_scgc.slice(0).fill(initial_temperatures[0]); + didWork = true; } // fill the second+ grid cells with the bottom temperature: @@ -208,6 +220,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..9f67e005 100644 --- a/src/neutrals_ics.cpp +++ b/src/neutrals_ics.cpp @@ -16,15 +16,18 @@ // file and fill with hydrostatic. // ----------------------------------------------------------------------------- -int Neutrals::initial_conditions(Grid grid, - Times time, - Indices indices) { +bool Neutrals::initial_conditions(Grid grid, + Times time, + Indices indices) { std::string function = "Neutrals::initial_conditions"; 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,73 @@ 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); + didWork = restart_file(input.get_restartin_dir(), DoRead); - if (!DidWork) - std::cout << "Reading Restart for Neutrals Failed!!!\n"; + 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)) { + 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 << " 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); - } + 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 - } // type = Msis + if (icsType == "planet") { - if (ics["type"] == "Planet") { + didWork = true; // --------------------------------------------------------------------- // This section assumes we want a hydrostatic solution given the @@ -165,8 +171,11 @@ int Neutrals::initial_conditions(Grid grid, } // type = planet } + if (!didWork) + report.error("Issue with initial conditions!"); + report.exit(function); - return iErr; + 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 f23e0cfe..01ea5f86 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 diff --git a/src/read_input_file.cpp b/src/read_input_file.cpp index 262400a1..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"]; @@ -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); } } @@ -68,12 +69,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()); diff --git a/src/solver_gradients.cpp b/src/solver_gradients.cpp index 6fb57b95..f03f4228 100644 --- a/src/solver_gradients.cpp +++ b/src/solver_gradients.cpp @@ -153,3 +153,77 @@ arma_cube calc_gradient_alt_4th(arma_cube value, Grid grid) { return gradient; } +// -------------------------------------------------------------------------- +// Calculate the gradient in cubesphere spatial discretization +// -------------------------------------------------------------------------- +std::vector calc_gradient_cubesphere(arma_cube value, Grid grid) { + // Must be used for cubesphere (Probably need a boolean check) + int64_t nXs = grid.get_nY(); + int64_t nYs = grid.get_nX(); + int64_t nGCs = grid.get_nGCs(); + int64_t nAlts = grid.get_nAlts(); + + // Initialize two arma cubes for return + arma_cube grad_lon(nXs, nYs, nAlts); + arma_cube grad_lat(nXs, nYs, nAlts); + + for (int64_t iAlt = 0; iAlt < nAlts; iAlt++) { + /** Extract Grid Features **/ + // Addition: Get a copy of dx dy + arma_mat curr_refx = grid.refx_scgc.slice(iAlt); + arma_mat curr_refy = grid.refy_scgc.slice(iAlt); + + // Get some dx dy metrics from the grid + precision_t dx = grid.drefx(iAlt); + precision_t dy = grid.drefy(iAlt); + + // Get values of current level + arma_mat curr_value = value.slice(iAlt); + + /** Calculate Gradient with Central Difference Scheme **/ + // Since Reference grid is orthogonal, we only need gradient along reference xy direction + // Then we convert gradient from xy direction ot lat-lon direction + + 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); + } + } + } 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); + } + } + } + + // 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); + } + + // Not initializing with array like procedure in case I get bugs + std::vector gradient; + gradient.push_back(grad_lon); + gradient.push_back(grad_lat); + + return gradient; +} \ No newline at end of file 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 << " " diff --git a/src/tools.cpp b/src/tools.cpp index 841ca678..92b1fb38 100644 --- a/src/tools.cpp +++ b/src/tools.cpp @@ -112,6 +112,20 @@ 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 // ----------------------------------------------------------------------------- diff --git a/src/transform.cpp b/src/transform.cpp index 54fac144..d4a657e0 100644 --- a/src/transform.cpp +++ b/src/transform.cpp @@ -6,6 +6,18 @@ #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 the ascii characters from a string into an int c-array // - this is so we can pass the array to Fortran. @@ -16,12 +28,10 @@ int* copy_string_to_int(std::string inString) { const int length = inString.length(); // declaring character array int* outArray = new int[400]; - for (int i = 0; i < length; i++) { + for (int i = 0; i < length; i++) outArray[i] = inString[i]; - } - for (int i = length; i < 400; i++) { + for (int i = length; i < 400; i++) outArray[i] = 0; - } return outArray; } 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() 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 ...') diff --git a/tests/gradient_cubesphere/aether.json.gradient_test b/tests/gradient_cubesphere/aether.json.gradient_test new file mode 100644 index 00000000..8f45120a --- /dev/null +++ b/tests/gradient_cubesphere/aether.json.gradient_test @@ -0,0 +1,18 @@ + +{ + + "GeoBlockSize" : { + "nLons" : 50, + "nLats" : 50, + "nAlts" : 10}, + + "CubeSphere" : { + "is" : true}, + + "GeoGrid" : { + "dAlt" : 0.25, + "IsUniformAlt" : false}, + + "PlanetFile" : "UA/inputs/earth.in" + +} diff --git a/tests/gradient_cubesphere/run_all.sh b/tests/gradient_cubesphere/run_all.sh new file mode 100755 index 00000000..54334e33 --- /dev/null +++ b/tests/gradient_cubesphere/run_all.sh @@ -0,0 +1,16 @@ +#!/bin/sh + +# remove old directories: +rm -rf build run.* + +# make code with gradient unit test on: +mkdir build +cd build +cmake -DTEST_GRADIENT=Y ../../.. +make -j4 +cd .. +cp -R ../../share/run ./run.gradient_test +cd run.gradient_test + +cp ../aether.json.gradient_test ./aether.json +mpirun -np 6 ./aether 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 - + +