diff --git a/include/advance.h b/include/advance.h index 0f906b8c..1f287c4b 100644 --- a/include/advance.h +++ b/include/advance.h @@ -33,20 +33,20 @@ bool advance(Planets &planet, - Grid &gGrid, - Grid &mGrid, - Times &time, - Euv &euv, - Neutrals &neutrals, - Neutrals &neutralsMag, - Ions &ions, - Ions &ionsMag, - Chemistry &chemistry, - Chemistry &chemistryMag, - Electrodynamics &electrodynamics, - Electrodynamics &electrodynamicsMag, - Indices &indices, - Logfile &logfile, - Logfile &logfileMag); + Grid &gGrid, + Grid &mGrid, + Times &time, + Euv &euv, + Neutrals &neutrals, + Neutrals &neutralsMag, + Ions &ions, + Ions &ionsMag, + Chemistry &chemistry, + Chemistry &chemistryMag, + Electrodynamics &electrodynamics, + Electrodynamics &electrodynamicsMag, + Indices &indices, + Logfile &logfile, + Logfile &logfileMag); #endif // INCLUDE_ADVANCE_H_ diff --git a/include/aurora.h b/include/aurora.h index 75889cc9..85cb0873 100644 --- a/include/aurora.h +++ b/include/aurora.h @@ -12,16 +12,16 @@ **/ void read_aurora(Neutrals &neutrals, - Ions &ions); + Ions &ions); arma_vec calculate_fang(float eflux, // in ergs/cm2/s - float avee, // in keV - float Ebin, // eV - arma_vec rhoH, - std::vector Ci, - float dE, // eV - arma_vec H, - bool DoDebug); + float avee, // in keV + float Ebin, // eV + arma_vec rhoH, + std::vector Ci, + float dE, // eV + arma_vec H, + bool DoDebug); /********************************************************************** * brief Read in a file containing information about splitting ionization @@ -32,7 +32,7 @@ arma_vec calculate_fang(float eflux, // in ergs/cm2/s **/ void calc_aurora(Grid grid, - Neutrals &neutrals, - Ions &ions); + Neutrals &neutrals, + Ions &ions); #endif // INCLUDE_AURORA_H_ diff --git a/include/bfield.h b/include/bfield.h index df64fac8..a7535f2c 100644 --- a/include/bfield.h +++ b/include/bfield.h @@ -11,18 +11,18 @@ struct bfield_info_type { }; arma_vec get_magnetic_pole(int IsNorth, - Planets planet); + Planets planet); bfield_info_type get_bfield(precision_t lon, precision_t lat, precision_t alt, - bool DoDebug, + bool DoDebug, Planets planet); bfield_info_type get_dipole(precision_t lon, precision_t lat, precision_t alt, - bool DoDebug, + bool DoDebug, Planets planet); #endif // INCLUDE_BFIELD_H_ diff --git a/include/calc_euv.h b/include/calc_euv.h index 561930a2..83b3509c 100644 --- a/include/calc_euv.h +++ b/include/calc_euv.h @@ -21,16 +21,16 @@ // ------------------------------------------------------------------------- bool calc_euv(Planets planet, - Grid grid, - Times time, - Euv &euv, - Neutrals &neutrals, - Ions &ions, - Indices indices); + Grid grid, + Times time, + Euv &euv, + Neutrals &neutrals, + Ions &ions, + Indices indices); void calc_ionization_heating(Euv euv, - Neutrals &neutrals, - Ions &ions); + Neutrals &neutrals, + Ions &ions); #endif // INCLUDE_CALC_EUV_H_ diff --git a/include/calc_grid_derived.h b/include/calc_grid_derived.h index 28d5237f..98c22be0 100644 --- a/include/calc_grid_derived.h +++ b/include/calc_grid_derived.h @@ -7,7 +7,7 @@ #include // ---------------------------------------------------------------------------- -// +// // ---------------------------------------------------------------------------- std::vector calc_bin_edges(std::vector centers); @@ -19,8 +19,8 @@ arma_vec calc_bin_widths(arma_vec centers); // ---------------------------------------------------------------------------- // A helper function for mapping grids // ---------------------------------------------------------------------------- -bool grid_match(Grid gGrid, - Grid mGrid, +bool grid_match(Grid gGrid, + Grid mGrid, Quadtree gQuadtree, Quadtree mQuadtree); diff --git a/include/chemistry.h b/include/chemistry.h index 118d0189..eef015c5 100644 --- a/include/chemistry.h +++ b/include/chemistry.h @@ -72,7 +72,7 @@ class Chemistry { /// type of formula to use for reaction rate: int type; /// name of the reaction - std::string name; + std::string name; }; @@ -93,12 +93,12 @@ class Chemistry { Ions &ions); private: - bool search(std::string name, - json &headers, + bool search(std::string name, + json &headers, std::vector &error); - bool check_chemistry_file(json &headers, - std::vector> csv, + bool check_chemistry_file(json &headers, + std::vector> csv, Report &report); int read_chemistry_file(Neutrals neutrals, @@ -107,7 +107,7 @@ class Chemistry { reaction_type interpret_reaction_line(const Neutrals &neutrals, const Ions &ions, const std::vector &line, - const json &headers); + const json &headers); void find_species_id(const std::string &name, const Neutrals &neutrals, diff --git a/include/collisions.h b/include/collisions.h index 1a59a19c..98b1173c 100644 --- a/include/collisions.h +++ b/include/collisions.h @@ -11,6 +11,6 @@ #include "ions.h" void calc_ion_neutral_coll_freq(Neutrals &neutrals, - Ions &ions); + Ions &ions); #endif // INCLUDE_COLLISIONS_H_ diff --git a/include/constants.h b/include/constants.h index 474eb8f4..16183e9d 100644 --- a/include/constants.h +++ b/include/constants.h @@ -6,6 +6,18 @@ #include +// ------------------------------------------------------------------------- +// Define some constants for the code so that all functions understand +// stuff +// These are not physical constants, but are useful references +// ------------------------------------------------------------------------- + +const int iSphere_ = 1; +const int iCubesphere_ = 2; +const int iDipole_ = 3; +const std::string neutralType_ = "neuGrid"; +const std::string ionType_ = "ionGrid"; + // ------------------------------------------------------------------------- // Physical Constants // - Naming standards: @@ -65,7 +77,7 @@ const double cJULIAN2000 = 2451545.0; // ------------------------------------------------------------------------- const precision_t cPI = 3.141592653589793; -const precision_t cTWOPI = 2*cPI; +const precision_t cTWOPI = 2 * cPI; // ------------------------------------------------------------------------- // Conversion Constants: @@ -75,8 +87,8 @@ const precision_t cTWOPI = 2*cPI; // - Names are all UPPER CASE otherwise // ------------------------------------------------------------------------- -const precision_t cDtoR = cPI/180.0; -const precision_t cRtoD = 180.0/cPI; +const precision_t cDtoR = cPI / 180.0; +const precision_t cRtoD = 180.0 / cPI; // ------------------------------------------------------------------------- // converting time between seconds and other units of time: @@ -99,7 +111,7 @@ const double cMtoS = 60.0; const double cStoM = 1.0 / cMtoS; // MilliSeconds <-> Seconds: -const double cMStoS = 1.0/1000.0; +const double cMStoS = 1.0 / 1000.0; const double cStoMS = 1000.0; // ------------------------------------------------------------------------- diff --git a/include/electrodynamics.h b/include/electrodynamics.h index 13a631ab..40cd5555 100644 --- a/include/electrodynamics.h +++ b/include/electrodynamics.h @@ -54,7 +54,7 @@ class Electrodynamics { This does the following: - initialize all variables to missing values - - read in file if it exists + - read in file if it exists **/ Electrodynamics(Times time); @@ -66,12 +66,12 @@ class Electrodynamics { \param time need current time \param ions Going to set the potential and aurora **/ - + bool update(Planets planet, - Grid gGrid, - Times time, - Indices &indices, - Ions &ions); + Grid gGrid, + Times time, + Indices &indices, + Ions &ions); /************************************************************** @@ -85,7 +85,7 @@ class Electrodynamics { **/ bool check_times(double inputStartTime, double inputEndTime); - + /************************************************************** \brief used in advance.cpp to get potential, eflux, avee @@ -95,9 +95,9 @@ class Electrodynamics { **/ std::tuple get_electrodynamics(arma_cube magLat, - arma_cube magLocalTime); + arma_mat, + arma_mat> get_electrodynamics(arma_cube magLat, + arma_cube magLocalTime); /************************************************************** \brief Gets interpolation indices @@ -105,7 +105,7 @@ class Electrodynamics { Performs 2d interpolation over search vector to get indices \param vals the 2d array that needs indices - \param search The vector of values to interpolate over + \param search The vector of values to interpolate over **/ arma_mat get_interpolation_indices(arma_mat vals, arma_vec search); @@ -121,7 +121,7 @@ class Electrodynamics { \param time the time requested. **/ - + void set_time(double time); /************************************************************** @@ -195,7 +195,7 @@ class Electrodynamics { \param value Value to assign to Kp index **/ void set_kp(precision_t value); - + /************************************************************** \brief Get 2D electric potential on specified grid @@ -209,7 +209,7 @@ class Electrodynamics { with the potentials in the grids **/ arma_cube get_potential(arma_cube magLat, - arma_cube magLocalTime); + arma_cube magLocalTime); /************************************************************** \brief Get 2D electron energy flux on specified grid @@ -266,13 +266,13 @@ class Electrodynamics { with the ion avee in the grids **/ arma_mat get_ion_avee(); - + /********************************************************************** \brief Check to see if internal state of class is ok **/ - + bool is_ok(); - + private: /// This is the interpolation method for time: @@ -284,7 +284,7 @@ class Electrodynamics { /// Use the next value: const int iNext_ = 2; // Use the closest value: - const int iClosest_ = 3; + const int iClosest_ = 3; /// Interpolate: const int iInterp_ = 4; @@ -307,7 +307,7 @@ class Electrodynamics { /// A 2d array of magnetic local times needed. Can set interpolation /// coefficients in all of the grids when this is called: arma_mat mlts_needed; - + /// These are all indices that may be needed by sub-models: precision_t imf_bx_needed; precision_t imf_by_needed; @@ -340,11 +340,11 @@ class Electrodynamics { /// If we don't read in an electrodynamics file, then this should be /// set to an auroral model to use. Need to add model types. std::string auroral_model_to_use; - + /// Set the interpolation indices as a float. For each interpolation index, - /// the integer portion is the current index, and the decimal part is the + /// the integer portion is the current index, and the decimal part is the /// percentage of the distance between the current index and the next - /// index. For example, a distance midway between index 45 and 46 + /// index. For example, a distance midway between index 45 and 46 /// would give an interpolation index of 45.5. /// For time, we are assuming that all grids have the same times or that /// there are no overlaps in time, I think. @@ -378,22 +378,22 @@ class Electrodynamics { /// Potential at current time: arma_mat potential_current; - + /// Vector of 2d electron energy flux (in ergs/cm2/s): std::vector energy_flux; /// Said energy flux at the current time: arma_mat energy_flux_current; - + /// Vector of 2d electron average energy (in keV): std::vector average_energy; /// Average energy at current time: arma_mat average_energy_current; - + /// Vector of 2d ion energy flux (in ergs/cm2/s): std::vector ion_energy_flux; /// ion energy flux at current time: arma_mat ion_energy_flux_current; - + /// Vector of 2d ion average energy (in keV): std::vector ion_average_energy; /// ion average energy at current time: @@ -401,7 +401,7 @@ class Electrodynamics { /// Set to 1 if ion precipitation is included, else set to 0: int DoesIncludeIonPrecip; - + /// This sets the priority of the grid. The higher the number, the /// more important it is, so it should overwrite any regions of /// a lower priority grid. For example, you could have a global @@ -419,28 +419,28 @@ class Electrodynamics { /// is outside of the mlt range of the grid, then the /// interpolation index should be set to -1: arma_mat mlts_indices; - + }; - + /// As described above, a structure containing the grid-based /// values of electrodynamics as a function of time. This is /// vector, because we can have nested grids, or, in theory, the /// grid could change as a function of time. You can then search /// for the apropriate grid in space and time. std::vector input_electrodynamics; - + /// Because each grid has a priority, we need to go through them in /// priority order, this is the sorted indices list, so that /// grid_order[0] points to the input_electrodynamics with the /// lowest priority, grid_order[1] points to the 2nd lowest, etc. std::vector grid_order; - + /// Number of input grids for electrodynamics: int nElectrodynamicsGrids; - + /// An internal variable to hold the state of the class bool IsOk; - + /************************************************************** \brief Reads a netcdf file that has the electrodynamics specification @@ -469,13 +469,13 @@ class Electrodynamics { grids, so that the values are overwritten. To keep it "functional", we pass in the last round of values and those are moved into the output values and then the overlapping region is - overwritten (e.g., in the get_potential function, the + overwritten (e.g., in the get_potential function, the grids need to be cycled through calling get_values with the potential on that grid and the interpolation indices for the grid. \param values_current the pot/eflux/avee/etc from input_electrodynamics grid - + \param lats_indices the interpolation indices for the current grid latitudes @@ -483,8 +483,8 @@ class Electrodynamics { grid mlts \param values_old the output of this function for the last grid - **/ - + **/ + arma_mat get_values(arma_mat matToInterpolateOn, int rows, int cols); void set_all_indices_for_ie(Times time, Indices &indices); diff --git a/include/euv.h b/include/euv.h index 138e74f3..a9dd826d 100644 --- a/include/euv.h +++ b/include/euv.h @@ -8,7 +8,7 @@ * \class Euv * * \brief Defines the Extreme Ultraviolet radiation above the atmosphere - * + * * The Euv class defines the EUV environment above the atmosphere. It * does this through the use of a CSV file that contains a bunch of * information. Namely: @@ -18,7 +18,7 @@ * * \author Aaron Ridley * - * \date 2021/03/28 + * \date 2021/03/28 * **************************************************************/ @@ -27,12 +27,12 @@ class Euv { -public: + public: /// whether to actuall use euv at all: bool doUse; - - /// number of wavelengths in spectrum: + + /// number of wavelengths in spectrum: int nWavelengths; // number of lines in the EUV CSV file: @@ -59,10 +59,10 @@ class Euv { /// EUV Spectrum, lower wavelength of the bins: std::vector wavelengths_short; - + /// EUV Spectrum, upper wavelength of the bins: std::vector wavelengths_long; - + /// EUV Spectrum, energy of bin: std::vector wavelengths_energy; @@ -81,7 +81,7 @@ class Euv { std::vector solomon_hfg_c1; std::vector solomon_hfg_c2; std::vector solomon_hfg_fref; - + /// NEUVAC model linear coefficients (1-3): std::vector neuvac_s1; std::vector neuvac_s2; @@ -93,7 +93,7 @@ class Euv { /// NEUVAC model intercept: std::vector neuvac_int; - + // -------------------------------------------------------------------- // Functions: @@ -109,13 +109,13 @@ class Euv { **/ 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 - **/ - bool solomon_hfg(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 + **/ + 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) @@ -135,7 +135,7 @@ class Euv { Reads through each row in the EUV CSV file and figures out whether the row is abs, ion, diss, and then figures out which neutral it is - acting on and which neutral or ion results from the action + acting on and which neutral or ion results from the action (e.g. O + photon -> O+, identifies O as ionization "loss" and O+ as an ionization "source") @@ -147,16 +147,16 @@ class Euv { /********************************************************************** \brief Check to see if internal state of class is ok **/ - + bool is_ok(); - -private: + + private: /********************************************************************** \brief Read in the EUV CSV file Read in the EUV CSV file that describes all of the wavelengths and - cross sections (and any other EUV - related things that are a + cross sections (and any other EUV - related things that are a function of wavelength) **/ bool read_file(); @@ -172,8 +172,8 @@ class Euv { \return values The values in the CSV row that matches the item (and item2) **/ bool slot_euv(std::string item, - std::string item2, - std::vector &values); + std::string item2, + std::vector &values); /// An internal variable to hold the state of the class bool IsOk; diff --git a/include/external_msis.h b/include/external_msis.h index d6285456..764557c5 100644 --- a/include/external_msis.h +++ b/include/external_msis.h @@ -8,13 +8,13 @@ * \class Msis * * \brief create an interface to the msis model - * + * * MSIS is a neutral model of the atmosphere, written in - * fortran and provided by NRL. + * fortran and provided by NRL. * * \author Aaron Ridley * - * \date 2023/04/30 + * \date 2023/04/30 * **************************************************************/ @@ -27,20 +27,20 @@ class Msis { bool set_f107(precision_t f107in, precision_t f107ain); bool set_ap(precision_t apin); bool set_locations(arma_vec longitude, - arma_vec latitude, - arma_vec altitude); + arma_vec latitude, + arma_vec altitude); bool set_locations(arma_mat longitude, - arma_mat latitude, - arma_mat altitude); + arma_mat latitude, + arma_mat altitude); bool set_locations(arma_cube longitude, - arma_cube latitude, - arma_cube altitude); + arma_cube latitude, + arma_cube altitude); arma_vec get_vec(std::string value); arma_mat get_mat(std::string value); arma_cube get_cube(std::string value); bool is_valid_species(std::string value); bool is_ok(); - + private: int iYear, iDay; @@ -54,10 +54,10 @@ class Msis { arma_cube altKm; std::vector msis_results; - bool didChange = true; + bool didChange = true; json value_lookup; bool isCompiled; - + bool reset_interface_variable_sizes(); bool reset_results(); }; diff --git a/include/file_input.h b/include/file_input.h index 460a6fe6..ee127f7e 100644 --- a/include/file_input.h +++ b/include/file_input.h @@ -31,9 +31,9 @@ std::vector> read_csv(std::ifstream &file_ptr); \param csvLines a matrix of strings **/ json put_csv_in_json_w_name(std::vector> - csvLines); + csvLines); json put_csv_in_json_wo_name(std::vector> - csvLines); + csvLines); /************************************************************** @@ -46,8 +46,8 @@ std::vector> read_ssv(std::ifstream &file_ptr); /************************************************************** \brief Reads either a comma-separated time or series of lines describing time - format is either - y, m, d, h, m, s, ms + format is either + y, m, d, h, m, s, ms or y m diff --git a/include/grid.h b/include/grid.h index c4187111..bc1cd15b 100644 --- a/include/grid.h +++ b/include/grid.h @@ -48,11 +48,7 @@ struct cubesphere_chars { class Grid { public: - const int iSphere_ = 1; - const int iCubesphere_ = 2; - const int iDipole_ = 3; int iGridShape_ = -1; - // Armidillo Cube Versions: // Cell Center Coordinates arma_cube geoLon_scgc, geoX_scgc; @@ -147,6 +143,7 @@ class Grid { // Matrices whose elements denote the altitude index of the interiormost ghost cell // in the k-up and k-down direction (altitude for geo grids, q for dipole). arma_mat first_lower_gc, first_upper_gc; + precision_t altitude_lower_bc; // Whether to close field lines on dipole grid (Always false for geo grids) bool IsClosed; diff --git a/include/indices.h b/include/indices.h index 0c58b216..1e1b6ba0 100644 --- a/include/indices.h +++ b/include/indices.h @@ -20,7 +20,7 @@ * \author Aaron Ridley * - * \date 2021/04/16 + * \date 2021/04/16 **************************************************************/ #include @@ -34,10 +34,10 @@ struct index_file_output_struct { /// number of times read in: int64_t nTimes; - + /// array of times that correspond to the values: std::vector times; - + /// number of variables read in: int nVars; @@ -62,9 +62,9 @@ void print_index_file_output_struct(index_file_output_struct contents); class Indices { -// ----------------------------------------------------------------------- -// Public functions and variables -// ----------------------------------------------------------------------- + // ----------------------------------------------------------------------- + // Public functions and variables + // ----------------------------------------------------------------------- public: @@ -88,10 +88,10 @@ class Indices { /************************************************************** \brief a series of functions that return the internal index number - In order to keep track of which index is which, the class uses + In order to keep track of which index is which, the class uses constants. These functions return these constants. The user doesn't really need to know about the constants, but they have to get the - constant (when reading the file, for example) and then provide that + constant (when reading the file, for example) and then provide that to the set index function. Conversely, we could create a bunch of set_ functions (such as the set_f107 function below). We figured that this minor inconvience is easier than making a bunch of set_ @@ -112,13 +112,17 @@ class Indices { int get_au_index_id(); int get_al_index_id(); + json get_all_indices(double time); + bool restart_file(std::string dir, bool DoRead, double time); + + /************************************************************** \brief Return the indices index of the variable name \param name the name of the variable to find the index for **/ - + int lookup_index_id(std::string name); - + /************************************************************** \brief This function sets the f107, does an 81 day ave, sets f107a too \param f107_contents contents from the f107 file (time, f107, etc.) @@ -134,9 +138,9 @@ class Indices { \param missing value for missing data **/ bool set_index(int index_id, - std::vector time, - std::vector values, - precision_t missing); + std::vector time, + std::vector values, + precision_t missing); /************************************************************** \brief set the index array into the indices class @@ -146,10 +150,10 @@ class Indices { \param missing value for missing data **/ bool set_index(std::string index_name, - std::vector timearray, - std::vector indexarray, - precision_t missing); - + std::vector timearray, + std::vector indexarray, + precision_t missing); + /************************************************************** \brief Perturbs the indices requested by user input **/ @@ -164,30 +168,44 @@ class Indices { **/ void perturb_index(int iIndex, int seed, json style, bool DoReport); + /************************************************************** + \brief Re-Perturbs the specific indices based on old values and the new value + \param iIndex which index to perturb + \param unperturbedValue unperturbed value (value read in at start)) + \param perturbedValue value that the code has now + \param newValue value that the restart index file contains + **/ + + void reperturb_index(int iIndex, + precision_t unperturbedValue, + precision_t perturbedValue, + precision_t newValue); + + /************************************************************** \brief The general function that returns the index value at the time \param time the time in seconds that the index is requested at \param the index to return (i.e., one of the constants defined above) **/ - precision_t get_index(double time, int index); + precision_t get_index(double time, int index, bool useNonperturbed = false); /************************************************************** * \brief Get the name of the indices at the specified index * \param iIndex which index to get name * \return The string of name if the function succeeds, empty string if iIndex is out of range **/ - std::string get_name(int iIndex); + std::string get_name(int iIndex); /************************************************************** \brief Return the number of the indices vector **/ int all_indices_array_size(); -// ----------------------------------------------------------------------- -// Private functions and variables -// ----------------------------------------------------------------------- + // ----------------------------------------------------------------------- + // Private functions and variables + // ----------------------------------------------------------------------- -private: + private: /// structure containing information about the specific index: struct index_time_pair { @@ -197,12 +215,17 @@ class Indices { /// a vector of values for the index: std::vector values; + std::vector originals; /// a vector of times for the values: std::vector times; /// the name of the index as a string: std::string name; + + bool didPerturb; + bool isAddPerturb; + bool isConstantPerturb; }; /// the vector that contains all of the indices vectors: diff --git a/include/init_mag_grid.h b/include/init_mag_grid.h index b1af1f4c..415978a1 100644 --- a/include/init_mag_grid.h +++ b/include/init_mag_grid.h @@ -13,12 +13,12 @@ bool init_dipole_grid(Grid &mGrid, Planets planet); // Analytic solution to get from q,p dipole coords to r,theta // q coordinate along b-field line -// p l-shell +// p l-shell // return (r,theta) std::pair qp_to_r_theta(precision_t q, precision_t p); // convert mag to geographic -std::vector mag_to_geo(arma_cube magLon, +std::vector mag_to_geo(arma_cube magLon, arma_cube magLat, arma_cube magAlt, Planets planet); diff --git a/include/inputs.h b/include/inputs.h index ce273f76..7c48dbf1 100644 --- a/include/inputs.h +++ b/include/inputs.h @@ -607,12 +607,12 @@ class Inputs { **/ std::vector get_satellite_dts(); - /********************************************************************** - \brief returns settings[" - \param - **/ + /********************************************************************** + \brief returns settings[" + \param + **/ json get_tests(); - + // General get_setting functions with error checks: /********************************************************************** diff --git a/include/logfile.h b/include/logfile.h index 1fe19d5d..f413a209 100644 --- a/include/logfile.h +++ b/include/logfile.h @@ -5,7 +5,7 @@ #define INCLUDE_LOGFILE_H_ /************************************************************** - * + * * logfile.h: * * Write the logfile @@ -19,7 +19,7 @@ /** * The class Satellite is used to track the satellites * Given any time, the user can obtain the geographic location of the satellite - * + * * ASSUMPTION : The satellite csv layout is the same as the following * year mon day hr min sec lon lat alt x y z vx vy vz * (int) (int) (int) (int) (int) (int) (degree) (degree) (km) (km) (km) (km) (km/s) (km/s) (km/s) @@ -27,11 +27,11 @@ class Satellite { -public: + public: /** * \brief Initialize the satellite class - * The name of the satellite is not allowed to have any characters which can + * The name of the satellite is not allowed to have any characters which can * terminate the read of a string including white space' ', endline'\n', and '\t' * Different satellites must have different names (not only input file names) * \param csv_in The path to the satellite csv file @@ -69,7 +69,7 @@ class Satellite { // DEBUG void print(); -private: + private: // The name of the satellite std::string name; @@ -93,7 +93,7 @@ class Satellite { class Logfile { -public: + public: /** * \brief Initialize the Logfile. @@ -101,7 +101,7 @@ class Logfile { * every dt time. */ Logfile(Indices &indices, int64_t iLog); - + /** * \brief Close the file stream if not append */ @@ -117,7 +117,7 @@ class Logfile { Grid &gGrid, Times &time); -private: + private: // The name of logfile std::string logfileName; @@ -133,7 +133,7 @@ class Logfile { bool doAppend; // A randomly chosen point for test - std::vector lla {2,2,2}; + std::vector lla {2, 2, 2}; }; #endif // INCLUDE_LOGFILE_H_ diff --git a/include/neutrals.h b/include/neutrals.h index 34f67e34..e51b3452 100644 --- a/include/neutrals.h +++ b/include/neutrals.h @@ -225,6 +225,7 @@ class Neutrals { std::vector initial_altitudes; std::vector initial_temperatures; int64_t nInitial_temps = 0; + precision_t altitude_of_bc; /// Number of species to advect: int nSpeciesAdvect; diff --git a/include/output.h b/include/output.h index d6503984..2dc3d7cd 100644 --- a/include/output.h +++ b/include/output.h @@ -9,7 +9,7 @@ /************************************************************** * \class Output * \brief A containing to allow storage of variables for output - * + * * Writing output is a multi-step process now: * 1. Create a container to store the variables you want to output * 2. Define the variables to output within the container @@ -17,13 +17,13 @@ * 4. Write the output * * \author Aaron Ridley - * \date 2021/10/21 + * \date 2021/10/21 **************************************************************/ class OutputContainer { public: - + /********************************************************************** \brief initialize the output container **/ @@ -63,8 +63,8 @@ class OutputContainer { \param value the array of the data to output **/ void store_variable(std::string name, - std::string unit, - arma_cube value); + std::string unit, + arma_cube value); /********************************************************************** \brief store a variable to the list of variables to output @@ -74,9 +74,9 @@ class OutputContainer { \param value the array of the data to output **/ void store_variable(std::string name, - std::string long_name, - std::string unit, - arma_cube value); + std::string long_name, + std::string unit, + arma_cube value); /********************************************************************** \brief Get an arma_cube from the Container @@ -129,12 +129,12 @@ class OutputContainer { \brief write a file with the information in the container **/ bool write(); - + /********************************************************************** \brief write a json header file with the information in the container **/ bool write_container_header(); - + /********************************************************************** \brief write a binary file with the information in the container **/ @@ -149,27 +149,27 @@ class OutputContainer { \brief write a netcdf file with the information in the container **/ bool write_container_netcdf(); - + /********************************************************************** \brief read from a file an load into the container **/ bool read(); - + /********************************************************************** \brief display information contained in the container **/ void display(); - + /********************************************************************** \brief read a netcdf file - put the information in the container **/ bool read_container_netcdf(); - + /********************************************************************** - \brief clears the vector of variables + \brief clears the vector of variables **/ void clear_variables(); - + private: /// User can set the directory for output @@ -206,7 +206,7 @@ class OutputContainer { /// The frequency of the output for this particular container: float dt_output; - + /// This is to allow the user to select different output formats int output_type; @@ -214,7 +214,7 @@ class OutputContainer { const int binary_type = 0; const int netcdf_type = 1; const int hdf5_type = 2; - + }; /********************************************************************** @@ -238,12 +238,12 @@ class OutputContainer { **/ bool output(const Neutrals &neutrals, - const Ions &ions, - Grid &grid, - Times time, - const Planets &planet); + const Ions &ions, + Grid &grid, + Times time, + const Planets &planet); void output_binary_3d(std::ofstream &binary, - arma_cube value); + arma_cube value); #endif // INCLUDE_OUTPUT_H_ diff --git a/include/parallel.h b/include/parallel.h index 11372be3..1a363569 100644 --- a/include/parallel.h +++ b/include/parallel.h @@ -4,7 +4,7 @@ #ifndef INCLUDE_PARALLEL_H_ #define INCLUDE_PARALLEL_H_ -/// Need MPI (message passing interface) to do parallel stuff: +/// Need MPI (message passing interface) to do parallel stuff: #include "mpi.h" /// number of processors in whole simulation @@ -51,10 +51,10 @@ bool init_parallel(Quadtree &quadtree, Quadtree &quadtree_ion); **/ bool pack_border(const arma_cube &value, - precision_t *packed, - int64_t *iCounter, - int64_t nG, - int iDir); + precision_t *packed, + int64_t *iCounter, + int64_t nG, + int iDir); /********************************************************************** \brief Unpack variable buffer after message pass @@ -72,16 +72,16 @@ bool pack_border(const arma_cube &value, **/ bool unpack_border(arma_cube &value, - precision_t *packed, - int64_t *iCounter, - int64_t nG, - int iDir, - bool DoReverseX, - bool DoReverseY, - bool XbecomesY); + precision_t *packed, + int64_t *iCounter, + int64_t nG, + int iDir, + bool DoReverseX, + bool DoReverseY, + bool XbecomesY); /********************************************************************** - \brief initialize the grid variables to set up ghostcell message passing + \brief initialize the grid variables to set up ghostcell message passing \param grid the grid to set up message passing on \param nVarsToPass how many variables to pass **/ @@ -96,8 +96,8 @@ bool exchange_sides_init(Grid &grid, int64_t nVarsToPass); **/ bool exchange_one_var(Grid &grid, - arma_cube &var_to_pass, - bool doReverseSignAcrossPole); + arma_cube &var_to_pass, + bool doReverseSignAcrossPole); /********************************************************************** \brief test the exchange messages one var function diff --git a/include/planets.h b/include/planets.h index a512daac..25ab1aa8 100644 --- a/include/planets.h +++ b/include/planets.h @@ -15,11 +15,11 @@ class Planets { -// ----------------------------------------------------------------------- -// Public functions and variables -// ----------------------------------------------------------------------- + // ----------------------------------------------------------------------- + // Public functions and variables + // ----------------------------------------------------------------------- -public: + public: // -------------------------------------------------------------------- // Functions: @@ -110,7 +110,7 @@ class Planets { precision_t get_dipole_strength(); /********************************************************************** - \brief Returns omega (rotation rate) of the planet + \brief Returns omega (rotation rate) of the planet **/ precision_t get_omega(); @@ -118,27 +118,32 @@ class Planets { \brief returns neutrals json for neutral density BCs **/ json get_neutrals(); - + /********************************************************************** \brief returns neutral temperature json for temperature ICs **/ json get_temperatures(); - + /********************************************************************** \brief returns ions json for ion density characteristics **/ json get_ions(); - + + /********************************************************************** + \brief returns altitude of the density boundary condition: + **/ + precision_t get_altitude_of_bc(); + /********************************************************************** \brief Check to see if internal state of class is ok **/ - + bool is_ok(); -// ----------------------------------------------------------------------- -// Private functions and variables -// ----------------------------------------------------------------------- - + // ----------------------------------------------------------------------- + // Private functions and variables + // ----------------------------------------------------------------------- + private: /// A structure to describe the planetary characteristics for each planet @@ -250,10 +255,14 @@ class Planets { /// Information about the initial temperature of the planet json temperatures; - + /// Information about the ions of the planet json ions; + /// This is needed to specify at what altitude the densities + /// are specified: + precision_t altitude_of_bc; + /// An internal variable to hold the state of the class bool IsOk; diff --git a/include/quadtree.h b/include/quadtree.h index 18b584c7..d098efcc 100644 --- a/include/quadtree.h +++ b/include/quadtree.h @@ -9,14 +9,14 @@ /************************************************************** * \class Quadtree * - * \brief Defines the quadtree for blocks - * + * \brief Defines the quadtree for blocks + * * Aether is logically an i, j, k grid structure. Aether does domain * decomposition on the (i, j) coordinates and each processor works on * the full domain of the k dimension. (e.g., nn spherical * coordinates, i = longitude, j = latitude, and k = altitude.) The * quadtree takes the (i, j) dimensions and makes blocks out of them. - * + * * In the quadtree there are a number of root nodes, which are then * divided into 2 x 2 blocks. Each of those can then be subdivided * into 2 x 2 blocks. Each block resides on a separate processor. @@ -32,13 +32,13 @@ * * \author Aaron Ridley * - * \date 2022/07/05 + * \date 2022/07/05 * **************************************************************/ class Quadtree { -public: + public: /// number of blocks in each direction: const uint64_t nLR = 2; @@ -86,7 +86,7 @@ class Quadtree { /// Number of root nodes: int64_t nRootNodes; - + /// The quadtree root nodes: std::vector root_nodes; @@ -98,7 +98,7 @@ class Quadtree { arma_vec limit_high = {0.0, 0.0, 0.0}; /// For the given processor, the side that it is on: uint64_t iSide = -1; - + /********************************************************************** \brief Initializes the quadtree **/ @@ -107,7 +107,7 @@ class Quadtree { /********************************************************************** \brief Builds the quadtree **/ - void build(std::string gridtype); + void build(std::string gridtype); /********************************************************************** \brief Makes a new node on the quadtree, recursively @@ -119,18 +119,18 @@ class Quadtree { \param iSide basically the root node, or the side of the cubesphere **/ qtnode new_node(arma_vec lower_left_norm_in, - arma_vec size_right_norm_in, - arma_vec size_up_norm_in, - uint64_t &iProc_in_out, - uint64_t depth_in, - uint64_t iSide); + arma_vec size_right_norm_in, + arma_vec size_up_norm_in, + uint64_t &iProc_in_out, + uint64_t depth_in, + uint64_t iSide); /********************************************************************** \brief Get different vectors from the node \param node which node to get the vector from \param which defines the vector to get: LL = lower left; - SR = size in the right/left direction; + SR = size in the right/left direction; SU = size in the up/down direction; MID = mid point of the node; **/ @@ -154,26 +154,26 @@ class Quadtree { int64_t find_root(arma_vec point); /********************************************************************** - \brief If the point is outside of the normalized limits of the + \brief If the point is outside of the normalized limits of the quadtree, this tries to put the point back into the domain \param point the x, y, z normalized coordinate of the point. **/ arma_vec wrap_point_sphere(arma_vec point); /********************************************************************** - \brief If the point is outside of the normalized limits of the + \brief If the point is outside of the normalized limits of the quadtree, this tries to put the point back into the domain \param point the x, y, z normalized coordinate of the point. **/ arma_vec wrap_point_cubesphere(arma_vec point); - + /********************************************************************** \brief Check to see if internal state of class is ok **/ - + bool is_ok(); -private: + private: /// Defines whether the quadtree state is ok: bool IsOk = true; @@ -183,7 +183,7 @@ class Quadtree { bool IsCubeSphere = false; /// Defines whether the quadtree is a dipole or not: bool IsDipole = false; - + }; #endif // INCLUDE_QUADTREE_H_ diff --git a/include/read_collision_file.h b/include/read_collision_file.h index 64a9fb91..a11d1c85 100644 --- a/include/read_collision_file.h +++ b/include/read_collision_file.h @@ -7,27 +7,27 @@ #include "../include/aether.h" void read_collision_file(Neutrals &neutrals, - Ions &ions); + Ions &ions); void parse_nu_in_table(std::vector> csv, - Neutrals &neutrals, - Ions &ions); + Neutrals &neutrals, + Ions &ions); void parse_resonant_nu_in_table(std::vector> csv, - Neutrals &neutrals, - Ions &ions); + Neutrals &neutrals, + Ions &ions); void parse_bst_in_table(std::vector> csv, - Neutrals &neutrals, - Ions &ions); + Neutrals &neutrals, + Ions &ions); void parse_diffexp_in_table(std::vector> csv, - Neutrals &neutrals); + Neutrals &neutrals); void parse_diff0_in_table(std::vector> csv, Neutrals &neutrals); void check_collision_frequncies(Ions ions, - Neutrals neutrals); + Neutrals neutrals); #endif // INCLUDE_COLLISION_FILE_H_ diff --git a/include/read_indices_files.h b/include/read_indices_files.h index bd20e2f1..52bd6449 100644 --- a/include/read_indices_files.h +++ b/include/read_indices_files.h @@ -17,7 +17,7 @@ /********************************************************************** \brief Reads in all of the indices files and stores them in Indices - This function goes through all of the input indices files and + This function goes through all of the input indices files and reads in the files, then stores the values into the Indices class. At this point, it can read in the following file types: 1. NGDC F10.7 files. @@ -32,7 +32,7 @@ bool read_and_store_indices(Indices &indices); \param f107_file the f10.7 file to read in **/ index_file_output_struct read_f107_file(std::string f107_file, - Indices indices); + Indices indices); /********************************************************************** \brief Read the OMNIWeb file format and store in the index_file struct @@ -40,7 +40,7 @@ index_file_output_struct read_f107_file(std::string f107_file, \param indices needed to get the indices index for each variable **/ index_file_output_struct read_omni_file(std::string omni_file, - Indices indices); + Indices indices); /********************************************************************** \brief This code compares a string to return the variable index diff --git a/include/report.h b/include/report.h index e29696eb..156a434d 100644 --- a/include/report.h +++ b/include/report.h @@ -28,11 +28,11 @@ class Report { -// ----------------------------------------------------------------------- -// Public functions and variables -// ----------------------------------------------------------------------- + // ----------------------------------------------------------------------- + // Public functions and variables + // ----------------------------------------------------------------------- -public: + public: // Functions: @@ -146,9 +146,9 @@ class Report { \param cFunctionName **/ void student_checker_function_name(bool isStudent, - std::string cStudentName, - int iFunctionNumber, - std::string cFunctionName); + std::string cStudentName, + int iFunctionNumber, + std::string cFunctionName); /************************************************************** \brief Starts timer and reports when entering a function, if applicable @@ -179,10 +179,10 @@ class Report { **/ void times(); -// ----------------------------------------------------------------------- -// Private functions and variables -// ----------------------------------------------------------------------- -private: + // ----------------------------------------------------------------------- + // Private functions and variables + // ----------------------------------------------------------------------- + private: /// global verbose level of the code int iVerbose; diff --git a/include/sizes.h b/include/sizes.h index ff10b6b2..5a02b9cb 100644 --- a/include/sizes.h +++ b/include/sizes.h @@ -7,7 +7,7 @@ // This is the file that defines the number of grid points in each // direction. The entire code is based on these numbers, so you need // to recompile if you change these numbers. -// +// // These are temporary and will eventually be removed. // This is for the geographic grid: diff --git a/include/solvers.h b/include/solvers.h index d7d2412a..53385833 100644 --- a/include/solvers.h +++ b/include/solvers.h @@ -64,10 +64,6 @@ arma_vec limiter_value(arma_vec projected, arma_vec values, int64_t nPts, } -void advect(Grid &grid, - Times &time, - Neutrals &neutrals); - arma_vec solver_conduction( arma_vec value, arma_vec lambda, @@ -84,6 +80,11 @@ arma_cube solver_chemistry(arma_cube density, arma_cube loss, precision_t dt); +arma_mat solver_chemistry(arma_mat density, + arma_mat source, + arma_mat loss, + precision_t dt); + std::vector coriolis(std::vector velocity, precision_t rotation_rate, arma_cube lat_scgc); diff --git a/include/sphere.h b/include/sphere.h index 32de5fdf..b2483752 100644 --- a/include/sphere.h +++ b/include/sphere.h @@ -12,20 +12,20 @@ *************************************************/ namespace Sphere { - /// The normalized origins of each face of the cube (i.e. corner) - static const arma_mat ORIGINS = { - { 0.0, -0.5, 0.0} - }; +/// The normalized origins of each face of the cube (i.e. corner) +static const arma_mat ORIGINS = { + { 0.0, -0.5, 0.0} +}; - /// Normalized right steps in cube - static const arma_mat RIGHTS = { - {2.0, 0.0, 0.0} - }; +/// Normalized right steps in cube +static const arma_mat RIGHTS = { + {2.0, 0.0, 0.0} +}; - /// Normalized right steps in cube - static const arma_mat UPS = { - {0.0, 1.0, 0.0} - }; +/// Normalized right steps in cube +static const arma_mat UPS = { + {0.0, 1.0, 0.0} +}; }; @@ -34,31 +34,31 @@ namespace Sphere { *************************************************/ namespace Sphere4 { - /// The normalized origins of each node (i.e. corner) - static const arma_mat ORIGINS = { - { 0.0, -0.5, 0.0}, - { 0.0, -0.25, 0.0}, - { 0.0, 0.0, 0.0}, - { 0.0, 0.25, 0.0} - }; - - /// Normalized right steps in node - static const arma_mat RIGHTS = { - {2.0, 0.0, 0.0}, - {2.0, 0.0, 0.0}, - {2.0, 0.0, 0.0}, - {2.0, 0.0, 0.0} - }; - - /// Normalized up steps in node - static const arma_mat UPS = { - {0.0, 0.25, 0.0}, - {0.0, 0.25, 0.0}, - {0.0, 0.25, 0.0}, - {0.0, 0.25, 0.0} - }; - - }; +/// The normalized origins of each node (i.e. corner) +static const arma_mat ORIGINS = { + { 0.0, -0.5, 0.0}, + { 0.0, -0.25, 0.0}, + { 0.0, 0.0, 0.0}, + { 0.0, 0.25, 0.0} +}; + +/// Normalized right steps in node +static const arma_mat RIGHTS = { + {2.0, 0.0, 0.0}, + {2.0, 0.0, 0.0}, + {2.0, 0.0, 0.0}, + {2.0, 0.0, 0.0} +}; + +/// Normalized up steps in node +static const arma_mat UPS = { + {0.0, 0.25, 0.0}, + {0.0, 0.25, 0.0}, + {0.0, 0.25, 0.0}, + {0.0, 0.25, 0.0} +}; + +}; /************************************************* * \brief A namespace with all (6-root) sphere grid logic. @@ -67,32 +67,32 @@ namespace Sphere6 { /// The normalized origins of each face of the cube (i.e. corner) static const arma_mat ORIGINS = { - { 0.0, -0.5, 0.0}, - {2.0/3.0, -0.5, 0.0}, - {4.0/3.0, -0.5, 0.0}, - { 0.0, 0.0, 0.0}, - {2.0/3.0, 0.0, 0.0}, - {4.0/3.0, 0.0, 0.0} + {0.0, -0.5, 0.0}, + {1.0, -0.5, 0.0}, + {0.0, -0.5 + 1.0 / 3.0, 0.0}, + {1.0, -0.5 + 1.0 / 3.0, 0.0}, + {0.0, -0.5 + 2.0 / 3.0, 0.0}, + {1.0, -0.5 + 2.0 / 3.0, 0.0} }; /// Normalized right steps in cube static const arma_mat RIGHTS = { - { 2.0/3.0, 0.0, 0.0}, - { 2.0/3.0, 0.0, 0.0}, - { 2.0/3.0, 0.0, 0.0}, - { 2.0/3.0, 0.0, 0.0}, - { 2.0/3.0, 0.0, 0.0}, - { 2.0/3.0, 0.0, 0.0} + { 1.0, 0.0, 0.0}, + { 1.0, 0.0, 0.0}, + { 1.0, 0.0, 0.0}, + { 1.0, 0.0, 0.0}, + { 1.0, 0.0, 0.0}, + { 1.0, 0.0, 0.0} }; /// Normalized right steps in cube static const arma_mat UPS = { - { 0.0, 0.5, 0.0}, - { 0.0, 0.5, 0.0}, - { 0.0, 0.5, 0.0}, - { 0.0, 0.5, 0.0}, - { 0.0, 0.5, 0.0}, - { 0.0, 0.5, 0.0} + { 0.0, 1.0 / 3.0, 0.0}, + { 0.0, 1.0 / 3.0, 0.0}, + { 0.0, 1.0 / 3.0, 0.0}, + { 0.0, 1.0 / 3.0, 0.0}, + { 0.0, 1.0 / 3.0, 0.0}, + { 0.0, 1.0 / 3.0, 0.0} }; } // CubeSphere:: diff --git a/include/test.h b/include/test.h index d6f9b287..3b021a68 100644 --- a/include/test.h +++ b/include/test.h @@ -9,7 +9,8 @@ // Gradient tests // Cubesphere is not done nor tested -bool test_gradient(Planets planet, Quadtree quadtree, json test_config, Grid gGrid, Grid mGrid); +bool test_gradient(Planets planet, Quadtree quadtree, json test_config, + Grid gGrid, Grid mGrid); bool test_gradient_cubesphere(Planets planet, Quadtree quadtree, Grid grid); bool test_gradient_ijk(Planets planet, Grid grid, bool debug); diff --git a/include/times.h b/include/times.h index f082f31e..12b891e0 100644 --- a/include/times.h +++ b/include/times.h @@ -5,7 +5,7 @@ #define INCLUDE_TIMES_H_ /************************************************************** - * + * * times.h: * * Functions that are assocuated with keeping track of time in @@ -20,7 +20,7 @@ class Times { -public: + public: /************************************************************** \brief Initialize the Times class @@ -61,7 +61,7 @@ class Times { \brief Sets the start, restart, and current times. This sets the start time, restart time, and current time to - the input time, and initializes iStep and dt, then calls + the input time, and initializes iStep and dt, then calls increment_time, which derives a bunch of other variables. \param itime year, month, day, hour, minute, second, millisecond vector @@ -161,21 +161,21 @@ class Times { \brief Get the current time as an array **/ std::vector get_iCurrent(); - + /************************************************************** \brief Get the current simulation time (sec since start) **/ double get_simulation_time(); - + /********************************************************************** \brief Read / Write restart files for time \param dir directory to write restart files \param DoRead read the restart files if true, write if false **/ - bool restart_file(std::string dir, bool DoRead); + bool restart_file(std::string dir, bool DoRead); + + private: -private: - // ------------------------------------------------------------- // These variables are for keeping track of the time. All in seconds // since reference time (except where noted). @@ -212,7 +212,7 @@ class Times { /// Universal time in hours precision_t ut; - + /// in weird JPL units precision_t orbittime; @@ -227,13 +227,13 @@ class Times { /// This is day of year (and NOT real Julian Day!) int jDay; - + /// This is Julian day double julian_day; /// represented as YYMMDD std::string sYMD; - + /// represented as HHMMSS std::string sHMS; diff --git a/include/transform.h b/include/transform.h index fe27f452..91b275cf 100644 --- a/include/transform.h +++ b/include/transform.h @@ -16,15 +16,15 @@ std::string mkupper(std::string inString); void copy_cube_to_array(arma_cube cube_in, precision_t *array_out); void copy_mat_to_array(arma_mat mat_in, - precision_t *array_out, - bool isFortran); + precision_t *array_out, + bool isFortran); void copy_array_to_mat(precision_t *array_in, arma_mat &mat_out, bool isFortran); void copy_vector_to_array(std::vector vector_in, - int64_t nElements, - precision_t *array_out); + int64_t nElements, + precision_t *array_out); // This is needed when sending strings to Fortran. // We do this by copying the ascii numbers into an integer array, @@ -34,13 +34,18 @@ int* copy_string_to_int(std::string inString); arma_cube calc_magnitude(std::vector xyz); std::vector transform_llr_to_xyz_3d(std::vector llr); std::vector transform_xyz_to_llr_3d(std::vector xyz); -std::vector rotate_around_x_3d(std::vector XYZ_in, precision_t angle); -std::vector rotate_around_y_3d(std::vector XYZ_in, precision_t angle); -std::vector rotate_around_z_3d(std::vector XYZ_in, precision_t angle); +std::vector rotate_around_x_3d(std::vector XYZ_in, + precision_t angle); +std::vector rotate_around_y_3d(std::vector XYZ_in, + precision_t angle); +std::vector rotate_around_z_3d(std::vector XYZ_in, + precision_t angle); void transform_llr_to_xyz(precision_t llr_in[3], precision_t xyz_out[3]); -void transform_rot_z(precision_t xyz_in[3], precision_t angle_in, precision_t xyz_out[3]); -void transform_rot_y(precision_t xyz_in[3], precision_t angle_in, precision_t xyz_out[3]); +void transform_rot_z(precision_t xyz_in[3], precision_t angle_in, + precision_t xyz_out[3]); +void transform_rot_y(precision_t xyz_in[3], precision_t angle_in, + precision_t xyz_out[3]); void transform_float_vector_to_array(std::vector input, precision_t output[3]); diff --git a/share/run/UA/inputs/earth.in b/share/run/UA/inputs/earth.in index 6a968640..812097d7 100644 --- a/share/run/UA/inputs/earth.in +++ b/share/run/UA/inputs/earth.in @@ -25,6 +25,8 @@ BC is a density that is used in the lowest boundary cell if you In this example file, the values are from 96.87 km Jan 1, 2013 O_1D is made up. +#ALTITUDE_OF_BC +96.87 #NEUTRALS name, mass, vibration, thermal_cond, thermal_exp, advect, BC diff --git a/share/run/aether.json b/share/run/aether.json index b192b5b5..6325cc69 100644 --- a/share/run/aether.json +++ b/share/run/aether.json @@ -13,7 +13,7 @@ "EndTime" : [2011, 3, 20, 0, 1, 0], "neuGrid" : { - "Shape": "sphere4", + "Shape": "sphere", "nLonsPerBlock" : 20, "nLatsPerBlock" : 18, "nAlts" : 40, @@ -21,13 +21,12 @@ "IsUniformAlt" : false}, "ionGrid": { - "Shape": "dipole4", + "Shape": "sphere", "nLonsPerBlock": 24, - "nLatsPerBlock": 18, + "nLatsPerBlock": 22, "nAlts": 50, - "LatRange": [10, 80], - "AltRange": [80.0, 1000], - "LonRange": [0.0, 360.0]}, + "dAltScale" : 0.3, + "IsUniformAlt" : false}, "OmniwebFiles" : ["UA/inputs/omni_20110319.txt"], diff --git a/src/advance.cpp b/src/advance.cpp index f56e6aea..5d39bb6e 100644 --- a/src/advance.cpp +++ b/src/advance.cpp @@ -46,6 +46,8 @@ bool advance(Planets &planet, didWork = neutralsMag.check_for_nonfinites("Top of Advance - ion grid"); } + json dummy = indices.get_all_indices(time.get_current()); + gGrid.calc_sza(planet, time); mGrid.calc_sza(planet, time); @@ -258,13 +260,16 @@ bool advance(Planets &planet, if (time.check_time_gate(input.get_dt_write_restarts())) { report.print(3, "Writing restart files"); - neutrals.restart_file(input.get_restartout_dir(), gGrid.get_gridtype(), + neutrals.restart_file(input.get_restartout_dir(), + gGrid.get_gridtype(), DoWrite); - neutralsMag.restart_file(input.get_restartout_dir(), mGrid.get_gridtype(), + neutralsMag.restart_file(input.get_restartout_dir(), + mGrid.get_gridtype(), DoWrite); ions.restart_file(input.get_restartout_dir(), gGrid.get_gridtype(), DoWrite); ionsMag.restart_file(input.get_restartout_dir(), mGrid.get_gridtype(), DoWrite); time.restart_file(input.get_restartout_dir(), DoWrite); + indices.restart_file(input.get_restartout_dir(), DoWrite, time.get_current()); } } diff --git a/src/calc_dt.cpp b/src/calc_dt.cpp index 96658471..1ab838bd 100644 --- a/src/calc_dt.cpp +++ b/src/calc_dt.cpp @@ -16,7 +16,7 @@ precision_t calc_dt(Grid grid, std::vector cMax_vcgc) { precision_t dt; - if (grid.iGridShape_ == grid.iCubesphere_) + if (grid.iGridShape_ == iCubesphere_) dt = calc_dt_cubesphere(grid, cMax_vcgc); else dt = calc_dt_sphere(grid, cMax_vcgc); @@ -110,13 +110,6 @@ precision_t calc_dt_cubesphere(Grid grid, std::vector cMax_vcgc) { //dty.slice(iAlt) = grid.drefy(iAlt) * dummy_1 / u2; } - //if (iProc == 0) - // display_cube("dtx : ", dtx); - - //if (iProc == 0) - // display_cube("dty : ", dty); - - // Take minimum dts in each direction: dta(0) = dtx.min(); dta(1) = dty.min(); diff --git a/src/calc_electron_temperature.cpp b/src/calc_electron_temperature.cpp index a5a00552..5b459982 100644 --- a/src/calc_electron_temperature.cpp +++ b/src/calc_electron_temperature.cpp @@ -616,7 +616,7 @@ arma_mat Ions::calc_thermoelectric_current(Grid &grid) { // with the dipole, the field-aligned current is in the k^ direction // But we do not solve for e- velocity (and exb is 0 parallel to B), so we cannot do this: - // if (grid.iGridShape_ == grid.iDipole_){ + // if (grid.iGridShape_ == iDipole_){ // for (int64_t iAlt = 0; iAlt < ions.density_scgc.n_slices; iAlt++){ // JParallel += (ions.density_scgc.slice(iAlt) * cE % (ions.velocity_vcgc[2].slice(iAlt) - ions.exb_vcgc[2].slice(iAlt))) // * grid.dalt_center_scgc[iAlt]; diff --git a/src/calc_ion_drift.cpp b/src/calc_ion_drift.cpp index 0cc9821a..19970f1a 100644 --- a/src/calc_ion_drift.cpp +++ b/src/calc_ion_drift.cpp @@ -26,11 +26,11 @@ void Ions::calc_efield(Grid grid) { // -------------------------------------------------------------------------- void Ions::calc_exb_drift(Grid grid) { - + std::string function = "Ions::calc_exb"; static int iFunction = -1; report.enter(function, iFunction); - + arma_cube bmag2 = (grid.bfield_mag_scgc) % (grid.bfield_mag_scgc); exb_vcgc = cross_product(efield_vcgc, grid.bfield_vcgc); diff --git a/src/calc_ion_temperature.cpp b/src/calc_ion_temperature.cpp index 42c5f5b1..72f6df5b 100644 --- a/src/calc_ion_temperature.cpp +++ b/src/calc_ion_temperature.cpp @@ -27,18 +27,17 @@ void Ions::init_ion_temperature(Neutrals neutrals, Grid grid) { temperature_scgc = neutrals.temperature_scgc; - // For electron temperature, we need to check if some species are present or not. + // For electron temperature, we need to check if some species are present or not. // Do this check now & warn if needed: if ((neutrals.get_species_id("O") == -1) - || (neutrals.get_species_id("O2") == -1) - ||(neutrals.get_species_id("N2") == -1)){ - if (input.get_do_photoelectron_heating() + || (neutrals.get_species_id("O2") == -1) + || (neutrals.get_species_id("N2") == -1)) { + if (input.get_do_photoelectron_heating() || input.get_do_ionization_heating() - || input.get_do_electron_neutral_elastic_collisional_heating()) { - report.error("Your electron temperature sources require neutral O, O2, and N2 to be present."); - } + || input.get_do_electron_neutral_elastic_collisional_heating()) + report.error("Your electron temperature sources require neutral O, O2, and N2 to be present."); } - + return; } diff --git a/src/calc_neutral_derived.cpp b/src/calc_neutral_derived.cpp index d6802d93..d4959b8b 100644 --- a/src/calc_neutral_derived.cpp +++ b/src/calc_neutral_derived.cpp @@ -407,7 +407,7 @@ precision_t Neutrals::calc_dt(Grid grid) { precision_t dt; - if (grid.iGridShape_ == grid.iCubesphere_) + if (grid.iGridShape_ == iCubesphere_) dt = calc_dt_cubesphere(grid); else { int iDir; diff --git a/src/chemistry.cpp b/src/chemistry.cpp index a7713c2f..9eb2ebf7 100644 --- a/src/chemistry.cpp +++ b/src/chemistry.cpp @@ -25,7 +25,8 @@ Chemistry::Chemistry(Neutrals neutrals, std::string function = "Chemistry::Chemistry"; //record current function static int iFunction = -1; //usually -1 for report function - report.enter(function, iFunction); //keeps track of functions for: verbose levels, etc. + report.enter(function, + iFunction); //keeps track of functions for: verbose levels, etc. if (read_chemistry_file(neutrals, ions) > 0) { //searching for valid chem file report.print(0, "Could not read chemistry file!"); @@ -529,10 +530,11 @@ int Chemistry::read_chemistry_file(Neutrals neutrals, // Interpret a comma separated line of the chemical reaction file // ----------------------------------------------------------------------------- -Chemistry::reaction_type Chemistry::interpret_reaction_line(const Neutrals &neutrals, - const Ions &ions, - const std::vector &line, - const json &headers) { +Chemistry::reaction_type Chemistry::interpret_reaction_line( + const Neutrals &neutrals, + const Ions &ions, + const std::vector &line, + const json &headers) { std::string function = "Chemistry::interpret_reaction_line"; static int iFunction = -1; @@ -647,13 +649,15 @@ void Chemistry::find_species_id(const std::string &name, int iSpecies; IsNeutral = false; - id_ = neutrals.get_species_id(name); //from earth.in, starts at 0 w/ first species under "#NEUTRALS",(neutrals.cpp) + id_ = neutrals.get_species_id( + name); //from earth.in, starts at 0 w/ first species under "#NEUTRALS",(neutrals.cpp) if (id_ > -1) IsNeutral = true; else - id_ = ions.get_species_id(name);//from earth.in, starts at 0 w/ first species under "#IONS",(ions.cpp) + id_ = ions.get_species_id( + name);//from earth.in, starts at 0 w/ first species under "#IONS",(ions.cpp) report.exit(function); return; @@ -671,23 +675,23 @@ void Chemistry::display_reaction(Chemistry::reaction_type reaction) { std::cout << "Number of Sources : " << reaction.nSources << "\n"; for (i = 0; i < reaction.nLosses; i++) // First line for reaction - if (i < reaction.nLosses - 1) {// + if (i < reaction.nLosses - 1) // std::cout << reaction.losses_names[i] << " + "; - } else {// + + else // std::cout << reaction.losses_names[i] << " -> "; - } for (i = 0; i < reaction.nSources; i++) - if (i < reaction.nSources - 1) {// + if (i < reaction.nSources - 1) // std::cout << reaction.sources_names[i] << " + "; - } else {// + + else // std::cout << reaction.sources_names[i] << " (RR : " << reaction.rate << ")\n"; - } for (i = 0; i < reaction.nLosses; i++)//Second line for reaction if (i < reaction.nLosses - 1) {// std::cout << reaction.losses_ids[i] - << "(" << reaction.losses_IsNeutral[i] << ")" << " + "; + << "(" << reaction.losses_IsNeutral[i] << ")" << " + "; } else {// std::cout << reaction.losses_ids[i] << "(" << reaction.losses_IsNeutral[i] << ")" << " -> "; @@ -696,11 +700,11 @@ void Chemistry::display_reaction(Chemistry::reaction_type reaction) { for (i = 0; i < reaction.nSources; i++) if (i < reaction.nSources - 1) {// std::cout << reaction.sources_ids[i] - << "(" << reaction.sources_IsNeutral[i] << ")" << " + "; + << "(" << reaction.sources_IsNeutral[i] << ")" << " + "; } else {// std::cout << reaction.sources_ids[i] - << "(" << reaction.sources_IsNeutral[i] - << ")" << " (RR : " << reaction.rate << ")\n"; + << "(" << reaction.sources_IsNeutral[i] + << ")" << " (RR : " << reaction.rate << ")\n"; } diff --git a/src/file_input.cpp b/src/file_input.cpp index 29e3d67a..667f8e03 100644 --- a/src/file_input.cpp +++ b/src/file_input.cpp @@ -170,7 +170,7 @@ precision_t read_float(std::ifstream &file_ptr, std::string hash) { line = strip_string_end(line); try { - output = stoi(line); + output = stof(line); } catch (...) { std::cout << "Issue in read_float!\n"; std::cout << "In hash: "; diff --git a/src/fill_grid.cpp b/src/fill_grid.cpp index e98e0603..3006b827 100644 --- a/src/fill_grid.cpp +++ b/src/fill_grid.cpp @@ -173,7 +173,7 @@ void Grid::fill_grid_bfield(Planets planet) { // Magnetic coordinates: // init_mag grid already initializes magLon & magInvLat - if (IsGeoGrid) { + if (iGridShape_ != iDipole_) { magInvLat_scgc(iLon, iLat, iAlt) = bfield_info.lat; magLon_scgc(iLon, iLat, iAlt) = bfield_info.lon; } @@ -192,7 +192,7 @@ void Grid::fill_grid_bfield(Planets planet) { // Now we modify the dipole's magnetic field to account for any imprecision. // Take the bfield_mag and put it into the third component (b-hat = k-hat) - if (IsDipole) { + if (iGridShape_ == iDipole_) { bfield_vcgc[2] = bfield_mag_scgc % sign(magInvLat_scgc * -1.0); bfield_vcgc[1].zeros(); bfield_vcgc[0].zeros(); @@ -228,6 +228,7 @@ void Grid::fill_grid_radius(Planets planet) { // This generalizes things so that radius could be a function of all // three dimensions. The Cubesphere has different latitudes in the first // and second dimensions. + for (iLon = 0; iLon < nLons; iLon++) for (iLat = 0; iLat < nLats; iLat++) for (iAlt = 0; iAlt < nAlts; iAlt++) diff --git a/src/grid.cpp b/src/grid.cpp index cbe9f409..974acf44 100644 --- a/src/grid.cpp +++ b/src/grid.cpp @@ -32,11 +32,15 @@ Grid::Grid(std::string gridtype) { if (iGridShape_ == iCubesphere_) { if (grid_input.nX > grid_input.nY) { report.error("Cubesphere grid: nX > nY, reducing nX"); + report.print(0, gridType + + ": Cubesphere selected, but nX /= nY, reducing nX"); grid_input.nX = grid_input.nY; } if (grid_input.nY > grid_input.nX) { report.error("Cubesphere grid: nY > nX, reducing nY"); + report.print(0, gridType + + ": Cubesphere selected, but nX /= nY, reducing nY"); grid_input.nY = grid_input.nX; } } @@ -292,6 +296,7 @@ Grid::Grid(std::string gridtype) { UseThisCell.fill(true); first_lower_gc.set_size(nX, nY); first_upper_gc.set_size(nX, nY); + altitude_lower_bc = 0.0; cent_acc_vcgc = make_cube_vector(nLons, nLats, nAlts, 3); diff --git a/src/grid_match.cpp b/src/grid_match.cpp index 1df5ce97..2497a13a 100644 --- a/src/grid_match.cpp +++ b/src/grid_match.cpp @@ -25,7 +25,7 @@ bool grid_match(Grid gGrid, lon = mGrid.geoLon_scgc(iX, iY, iZ); lat = mGrid.geoLat_scgc(iX, iY, iZ); - if (gGrid.iGridShape_ == gGrid.iSphere_) { + if (gGrid.iGridShape_ == iSphere_) { norms(0) = lon / cPI; norms(1) = lat / cPI; norms(2) = 0.0; @@ -34,6 +34,7 @@ bool grid_match(Grid gGrid, norms = sphere_to_cube(lon, lat); iNode = gQuadtree.find_point(norms); } + if (report.test_verbose(6)) std::cout << "lon, lat, node: " << lon*cRtoD << " " << lat*cRtoD << " " diff --git a/src/grid_spacing.cpp b/src/grid_spacing.cpp index ad6d5b45..f1a55c92 100644 --- a/src/grid_spacing.cpp +++ b/src/grid_spacing.cpp @@ -142,12 +142,13 @@ void Grid::calc_k_grid_spacing() { radius_scgc.slice(iZ) - radius_scgc.slice(iZ - 1); // For the sphere & cubesphere, k is in meters: - if (iGridShape_ == iSphere_ || iGridShape_ == iCubesphere_){ + if (iGridShape_ == iSphere_ || iGridShape_ == iCubesphere_) { dk_center_m_scgc = dk_center_scgc; dk_edge_m = dk_edge; } + // This needs to be turned into a distance for the dipole: - if (iGridShape_ == iDipole_){ + if (iGridShape_ == iDipole_) { // the dk's may be negative (not allowed). make sure they are positive // this gets rid of SO many errors... dk_center_scgc = abs(dk_center_scgc); @@ -243,14 +244,15 @@ void Grid::calc_i_grid_spacing() { di_center_m_scgc = di_center_scgc % radius_scgc; di_edge_m = di_edge % radius_scgc; - // If the shape is a sphere or dipole, then the first coordinate is longitude. - // The physical distance needs to be changed by the cos of the latitude, + // If the shape is a sphere or dipole, then the first coordinate is longitude. + // The physical distance needs to be changed by the cos of the latitude, // which is the j coordinate in the sphere (different for dipole). if (iGridShape_ == iSphere_) { di_center_m_scgc = di_center_m_scgc % abs(cos(j_center_scgc)); // edge is in-line with the j center di_edge_m = di_edge_m % abs(cos(j_center_scgc)); } + // Dipole will use cos(magLat) if (iGridShape_ == iDipole_) { di_center_m_scgc = di_center_m_scgc % abs(cos(magLat_scgc)); @@ -343,8 +345,10 @@ void Grid::calc_j_grid_spacing() { // Dipole will have different scaling... if (iGridShape_ == iDipole_) { - dj_center_m_scgc = radius_scgc % dj_center_scgc % pow(cos(magLat_scgc), 3) / delTheta(magLat_scgc) % sign(magLat_scgc); - dj_edge_m = radius_scgc % dj_edge % pow(cos(magLat_scgc), 3) / delTheta(magLat_scgc) % sign(magLat_scgc); + dj_center_m_scgc = radius_scgc % dj_center_scgc % pow(cos(magLat_scgc), + 3) / delTheta(magLat_scgc) % sign(magLat_scgc); + dj_edge_m = radius_scgc % dj_edge % pow(cos(magLat_scgc), + 3) / delTheta(magLat_scgc) % sign(magLat_scgc); } // For a stretched grid, calculate some useful quantities: diff --git a/src/grid_sphere.cpp b/src/grid_sphere.cpp index 8f940ff9..d3bb2120 100644 --- a/src/grid_sphere.cpp +++ b/src/grid_sphere.cpp @@ -124,6 +124,11 @@ void Grid::create_sphere_grid(Quadtree quadtree) { for (iLon = 0; iLon < nLons; iLon++) lon1d(iLon) = lon0 + (iLon - nGCs + 0.5) * dlon; + if (report.test_verbose(1)) { + std::cout << function << ": " << lon0 << " " << dlon << "\n"; + display_vector("in function " + function + " lon1d : ", lon1d * cRtoD); + } + for (iLat = 0; iLat < nLats; iLat++) { for (iAlt = 0; iAlt < nAlts; iAlt++) { geoLon_scgc.subcube(0, iLat, iAlt, nLons - 1, iLat, iAlt) = lon1d; @@ -146,6 +151,12 @@ void Grid::create_sphere_grid(Quadtree quadtree) { for (iLat = 0; iLat < nLats; iLat++) lat1d(iLat) = lat0 + (iLat - nGCs + 0.5) * dlat; + if (report.test_verbose(1)) { + std::cout << function << ": " << lat0 << " " << dlat << "\n"; + + display_vector("in function " + function + " lat1d : ", lat1d * cRtoD); + } + for (iLon = 0; iLon < nLons; iLon++) { for (iAlt = 0; iAlt < nAlts; iAlt++) { geoLat_scgc.subcube(iLon, 0, iAlt, iLon, nLats - 1, iAlt) = lat1d; diff --git a/src/indices.cpp b/src/indices.cpp index 4dedf81d..b26e2b5a 100644 --- a/src/indices.cpp +++ b/src/indices.cpp @@ -20,6 +20,9 @@ Indices::Indices() { index_time_pair single_index; single_index.nValues = 0; single_index.name = ""; + single_index.didPerturb = false; + single_index.isAddPerturb = false; + single_index.isConstantPerturb = false; std::string lookup_file = input.get_indices_lookup_file(); indices_lookup = read_json(lookup_file); @@ -143,7 +146,7 @@ bool read_and_store_indices(Indices &indices) { bool Indices::perturb() { bool DidWork = true; bool DoReport = false; - int64_t iDebug = 2; + int64_t iDebug = 0; json perturb_values = input.get_perturb_values(); @@ -152,7 +155,7 @@ bool Indices::perturb() { for (auto it = perturb_values.begin(); it != perturb_values.end(); ++it) { std::string name = it.key(); - if (name != "Chemistry") { + if (name != "Chemistry" && name != "restart_control") { if (report.test_verbose(iDebug)) { std::cout << "Perturbing Index : " << name << "\n"; @@ -181,6 +184,76 @@ bool Indices::perturb() { // Perturb a specific index in the way the user requested // ---------------------------------------------------------------------- +void Indices::reperturb_index(int iIndex, + precision_t unperturbedValue, + precision_t perturbedValue, + precision_t newValue) { + + int64_t nValues = all_indices_arrays[iIndex].nValues; + + if (all_indices_arrays[iIndex].didPerturb && + all_indices_arrays[iIndex].isConstantPerturb) { + precision_t perturb; + + if (all_indices_arrays[iIndex].isAddPerturb) { + // constant, non-normalized value: + perturb = newValue - unperturbedValue; + + if (iGrid == 0) + std::cout << " -> New Added Perturb : " << perturb << "\n"; + + for (int64_t iValue = 0; iValue < nValues; iValue++) { + all_indices_arrays[iIndex].values[iValue] = + all_indices_arrays[iIndex].originals[iValue] + perturb; + } + + } else { + // constant, normalized value: + perturb = newValue / unperturbedValue; + + if (iGrid == 0) + std::cout << " -> New (normalized) Multiplied Perturb : " + << perturb << "\n"; + + for (int64_t iValue = 0; iValue < nValues; iValue++) { + + all_indices_arrays[iIndex].values[iValue] = + all_indices_arrays[iIndex].originals[iValue] * perturb; + } + } + + } else { + std::string mess = "Reperturb index: don't know how to "; + mess = mess + "handle non-perturb or nonconstant perturb"; + report.error(mess); + } + +} + +// ---------------------------------------------------------------------- +// Perturb a specific index in the way the user requested +/* +The way this code works is that you can perturb things in different ways. +Multiply by a constant value: + - if the mean is 1.0, then the perturbed value will be unbiased + - if the mean is above or below 1.0, it will be biased. + - the standard deviation is normalized to 1, so it is a percentage + of the value. + - This will come up with a value that you multiply all of the values by, + like 0.843 or 1.203. +Multiply by a non-constant value: + - same as above, but each value will have a different random number + instead of a single (constant) value +Add a constant value: + - the mean and standard deviation are NOT normalized. + - a single value then derived given the mean and the standard dev. + - an unbiased value would have a mean = 0 +Add a non-constant value: + - same as above, but each value will have a different random number + instead of a single (constant) value +*/ +// ---------------------------------------------------------------------- + void Indices::perturb_index(int iIndex, int seed, json style, bool DoReport) { @@ -191,6 +264,8 @@ void Indices::perturb_index(int iIndex, int seed, bool add = true; bool constant = false; + all_indices_arrays[iIndex].didPerturb = true; + if (style.contains("Mean")) mean = style["Mean"]; @@ -200,15 +275,21 @@ void Indices::perturb_index(int iIndex, int seed, std = standard_deviation(all_indices_arrays[iIndex].values); // Add or Multiply the random values - if (style.contains("Add")) + if (style.contains("Add")) { add = style["Add"]; + if (add) + all_indices_arrays[iIndex].isAddPerturb = true; + } + // Only one value for all elements or individual values for elements if (style.contains("Constant")) constant = style["Constant"]; - if (constant) + if (constant) { nV = 1; + all_indices_arrays[iIndex].isConstantPerturb = true; + } std::vector perturbations = get_normal_random_vect(mean, std, @@ -220,6 +301,9 @@ void Indices::perturb_index(int iIndex, int seed, if (!constant) iV = iValue; + all_indices_arrays[iIndex].originals.push_back( + all_indices_arrays[iIndex].values[iValue]); + if (add) { if (DoReport && iValue == 0) std::cout << " ==> Adding " << perturbations[iV] << "\n"; @@ -316,7 +400,8 @@ precision_t Indices:: get_f107a(double time) { // This is the general function for getting an index // ---------------------------------------------------------------------- -precision_t Indices::get_index(double time, int index) { +precision_t Indices::get_index(double time, int index, + bool useNonperturbed /* = false */) { int64_t iLow, iMid, iHigh; @@ -353,8 +438,14 @@ precision_t Indices::get_index(double time, int index) { all_indices_arrays[index].times[iMid]); precision_t x = (time - all_indices_arrays[index].times[iMid]) / dt; - precision_t value = (1.0 - x) * all_indices_arrays[index].values[iMid] + - x * all_indices_arrays[index].values[iMid + 1]; + precision_t value; + + if (useNonperturbed) + value = (1.0 - x) * all_indices_arrays[index].originals[iMid] + + x * all_indices_arrays[index].originals[iMid + 1]; + else + value = (1.0 - x) * all_indices_arrays[index].values[iMid] + + x * all_indices_arrays[index].values[iMid + 1]; return value; } @@ -417,6 +508,81 @@ bool Indices::set_index(int index, return DidWork; } +json Indices::get_all_indices(double time) { + json outputJson; + + int64_t iIndex; + precision_t value; + + for (iIndex = 0; iIndex < nIndices; iIndex++) { + if (all_indices_arrays[iIndex].nValues > 0) { + value = get_index(time, iIndex); + outputJson[all_indices_arrays[iIndex].name] = value; + } + } + + return outputJson; +} + +// ----------------------------------------------------------------------------- +// This is for restarting the code. Either write or read the time. +// ----------------------------------------------------------------------------- + +bool Indices::restart_file(std::string dir, bool DoRead, double time) { + + std::string filename; + bool DidWork = true; + filename = dir + "/indices_" + cMember + ".json"; + + json restart_indices_json, original_indices_json; + + if (DoRead) { + restart_indices_json = read_json(filename); + + if (report.test_verbose(1)) { + std::cout << "Restarted indices, Current time : "; + std::cout << std::setw(2) << restart_indices_json << "\n"; + } + + original_indices_json = get_all_indices(time); + precision_t orig, rest, unperturbed; + + for (auto it = original_indices_json.begin(); + it != original_indices_json.end(); ++it) { + std::string name = it.key(); + orig = original_indices_json[name]; + rest = restart_indices_json[name]; + + if (abs(orig - rest) > cSmall) { + int iIndex = lookup_index_id(name); + unperturbed = get_index(time, iIndex, true); + + if (iGrid == 0) + std::cout << " -> Index was altered during restart : " + << name << " -> index number: " + << iIndex << "; orig, rest, un " + << orig << " " + << rest << " " + << unperturbed << " " + << all_indices_arrays[iIndex].isAddPerturb << "\n"; + + reperturb_index(iIndex, unperturbed, orig, rest); + + } + } + + + } else { + restart_indices_json = get_all_indices(time); + + if (iGrid == 0) + DidWork = write_json(filename, restart_indices_json); + } + + return DidWork; +} + + // ---------------------------------------------------------------------- // Dump the contents of an index_file_output_struct // ---------------------------------------------------------------------- diff --git a/src/init_geo_grid.cpp b/src/init_geo_grid.cpp index ed333537..ddb19d2c 100644 --- a/src/init_geo_grid.cpp +++ b/src/init_geo_grid.cpp @@ -20,7 +20,9 @@ void Grid::create_altitudes(Planets planet) { arma_vec alt1d(nAlts); - Inputs::grid_input_struct grid_input = input.get_grid_inputs("neuGrid"); + Inputs::grid_input_struct grid_input; + + grid_input = input.get_grid_inputs(gridType); if (grid_input.IsUniformAlt) { for (iAlt = 0; iAlt < nAlts; iAlt++) @@ -186,17 +188,17 @@ bool Grid::init_geo_grid(Quadtree quadtree, report.enter(function, iFunction); bool DidWork = true; - IsGeoGrid = 1; + IsGeoGrid = true; if (iGridShape_ == iCubesphere_) { - report.print(0, "Creating Cubesphere Grid"); + report.print(0, "Creating Cubesphere Grid for : " + gridType); if (!Is0D & !Is1Dz) create_cubesphere_connection(quadtree); IsCubeSphereGrid = true; } else { - report.print(0, "Creating Spherical Grid"); + report.print(0, "Creating Spherical Grid for : " + gridType); if (!Is0D & !Is1Dz) create_sphere_connection(quadtree); @@ -217,6 +219,9 @@ bool Grid::init_geo_grid(Quadtree quadtree, //MPI_Barrier(aether_comm); create_altitudes(planet); + // set the altitude of the lower boundary values: + altitude_lower_bc = planet.get_altitude_of_bc(); + init_connection(); //DidWork = write_restart(input.get_restartout_dir()); @@ -230,14 +235,14 @@ bool Grid::init_geo_grid(Quadtree quadtree, if (iGridShape_ == iCubesphere_) { correct_xy_grid(planet); // New functions for equal-angular grid (center, left, down): - report.print(3, "Scaling Cube by Radius"); + report.print(2, "Scaling Cube by Radius"); scale_cube_by_radius(cubeC); scale_cube_by_radius(cubeL); scale_cube_by_radius(cubeD); - report.print(3, "Done Scaling Cube by Radius"); + report.print(2, "Done Scaling Cube by Radius"); } - if (IsMagGrid) { + if (gridType == ionType_) { report.print(0, "--> Grid is Magnetic, so rotating"); std::vector llr, xyz, xyzRot1, xyzRot2; llr.push_back(geoLon_scgc); @@ -264,17 +269,17 @@ bool Grid::init_geo_grid(Quadtree quadtree, // Calculate PFPC coordinates (i.e., XYZ from LLR) calc_xyz(planet); - // Calculate grid spacing calc_grid_spacing(planet); //calculate radial unit vector (for spherical or oblate planet) calc_rad_unit(planet); // Calculate gravity (including J2 term, if desired) calc_gravity(planet); - // Calculate magnetic field and magnetic coordinates: fill_grid_bfield(planet); + write_restart(input.get_restartout_dir()); + // Throw a little message for students: report.student_checker_function_name(input.get_is_student(), input.get_student_name(), diff --git a/src/init_mag_grid.cpp b/src/init_mag_grid.cpp index 8fc17d2d..be4bf6c1 100644 --- a/src/init_mag_grid.cpp +++ b/src/init_mag_grid.cpp @@ -186,16 +186,16 @@ bool Grid::init_dipole_grid(Quadtree quadtree_ion, Planets planet) { IsCubeSphereGrid = false; IsDipole = true; - report.print(0, "Creating inter-node dipole connections"); + report.print(0, "Creating inter-node dipole connections for: " + gridType); if (!Is0D & !Is1Dz) create_dipole_connection(quadtree_ion); - report.print(0, "Creating Dipole Grid"); + report.print(0, "Creating Dipole Grid for: " + gridType); report.print(3, "Getting grid inputs for dipole grid"); - Inputs::grid_input_struct grid_input = input.get_grid_inputs("ionGrid"); + Inputs::grid_input_struct grid_input = input.get_grid_inputs(gridType); // Number of ghost cells: int64_t nGCs = get_nGCs(); @@ -215,6 +215,10 @@ bool Grid::init_dipole_grid(Quadtree quadtree_ion, Planets planet) { precision_t min_alt_re = (min_alt + planetRadius) / planetRadius; precision_t max_alt_re = (max_alt + planetRadius) / planetRadius; + // set the altitude of the lower boundary from the planet file + // -- this is used for setting densities hydrostatically. + altitude_lower_bc = planet.get_altitude_of_bc(); + if (nAlts % 2 != 0) { report.error("nAlts must be even!"); DidWork = false; @@ -625,6 +629,8 @@ bool Grid::init_dipole_grid(Quadtree quadtree_ion, Planets planet) { fill_grid_bfield(planet); report.print(4, "Done filling dipole grid with b-field!"); + write_restart(input.get_restartout_dir()); + report.exit(function); return DidWork; } diff --git a/src/init_parallel.cpp b/src/init_parallel.cpp index e89751d2..4da0c6a1 100644 --- a/src/init_parallel.cpp +++ b/src/init_parallel.cpp @@ -105,9 +105,8 @@ bool init_parallel(Quadtree &quadtree, Quadtree &quadtree_ion) { if (report.test_verbose(2)) std::cout << "seed : " << seed << "\n"; - quadtree.build("neuGrid"); - // #TODO - quadtree_ion.build("ionGrid"); + quadtree.build(neutralType_); + quadtree_ion.build(ionType_); } else { if (iProc == 0) { diff --git a/src/inputs.cpp b/src/inputs.cpp index d4303538..21dc05e0 100644 --- a/src/inputs.cpp +++ b/src/inputs.cpp @@ -107,7 +107,7 @@ std::string dummy_string = "unknown"; bool Inputs::check_settings(std::string key1, std::string key2) { - if (report.test_verbose(5)) + if (report.test_verbose(10)) std::cout << "checking setting : " << key1 << " and " << key2 << "\n"; @@ -132,7 +132,7 @@ bool Inputs::check_settings(std::string key1, // 1 key: bool Inputs::check_settings(std::string key1) { - if (report.test_verbose(5)) + if (report.test_verbose(10)) std::cout << "checking setting : " << key1 << "\n"; // try to find the keys first @@ -466,6 +466,25 @@ Inputs::grid_input_struct Inputs::get_grid_inputs(std::string gridtype) { grid_specs.alt_min = min_max[0]; grid_specs.alt_max = min_max[1]; + precision_t minDipoleLat = 10.0 * cDtoR; + precision_t maxDipoleLat = 80.0 * cDtoR; + + if (grid_specs.lat_min < minDipoleLat) { + grid_specs.lat_min = minDipoleLat; + report.print(0, "Error in setting min lat for " + + grid_specs.shape + + " - moving to 10 deg"); + report.error("Setting min dipole lat to 10.0"); + } + + if (grid_specs.lat_max > maxDipoleLat) { + grid_specs.lat_max = maxDipoleLat; + report.print(0, "Error in setting max lat for " + + grid_specs.shape + + " - moving to 80 deg"); + report.error("Setting max dipole lat to 80.0"); + } + } else { grid_specs.alt_min = check_settings_pt(gridtype, "MinAlt"); grid_specs.alt_file = check_settings_str(gridtype, "AltFile"); diff --git a/src/main/main.cpp b/src/main/main.cpp index 54e8316f..e3405588 100644 --- a/src/main/main.cpp +++ b/src/main/main.cpp @@ -38,12 +38,12 @@ int main() { // For now, the number of processors and blocks are set by the // neutral grid shape, since this could be sphere (1 root) or // cubesphere (6 root) - Quadtree quadtree(input.get_grid_shape("neuGrid")); + Quadtree quadtree(input.get_grid_shape(neutralType_)); if (!quadtree.is_ok()) throw std::string("quadtree for neutrals initialization failed!"); - Quadtree quadtree_ion(input.get_grid_shape("ionGrid")); + Quadtree quadtree_ion(input.get_grid_shape(ionType_)); if (!quadtree_ion.is_ok()) throw std::string("quadtree for ions initialization failed!"); @@ -84,8 +84,8 @@ int main() { // Perturb the inputs if user has asked for this indices.perturb(); - // Initialize Geographic grid: - Grid gGrid("neuGrid"); + // Initialize neutral grid: + Grid gGrid(neutralType_); didWork = gGrid.init_geo_grid(quadtree, planet); MPI_Barrier(aether_comm); @@ -102,25 +102,24 @@ int main() { if (input.get_cent_acc()) gGrid.calc_cent_acc(planet); - // Initialize Magnetic grid: - Grid mGrid("ionGrid"); + // Initialize ion grid: + Grid mGrid(ionType_); - if (mGrid.iGridShape_ == mGrid.iDipole_) { + if (mGrid.iGridShape_ == iDipole_) { didWork = mGrid.init_dipole_grid(quadtree_ion, planet); if (!didWork) throw std::string("init_dipole_grid failed!"); } else { - std::cout << "Making Spherical Magnetic Grid\n"; + report.print(1, "Making Spherical Magnetic Grid\n"); mGrid.set_IsDipole(false); - didWork = mGrid.init_geo_grid(quadtree, planet); + didWork = mGrid.init_geo_grid(quadtree_ion, planet); mGrid.set_IsGeoGrid(false); } if (!didWork) throw std::string("Initializing magneitic grid failed!"); - didWork = grid_match(gGrid, mGrid, quadtree, quadtree_ion); // Initialize Neutrals on geographic grid: @@ -197,6 +196,14 @@ int main() { if (!didWork) throw std::string("Reading Restart for time Failed!!!\n"); + + didWork = indices.restart_file(input.get_restartin_dir(), + true, + time.get_current()); + + if (!didWork) + throw std::string("Reading Restart for Indices Failed!!!\n"); + } // This is for the initial output. If it is not a restart, this will go: diff --git a/src/neutrals_advect.cpp b/src/neutrals_advect.cpp index dcdd67cd..613d4130 100644 --- a/src/neutrals_advect.cpp +++ b/src/neutrals_advect.cpp @@ -42,7 +42,7 @@ bool Neutrals::advect_horizontal(Grid & grid, Times & time) { static int iFunction = -1; report.enter(function, iFunction); - if (grid.iGridShape_ == grid.iCubesphere_) { + if (grid.iGridShape_ == iCubesphere_) { if (input.get_advection_neutrals_horizontal() == "advect_test") solver_horizontal_RK4_advection(grid, time); else if (input.get_advection_neutrals_horizontal() == "fv") diff --git a/src/neutrals_bcs.cpp b/src/neutrals_bcs.cpp index 59990c51..ecb87da2 100644 --- a/src/neutrals_bcs.cpp +++ b/src/neutrals_bcs.cpp @@ -207,8 +207,23 @@ bool Neutrals::set_lower_bcs(Grid grid, temperature_scgc.subcube(iLon, iLat, 0, iLon, iLat, iAlt - 1).fill( temperature_scgc(iLon, iLat, iAlt)); + precision_t t = temperature_scgc(iLon, iLat, 0); + precision_t g = abs(grid.gravity_vcgc[2](iLon, iLat, iAlt)); + + precision_t alt1 = grid.geoAlt_scgc(iLon, iLat, iAlt); + precision_t alt0 = grid.altitude_lower_bc; + precision_t dz = alt1 - alt0; + for (iSpecies = 0; iSpecies < nSpecies; iSpecies++) { + precision_t m = mean_major_mass_scgc(iLon, iLat, iAlt); + + //if (m == 0) + m = species[iSpecies].mass; + + precision_t h = cKB * t / (m * g); + precision_t factor = exp(-dz / h); + //----------------------------------------------- // Planet BCs - set to fixed constant values. //----------------------------------------------- @@ -217,6 +232,7 @@ bool Neutrals::set_lower_bcs(Grid grid, // Fill all lower ghost cells density with lower boundary condition: species[iSpecies].density_scgc.subcube(iLon, iLat, 0, iLon, iLat, iAlt - 1).fill( + factor * species[iSpecies].lower_bc_density); } // planet bc type diff --git a/src/output.cpp b/src/output.cpp index f4adf0d9..08ba1dd4 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -362,7 +362,7 @@ bool output(const Neutrals &neutrals, report.error("File output type not found!"); didWork = false; } else { - if (grid.get_IsGeoGrid()) + if (grid.get_gridtype() == neutralType_) filename = filename + "G_"; else filename = filename + "M_"; diff --git a/src/output_netcdf.cpp b/src/output_netcdf.cpp index 4a30ef9e..0e2501bf 100644 --- a/src/output_netcdf.cpp +++ b/src/output_netcdf.cpp @@ -72,8 +72,10 @@ bool OutputContainer::read_container_netcdf() { std::string UNITS = "units"; try { - std::cout << "Reading NetCDF file into container : " - << whole_filename << "\n"; + if (report.test_verbose(0)) + std::cout << "Reading NetCDF file into container : " + << whole_filename << "\n"; + NcFile ncdf_file_in(whole_filename, NcFile::read); std::multimap variables_in_file; diff --git a/src/planets.cpp b/src/planets.cpp index 0528cfd7..200786b7 100644 --- a/src/planets.cpp +++ b/src/planets.cpp @@ -502,6 +502,14 @@ json Planets::get_ions() { return ions; } +// -------------------------------------------------------------------------- +// returns altitude of the density boundary condition +// -------------------------------------------------------------------------- + +precision_t Planets::get_altitude_of_bc() { + return altitude_of_bc; +} + // ----------------------------------------------------------------------------- // Read in the planet specific file that describes the species // ----------------------------------------------------------------------------- @@ -544,6 +552,17 @@ bool Planets::read_planet_specific_file() { std::cout << neutrals << "\n"; } // #neutrals + if (hash == "#altitude_of_bc") { + report.print(iDebug, "Found #altitude_of_bc!"); + altitude_of_bc = read_float(infile_ptr, "#altitude_of_bc"); + // Units read in = km + // Units needed in code = m: + altitude_of_bc = altitude_of_bc * 1000.0; + + if (report.test_verbose(iDebug)) + std::cout << altitude_of_bc << "\n"; + } // #altitude_of_bc + if (hash == "#temperature") { report.print(iDebug, "Found #temperatures!"); std::vector> lines = read_csv(infile_ptr); diff --git a/src/quadtree.cpp b/src/quadtree.cpp index aa61b98e..832e5a78 100644 --- a/src/quadtree.cpp +++ b/src/quadtree.cpp @@ -1,12 +1,17 @@ // Copyright 2024, the Aether Development Team (see doc/dev_team.md for members) // Full license can be found in License.md -// Need to allow more types of grids. We have two axes of grids, really: -// - Neutral -// - Ion -// Within each of those, we can have several types of grids: +// The way that the quadtree works is that you start with a given number of +// root nodes. When the user wants higher resolution, they as for 4 times +// more processors, and then each root node is broken into 4. If 16 times +// more processors are asked for, then those 4 nodes are each broken into +// 4 more nodes. This goes on for as many processors as the user would like, +// but the processors has to equal nRootNodes * 4^depth +// +// The following grid shapes are suppored at this time: // - Cubesphere, this has 6 root nodes (2 polar, 4 equatorial) // - Sphere, this has 1 root node (whole grid) +// - Sphere4, this is a spherical grid, but has 4 root nodes (2 lats, 2 lons) // - Sphere6, this is a spherical grid, but has 6 root nodes (2 lats, 3 lons) // - Dipole4, which has 4 root nodes (4 lats, 1 lon) // - Dipole6, which has 4 root nodes (6 lats, 1 lon) @@ -34,6 +39,11 @@ Quadtree::Quadtree(std::string shapeInput) { IsOk = true; } + if (shape == "sphere6") { + nRootNodes = 6; + IsOk = true; + } + if (shape == "dipole4") { nRootNodes = 4; IsOk = true; @@ -63,17 +73,28 @@ bool Quadtree::is_ok() { void Quadtree::build(std::string gridtype) { + std::string function = "Quadtree::build"; + static int iFunction = -1; + report.enter(function, iFunction); + + IsOk = false; + arma_mat origins; arma_mat rights; arma_mat ups; Inputs::grid_input_struct grid_input = input.get_grid_inputs(gridtype); + // Here we are taking the shape and getting the sizes and positions + // of the root nodes. These are defined in the different header files + // such as sphere.h, cubesphere.h, and dipole.h + if (grid_input.shape == "cubesphere") { origins = CubeSphere::ORIGINS; rights = CubeSphere::RIGHTS; ups = CubeSphere::UPS; IsCubeSphere = true; + IsOk = true; } if (grid_input.shape == "sphere") { @@ -81,6 +102,7 @@ void Quadtree::build(std::string gridtype) { rights = Sphere::RIGHTS; ups = Sphere::UPS; IsSphere = true; + IsOk = true; } if (grid_input.shape == "sphere4") { @@ -88,6 +110,15 @@ void Quadtree::build(std::string gridtype) { rights = Sphere4::RIGHTS; ups = Sphere4::UPS; IsSphere = true; + IsOk = true; + } + + if (grid_input.shape == "sphere6") { + origins = Sphere6::ORIGINS; + rights = Sphere6::RIGHTS; + ups = Sphere6::UPS; + IsSphere = true; + IsOk = true; } if (grid_input.shape == "dipole4") { @@ -95,6 +126,7 @@ void Quadtree::build(std::string gridtype) { rights = Dipole4::RIGHTS; ups = Dipole4::UPS; IsDipole = true; + IsOk = true; } if (grid_input.shape == "dipole6") { @@ -102,7 +134,14 @@ void Quadtree::build(std::string gridtype) { rights = Dipole6::RIGHTS; ups = Dipole6::UPS; IsDipole = true; + IsOk = true; + } + // If we can't find the shape, then there is a big problem + if (!IsOk) { + report.error("quadtree shape not found (in build): " + grid_input.shape); + report.exit(function); + return; } arma_vec o(3), r(3), u(3); @@ -129,6 +168,9 @@ void Quadtree::build(std::string gridtype) { // Before we build the quadtree, we need to allow the user to // restrict the domain. This will only work for the spherical // grid so far: + // (as programmed, this should work ok for the sphere and + // dipole shapes, but will never work for cubesphere. For the + // cubesphere grid, it is much more complicated.) if ((grid_input.lon_min > 0.0 || grid_input.lon_max < 2.0 * cPI || @@ -163,6 +205,9 @@ void Quadtree::build(std::string gridtype) { tmp = new_node(o, r, u, iP, iDepth, iNode); root_nodes.push_back(tmp); } + + report.exit(function); + return; } // -------------------------------------------------------------------------- @@ -382,7 +427,14 @@ int64_t Quadtree::find_point(arma_vec point, Quadtree::qtnode node) { } // -------------------------------------------------------------------------- -// +// This takes a normalized point, figures out if it is beyond the limits +// of the root node, and if it is, then moves the coordinates onto the other +// node. +// This is pretty much useful for the CubeSphere, since when you go over +// the edge of one side, you are technically then on another side. This can +// happen on the sphere also, when you go across the 0/360 line (or 0/2 line +// in normalized coordinates). If can aslo happen at the poles when you go +// over the pole. // -------------------------------------------------------------------------- arma_vec Quadtree::wrap_point_sphere(arma_vec point) { @@ -435,7 +487,8 @@ arma_vec Quadtree::wrap_point_sphere(arma_vec point) { } // -------------------------------------------------------------------------- -// +// Well, ok - the above wrap_point seems to only work for the sphere and +// dipole shape, which this is specially designed for the cubesphere // -------------------------------------------------------------------------- arma_vec Quadtree::wrap_point_cubesphere(arma_vec point) { diff --git a/src/read_input_file.cpp b/src/read_input_file.cpp index 389e0878..66a5e883 100644 --- a/src/read_input_file.cpp +++ b/src/read_input_file.cpp @@ -24,6 +24,9 @@ bool Inputs::read_inputs_json(Times &time) { json defaults; json user_inputs; + // allow changing of perturbations during the restart process: + json perturbations; + isOk = true; // Set the default values first: @@ -54,6 +57,9 @@ bool Inputs::read_inputs_json(Times &time) { // if they really want: restart_inputs["Logfile"]["append"] = true; settings.merge_patch(restart_inputs); + + if (restart_inputs.contains("Perturb")) + perturbations["Perturb"] = restart_inputs["Perturb"]; } } } @@ -62,6 +68,21 @@ bool Inputs::read_inputs_json(Times &time) { // settings, with the default/restart settings being the default: settings.merge_patch(user_inputs); + // There are perturbations in the restart files: + if (perturbations.contains("Perturb")) + + // If the user wants the restart perturbations to overwrite the + // aether.json perturbations, then do it: + if (user_inputs.contains("Perturb")) { + if (user_inputs["Perturb"].contains("restart_control")) + if (user_inputs["Perturb"]["restart_control"]) + settings.merge_patch(perturbations); + } else + // if there are perturbations in the restart files, but none + // in the user files, then push the restart perturbations into + // the settings to make them consistent + settings.merge_patch(perturbations); + //change planet file to the one specified on aether.json: if (isOk) settings["PlanetSpeciesFile"] = get_setting_str("Planet", "file"); diff --git a/src/solver_advection.cpp b/src/solver_advection.cpp index 48e88688..8da03682 100644 --- a/src/solver_advection.cpp +++ b/src/solver_advection.cpp @@ -284,7 +284,7 @@ precision_t calc_dt(arma_mat &xWidth, void Neutrals::advect_sphere(Grid &grid, Times &time) { - std::string function = "advect"; + std::string function = "advect_sphere"; static int iFunction = -1; report.enter(function, iFunction); @@ -369,10 +369,12 @@ void Neutrals::advect_sphere(Grid &grid, Times &time) { xMomentum = rho % xVel; yMomentum = rho % yVel; - x = grid.x_Center.slice(iAlt) * grid.radius_scgc(1, 1, iAlt); - y = grid.y_Center.slice(iAlt) * grid.radius_scgc(1, 1, iAlt); - xEdges = grid.x_Left.slice(iAlt) * grid.radius_scgc(1, 1, iAlt); - yEdges = grid.y_Down.slice(iAlt) * grid.radius_scgc(1, 1, iAlt); + precision_t radius = grid.radius_scgc(1, 1, iAlt); + + x = grid.x_Center.slice(iAlt) * radius; + y = grid.y_Center.slice(iAlt) * radius; + xEdges = grid.x_Left.slice(iAlt) * radius; + yEdges = grid.y_Down.slice(iAlt) * radius; rhoP = project_to_edges(rho, x, xEdges, y, yEdges, nGCs); xVelP = project_to_edges(xVel, x, xEdges, y, yEdges, nGCs); diff --git a/src/solver_chemistry.cpp b/src/solver_chemistry.cpp index 3c9fd21d..ce40ca00 100644 --- a/src/solver_chemistry.cpp +++ b/src/solver_chemistry.cpp @@ -12,8 +12,20 @@ arma_cube solver_chemistry(arma_cube density, arma_cube source, arma_cube loss, precision_t dt) { - arma_cube normalized_loss = loss / (density + 1e-6); + arma_cube normalized_loss; + normalized_loss = loss / (density + 1e-6); arma_cube new_density = (density + dt * source) / (1.0 + dt * normalized_loss); return new_density; } + +arma_mat solver_chemistry(arma_mat density, + arma_mat source, + arma_mat loss, + precision_t dt) { + arma_mat normalized_loss; + normalized_loss = loss / (density + 1e-6); + arma_mat new_density = (density + dt * source) / + (1.0 + dt * normalized_loss); + return new_density; +} diff --git a/src/solver_conduction.cpp b/src/solver_conduction.cpp index ff353a11..e6de684a 100644 --- a/src/solver_conduction.cpp +++ b/src/solver_conduction.cpp @@ -29,7 +29,7 @@ arma_vec solver_conduction(arma_vec value, int64_t nGCs, bool return_diff, // (optional) False by default (return new `value`) arma_vec source2 // (optional) Sources dependent on `value` - ) { + ) { int64_t nPts = value.n_elem; @@ -52,7 +52,7 @@ arma_vec solver_conduction(arma_vec value, conduction.zeros(); // If source2 is not given, set it to zero: - if (source2.n_elem == 0){ + if (source2.n_elem == 0) { source2.set_size(source.n_elem); source2.zeros(); } @@ -64,7 +64,8 @@ arma_vec solver_conduction(arma_vec value, arma_vec a = di / du22 % r - dl / du12 % r % r; arma_vec c = di / du22 + dl / du12; - arma_vec b = -1.0 / m - di / du22 % (1.0 + r) - dl / du12 % (1.0 - r % r) + source2 % front * dt; + arma_vec b = -1.0 / m - di / du22 % (1.0 + r) - dl / du12 % + (1.0 - r % r) + source2 % front * dt; arma_vec d = -1.0 * (value / m + source % front * dt); // Lower BCs (fixed value): diff --git a/src/solver_gradients.cpp b/src/solver_gradients.cpp index 73754c10..6a0d0ac2 100644 --- a/src/solver_gradients.cpp +++ b/src/solver_gradients.cpp @@ -18,9 +18,9 @@ std::vector calc_gradient_vector(arma_cube value_scgc, Grid grid) { display_vector("gradient, value : ", value_scgc.tube(9, 9)); } - if (grid.iGridShape_ == grid.iCubesphere_) + if (grid.iGridShape_ == iCubesphere_) gradient_vcgc = calc_gradient_cubesphere(value_scgc, grid); - else if (grid.iGridShape_ == grid.iDipole_) + else if (grid.iGridShape_ == iDipole_) gradient_vcgc = calc_gradient_dipole(value_scgc, grid); else { diff --git a/src/solver_grid_interpolation.cpp b/src/solver_grid_interpolation.cpp index 409ffd3d..e5e294bb 100644 --- a/src/solver_grid_interpolation.cpp +++ b/src/solver_grid_interpolation.cpp @@ -382,7 +382,8 @@ void Grid::set_interp_coef_dipole(const dipole_range &dr, abs(j_center_scgc.tube(coef.iRow, coef.iCol)), nGCs); // need alt index to find lat coef - coef.iAlt = bisect_search_array(alt_in, k_center_scgc.tube(coef.iRow, coef.iCol), + coef.iAlt = bisect_search_array(alt_in, k_center_scgc.tube(coef.iRow, + coef.iCol), nGCs); // then we can do the ratios: @@ -425,6 +426,8 @@ bool Grid::set_interpolation_coefs(const std::vector &i_coords, static int iFunction = -1; report.enter(function, iFunction); + report.print(1, "interpolation gridtype : " + gridType); + // If the size of Lons, Lats and Alts are not the same, return false if (i_coords.size() != j_coords.size() || j_coords.size() != k_coords.size()) { report.error("Length of i,j,k vectors do not match!"); @@ -434,28 +437,31 @@ bool Grid::set_interpolation_coefs(const std::vector &i_coords, // Clear the previous interpolation coefficients interp_coefs.clear(); - if (IsGeoGrid) { - // Handle according to whether it is cubesphere or not - if (IsCubeSphereGrid) { - // Calculate the range of the grid - struct cubesphere_range cr; - get_cubesphere_grid_range(cr); - - // Calculate the index and coefficients for each point - for (size_t i = 0; i < i_coords.size(); ++i) - set_interp_coef_cubesphere(cr, i_coords[i], j_coords[i], k_coords[i]); - } else if (IsLatLonGrid) { - // Calculate the range of the grid - struct sphere_range sr; - get_sphere_grid_range(sr); - - // Calculate the index and coefficients for each point - for (size_t i = 0; i < i_coords.size(); ++i) - set_interp_coef_sphere(sr, i_coords[i], j_coords[i], k_coords[i]); - } + if (iGridShape_ == iCubesphere_) { + report.print(1, "interpolation grid is cubesphere"); + + // Calculate the range of the grid + struct cubesphere_range cr; + get_cubesphere_grid_range(cr); + + // Calculate the index and coefficients for each point + for (size_t i = 0; i < i_coords.size(); ++i) + set_interp_coef_cubesphere(cr, i_coords[i], j_coords[i], k_coords[i]); } - else { // IsDipole + if (iGridShape_ == iSphere_) { + report.print(1, "interpolation grid is sphere"); + + struct sphere_range sr; + get_sphere_grid_range(sr); + + // Calculate the index and coefficients for each point + for (size_t i = 0; i < i_coords.size(); ++i) + set_interp_coef_sphere(sr, i_coords[i], j_coords[i], k_coords[i]); + } + + if (iGridShape_ == iDipole_) { // IsDipole + report.print(1, "interpolation grid is dipole"); // Calculate the range of the grid struct dipole_range dr; @@ -489,22 +495,27 @@ bool Grid::set_interpolation_coefs(const std::vector &i_coords, } else { - magCoords = {vec2cube(i_coords),vec2cube(j_coords), vec2cube(k_coords)}; + magCoords = {vec2cube(i_coords), vec2cube(j_coords), vec2cube(k_coords)}; } // std::vector dipcoords = geo_to_mag(i_coords[0], j_coords[0], k_coords[0], planet); std::vector dipCoords; + if (!areLocsIJK) { std::vector planet_radii(nPts); + for (iLoc = 0; iLoc < nPts; iLoc++) { // Convert from mag->dipole coordinates. - if (areLocsGeo){ // we were given the geo-latitude + if (areLocsGeo) // we were given the geo-latitude planet_radii[iLoc] = planet.get_radius(j_coords[iLoc]); - } else{ + + else { // equatorial radius :( planet_radii[iLoc] = planet.get_radius(0.0); } - dipCoords = mag_to_ijk(i_coords[iLoc], j_coords[iLoc], k_coords[iLoc], planet_radii[iLoc]); + + dipCoords = mag_to_ijk(i_coords[iLoc], j_coords[iLoc], k_coords[iLoc], + planet_radii[iLoc]); mlon[iLoc] = dipCoords[0]; p_coord[iLoc] = dipCoords[1]; q_coord[iLoc] = dipCoords[2]; diff --git a/src/time.cpp b/src/time.cpp index 6accc347..e15d4a0d 100644 --- a/src/time.cpp +++ b/src/time.cpp @@ -58,11 +58,22 @@ bool Times::restart_file(std::string dir, bool DoRead) { iStep--; dt = 0; increment_time(); - std::cout << "Restarted time, Current time : "; - display_itime(iCurrent); + + if (report.test_verbose(0)) { + std::cout << "Restarted time, Current time : "; + display_itime(iCurrent); + } } else { - restart_time_json = { {"currenttime", current}, - {"istep", iStep} + restart_time_json = { + {"currenttime", current}, + {"istep", iStep}, + {"year", year}, + {"month", month}, + {"day", day}, + {"hour", hour}, + {"minute", minute}, + {"second", second}, + {"millisecond", milli}, }; DidWork = write_json(filename, restart_time_json); } diff --git a/srcPython/postAether.py b/srcPython/postAether.py index 694b2ad4..a938ffb2 100755 --- a/srcPython/postAether.py +++ b/srcPython/postAether.py @@ -1024,9 +1024,11 @@ def calc_blocks(dataToWrite, iLon_ = 0, iLat_ = 1): iBlock = 1 while (iBlock < nBlocksTotal): - if (np.abs(dataToWrite[iBlock][iLon_][0, 0, 0] - testLon) > 0.001): + if ((np.abs(dataToWrite[iBlock][iLon_][0, 0, 0] - testLon) > 0.1) and + (np.abs(dataToWrite[iBlock][iLat_][0, 0, 0] - testLat) < 0.1)): nBlocksLon = nBlocksLon + 1 - if (np.abs(dataToWrite[iBlock][iLat_][0, 0, 0] - testLat) > 0.001): + if ((np.abs(dataToWrite[iBlock][iLat_][0, 0, 0] - testLat) > 0.1) and + (np.abs(dataToWrite[iBlock][iLon_][0, 0, 0] - testLon) < 0.1)): nBlocksLat = nBlocksLat + 1 iBlock = iBlock + 1 @@ -1064,20 +1066,20 @@ def consolidate_blocks(originalData, iLon_ = 0, iLat_ = 1): consolidatedData[key] = originalData[0][key] else: # need to move data over + #print('variable : ', key) data = np.zeros((nLonsTotal, nLatsTotal, nAlts)) for iBlock in range(nBlocks): # interior points: - iLatS = int((originalData[iBlock][iLat_][nGCs, nGCs, nGCs] - Lat0)/dLat) + nGCs + iLatS = int(round((originalData[iBlock][iLat_][nGCs, nGCs, nGCs] - Lat0)/dLat)) + nGCs iLatE = iLatS + nLats - iLonS = int((originalData[iBlock][iLon_][nGCs, nGCs, nGCs] - Lon0)/dLon) + nGCs + iLonS = int(round((originalData[iBlock][iLon_][nGCs, nGCs, nGCs] - Lon0)/dLon)) + nGCs iLonE = iLonS + nLons iLonSO = nGCs iLonEO = nGCs + nLons iLatSO = nGCs iLatEO = iLatSO + nLats - #print(iLonS, iLonE, nLonsTotal, ' -> ', iLonSO, iLonEO, nLons) - #print(iLatS, iLatE, nLatsTotal, ' -> ', iLatSO, iLatEO, nLats) - #print(nGCs) + #print('lons : ', iLonS, iLonE, nLonsTotal, ' -> ', iLonSO, iLonEO, nLons) + #print('lats : ', iLatS, iLatE, nLatsTotal, ' -> ', iLatSO, iLatEO, nLats) data[iLonS:iLonE, iLatS:iLatE, 0:nAlts] = \ originalData[iBlock][key][iLonSO:iLonEO, iLatSO:iLatEO, 0:nAlts] diff --git a/tests/grid_shapes/aether.json.whole b/tests/grid_shapes/aether.json.whole new file mode 100644 index 00000000..84033c7b --- /dev/null +++ b/tests/grid_shapes/aether.json.whole @@ -0,0 +1,42 @@ + +{ + "Debug" : { + "dt": 1.0, + "TimingPercent": 10.0, + "iVerbose" : 0}, + + "EndTime" : [2011, 3, 20, 0, 1, 0], + + "neuGrid" : { + "Shape": "cubesphere", + "nLonsPerBlock" : 22, + "nLatsPerBlock" : 22, + "nAlts" : 30, + "MinAlt": 95, + "LatRange": [-90, 90], + "LonRange": [0.0, 360.0], + "dAltScale" : 0.3, + "IsUniformAlt" : false}, + + "ionGrid": { + "Shape": "sphere6", + "nLonsPerBlock": 32, + "nLatsPerBlock": 16, + "nAlts": 40, + "MinAlt": 80, + "LatRange": [-90.0, 90.0], + "LonRange": [0.0, 360.0], + "dAltScale" : 0.3, + "IsUniformAlt" : false}, + + "OmniwebFiles" : ["UA/inputs/omni_20110319.txt"], + + "Electrodynamics" : { + "Potential" : "Weimer05", + "DiffuseAurora" : "fta"}, + + "Outputs" : { + "type" : ["states"], + "dt" : [60] } + +} diff --git a/tests/grid_shapes/aether_cube_cube.json b/tests/grid_shapes/aether_cube_cube.json new file mode 100644 index 00000000..b8e19691 --- /dev/null +++ b/tests/grid_shapes/aether_cube_cube.json @@ -0,0 +1,42 @@ + +{ + "Debug" : { + "dt": 1.0, + "TimingPercent": 10.0, + "iVerbose" : 0}, + + "EndTime" : [2011, 3, 20, 0, 1, 0], + + "neuGrid" : { + "Shape": "cubesphere", + "nLonsPerBlock" : 22, + "nLatsPerBlock" : 22, + "nAlts" : 30, + "MinAlt": 95, + "LatRange": [-90, 90], + "LonRange": [0.0, 360.0], + "dAltScale" : 0.3, + "IsUniformAlt" : false}, + + "ionGrid": { + "Shape": "cubesphere", + "nLonsPerBlock": 32, + "nLatsPerBlock": 32, + "nAlts": 40, + "MinAlt": 80, + "LatRange": [-90, 90], + "LonRange": [0.0, 360.0], + "dAltScale" : 0.3, + "IsUniformAlt" : false}, + + "OmniwebFiles" : ["UA/inputs/omni_20110319.txt"], + + "Electrodynamics" : { + "Potential" : "Weimer05", + "DiffuseAurora" : "fta"}, + + "Outputs" : { + "type" : ["states"], + "dt" : [60] } + +} diff --git a/tests/grid_shapes/aether_cube_dipole6.json b/tests/grid_shapes/aether_cube_dipole6.json new file mode 100644 index 00000000..696d0352 --- /dev/null +++ b/tests/grid_shapes/aether_cube_dipole6.json @@ -0,0 +1,42 @@ + +{ + "Debug" : { + "dt": 1.0, + "TimingPercent": 10.0, + "iVerbose" : 0}, + + "EndTime" : [2011, 3, 20, 0, 1, 0], + + "neuGrid" : { + "Shape": "cubesphere", + "nLonsPerBlock" : 22, + "nLatsPerBlock" : 22, + "nAlts" : 30, + "MinAlt": 95, + "LatRange": [-90, 90], + "LonRange": [0.0, 360.0], + "dAltScale" : 0.3, + "IsUniformAlt" : false}, + + "ionGrid": { + "Shape": "dipole6", + "nLonsPerBlock": 32, + "nLatsPerBlock": 16, + "nAlts": 40, + "MinAlt": 80, + "LatRange": [10.0, 80.0], + "LonRange": [0.0, 360.0], + "dAltScale" : 0.3, + "IsUniformAlt" : false}, + + "OmniwebFiles" : ["UA/inputs/omni_20110319.txt"], + + "Electrodynamics" : { + "Potential" : "Weimer05", + "DiffuseAurora" : "fta"}, + + "Outputs" : { + "type" : ["states"], + "dt" : [60] } + +} diff --git a/tests/grid_shapes/aether_cube_sphere6.json b/tests/grid_shapes/aether_cube_sphere6.json new file mode 100644 index 00000000..84033c7b --- /dev/null +++ b/tests/grid_shapes/aether_cube_sphere6.json @@ -0,0 +1,42 @@ + +{ + "Debug" : { + "dt": 1.0, + "TimingPercent": 10.0, + "iVerbose" : 0}, + + "EndTime" : [2011, 3, 20, 0, 1, 0], + + "neuGrid" : { + "Shape": "cubesphere", + "nLonsPerBlock" : 22, + "nLatsPerBlock" : 22, + "nAlts" : 30, + "MinAlt": 95, + "LatRange": [-90, 90], + "LonRange": [0.0, 360.0], + "dAltScale" : 0.3, + "IsUniformAlt" : false}, + + "ionGrid": { + "Shape": "sphere6", + "nLonsPerBlock": 32, + "nLatsPerBlock": 16, + "nAlts": 40, + "MinAlt": 80, + "LatRange": [-90.0, 90.0], + "LonRange": [0.0, 360.0], + "dAltScale" : 0.3, + "IsUniformAlt" : false}, + + "OmniwebFiles" : ["UA/inputs/omni_20110319.txt"], + + "Electrodynamics" : { + "Potential" : "Weimer05", + "DiffuseAurora" : "fta"}, + + "Outputs" : { + "type" : ["states"], + "dt" : [60] } + +} diff --git a/tests/grid_shapes/aether_sphere4_dipole4.json b/tests/grid_shapes/aether_sphere4_dipole4.json new file mode 100644 index 00000000..19b55f43 --- /dev/null +++ b/tests/grid_shapes/aether_sphere4_dipole4.json @@ -0,0 +1,42 @@ + +{ + "Debug" : { + "dt": 1.0, + "TimingPercent": 10.0, + "iVerbose" : 0}, + + "EndTime" : [2011, 3, 20, 0, 1, 0], + + "neuGrid" : { + "Shape": "sphere4", + "nLonsPerBlock" : 22, + "nLatsPerBlock" : 22, + "nAlts" : 30, + "MinAlt": 95, + "LatRange": [-90, 90], + "LonRange": [0.0, 360.0], + "dAltScale" : 0.3, + "IsUniformAlt" : false}, + + "ionGrid": { + "Shape": "dipole4", + "nLonsPerBlock": 32, + "nLatsPerBlock": 16, + "nAlts": 40, + "MinAlt": 80, + "LatRange": [10.0, 80.0], + "LonRange": [0.0, 360.0], + "dAltScale" : 0.3, + "IsUniformAlt" : false}, + + "OmniwebFiles" : ["UA/inputs/omni_20110319.txt"], + + "Electrodynamics" : { + "Potential" : "Weimer05", + "DiffuseAurora" : "fta"}, + + "Outputs" : { + "type" : ["states"], + "dt" : [60] } + +} diff --git a/tests/grid_shapes/aether_sphere4_sphere4.json b/tests/grid_shapes/aether_sphere4_sphere4.json new file mode 100644 index 00000000..b402b0f5 --- /dev/null +++ b/tests/grid_shapes/aether_sphere4_sphere4.json @@ -0,0 +1,42 @@ + +{ + "Debug" : { + "dt": 1.0, + "TimingPercent": 10.0, + "iVerbose" : 0}, + + "EndTime" : [2011, 3, 20, 0, 1, 0], + + "neuGrid" : { + "Shape": "sphere4", + "nLonsPerBlock" : 22, + "nLatsPerBlock" : 22, + "nAlts" : 30, + "MinAlt": 95, + "LatRange": [-90, 90], + "LonRange": [0.0, 360.0], + "dAltScale" : 0.3, + "IsUniformAlt" : false}, + + "ionGrid": { + "Shape": "sphere4", + "nLonsPerBlock": 32, + "nLatsPerBlock": 16, + "nAlts": 40, + "MinAlt": 80, + "LatRange": [-90, 90], + "LonRange": [0.0, 360.0], + "dAltScale" : 0.3, + "IsUniformAlt" : false}, + + "OmniwebFiles" : ["UA/inputs/omni_20110319.txt"], + + "Electrodynamics" : { + "Potential" : "Weimer05", + "DiffuseAurora" : "fta"}, + + "Outputs" : { + "type" : ["states"], + "dt" : [60] } + +} diff --git a/tests/grid_shapes/aether_sphere6_dipole6.json b/tests/grid_shapes/aether_sphere6_dipole6.json new file mode 100644 index 00000000..c0bac693 --- /dev/null +++ b/tests/grid_shapes/aether_sphere6_dipole6.json @@ -0,0 +1,42 @@ + +{ + "Debug" : { + "dt": 1.0, + "TimingPercent": 10.0, + "iVerbose" : 0}, + + "EndTime" : [2011, 3, 20, 0, 1, 0], + + "neuGrid" : { + "Shape": "sphere6", + "nLonsPerBlock" : 26, + "nLatsPerBlock" : 22, + "nAlts" : 30, + "MinAlt": 95, + "LatRange": [-90, 90], + "LonRange": [0.0, 360.0], + "dAltScale" : 0.3, + "IsUniformAlt" : false}, + + "ionGrid": { + "Shape": "dipole6", + "nLonsPerBlock": 32, + "nLatsPerBlock": 16, + "nAlts": 40, + "MinAlt": 80, + "LatRange": [10.0, 80.0], + "LonRange": [0.0, 360.0], + "dAltScale" : 0.3, + "IsUniformAlt" : false}, + + "OmniwebFiles" : ["UA/inputs/omni_20110319.txt"], + + "Electrodynamics" : { + "Potential" : "Weimer05", + "DiffuseAurora" : "fta"}, + + "Outputs" : { + "type" : ["states"], + "dt" : [60] } + +} diff --git a/tests/grid_shapes/aether_sphere6_sphere6.json b/tests/grid_shapes/aether_sphere6_sphere6.json new file mode 100644 index 00000000..3b8a1c68 --- /dev/null +++ b/tests/grid_shapes/aether_sphere6_sphere6.json @@ -0,0 +1,42 @@ + +{ + "Debug" : { + "dt": 1.0, + "TimingPercent": 10.0, + "iVerbose" : 0}, + + "EndTime" : [2011, 3, 20, 0, 1, 0], + + "neuGrid" : { + "Shape": "sphere6", + "nLonsPerBlock" : 22, + "nLatsPerBlock" : 22, + "nAlts" : 30, + "MinAlt": 95, + "LatRange": [-90, 90], + "LonRange": [0.0, 360.0], + "dAltScale" : 0.3, + "IsUniformAlt" : false}, + + "ionGrid": { + "Shape": "sphere6", + "nLonsPerBlock": 32, + "nLatsPerBlock": 16, + "nAlts": 40, + "MinAlt": 80, + "LatRange": [-90.0, 90.0], + "LonRange": [0.0, 360.0], + "dAltScale" : 0.3, + "IsUniformAlt" : false}, + + "OmniwebFiles" : ["UA/inputs/omni_20110319.txt"], + + "Electrodynamics" : { + "Potential" : "Weimer05", + "DiffuseAurora" : "fta"}, + + "Outputs" : { + "type" : ["states"], + "dt" : [60] } + +} diff --git a/tests/grid_shapes/aether_sphere_sphere.json b/tests/grid_shapes/aether_sphere_sphere.json new file mode 100644 index 00000000..935880b3 --- /dev/null +++ b/tests/grid_shapes/aether_sphere_sphere.json @@ -0,0 +1,42 @@ + +{ + "Debug" : { + "dt": 1.0, + "TimingPercent": 10.0, + "iVerbose" : 0}, + + "EndTime" : [2011, 3, 20, 0, 1, 0], + + "neuGrid" : { + "Shape": "sphere", + "nLonsPerBlock" : 32, + "nLatsPerBlock" : 16, + "nAlts" : 30, + "MinAlt": 95, + "LatRange": [-90, 90], + "LonRange": [0.0, 360.0], + "dAltScale" : 0.3, + "IsUniformAlt" : false}, + + "ionGrid": { + "Shape": "sphere", + "nLonsPerBlock": 42, + "nLatsPerBlock": 22, + "nAlts": 40, + "MinAlt": 80, + "LatRange": [-90, 90], + "LonRange": [0.0, 360.0], + "dAltScale" : 0.3, + "IsUniformAlt" : false}, + + "OmniwebFiles" : ["UA/inputs/omni_20110319.txt"], + + "Electrodynamics" : { + "Potential" : "Weimer05", + "DiffuseAurora" : "fta"}, + + "Outputs" : { + "type" : ["states"], + "dt" : [60] } + +} diff --git a/tests/grid_shapes/run_test.sh b/tests/grid_shapes/run_test.sh new file mode 100755 index 00000000..4c87f3d0 --- /dev/null +++ b/tests/grid_shapes/run_test.sh @@ -0,0 +1,92 @@ +#!/bin/sh + +RUN=sphere_sphere +PE=1 +rm -rf ./run.${RUN} +cp -R ../../share/run ./run.${RUN} +cd run.${RUN} +cp ../aether_${RUN}.json ./aether.json +mpirun -np ${PE} ./aether +../../../srcPython/postAether.py -rm +cd .. + +RUN=sphere_sphere +PE=4 +rm -rf ./run.${RUN} +cp -R ../../share/run ./run.${RUN} +cd run.${RUN} +cp ../aether_${RUN}.json ./aether.json +mpirun -np ${PE} ./aether +../../../srcPython/postAether.py -rm +cd .. + +RUN=sphere4_sphere4 +PE=4 +rm -rf ./run.${RUN} +cp -R ../../share/run ./run.${RUN} +cd run.${RUN} +cp ../aether_${RUN}.json ./aether.json +mpirun -np ${PE} ./aether +../../../srcPython/postAether.py -rm +cd .. + +RUN=sphere6_sphere6 +PE=6 +rm -rf ./run.${RUN} +cp -R ../../share/run ./run.${RUN} +cd run.${RUN} +cp ../aether_${RUN}.json ./aether.json +mpirun -np ${PE} ./aether +../../../srcPython/postAether.py -rm +cd .. + +RUN=cube_cube +PE=6 +rm -rf ./run.${RUN} +cp -R ../../share/run ./run.${RUN} +cd run.${RUN} +cp ../aether_${RUN}.json ./aether.json +mpirun -np ${PE} ./aether +../../../srcPython/postAether.py -rm +cd .. + +RUN=cube_sphere6 +PE=6 +rm -rf ./run.${RUN} +cp -R ../../share/run ./run.${RUN} +cd run.${RUN} +cp ../aether_${RUN}.json ./aether.json +mpirun -np ${PE} ./aether +../../../srcPython/postAether.py -rm +cd .. + +RUN=sphere4_dipole4 +PE=4 +rm -rf ./run.${RUN} +cp -R ../../share/run ./run.${RUN} +cd run.${RUN} +cp ../aether_${RUN}.json ./aether.json +mpirun -np ${PE} ./aether +../../../srcPython/postAether.py -rm +cd .. + +RUN=cube_dipole6 +PE=6 +rm -rf ./run.${RUN} +cp -R ../../share/run ./run.${RUN} +cd run.${RUN} +cp ../aether_${RUN}.json ./aether.json +mpirun -np ${PE} ./aether +../../../srcPython/postAether.py -rm +cd .. + +RUN=sphere6_dipole6 +PE=6 +rm -rf ./run.${RUN} +cp -R ../../share/run ./run.${RUN} +cd run.${RUN} +cp ../aether_${RUN}.json ./aether.json +mpirun -np ${PE} ./aether +../../../srcPython/postAether.py -rm +cd .. + diff --git a/tests/restart_cubesphere/aether.whole.json b/tests/restart_cubesphere/aether.whole.json index 8a70fef8..7e683f29 100644 --- a/tests/restart_cubesphere/aether.whole.json +++ b/tests/restart_cubesphere/aether.whole.json @@ -15,8 +15,6 @@ "Electrodynamics" : { "File" : "UA/inputs/b20110320n_omni.bin"}, - "OmniwebFiles" : ["UA/inputs/omni_20110319.txt"], - "CubeSphere" : { "is" : true}, diff --git a/tests/restart_ensembles/aether.json.whole b/tests/restart_ensembles/aether.json.whole index 5cb6453a..797526d5 100644 --- a/tests/restart_ensembles/aether.json.whole +++ b/tests/restart_ensembles/aether.json.whole @@ -1,7 +1,7 @@ { "Ensembles" : { - "nMembers" : 5}, + "nMembers" : 3}, "Perturb": { "f107" : { "Mean" : 1.0, @@ -17,23 +17,32 @@ "iVerbose" : 0}, "StartTime" : [2011, 3, 20, 0, 0, 0], - "EndTime" : [2011, 3, 20, 0, 10, 0], + "EndTime" : [2011, 3, 20, 0, 6, 0], + "Perturb": { + "f107" : { "Mean" : 1.0, + "Std" : 0.10, + "Add" : false, + "Constant" : true}}, + "neuGrid" : { - "nLons" : 12, - "nLats" : 12, - "nAlts" : 30}, - + "Shape": "cubesphere", + "nLonsPerBlock" : 20, + "nLatsPerBlock" : 20, + "nAlts" : 30, + "dAltScale" : 0.3, + "IsUniformAlt" : false}, + "ionGrid": { - "dAltStretch": 0.2, - "LatStretch": 1, - "Shape": "dipole", - "nLonsPerBlock": 18, - "nLatsPerBlock" : 18, - "nAlts":36, - "LatMax":88, - "MinAlt": 80.0, - "MinApex": 125.0}, + "Shape": "cubesphere", + "nLonsPerBlock": 16, + "nLatsPerBlock": 16, + "nAlts": 30, + "LatRange": [-90, 90], + "AltRange": [100.0, 1000], + "LonRange": [0.0, 360.0], + "dAltScale" : 0.3, + "IsUniformAlt" : false}, "InitialConditions" : { "type" : "msis"}, diff --git a/tests/restart_ensembles/run_all.sh b/tests/restart_ensembles/run_all.sh index 6d2c6d49..636549b8 100755 --- a/tests/restart_ensembles/run_all.sh +++ b/tests/restart_ensembles/run_all.sh @@ -1,11 +1,17 @@ #!/bin/sh -NPROC=5 +# cubesphere has 6 blocks: +NBLOCKS=6 +# run with 3 members: +NMEMBERS=2 +# run for a total of 180s TOTALTIME=180 +# this is the mpi command MPI=/usr/bin/mpirun +# stop this many times NTIMES=2 # include -dowhole to run the whole simulation as comparison: -../../srcPython/run_restarts.py -totaltime=${TOTALTIME} -mpi=${MPI} -rundir=../../share/run -ensembles=${NPROC} -restarts=${NTIMES} +../../srcPython/run_restarts.py -totaltime=${TOTALTIME} -mpi=${MPI} -rundir=../../share/run -ensembles=${NMEMBERS} -restarts=${NTIMES} -blocks=${NBLOCKS} diff --git a/tests/restarts/aether.first.json b/tests/restarts/aether.first.json index 9fec0780..d8ae5f16 100644 --- a/tests/restarts/aether.first.json +++ b/tests/restarts/aether.first.json @@ -6,21 +6,23 @@ "dt" : 10.0}, "neuGrid" : { + "Shape" : "sphere", "nLons" : 12, "nLats" : 12, - "nAlts" : 30}, - + "nAlts" : 30, + "dAltScale" : 0.3, + "IsUniformAlt" : false}, + "ionGrid": { - "dAltStretch": 0.2, - "LatStretch": 1, - "Shape": "dipole", + "Shape": "sphere", "nLonsPerBlock": 18, "nLatsPerBlock" : 18, "nAlts":36, - "LatMax":88, - "MinAlt": 80.0, - "MinApex": 125.0 - }, + "LatRange": [-90, 90], + "AltRange": [100.0, 1000], + "LonRange": [0.0, 360.0], + "dAltScale" : 0.3, + "IsUniformAlt" : false}, "StartTime" : [2011, 3, 20, 0, 0, 0], "EndTime" : [2011, 3, 20, 0, 5, 0], @@ -30,6 +32,12 @@ "DiffuseAurora" : "fta", "File": "UA/inputs/b20110320n_omni.bin"}, + "InitialConditions" : { + "type" : "msis"}, + + "BoundaryConditions" : { + "type" : "msis"}, + "OmniwebFiles" : ["UA/inputs/omni_20110319.txt"], "Outputs" : { @@ -38,6 +46,8 @@ "Restart" : { "do" : false, - "dt" : 900.0} + "dt" : 900.0}, + + "PlanetFile" : "UA/inputs/earth.in" } diff --git a/tests/restarts/aether.whole.json b/tests/restarts/aether.whole.json index 53b964d9..f421033d 100644 --- a/tests/restarts/aether.whole.json +++ b/tests/restarts/aether.whole.json @@ -6,21 +6,23 @@ "dt" : 10.0}, "neuGrid" : { + "Shape" : "sphere", "nLons" : 12, "nLats" : 12, - "nAlts" : 30}, - + "nAlts" : 30, + "dAltScale" : 0.3, + "IsUniformAlt" : false}, + "ionGrid": { - "dAltStretch": 0.2, - "LatStretch": 1, - "Shape": "dipole", + "Shape": "sphere", "nLonsPerBlock": 18, "nLatsPerBlock" : 18, "nAlts":36, - "LatMax":88, - "MinAlt": 80.0, - "MinApex": 125.0 - }, + "LatRange": [-90, 90], + "AltRange": [100.0, 1000], + "LonRange": [0.0, 360.0], + "dAltScale" : 0.3, + "IsUniformAlt" : false}, "StartTime" : [2011, 3, 20, 0, 0, 0], "EndTime" : [2011, 3, 20, 0, 10, 0], @@ -30,6 +32,12 @@ "DiffuseAurora" : "fta", "File": "UA/inputs/b20110320n_omni.bin"}, + "InitialConditions" : { + "type" : "msis"}, + + "BoundaryConditions" : { + "type" : "msis"}, + "OmniwebFiles" : ["UA/inputs/omni_20110319.txt"], "Outputs" : { @@ -38,6 +46,8 @@ "Restart" : { "do" : false, - "dt" : 900.0} + "dt" : 900.0}, + + "PlanetFile" : "UA/inputs/earth.in" }