From 0003b19b6bd15e0b5d05a0951af3281c0d9bf32f Mon Sep 17 00:00:00 2001 From: Aaron Ridley Date: Mon, 20 Nov 2023 07:38:12 -0500 Subject: [PATCH] STY: astyle --- src/advance.cpp | 1 + src/calc_neutral_derived.cpp | 22 ++--- src/electrodynamics.cpp | 151 ++++++++++++++++++----------------- src/exchange_messages_v2.cpp | 6 +- src/init_geo_grid.cpp | 2 + src/inputs.cpp | 8 +- src/neutrals.cpp | 7 +- src/solver_gradients.cpp | 7 +- src/tools.cpp | 17 ++-- src/transform.cpp | 19 +++-- 10 files changed, 135 insertions(+), 105 deletions(-) diff --git a/src/advance.cpp b/src/advance.cpp index 9abc23b2..c1269cd1 100644 --- a/src/advance.cpp +++ b/src/advance.cpp @@ -59,6 +59,7 @@ bool advance(Planets &planet, // first neutrals.calc_scale_height(gGrid); + if (didWork) didWork = neutrals.set_bcs(gGrid, time, indices); diff --git a/src/calc_neutral_derived.cpp b/src/calc_neutral_derived.cpp index da069cc3..32ca53fe 100644 --- a/src/calc_neutral_derived.cpp +++ b/src/calc_neutral_derived.cpp @@ -393,7 +393,7 @@ precision_t Neutrals::calc_dt_cubesphere(Grid grid) { int64_t nAlts = grid.get_nAlts(); int64_t nXs = grid.get_nLons(); int64_t nYs = grid.get_nLats(); - + // dtx dty for reference coordinate system arma_cube dtx(size(cMax_vcgc[0])); arma_cube dty(size(cMax_vcgc[0])); @@ -405,17 +405,17 @@ precision_t Neutrals::calc_dt_cubesphere(Grid grid) { for (int iAlt = 0; iAlt < nAlts; iAlt++) { // Conver cMax to contravariant velocity first arma_mat u1 = sqrt( - cMax_vcgc[0].slice(iAlt) % grid.A11_inv_scgc.slice(iAlt) % - cMax_vcgc[0].slice(iAlt) % grid.A11_inv_scgc.slice(iAlt) + - cMax_vcgc[1].slice(iAlt) % grid.A12_inv_scgc.slice(iAlt) % - cMax_vcgc[1].slice(iAlt) % grid.A12_inv_scgc.slice(iAlt)); + cMax_vcgc[0].slice(iAlt) % grid.A11_inv_scgc.slice(iAlt) % + cMax_vcgc[0].slice(iAlt) % grid.A11_inv_scgc.slice(iAlt) + + cMax_vcgc[1].slice(iAlt) % grid.A12_inv_scgc.slice(iAlt) % + cMax_vcgc[1].slice(iAlt) % grid.A12_inv_scgc.slice(iAlt)); arma_mat u2 = sqrt( - cMax_vcgc[0].slice(iAlt) % grid.A21_inv_scgc.slice(iAlt) % - cMax_vcgc[0].slice(iAlt) % grid.A21_inv_scgc.slice(iAlt) + - cMax_vcgc[1].slice(iAlt) % grid.A22_inv_scgc.slice(iAlt) % - cMax_vcgc[1].slice(iAlt) % grid.A22_inv_scgc.slice(iAlt)); - dtx.slice(iAlt) = grid.drefx(iAlt) * dummy_1 / u1; - dty.slice(iAlt) = grid.drefy(iAlt) * dummy_1 / u2; + cMax_vcgc[0].slice(iAlt) % grid.A21_inv_scgc.slice(iAlt) % + cMax_vcgc[0].slice(iAlt) % grid.A21_inv_scgc.slice(iAlt) + + cMax_vcgc[1].slice(iAlt) % grid.A22_inv_scgc.slice(iAlt) % + cMax_vcgc[1].slice(iAlt) % grid.A22_inv_scgc.slice(iAlt)); + dtx.slice(iAlt) = grid.drefx(iAlt) * dummy_1 / u1; + dty.slice(iAlt) = grid.drefy(iAlt) * dummy_1 / u2; } // simply some things, and just take the bulk value for now: diff --git a/src/electrodynamics.cpp b/src/electrodynamics.cpp index c52eee44..5b930da4 100644 --- a/src/electrodynamics.cpp +++ b/src/electrodynamics.cpp @@ -5,7 +5,7 @@ // ----------------------------------------------------------------------------- -// Call (fortran) ionospheric electrodynamics if we enable this +// Call (fortran) ionospheric electrodynamics if we enable this // in the cmake file // ----------------------------------------------------------------------------- @@ -27,36 +27,40 @@ Electrodynamics::Electrodynamics(Times time) { HaveFortranIe = false; isSet = false; - #ifdef FORTRAN - // Initialize IE components (reading in the data files): - std::string ieDir = input.get_electrodynamics_dir(); - int* ieDir_iArray = copy_string_to_int(ieDir); - std::string efield = input.get_potential_model(); - std::string aurora = input.get_diffuse_auroral_model(); +#ifdef FORTRAN + // Initialize IE components (reading in the data files): + std::string ieDir = input.get_electrodynamics_dir(); + int* ieDir_iArray = copy_string_to_int(ieDir); + std::string efield = input.get_potential_model(); + std::string aurora = input.get_diffuse_auroral_model(); - if (efield.length() == 0 & aurora.length() == 0) { - HaveFortranIe = false; + if (efield.length() == 0 & aurora.length() == 0) + HaveFortranIe = false; + + else { + int* efield_iArray = copy_string_to_int(efield); + int* aurora_iArray = copy_string_to_int(aurora); + ie_init_library(ieDir_iArray, efield_iArray, aurora_iArray, &iError); + + if (iError > 0) { + report.print(0, "Error in setting fortran IE!"); + IsOk = false; } else { - int* efield_iArray = copy_string_to_int(efield); - int* aurora_iArray = copy_string_to_int(aurora); - ie_init_library(ieDir_iArray, efield_iArray, aurora_iArray, &iError); - if (iError > 0) { - report.print(0,"Error in setting fortran IE!"); - IsOk = false; - } else { - HaveFortranIe = true; - isSet = true; - } + HaveFortranIe = true; + isSet = true; } - isCompiled = true; - #else - isCompiled = false; - #endif + } + + isCompiled = true; +#else + isCompiled = false; +#endif // If we don't set IE through Fortran, then try to set it through a file: if (!HaveFortranIe) { std::string electrodynamics_file = input.get_electrodynamics_file(); + if (electrodynamics_file.length() > 0 & electrodynamics_file != "none") { // This function sets HaveElectrodynamicsFile = true. @@ -83,7 +87,7 @@ Electrodynamics::Electrodynamics(Times time) { // ----------------------------------------------------------------------------- void Electrodynamics::set_all_indices_for_ie(Times time, - Indices &indices) { + Indices &indices) { std::string function = "Electrodynamics::set_all_indices_for_ie"; static int iFunction = -1; report.enter(function, iFunction); @@ -111,17 +115,18 @@ void Electrodynamics::set_all_indices_for_ie(Times time, std::cout << "sw v : " << iVx_ << " " << swv << "\n"; std::cout << "sw n : " << iN_ << " " << swn << "\n"; } - #ifdef FORTRAN - ie_set_time(&time_now); - ie_set_imfby(&imfby); - ie_set_imfbz(&imfbz); - ie_set_swv(&swv); - ie_set_swn(&swn); - ie_set_ae(&ae); - ie_set_au(&au); - ie_set_al(&al); - ie_set_hp_from_ae(&ae); - #endif + +#ifdef FORTRAN + ie_set_time(&time_now); + ie_set_imfby(&imfby); + ie_set_imfbz(&imfbz); + ie_set_swv(&swv); + ie_set_swn(&swn); + ie_set_ae(&ae); + ie_set_au(&au); + ie_set_al(&al); + ie_set_hp_from_ae(&ae); +#endif report.exit(function); return; @@ -154,60 +159,62 @@ bool Electrodynamics::update(Planets planet, ions.eflux.zeros(); ions.avee.ones(); - #ifdef FORTRAN - if (HaveFortranIe) { - report.print(3, "Using Fortran Electrodynamics!"); - set_all_indices_for_ie(time, indices); - - if (!IsAllocated) { - int nXs = gGrid.get_nX(); - ie_set_nxs(&nXs); - int nYs = gGrid.get_nY(); - ie_set_nys(&nYs); - int64_t iTotal = nXs * nYs; - mlt2d = static_cast(malloc(iTotal * sizeof(float))); - lat2d = static_cast(malloc(iTotal * sizeof(float))); - pot2d = static_cast(malloc(iTotal * sizeof(float))); - eflux2d = static_cast(malloc(iTotal * sizeof(float))); - avee2d = static_cast(malloc(iTotal * sizeof(float))); - IsAllocated = true; - } +#ifdef FORTRAN + + if (HaveFortranIe) { + report.print(3, "Using Fortran Electrodynamics!"); + set_all_indices_for_ie(time, indices); + + if (!IsAllocated) { + int nXs = gGrid.get_nX(); + ie_set_nxs(&nXs); + int nYs = gGrid.get_nY(); + ie_set_nys(&nYs); + int64_t iTotal = nXs * nYs; + mlt2d = static_cast(malloc(iTotal * sizeof(float))); + lat2d = static_cast(malloc(iTotal * sizeof(float))); + pot2d = static_cast(malloc(iTotal * sizeof(float))); + eflux2d = static_cast(malloc(iTotal * sizeof(float))); + avee2d = static_cast(malloc(iTotal * sizeof(float))); + IsAllocated = true; + } - int64_t nZs = gGrid.get_nZ(); - int64_t iZ; + int64_t nZs = gGrid.get_nZ(); + int64_t iZ; - int iError; + int iError; - for (iZ = 0; iZ < nZs; iZ++) { - copy_mat_to_array(gGrid.magLocalTime_scgc.slice(iZ), mlt2d, true); - copy_mat_to_array(gGrid.magLat_scgc.slice(iZ), lat2d, true); + for (iZ = 0; iZ < nZs; iZ++) { + copy_mat_to_array(gGrid.magLocalTime_scgc.slice(iZ), mlt2d, true); + copy_mat_to_array(gGrid.magLat_scgc.slice(iZ), lat2d, true); - ie_set_mlts(mlt2d, &iError); - ie_set_lats(lat2d, &iError); - ie_update_grid(&iError); + ie_set_mlts(mlt2d, &iError); + ie_set_lats(lat2d, &iError); + ie_update_grid(&iError); - ie_get_potential(pot2d, &iError); - copy_array_to_mat(pot2d, ions.potential_scgc.slice(iZ), true); + ie_get_potential(pot2d, &iError); + copy_array_to_mat(pot2d, ions.potential_scgc.slice(iZ), true); - if (iZ == nZs-1) { - ie_get_electron_diffuse_aurora(eflux2d, avee2d, &iError); - copy_array_to_mat(avee2d, ions.avee, true); - copy_array_to_mat(eflux2d, ions.eflux, true); - } + if (iZ == nZs - 1) { + ie_get_electron_diffuse_aurora(eflux2d, avee2d, &iError); + copy_array_to_mat(avee2d, ions.avee, true); + copy_array_to_mat(eflux2d, ions.eflux, true); } } - #endif + } + +#endif if (HaveElectrodynamicsFile) { report.print(3, "Setting electrodynamics from file!"); auto electrodynamics_values = get_electrodynamics(gGrid.magLat_scgc, - gGrid.magLocalTime_scgc); + gGrid.magLocalTime_scgc); ions.potential_scgc = std::get<0>(electrodynamics_values); ions.eflux = std::get<1>(electrodynamics_values); ions.avee = std::get<2>(electrodynamics_values); } - } + } report.exit(function); return true; diff --git a/src/exchange_messages_v2.cpp b/src/exchange_messages_v2.cpp index 39104f29..b6ceeaff 100644 --- a/src/exchange_messages_v2.cpp +++ b/src/exchange_messages_v2.cpp @@ -602,10 +602,12 @@ void Grid::exchange(arma_cube &data, const bool pole_inverse) { bool Neutrals::exchange(Grid &grid) { // For each species, exchange if its DoAdvect is true int64_t nGCs = grid.get_nGCs(); + for (int i = 0; i < nSpecies; ++i) { if (species[i].DoAdvect) grid.exchange(species[i].density_scgc, false); - fill_corners(species[i].density_scgc, nGCs); + + fill_corners(species[i].density_scgc, nGCs); } // Exchange temperature @@ -617,7 +619,9 @@ bool Neutrals::exchange(Grid &grid) { grid.exchange(velocity_vcgc[0], true); grid.exchange(velocity_vcgc[1], true); grid.exchange(velocity_vcgc[2], false); + for (int iDir = 0; iDir < 3; iDir++) fill_corners(velocity_vcgc[iDir], nGCs); + return true; } diff --git a/src/init_geo_grid.cpp b/src/init_geo_grid.cpp index 0a379f28..a904366d 100644 --- a/src/init_geo_grid.cpp +++ b/src/init_geo_grid.cpp @@ -984,10 +984,12 @@ bool Grid::init_geo_grid(Quadtree quadtree, // Calculate the radius (for spherical or non-spherical) fill_grid_radius(planet); + // Correct the reference grid with correct length scale: // (with R = actual radius) if (input.get_is_cubesphere()) correct_xy_grid(planet); + // Calculate grid spacing calc_grid_spacing(planet); //calculate radial unit vector (for spherical or oblate planet) diff --git a/src/inputs.cpp b/src/inputs.cpp index 186bed3a..3dc42517 100644 --- a/src/inputs.cpp +++ b/src/inputs.cpp @@ -461,7 +461,7 @@ bool Inputs::get_euv_douse() { } // ----------------------------------------------------------------------- -// Return the Electrodynamics Dir - this is where all of the +// Return the Electrodynamics Dir - this is where all of the // files that are for the empirical models reside // ----------------------------------------------------------------------- @@ -470,7 +470,7 @@ std::string Inputs::get_electrodynamics_north_file() { } // ----------------------------------------------------------------------- -// Return the Electrodynamics Dir - this is where all of the +// Return the Electrodynamics Dir - this is where all of the // files that are for the empirical models reside // ----------------------------------------------------------------------- @@ -479,7 +479,7 @@ std::string Inputs::get_electrodynamics_south_file() { } // ----------------------------------------------------------------------- -// Return the Electrodynamics Dir - this is where all of the +// Return the Electrodynamics Dir - this is where all of the // files that are for the empirical models reside // ----------------------------------------------------------------------- @@ -488,7 +488,7 @@ std::string Inputs::get_electrodynamics_file() { } // ----------------------------------------------------------------------- -// Return the Electrodynamics Dir - this is where all of the +// Return the Electrodynamics Dir - this is where all of the // files that are for the empirical models reside // ----------------------------------------------------------------------- diff --git a/src/neutrals.cpp b/src/neutrals.cpp index d9ff0b69..0525b237 100644 --- a/src/neutrals.cpp +++ b/src/neutrals.cpp @@ -290,22 +290,25 @@ bool Neutrals::check_for_nonfinites() { bool didWork = true; isBad = !all_finite(density_scgc, "density_scgc"); + if (isBad) { report.error("non-finite found in neutral density!"); didWork = false; } isBad = !all_finite(temperature_scgc, "temperature_scgc"); + if (isBad) { report.error("non-finite found in neutral temperature!"); didWork = false; - } + } isBad = !all_finite(velocity_vcgc, "velocity_vcgc"); + if (isBad) { report.error("non-finite found in neutral velocity!"); didWork = false; - } + } return didWork; } diff --git a/src/solver_gradients.cpp b/src/solver_gradients.cpp index 0ba729fd..f2bf0cc9 100644 --- a/src/solver_gradients.cpp +++ b/src/solver_gradients.cpp @@ -11,13 +11,15 @@ std::vector calc_gradient_vector(arma_cube value_scgc, Grid grid) { std::vector gradient_vcgc; - if (input.get_is_cubesphere()) { + if (input.get_is_cubesphere()) gradient_vcgc = calc_gradient_cubesphere(value_scgc, grid); - } else { + + else { gradient_vcgc.push_back(calc_gradient_lon(value_scgc, grid)); gradient_vcgc.push_back(calc_gradient_lat(value_scgc, grid)); gradient_vcgc.push_back(calc_gradient_alt(value_scgc, grid)); } + return gradient_vcgc; } @@ -215,6 +217,7 @@ std::vector calc_gradient_cubesphere(arma_cube value, Grid grid) { } } } + // We then use A transformation matrices to convert grad_xy to grad_latlon // Ref -> Physical, we use A matrix grad_lon.slice(iAlt) = grad_x_curr % grid.A11_inv_scgc.slice( diff --git a/src/tools.cpp b/src/tools.cpp index 70610262..973cd234 100644 --- a/src/tools.cpp +++ b/src/tools.cpp @@ -18,22 +18,23 @@ void fill_corners(arma_cube &values, int64_t nGCs) { for (iGCy = 0; iGCy < nGCs; iGCy++) { // lower left: values.tube(iGCx, iGCy) = 0.5 * ( - values.tube(iGCx, nGCs) + - values.tube(nGCs, iGCy)); + values.tube(iGCx, nGCs) + + values.tube(nGCs, iGCy)); // lower right: values.tube(nXs - iGCx - 1, iGCy) = 0.5 * ( - values.tube(nXs - iGCx - 1, nGCs) + - values.tube(nXs - nGCs - 1, iGCy)); + values.tube(nXs - iGCx - 1, nGCs) + + values.tube(nXs - nGCs - 1, iGCy)); // upper left: values.tube(iGCx, nYs - iGCy - 1) = 0.5 * ( - values.tube(iGCx, nYs - nGCs - 1) + - values.tube(nGCs, nYs - iGCy - 1)); + values.tube(iGCx, nYs - nGCs - 1) + + values.tube(nGCs, nYs - iGCy - 1)); // upper right: values.tube(nXs - iGCx - 1, nYs - iGCy - 1) = 0.5 * ( - values.tube(nXs - iGCx - 1, nYs - nGCs - 1) + - values.tube(nXs - nGCs - 1, nYs - iGCy - 1)); + values.tube(nXs - iGCx - 1, nYs - nGCs - 1) + + values.tube(nXs - nGCs - 1, nYs - iGCy - 1)); } } + return; } diff --git a/src/transform.cpp b/src/transform.cpp index d4a657e0..b18a6fab 100644 --- a/src/transform.cpp +++ b/src/transform.cpp @@ -13,8 +13,10 @@ std::string mklower(std::string inString) { std::string outString = inString; int64_t nChars = outString.length(); + for (int64_t iChar = 0; iChar < nChars; iChar++) outString[iChar] = tolower(outString[iChar]); + return outString; } @@ -25,13 +27,16 @@ std::string mklower(std::string inString) { // ----------------------------------------------------------------------- int* copy_string_to_int(std::string inString) { - const int length = inString.length(); + const int length = inString.length(); // declaring character array - int* outArray = new int[400]; + int* outArray = new int[400]; + for (int i = 0; i < length; i++) outArray[i] = inString[i]; + for (int i = length; i < 400; i++) outArray[i] = 0; + return outArray; } @@ -75,8 +80,8 @@ void copy_cube_to_array(arma_cube cube_in, // ----------------------------------------------------------------------- void copy_mat_to_array(arma_mat mat_in, - float *array_out, - bool isFortran) { + float *array_out, + bool isFortran) { int64_t nX = mat_in.n_rows; int64_t nY = mat_in.n_cols; @@ -88,14 +93,16 @@ void copy_mat_to_array(arma_mat mat_in, index = iY * nX + iX; else index = iX * nY + iY; + array_out[index] = mat_in(iX, iY); } } + return; } // ----------------------------------------------------------------------- -// copy from a 2d c-native array to +// copy from a 2d c-native array to // an already defined (!!) armidillo matrix // - the code uses the dimension of the matrix to figure out // the size of the array, so be careful!!! There is no checking! @@ -116,9 +123,11 @@ void copy_array_to_mat(float *array_in, index = iY * nX + iX; else index = iX * nY + iY; + mat_out(iX, iY) = array_in[index]; } } + return; }