diff --git a/src/advance.cpp b/src/advance.cpp index 5d59c03f..ff78c827 100644 --- a/src/advance.cpp +++ b/src/advance.cpp @@ -11,18 +11,18 @@ // ----------------------------------------------------------------------------- bool advance(Planets &planet, - Grid &gGrid, - Times &time, - Euv &euv, - Neutrals &neutrals, - Ions &ions, - Chemistry &chemistry, - Electrodynamics &electrodynamics, - Indices &indices, - Logfile &logfile) { + Grid &gGrid, + Times &time, + Euv &euv, + Neutrals &neutrals, + Ions &ions, + Chemistry &chemistry, + Electrodynamics &electrodynamics, + Indices &indices, + Logfile &logfile) { bool didWork = true; - + std::string function = "advance"; static int iFunction = -1; report.enter(function, iFunction); @@ -66,19 +66,19 @@ bool advance(Planets &planet, if (didWork) didWork = calc_euv(planet, - gGrid, - time, - euv, - neutrals, - ions, - indices); + gGrid, + time, + euv, + neutrals, + ions, + indices); if (didWork) didWork = electrodynamics.update(planet, - gGrid, - time, - ions); - + gGrid, + time, + ions); + if (didWork) { calc_ion_neutral_coll_freq(neutrals, ions); ions.calc_ion_drift(neutrals, gGrid, time.get_dt()); @@ -111,7 +111,7 @@ bool advance(Planets &planet, neutrals.exchange(gGrid); else neutrals.exchange_old(gGrid); - + time.increment_time(); if (time.check_time_gate(input.get_dt_write_restarts())) { @@ -130,7 +130,7 @@ bool advance(Planets &planet, if (!didWork) report.error("Error in Advance!"); - + report.exit(function); return didWork; } diff --git a/src/calc_euv.cpp b/src/calc_euv.cpp index c4e400ed..09be872e 100644 --- a/src/calc_euv.cpp +++ b/src/calc_euv.cpp @@ -14,25 +14,25 @@ // ----------------------------------------------------------------------------- bool calc_euv(Planets planet, - Grid grid, - Times time, - Euv &euv, - Neutrals &neutrals, - Ions &ions, - Indices indices) { + Grid grid, + Times time, + Euv &euv, + Neutrals &neutrals, + Ions &ions, + Indices indices) { bool didWork = true; if (!time.check_time_gate(input.get_dt_euv())) return true; - + std::string function = "Euv::calc_euv"; static int iFunction = -1; report.enter(function, iFunction); if (input.get_is_student()) report.print(-1, "(2) What function is this " + - input.get_student_name() + "?"); + input.get_student_name() + "?"); // Chapman integrals for EUV energy deposition: // Need chapman integrals for aurora too! @@ -42,6 +42,7 @@ bool calc_euv(Planets planet, // set didWork to false in order to catch bad euv models: didWork = false; std::string euvModel = mklower(input.get_euv_model()); + if (euvModel == "euvac") didWork = euv.euvac(time, indices); else if (euvModel == "neuvac") @@ -51,13 +52,14 @@ bool calc_euv(Planets planet, if (didWork) euv.scale_from_1au(planet, time); - + calc_ionization_heating(euv, neutrals, ions); } else neutrals.heating_euv_scgc.zeros(); if (!didWork) - report.error("Error in calc_euv! Check euv models."); + report.error("Error in calc_euv! Check euv models."); + report.exit(function); return didWork; } diff --git a/src/calc_neutral_derived.cpp b/src/calc_neutral_derived.cpp index b8d28423..9e1c64a4 100644 --- a/src/calc_neutral_derived.cpp +++ b/src/calc_neutral_derived.cpp @@ -349,7 +349,7 @@ precision_t Neutrals::calc_dt(Grid grid) { int64_t nAlts = grid.get_nAlts(); int64_t nXs = grid.get_nLons(); int64_t nYs = grid.get_nLats(); - + // dtx dty for reference coordinate system arma_cube dtx(size(cMax_vcgc[0])); arma_cube dty(size(cMax_vcgc[0])); @@ -360,15 +360,16 @@ precision_t Neutrals::calc_dt(Grid grid) { // Loop through altitudes for (int iAlt = 0; iAlt < nAlts; iAlt++) { // Conver cMax to contravariant velocity first - arma_mat u1 = cMax_vcgc[0].slice(iAlt) % grid.A11_inv_scgc.slice(iAlt) - + cMax_vcgc[1].slice(iAlt) % grid.A12_inv_scgc.slice(iAlt); - arma_mat u2 = cMax_vcgc[0].slice(iAlt) % grid.A21_inv_scgc.slice(iAlt) - + cMax_vcgc[1].slice(iAlt) % grid.A22_inv_scgc.slice(iAlt); + arma_mat u1 = cMax_vcgc[0].slice(iAlt) % grid.A11_inv_scgc.slice(iAlt) + + cMax_vcgc[1].slice(iAlt) % grid.A12_inv_scgc.slice(iAlt); + arma_mat u2 = cMax_vcgc[0].slice(iAlt) % grid.A21_inv_scgc.slice(iAlt) + + cMax_vcgc[1].slice(iAlt) % grid.A22_inv_scgc.slice(iAlt); // Extract a scalar dx and divide by contravariant velocity - dtx.slice(iAlt) = grid.drefx(iAlt)*dummy_1 / u1; - dty.slice(iAlt) = grid.drefy(iAlt)*dummy_1 / u2; + dtx.slice(iAlt) = grid.drefx(iAlt) * dummy_1 / u1; + dty.slice(iAlt) = grid.drefy(iAlt) * dummy_1 / u2; } + dta(0) = dtx.min(); dta(1) = dty.min(); } else { @@ -378,7 +379,7 @@ precision_t Neutrals::calc_dt(Grid grid) { arma_cube dty = grid.dlat_center_dist_scgc / cMax_vcgc[1]; dta(1) = dty.min(); } - + if (input.get_nAltsGeo() > 1) { arma_cube dtz = grid.dalt_center_scgc / cMax_vcgc[2]; dta(2) = dtz.min(); diff --git a/src/electrodynamics.cpp b/src/electrodynamics.cpp index 14d757fa..476a34f0 100644 --- a/src/electrodynamics.cpp +++ b/src/electrodynamics.cpp @@ -34,9 +34,9 @@ Electrodynamics::Electrodynamics(Times time) { // ----------------------------------------------------------------------------- bool Electrodynamics::update(Planets planet, - Grid gGrid, - Times time, - Ions &ions) { + Grid gGrid, + Times time, + Ions &ions) { std::string function = "Electrodynamics::update"; static int iFunction = -1; @@ -183,9 +183,9 @@ arma_mat Electrodynamics::get_values(arma_mat matToInterpolateOn, // ----------------------------------------------------------------------------- std::tuple Electrodynamics::get_electrodynamics(arma_cube magLat, - arma_cube magLocalTime) { + arma_mat, + arma_mat> Electrodynamics::get_electrodynamics(arma_cube magLat, +arma_cube magLocalTime) { arma_cube pot; arma_mat eflux; arma_mat avee; diff --git a/src/euv.cpp b/src/euv.cpp index 0084eaea..c4a4069a 100644 --- a/src/euv.cpp +++ b/src/euv.cpp @@ -313,7 +313,7 @@ bool Euv::pair_euv(Neutrals &neutrals, // -------------------------------------------------------------------------- void Euv::scale_from_1au(Planets planet, - Times time) { + Times time) { precision_t d = planet.get_star_to_planet_dist(time); precision_t scale = 1.0 / (d * d); @@ -331,7 +331,7 @@ void Euv::scale_from_1au(Planets planet, // -------------------------------------------------------------------------- bool Euv::euvac(Times time, - Indices indices) { + Indices indices) { bool didWork = true; precision_t slope; @@ -375,7 +375,7 @@ bool Euv::euvac(Times time, // -------------------------------------------------------------------------- bool Euv::neuvac(Times time, - Indices indices) { + Indices indices) { bool didWork = true; precision_t slope; @@ -419,7 +419,7 @@ bool Euv::neuvac(Times time, // -------------------------------------------------------------------------- bool Euv::solomon_hfg(Times time, - Indices indices) { + Indices indices) { std::string function = "Euv::solomon_hfg"; static int iFunction = -1; diff --git a/src/exchange_messages_v2.cpp b/src/exchange_messages_v2.cpp index d64e9c73..220b4bf5 100644 --- a/src/exchange_messages_v2.cpp +++ b/src/exchange_messages_v2.cpp @@ -20,35 +20,35 @@ int pos_iProc_surface(precision_t oLR, precision_t LR, precision_t DU, int numProcs) { - // The iProc of the lower left corner - int retProc = 0; - - // Repeat until there is only one processor in the described area - while (numProcs != 1) { - // Divide into 4 parts - dLR /= 2; - dDU /= 2; - - if (LR >= oLR + dLR) { - // The point is on the right, and the iProc of the lower left - // corner where the point is located is [1/4] or [3/4] - // * remainProc, depending on DU - retProc += numProcs / 4; - // Update the origin - oLR += dLR; - } - - if (DU >= oDU + dDU) { - // The point is above - retProc += numProcs / 2; - // Update the origin - oDU += dDU; - } + // The iProc of the lower left corner + int retProc = 0; + + // Repeat until there is only one processor in the described area + while (numProcs != 1) { + // Divide into 4 parts + dLR /= 2; + dDU /= 2; + + if (LR >= oLR + dLR) { + // The point is on the right, and the iProc of the lower left + // corner where the point is located is [1/4] or [3/4] + // * remainProc, depending on DU + retProc += numProcs / 4; + // Update the origin + oLR += dLR; + } - numProcs /= 4; + if (DU >= oDU + dDU) { + // The point is above + retProc += numProcs / 2; + // Update the origin + oDU += dDU; } - return retProc; + numProcs /= 4; + } + + return retProc; } // -------------------------------------------------------------------------- @@ -56,13 +56,13 @@ int pos_iProc_surface(precision_t oLR, // -------------------------------------------------------------------------- int pos_iProc_sphere(precision_t lon_in, precision_t lat_in) { - return pos_iProc_surface(Sphere::ORIGINS(0) * cPI, - Sphere::ORIGINS(1) * cPI, - Sphere::RIGHTS(0) * cPI, - Sphere::UPS(1) * cPI, - lon_in, - lat_in, - nGrids); + return pos_iProc_surface(Sphere::ORIGINS(0) * cPI, + Sphere::ORIGINS(1) * cPI, + Sphere::RIGHTS(0) * cPI, + Sphere::UPS(1) * cPI, + lon_in, + lat_in, + nGrids); } // -------------------------------------------------------------------------- @@ -70,11 +70,11 @@ int pos_iProc_sphere(precision_t lon_in, precision_t lat_in) { // -------------------------------------------------------------------------- bool inRange(precision_t a, precision_t b, precision_t c) { - if (b > 0) { - return c >= a && c < a + b; - } else { - return c <= a && c > a + b; - } + if (b > 0) + return c >= a && c < a + b; + + else + return c <= a && c > a + b; } // -------------------------------------------------------------------------- @@ -82,60 +82,61 @@ bool inRange(precision_t a, precision_t b, precision_t c) { // -------------------------------------------------------------------------- int pos_iProc_cube(precision_t lon_in, precision_t lat_in) { - // Transfer polar coordinate to cartesian coordinate - arma_vec point = sphere_to_cube(lon_in, lat_in); - - // Store the index of surface, LR and DU - // INITIAL VALUES ARE HARD-CODED BASED ON SURFACE 5(6) TO - // SOLVE BOUNDARY PROBLEMS - int iFace = 0, iLR = 1, iDU = 0; - bool on_surface = false; - // Find the surface where the point is located - for (;iFace < 6 && !on_surface; ++iFace) { - // Assume that the point is on current surface - on_surface = true; - int sLR = -1, sDU = -1; - - // Check three dimensions - for (int j = 0; j < 3 && on_surface; ++j) { - if (CubeSphere::RIGHTS(iFace, j) != 0) { - sLR = j; - on_surface = inRange(CubeSphere::ORIGINS(iFace, j), - CubeSphere::RIGHTS(iFace, j), - point(j)); - } else if (CubeSphere::UPS(iFace, j) != 0) { - sDU = j; - on_surface = inRange(CubeSphere::ORIGINS(iFace, j), - CubeSphere::UPS(iFace, j), - point(j)); - } else { - on_surface = (CubeSphere::ORIGINS(iFace, j) == point(j)); - } - } + // Transfer polar coordinate to cartesian coordinate + arma_vec point = sphere_to_cube(lon_in, lat_in); + + // Store the index of surface, LR and DU + // INITIAL VALUES ARE HARD-CODED BASED ON SURFACE 5(6) TO + // SOLVE BOUNDARY PROBLEMS + int iFace = 0, iLR = 1, iDU = 0; + bool on_surface = false; + + // Find the surface where the point is located + for (; iFace < 6 && !on_surface; ++iFace) { + // Assume that the point is on current surface + on_surface = true; + int sLR = -1, sDU = -1; + + // Check three dimensions + for (int j = 0; j < 3 && on_surface; ++j) { + if (CubeSphere::RIGHTS(iFace, j) != 0) { + sLR = j; + on_surface = inRange(CubeSphere::ORIGINS(iFace, j), + CubeSphere::RIGHTS(iFace, j), + point(j)); + } else if (CubeSphere::UPS(iFace, j) != 0) { + sDU = j; + on_surface = inRange(CubeSphere::ORIGINS(iFace, j), + CubeSphere::UPS(iFace, j), + point(j)); + } else + on_surface = (CubeSphere::ORIGINS(iFace, j) == point(j)); + } - // Assign iLR and iDU if the point is indeed on current surface - if (on_surface) { - iLR = sLR; - iDU = sDU; - } + // Assign iLR and iDU if the point is indeed on current surface + if (on_surface) { + iLR = sLR; + iDU = sDU; } - // Minus 1 due to for loop mechanism - --iFace; - - // Change the sign if delta is negative - precision_t coefLR, coefDU; - coefLR = CubeSphere::RIGHTS(iFace, iLR) > 0 ? 1.0 : -1.0; - coefDU = CubeSphere::UPS(iFace, iDU) > 0 ? 1.0 : -1.0; - // Compute iProc in the surface - int ip = pos_iProc_surface(CubeSphere::ORIGINS(iFace, iLR) * coefLR, - CubeSphere::ORIGINS(iFace, iDU) * coefDU, - CubeSphere::RIGHTS(iFace, iLR) * coefLR, - CubeSphere::UPS(iFace, iDU) * coefDU, - point(iLR) * coefLR, - point(iDU) * coefDU, - nGrids / 6); - // Consider surface number - return ip + iFace * nGrids / 6; + } + + // Minus 1 due to for loop mechanism + --iFace; + + // Change the sign if delta is negative + precision_t coefLR, coefDU; + coefLR = CubeSphere::RIGHTS(iFace, iLR) > 0 ? 1.0 : -1.0; + coefDU = CubeSphere::UPS(iFace, iDU) > 0 ? 1.0 : -1.0; + // Compute iProc in the surface + int ip = pos_iProc_surface(CubeSphere::ORIGINS(iFace, iLR) * coefLR, + CubeSphere::ORIGINS(iFace, iDU) * coefDU, + CubeSphere::RIGHTS(iFace, iLR) * coefLR, + CubeSphere::UPS(iFace, iDU) * coefDU, + point(iLR) * coefLR, + point(iDU) * coefDU, + nGrids / 6); + // Consider surface number + return ip + iFace * nGrids / 6; } /******** Functions to find iproc ends ********/ @@ -149,16 +150,16 @@ int pos_iProc_cube(precision_t lon_in, precision_t lat_in) { template struct mpi_msg_t { - // The message to send - std::vector buf; - // The processor to send - int rank; - // Default constructor - mpi_msg_t() - : buf(), rank() {} - // Move constructor - mpi_msg_t(std::vector &&_buf, const int _rank) - : buf(std::move(_buf)), rank(_rank) {} + // The message to send + std::vector buf; + // The processor to send + int rank; + // Default constructor + mpi_msg_t() + : buf(), rank() {} + // Move constructor + mpi_msg_t(std::vector &&_buf, const int _rank) + : buf(std::move(_buf)), rank(_rank) {} }; // -------------------------------------------------------------------------- @@ -168,15 +169,17 @@ struct mpi_msg_t { template std::ostream& operator<<(std::ostream &out, const std::vector> &msg) { - for (auto& it : msg) { - out << "\tTo " << it.rank << ":\n"; - for (auto& it2 : it.buf) { - out << "\t\t" << it2 << '\n'; - } - out << '\n'; - } + for (auto& it : msg) { + out << "\tTo " << it.rank << ":\n"; + + for (auto& it2 : it.buf) + out << "\t\t" << it2 << '\n'; + out << '\n'; - return out; + } + + out << '\n'; + return out; } /** @@ -194,146 +197,158 @@ template bool CUSTOM_MPI_SENDRECV(const std::vector> &msg_in, std::vector> &msg_out, const MPI_Comm &comm) { - // Promise previous send and receive have all ended - MPI_Barrier(comm); - - // Clear the output variables - msg_out.clear(); - - // Get the rank the size of communicator - // Processor with rank 0 checks whether all ends - const int num = msg_in.size(); - int comm_rank, comm_size; - MPI_Comm_rank(comm, &comm_rank); - MPI_Comm_size(comm, &comm_size); - - // Define constants - const int TAG_RESERVED = 3; - const char ALL_FINISH_TAG = 0; - const char PROC_FINISH_TAG = 1; - const char ACK_RECV_TAG = 2; - - // Only used to check correctness - std::vector request; - - // Send - for (int i = 0; i < num; ++i) { - request.emplace_back(); - MPI_Isend(msg_in[i].buf.data(), - msg_in[i].buf.size() * sizeof(T), - MPI_BYTE, - msg_in[i].rank, - i + TAG_RESERVED, - comm, - &request.back()); + // Promise previous send and receive have all ended + MPI_Barrier(comm); + + // Clear the output variables + msg_out.clear(); + + // Get the rank the size of communicator + // Processor with rank 0 checks whether all ends + const int num = msg_in.size(); + int comm_rank, comm_size; + MPI_Comm_rank(comm, &comm_rank); + MPI_Comm_size(comm, &comm_size); + + // Define constants + const int TAG_RESERVED = 3; + const char ALL_FINISH_TAG = 0; + const char PROC_FINISH_TAG = 1; + const char ACK_RECV_TAG = 2; + + // Only used to check correctness + std::vector request; + + // Send + for (int i = 0; i < num; ++i) { + request.emplace_back(); + MPI_Isend(msg_in[i].buf.data(), + msg_in[i].buf.size() * sizeof(T), + MPI_BYTE, + msg_in[i].rank, + i + TAG_RESERVED, + comm, + &request.back()); + } + + // Receive + // Record how many messages are sent successfully and + // how many processors finish sending messages + // finished_proc is only used by iProc 0 + int success_sent = 0, finished_proc = 0; + bool done = false; + + // If nothing needs to send, then the processor already finish sending + if (!num) { + if (!comm_rank) { + // Increase finished_proc if iProc = 0 + ++finished_proc; + } else { + // Send signal to iProc 0 otherwise + request.emplace_back(); + MPI_Isend(NULL, 0, MPI_BYTE, 0, PROC_FINISH_TAG, + comm, &request.back()); } - - // Receive - // Record how many messages are sent successfully and - // how many processors finish sending messages - // finished_proc is only used by iProc 0 - int success_sent = 0, finished_proc = 0; - bool done = false; - // If nothing needs to send, then the processor already finish sending - if (!num) { - if (!comm_rank) { - // Increase finished_proc if iProc = 0 + } + + while (!done) { + MPI_Status status; + int count; + // Probe the size of next message + MPI_Probe(MPI_ANY_SOURCE, MPI_ANY_TAG, comm, &status); + MPI_Get_count(&status, MPI_BYTE, &count); + + if (count) { + // Receiving actual messages + mpi_msg_t msg; + // Set rank and count, create space to receive message + msg.buf.resize(count / sizeof(T)); + msg.rank = status.MPI_SOURCE; + MPI_Recv(msg.buf.data(), + count, + MPI_BYTE, + status.MPI_SOURCE, + status.MPI_TAG, + comm, + MPI_STATUS_IGNORE); + msg_out.push_back(msg); + // Send ack message to the origin of this message + request.emplace_back(); + MPI_Isend(NULL, 0, MPI_BYTE, status.MPI_SOURCE, + ACK_RECV_TAG, comm, &request.back()); + } else { + // The message is a control signal with size = 0 + MPI_Recv(NULL, 0, MPI_BYTE, status.MPI_SOURCE, + status.MPI_TAG, comm, &status); + + switch (status.MPI_TAG) { + case ACK_RECV_TAG: { + // Increase the success_sent + ++success_sent; + + if (success_sent == num) { + // The processor finish sending + if (!comm_rank) { + // Increase finished_proc for rank 0 ++finished_proc; - } else { - // Send signal to iProc 0 otherwise + } else { + // Send proc finish to rank 0 for other processors request.emplace_back(); MPI_Isend(NULL, 0, MPI_BYTE, 0, PROC_FINISH_TAG, comm, &request.back()); + } } - } - while (!done) { - MPI_Status status; - int count; - // Probe the size of next message - MPI_Probe(MPI_ANY_SOURCE, MPI_ANY_TAG, comm, &status); - MPI_Get_count(&status, MPI_BYTE, &count); - if (count) { - // Receiving actual messages - mpi_msg_t msg; - // Set rank and count, create space to receive message - msg.buf.resize(count / sizeof(T)); - msg.rank = status.MPI_SOURCE; - MPI_Recv(msg.buf.data(), - count, - MPI_BYTE, - status.MPI_SOURCE, - status.MPI_TAG, - comm, - MPI_STATUS_IGNORE); - msg_out.push_back(msg); - // Send ack message to the origin of this message - request.emplace_back(); - MPI_Isend(NULL, 0, MPI_BYTE, status.MPI_SOURCE, - ACK_RECV_TAG, comm, &request.back()); - } else { - // The message is a control signal with size = 0 - MPI_Recv(NULL, 0, MPI_BYTE, status.MPI_SOURCE, - status.MPI_TAG, comm, &status); - switch (status.MPI_TAG) { - case ACK_RECV_TAG: { - // Increase the success_sent - ++success_sent; - if (success_sent == num) { - // The processor finish sending - if (!comm_rank) { - // Increase finished_proc for rank 0 - ++finished_proc; - } else { - // Send proc finish to rank 0 for other processors - request.emplace_back(); - MPI_Isend(NULL, 0, MPI_BYTE, 0, PROC_FINISH_TAG, - comm, &request.back()); - } - } - break; - } - case PROC_FINISH_TAG: { - // Only processor with rank 0 will get this tag - ++finished_proc; - break; - } - case ALL_FINISH_TAG: { - // Exit the while loop - done = true; - break; - } - default: { - // This should never be executed - return false; - } - } - // Check whether everything is done by processor with rank 0 - if (finished_proc == comm_size) { - // Send all finish to all other processors - for (int i = 1; i < comm_size; ++i) { - request.emplace_back(); - MPI_Isend(NULL, 0, MPI_BYTE, i, ALL_FINISH_TAG, - comm, &request.back()); - } - done = true; - } + + break; + } + + case PROC_FINISH_TAG: { + // Only processor with rank 0 will get this tag + ++finished_proc; + break; + } + + case ALL_FINISH_TAG: { + // Exit the while loop + done = true; + break; + } + + default: { + // This should never be executed + return false; + } + } + + // Check whether everything is done by processor with rank 0 + if (finished_proc == comm_size) { + // Send all finish to all other processors + for (int i = 1; i < comm_size; ++i) { + request.emplace_back(); + MPI_Isend(NULL, 0, MPI_BYTE, i, ALL_FINISH_TAG, + comm, &request.back()); } + + done = true; + } } + } - // Correctness check using request - for (auto& it : request) { - int flag; - MPI_Test(&it, &flag, MPI_STATUS_IGNORE); - if (!flag) { - // This should never be executed - return false; - } + // Correctness check using request + for (auto& it : request) { + int flag; + MPI_Test(&it, &flag, MPI_STATUS_IGNORE); + + if (!flag) { + // This should never be executed + return false; } + } - // Promise all send and receive end - MPI_Barrier(comm); + // Promise all send and receive end + MPI_Barrier(comm); - return true; + return true; } /** @@ -350,49 +365,50 @@ template bool CUSTOM_MPI_SENDRECV2(const std::vector> &msg_in, std::vector> &msg_out, const MPI_Comm &comm) { - // Promise previous send and receive have all ended - MPI_Barrier(comm); - - // Only used to check correctness - std::vector request; - - // Send - for (auto& it : msg_in) { - request.emplace_back(); - MPI_Isend(it.buf.data(), - it.buf.size() * sizeof(T), - MPI_BYTE, - it.rank, - 0, - comm, - &request.back()); + // Promise previous send and receive have all ended + MPI_Barrier(comm); + + // Only used to check correctness + std::vector request; + + // Send + for (auto& it : msg_in) { + request.emplace_back(); + MPI_Isend(it.buf.data(), + it.buf.size() * sizeof(T), + MPI_BYTE, + it.rank, + 0, + comm, + &request.back()); + } + + // Receive + for (auto& it : msg_out) { + MPI_Recv(it.buf.data(), + it.buf.size() * sizeof(T), + MPI_BYTE, + it.rank, + 0, + comm, + MPI_STATUS_IGNORE); + } + + // Promise all send and receive end + MPI_Barrier(comm); + + // Correctness check using request + for (auto& it : request) { + int flag; + MPI_Test(&it, &flag, MPI_STATUS_IGNORE); + + if (!flag) { + // This should never be executed + return false; } + } - // Receive - for (auto& it : msg_out) { - MPI_Recv(it.buf.data(), - it.buf.size() * sizeof(T), - MPI_BYTE, - it.rank, - 0, - comm, - MPI_STATUS_IGNORE); - } - - // Promise all send and receive end - MPI_Barrier(comm); - - // Correctness check using request - for (auto& it : request) { - int flag; - MPI_Test(&it, &flag, MPI_STATUS_IGNORE); - if (!flag) { - // This should never be executed - return false; - } - } - - return true; + return true; } /******** Functions related to MPI ends ********/ @@ -402,8 +418,8 @@ bool CUSTOM_MPI_SENDRECV2(const std::vector> &msg_in, // -------------------------------------------------------------------------- struct lon_lat_t { - precision_t lon; - precision_t lat; + precision_t lon; + precision_t lat; }; // -------------------------------------------------------------------------- @@ -411,8 +427,8 @@ struct lon_lat_t { // -------------------------------------------------------------------------- std::ostream& operator<<(std::ostream &out, const lon_lat_t &lonlat) { - out << "lon = " << lonlat.lon << "\tlat = " << lonlat.lat << '\n'; - return out; + out << "lon = " << lonlat.lon << "\tlat = " << lonlat.lat << '\n'; + return out; } // -------------------------------------------------------------------------- @@ -421,23 +437,25 @@ std::ostream& operator<<(std::ostream &out, const lon_lat_t &lonlat) { // -------------------------------------------------------------------------- int wrap_point_sphere(precision_t &lon, precision_t &lat) { - int inverse = 1; - if (lat < -0.5 * cPI) { - lat = -cPI - lat; - inverse = -1; - lon += cPI; - } else if (lat > 0.5 * cPI) { - lat = cPI - lat; - inverse = -1; - lon += cPI; - } - while (lon < 0) { - lon += cTWOPI; - } - while (lon > cTWOPI) { - lon -= cTWOPI; - } - return inverse; + int inverse = 1; + + if (lat < -0.5 * cPI) { + lat = -cPI - lat; + inverse = -1; + lon += cPI; + } else if (lat > 0.5 * cPI) { + lat = cPI - lat; + inverse = -1; + lon += cPI; + } + + while (lon < 0) + lon += cTWOPI; + + while (lon > cTWOPI) + lon -= cTWOPI; + + return inverse; } // -------------------------------------------------------------------------- @@ -445,85 +463,92 @@ int wrap_point_sphere(precision_t &lon, precision_t &lat) { // -------------------------------------------------------------------------- void Grid::init_connection() { - // Initialize communicator. Message exchange is limited in each - // ensemble and rank 0 in each ensemble has special purpose - MPI_Comm_split(aether_comm, iMember, iGrid, &grid_comm); - - // Classify the {lon, lat} for all ghost cells based on destination - std::unordered_map> lonlatmap; - for (int64_t iLon = 0; iLon < nLons; ++iLon) { - for (int64_t iLat = 0; iLat < nLats; ++iLat) { - // Read the lon, lat of ghost cell - precision_t lon = geoLon_scgc(iLon, iLat, 0); - precision_t lat = geoLat_scgc(iLon, iLat, 0); - // Find which processor can handle it and whether - // the message crosses the north/south pole - int dest, inverse = 1; - if (IsCubeSphereGrid) { - dest = pos_iProc_cube(lon, lat); - } else { - // Spherical grid needs wrap - inverse = wrap_point_sphere(lon, lat); - dest = pos_iProc_sphere(lon, lat); - } - - // Prepare to send the lon, lat to dest processor - lonlatmap[dest].push_back({lon, lat}); - // Record where to place data when receiving from that processor - exch_recv[dest].push_back({iLon, iLat, inverse}); - // Skip non ghost cells - if (iLat == nGCs - 1 && iLon >= nGCs && iLon < nLons - nGCs) { - iLat = nLats - nGCs - 1; - } - } + // Initialize communicator. Message exchange is limited in each + // ensemble and rank 0 in each ensemble has special purpose + MPI_Comm_split(aether_comm, iMember, iGrid, &grid_comm); + + // Classify the {lon, lat} for all ghost cells based on destination + std::unordered_map> lonlatmap; + + for (int64_t iLon = 0; iLon < nLons; ++iLon) { + for (int64_t iLat = 0; iLat < nLats; ++iLat) { + // Read the lon, lat of ghost cell + precision_t lon = geoLon_scgc(iLon, iLat, 0); + precision_t lat = geoLat_scgc(iLon, iLat, 0); + // Find which processor can handle it and whether + // the message crosses the north/south pole + int dest, inverse = 1; + + if (IsCubeSphereGrid) + dest = pos_iProc_cube(lon, lat); + + else { + // Spherical grid needs wrap + inverse = wrap_point_sphere(lon, lat); + dest = pos_iProc_sphere(lon, lat); + } + + // Prepare to send the lon, lat to dest processor + lonlatmap[dest].push_back({lon, lat}); + // Record where to place data when receiving from that processor + exch_recv[dest].push_back({iLon, iLat, inverse}); + + // Skip non ghost cells + if (iLat == nGCs - 1 && iLon >= nGCs && iLon < nLons - nGCs) + iLat = nLats - nGCs - 1; } - - // Construct send messages - std::vector> msg_in; - for (auto& it : lonlatmap) { - msg_in.emplace_back(std::move(it.second), it.first); + } + + // Construct send messages + std::vector> msg_in; + + for (auto& it : lonlatmap) + msg_in.emplace_back(std::move(it.second), it.first); + + // Create space to receive messages + std::vector> msg_out; + + // Send all longitude and latitude to where it can be handled + // Receive all longitude and latitude this processor needs to handle + CUSTOM_MPI_SENDRECV(msg_in, msg_out, grid_comm); + + // Build interpolation coefficients based on lon and lat + for (auto& it : msg_out) { + // Initialize the inputs to set_interp_coefs + std::vector Lons; + std::vector Lats; + std::vector Alts; + + // Combine each altitude with lat_lon into points + for (int64_t iAlt = nGCs; iAlt < nAlts - nGCs; ++iAlt) { + for (auto& it2 : it.buf) { + Lons.push_back(it2.lon); + Lats.push_back(it2.lat); + Alts.push_back(geoAlt_scgc(0, 0, iAlt)); + } } - // Create space to receive messages - std::vector> msg_out; - - // Send all longitude and latitude to where it can be handled - // Receive all longitude and latitude this processor needs to handle - CUSTOM_MPI_SENDRECV(msg_in, msg_out, grid_comm); - - // Build interpolation coefficients based on lon and lat - for (auto& it : msg_out) { - // Initialize the inputs to set_interp_coefs - std::vector Lons; - std::vector Lats; - std::vector Alts; - - // Combine each altitude with lat_lon into points - for (int64_t iAlt = nGCs; iAlt < nAlts - nGCs; ++iAlt) { - for (auto& it2 : it.buf) { - Lons.push_back(it2.lon); - Lats.push_back(it2.lat); - Alts.push_back(geoAlt_scgc(0, 0, iAlt)); - } - } - // Compute interpolation coefficients - set_interpolation_coefs(Lons, Lats, Alts); - - // Fix precision issues - for (auto& it : interp_coefs) { - if (it.rRow < 0.0001) { - it.rRow = 0; - } else if (it.rRow > 0.9999) { - it.rRow = 1; - } - if (it.rCol < 0.0001) { - it.rCol = 0; - } else if (it.rCol > 0.9999) { - it.rCol = 1; - } - } - // Store the coefficients - exch_send[it.rank] = std::move(interp_coefs); + + // Compute interpolation coefficients + set_interpolation_coefs(Lons, Lats, Alts); + + // Fix precision issues + for (auto& it : interp_coefs) { + if (it.rRow < 0.0001) + it.rRow = 0; + + else if (it.rRow > 0.9999) + it.rRow = 1; + + if (it.rCol < 0.0001) + it.rCol = 0; + + else if (it.rCol > 0.9999) + it.rCol = 1; } + + // Store the coefficients + exch_send[it.rank] = std::move(interp_coefs); + } } // -------------------------------------------------------------------------- @@ -531,41 +556,43 @@ void Grid::init_connection() { // -------------------------------------------------------------------------- void Grid::exchange(arma_cube &data, const bool pole_inverse) { - // Construct send and receive buffer - std::vector> msg_send; - std::vector> msg_recv; - - // Do the interpolation for other processors - for (auto& it : exch_send) { - std::swap(interp_coefs, it.second); - msg_send.emplace_back(get_interpolation_values(data), - it.first); - std::swap(interp_coefs, it.second); - } - - // Initialize the size of receive buffer - const int64_t numalt = nAlts - 2 * nGCs; - for (auto& it : exch_recv) { - msg_recv.emplace_back(std::vector(it.second.size()*numalt), - it.first); - } - - // Exchange message - CUSTOM_MPI_SENDRECV2(msg_send, msg_recv, grid_comm); - - // Fill the ghost cell with received message - for (auto& it : msg_recv) { - const std::vector &idx = exch_recv[it.rank]; - for (int64_t iAlt = nGCs; iAlt < nAlts - nGCs; ++iAlt) { - for (size_t i = 0; i < idx.size(); ++i) { - if (pole_inverse) { - it.buf[(iAlt - nGCs) * idx.size() + i] *= idx[i].inverse; - } - data(idx[i].ilon, idx[i].ilat, iAlt) = - it.buf[(iAlt - nGCs) * idx.size() + i]; - } - } + // Construct send and receive buffer + std::vector> msg_send; + std::vector> msg_recv; + + // Do the interpolation for other processors + for (auto& it : exch_send) { + std::swap(interp_coefs, it.second); + msg_send.emplace_back(get_interpolation_values(data), + it.first); + std::swap(interp_coefs, it.second); + } + + // Initialize the size of receive buffer + const int64_t numalt = nAlts - 2 * nGCs; + + for (auto& it : exch_recv) { + msg_recv.emplace_back(std::vector(it.second.size()*numalt), + it.first); + } + + // Exchange message + CUSTOM_MPI_SENDRECV2(msg_send, msg_recv, grid_comm); + + // Fill the ghost cell with received message + for (auto& it : msg_recv) { + const std::vector &idx = exch_recv[it.rank]; + + for (int64_t iAlt = nGCs; iAlt < nAlts - nGCs; ++iAlt) { + for (size_t i = 0; i < idx.size(); ++i) { + if (pole_inverse) + it.buf[(iAlt - nGCs) * idx.size() + i] *= idx[i].inverse; + + data(idx[i].ilon, idx[i].ilat, iAlt) = + it.buf[(iAlt - nGCs) * idx.size() + i]; + } } + } } // ----------------------------------------------------------------------------- @@ -573,17 +600,17 @@ void Grid::exchange(arma_cube &data, const bool pole_inverse) { // ----------------------------------------------------------------------------- bool Neutrals::exchange(Grid &grid) { - // For each species, exchange if its DoAdvect is true - for (int i = 0; i < nSpecies; ++i) { - if (species[i].DoAdvect) { - grid.exchange(species[i].density_scgc, false); - } - } - // Exchange temperature - grid.exchange(temperature_scgc, false); - // Exchange velocity - grid.exchange(velocity_vcgc[0], true); - grid.exchange(velocity_vcgc[1], true); - grid.exchange(velocity_vcgc[2], false); - return true; + // For each species, exchange if its DoAdvect is true + for (int i = 0; i < nSpecies; ++i) { + if (species[i].DoAdvect) + grid.exchange(species[i].density_scgc, false); + } + + // Exchange temperature + grid.exchange(temperature_scgc, false); + // Exchange velocity + grid.exchange(velocity_vcgc[0], true); + grid.exchange(velocity_vcgc[1], true); + grid.exchange(velocity_vcgc[2], false); + return true; } diff --git a/src/init_geo_grid.cpp b/src/init_geo_grid.cpp index c7e317a2..e50c1e7f 100644 --- a/src/init_geo_grid.cpp +++ b/src/init_geo_grid.cpp @@ -279,7 +279,8 @@ void transformation_metrics(Quadtree quadtree, A12_inv(i, j) = 0; A21_inv(i, j) = p2 * tan(latp) * tan(lonp - cPI / 4. - 3 * cPI / 2.); A22_inv(i, j) = p2 / cos(latp); - } else if (quadtree.iSide == 6 - 1) { // Face 5 and 6 are flipped than Nair's formulation + } else if (quadtree.iSide == 6 - + 1) { // Face 5 and 6 are flipped than Nair's formulation double p1 = R * sin(latp) / a; A11(i, j) = p1 * cos(lonp - 3 * cPI / 4.); A12(i, j) = p1 * sin(lonp - 3 * cPI / 4.); @@ -291,7 +292,8 @@ void transformation_metrics(Quadtree quadtree, A12_inv(i, j) = -p2 * sin(lonp - 3 * cPI / 4.); A21_inv(i, j) = p2 * sin(latp) * sin(lonp - 3 * cPI / 4.); A22_inv(i, j) = p2 * cos(lonp - 3 * cPI / 4.); - } else if (quadtree.iSide == 5 - 1) { // Face 5 and 6 are flipped than Nair's formulation + } else if (quadtree.iSide == 5 - + 1) { // Face 5 and 6 are flipped than Nair's formulation double p1 = R * sin(latp) / a; A11(i, j) = -p1 * cos(lonp - 3 * cPI / 4.); A12(i, j) = p1 * sin(lonp - 3 * cPI / 4.); @@ -862,7 +864,7 @@ void Grid::correct_xy_grid(Planets planet) { int64_t iAlt; // initialize grid drefx drefy - drefx = arma_vec(nAlts); + drefx = arma_vec(nAlts); drefy = arma_vec(nAlts); // Planet.get_radius() takes in latitude diff --git a/src/inputs.cpp b/src/inputs.cpp index 7385764d..ac912957 100644 --- a/src/inputs.cpp +++ b/src/inputs.cpp @@ -91,8 +91,10 @@ bool Inputs::write_restart() { std::string Inputs::get_logfile() { std::string logfile = check_settings_str("Logfile", "name"); + if (nMembers > 1) logfile = add_cmember(logfile); + return logfile; } @@ -230,7 +232,7 @@ std::string Inputs::get_settings_str(std::string key1, } int Inputs::get_settings(std::string key1, - std::string key2) { + std::string key2) { int value = -1; if (settings.find(key1) != settings.end()) @@ -255,8 +257,8 @@ bool Inputs::check_settings(std::string key1, if (report.test_verbose(1)) std::cout << "checking setting : " - << key1 << " and " - << key2 << "\n"; + << key1 << " and " + << key2 << "\n"; //try to find the keys first if (settings.find(key1) != settings.end()) { @@ -516,6 +518,7 @@ int Inputs::get_original_seed() { std::cout << "Error in getting seed!\n"; return 0; } + return settings.at("Seed"); } @@ -582,20 +585,24 @@ bool Inputs::set_verbose(json in) { bool DidWork = true; int iVerbose = -1; + // Want to set verbose level ASAP: if (in.contains("Debug")) { if (in.at("Debug").contains("iVerbose")) { iVerbose = in.at("Debug").at("iVerbose"); + if (in.at("Debug").contains("iProc")) { - if (iProc != in.at("Debug").at("iProc")) - iVerbose = -1; + if (iProc != in.at("Debug").at("iProc")) + iVerbose = -1; } - } + } } + if (iVerbose > 0) { std::cout << "Setting iVerbose : " << iVerbose << "\n"; report.set_verbose(iVerbose); } + return DidWork; } diff --git a/src/logfile.cpp b/src/logfile.cpp index 1d43d326..7d5db759 100644 --- a/src/logfile.cpp +++ b/src/logfile.cpp @@ -253,7 +253,7 @@ Logfile::Logfile(Indices &indices) { std::vector sat_dts = input.get_satellite_dts(); bool isRestart = input.get_do_restart(); - + // Only the 0th processor in each member needs to open the logfile if (iGrid == 0) { @@ -264,6 +264,7 @@ Logfile::Logfile(Indices &indices) { logfilestream.open(logfileName, std::ofstream::trunc); logfilestream.precision(4); } + // Report error if can not open the log file stream if (!logfilestream.is_open() & !isRestart) { // TRY TO EXIT GRACEFULLY HERE. ALL THE FOLLOWING CODE SHOULD @@ -297,6 +298,7 @@ Logfile::Logfile(Indices &indices) { // Write the header to log if (!isRestart) logfilestream << header_time << header_log << '\n'; + // Close the file stream if append if (doAppend) logfilestream.close(); diff --git a/src/neutrals.cpp b/src/neutrals.cpp index 0b5c73c7..519a9711 100644 --- a/src/neutrals.cpp +++ b/src/neutrals.cpp @@ -140,6 +140,7 @@ Neutrals::Neutrals(Grid grid, // This gets a bunch of the species-dependent characteristics: iErr = read_planet_file(planet); + if (iErr > 0) report.error("Error reading planet file!"); diff --git a/src/neutrals_bcs.cpp b/src/neutrals_bcs.cpp index b8d497e0..0c7c4017 100644 --- a/src/neutrals_bcs.cpp +++ b/src/neutrals_bcs.cpp @@ -32,8 +32,10 @@ bool Neutrals::set_bcs(Grid grid, if (input.get_nAltsGeo() > 1) { didWork = set_lower_bcs(grid, time, indices); + if (didWork) didWork = set_upper_bcs(grid); + if (didWork) calc_mass_density(); } @@ -132,6 +134,7 @@ bool Neutrals::set_lower_bcs(Grid grid, if (!msis.is_ok()) { didWork = false; report.error("MSIS initialization not ok"); + if (report.test_verbose(0)) { std::cout << "MSIS Boundary Conditions asked for, "; std::cout << "but MSIS is not compiled! Yikes!\n"; @@ -191,6 +194,7 @@ bool Neutrals::set_lower_bcs(Grid grid, species[iSpecies].density_scgc.slice(0). fill(species[iSpecies].lower_bc_density); } + temperature_scgc.slice(0).fill(initial_temperatures[0]); didWork = true; } @@ -220,7 +224,7 @@ bool Neutrals::set_lower_bcs(Grid grid, report.error("issue with lower BCs!"); report.error("maybe check boundaryconditions type : " + bcsType); } - + report.exit(function); return didWork; } diff --git a/src/neutrals_ics.cpp b/src/neutrals_ics.cpp index 558a5b07..9f67e005 100644 --- a/src/neutrals_ics.cpp +++ b/src/neutrals_ics.cpp @@ -17,8 +17,8 @@ // ----------------------------------------------------------------------------- bool Neutrals::initial_conditions(Grid grid, - Times time, - Indices indices) { + Times time, + Indices indices) { std::string function = "Neutrals::initial_conditions"; static int iFunction = -1; @@ -27,7 +27,7 @@ bool Neutrals::initial_conditions(Grid grid, // initialize didWork to false, so that we can catch if the initial // conditions are not actually set bool didWork = false; - + int64_t iLon, iLat, iAlt, iA; precision_t alt, r; int64_t nAlts = grid.get_nZ(); @@ -37,8 +37,10 @@ bool Neutrals::initial_conditions(Grid grid, if (input.get_do_restart()) { report.print(1, "Restarting! Reading neutral files!"); didWork = restart_file(input.get_restartin_dir(), DoRead); + if (!didWork) { - report.error("Reading Restart for Neutrals Failed!!!");\ + report.error("Reading Restart for Neutrals Failed!!!"); + \ } } else { @@ -53,48 +55,48 @@ bool Neutrals::initial_conditions(Grid grid, if (!msis.is_ok()) { didWork = false; - report.error("MSIS Initial Conditions asked for, "); + report.error("MSIS Initial Conditions asked for, "); report.error("but MSIS is not compiled/working!"); } else { - didWork = true; - msis.set_time(time); - precision_t f107 = indices.get_f107(time.get_current()); - precision_t f107a = indices.get_f107a(time.get_current()); - msis.set_f107(f107, f107a); - msis.set_ap(10.0); - msis.set_locations(grid.geoLon_scgc, - grid.geoLat_scgc, - grid.geoAlt_scgc); - - // This is just to check if MSIS is actually working: - if (msis.is_valid_species("Tn")) - // if it is, fill will temperature: - temperature_scgc = msis.get_cube("Tn"); - else - // if it is not, then fill with a value: - temperature_scgc.fill(initial_temperatures[0]); - - for (int iSpecies = 0; iSpecies < nSpecies; iSpecies++) { - if (report.test_verbose(3)) - std::cout << "Setting Species : " - << species[iSpecies].cName << "\n"; - - if (msis.is_valid_species(species[iSpecies].cName)) { - if (report.test_verbose(3)) - std::cout << " Found in MSIS!\n"; - - species[iSpecies].density_scgc = - msis.get_cube(species[iSpecies].cName); - } else { - if (report.test_verbose(3)) - std::cout << " NOT Found in MSIS - setting constant\n"; - - species[iSpecies].density_scgc.slice(0). - fill(species[iSpecies].lower_bc_density); - fill_with_hydrostatic(iSpecies, 1, nAlts, grid); - } - - } // for species + didWork = true; + msis.set_time(time); + precision_t f107 = indices.get_f107(time.get_current()); + precision_t f107a = indices.get_f107a(time.get_current()); + msis.set_f107(f107, f107a); + msis.set_ap(10.0); + msis.set_locations(grid.geoLon_scgc, + grid.geoLat_scgc, + grid.geoAlt_scgc); + + // This is just to check if MSIS is actually working: + if (msis.is_valid_species("Tn")) + // if it is, fill will temperature: + temperature_scgc = msis.get_cube("Tn"); + else + // if it is not, then fill with a value: + temperature_scgc.fill(initial_temperatures[0]); + + for (int iSpecies = 0; iSpecies < nSpecies; iSpecies++) { + if (report.test_verbose(3)) + std::cout << "Setting Species : " + << species[iSpecies].cName << "\n"; + + if (msis.is_valid_species(species[iSpecies].cName)) { + if (report.test_verbose(3)) + std::cout << " Found in MSIS!\n"; + + species[iSpecies].density_scgc = + msis.get_cube(species[iSpecies].cName); + } else { + if (report.test_verbose(3)) + std::cout << " NOT Found in MSIS - setting constant\n"; + + species[iSpecies].density_scgc.slice(0). + fill(species[iSpecies].lower_bc_density); + fill_with_hydrostatic(iSpecies, 1, nAlts, grid); + } + + } // for species } // msis init worked ok } // type = msis @@ -171,7 +173,7 @@ bool Neutrals::initial_conditions(Grid grid, if (!didWork) report.error("Issue with initial conditions!"); - + report.exit(function); return didWork; } diff --git a/src/neutrals_momentum_friction.cpp b/src/neutrals_momentum_friction.cpp index 78a738e7..30493882 100644 --- a/src/neutrals_momentum_friction.cpp +++ b/src/neutrals_momentum_friction.cpp @@ -94,7 +94,7 @@ void calc_neutral_friction(Neutrals &neutrals) { for (iSpecies = 0; iSpecies < neutrals.nSpecies; iSpecies++) for (iDir = 0; iDir < 3; iDir++) neutrals.species[iSpecies].acc_neutral_friction[iDir].zeros(); - + if (input.get_advection_neutrals_vertical() != "hydro") { arma_vec vels(neutrals.nSpeciesAdvect, fill::zeros); diff --git a/src/output.cpp b/src/output.cpp index f35b7ee5..01ea5f86 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -20,10 +20,10 @@ // ----------------------------------------------------------------------------- bool output(const Neutrals &neutrals, - const Ions &ions, - const Grid &grid, - Times time, - const Planets &planet) { + const Ions &ions, + const Grid &grid, + Times time, + const Planets &planet) { std::string function = "output"; static int iFunction = -1; @@ -311,8 +311,8 @@ bool output(const Neutrals &neutrals, // write output container if (!AllOutputContainers[iOutput].write()) { - report.error("Error in writing output container!"); - didWork = false; + report.error("Error in writing output container!"); + didWork = false; } // ------------------------------------------------------------ diff --git a/src/read_input_file.cpp b/src/read_input_file.cpp index 6acf4990..fbc3786d 100644 --- a/src/read_input_file.cpp +++ b/src/read_input_file.cpp @@ -29,7 +29,7 @@ bool Inputs::read_inputs_json(Times &time) { // Set the default values first: settings = read_json("UA/inputs/defaults.json"); DidWork = set_verbose(settings); - + // Set the planet-specific file (user can change this in aether.in file!): settings["PlanetSpeciesFile"] = settings["Planet"]["file"]; @@ -50,9 +50,9 @@ bool Inputs::read_inputs_json(Times &time) { restart_file = restart_file + "/settings.json"; json restart_inputs; restart_inputs = read_json(restart_file); - // This forces the logfile to append. User can override - // if they really want: - restart_inputs["Logfile"]["append"] = true; + // This forces the logfile to append. User can override + // if they really want: + restart_inputs["Logfile"]["append"] = true; settings.merge_patch(restart_inputs); } } diff --git a/src/solver_gradients.cpp b/src/solver_gradients.cpp index 0b587a98..f03f4228 100644 --- a/src/solver_gradients.cpp +++ b/src/solver_gradients.cpp @@ -186,35 +186,38 @@ std::vector calc_gradient_cubesphere(arma_cube value, Grid grid) { arma_mat grad_x_curr(nXs, nYs); arma_mat grad_y_curr(nXs, nYs); - + // Calc gradient on x and y direction (since reference grid) // Only update interior cells // May vectorize for future improvements - if (nGCs >=2) { // if more than 1 nGCs, we do fourth order, some foolproofing in case we go into debug hell - for (int j = nGCs; j < nYs - nGCs; j++) - { - for (int i = nGCs; i < nXs - nGCs; i++) - { - grad_x_curr(i, j) = (-curr_value(i+2,j) + 8*curr_value(i+1,j) - 8*curr_value(i-1,j) + curr_value(i-2,j))*(1./12./dx); - grad_y_curr(i, j) = (-curr_value(i,j+2) + 8*curr_value(i,j+1) - 8*curr_value(i,j-1) + curr_value(i,j-2))*(1./12./dy); + if (nGCs >= + 2) { // if more than 1 nGCs, we do fourth order, some foolproofing in case we go into debug hell + for (int j = nGCs; j < nYs - nGCs; j++) { + for (int i = nGCs; i < nXs - nGCs; i++) { + grad_x_curr(i, j) = (-curr_value(i + 2, j) + 8 * curr_value(i + 1, + j) - 8 * curr_value(i - 1, j) + curr_value(i - 2, j)) * (1. / 12. / dx); + grad_y_curr(i, j) = (-curr_value(i, j + 2) + 8 * curr_value(i, + j + 1) - 8 * curr_value(i, j - 1) + curr_value(i, j - 2)) * (1. / 12. / dy); } } } else { // otherwise we do second order - for (int j = nGCs; j < nYs - nGCs; j++) - { - for (int i = nGCs; i < nXs - nGCs; i++) - { - grad_x_curr(i, j) = (curr_value(i+1,j)-curr_value(i-1,j))*(1./2./dx); - grad_y_curr(i, j) = (curr_value(i,j+1)-curr_value(i,j-1))*(1./2./dy); + for (int j = nGCs; j < nYs - nGCs; j++) { + for (int i = nGCs; i < nXs - nGCs; i++) { + grad_x_curr(i, j) = (curr_value(i + 1, j) - curr_value(i - 1, + j)) * (1. / 2. / dx); + grad_y_curr(i, j) = (curr_value(i, j + 1) - curr_value(i, + j - 1)) * (1. / 2. / dy); } } } - + // We then use A transformation matrices to convert grad_xy to grad_latlon // Ref -> Physical, we use A matrix - grad_lon.slice(iAlt) = grad_x_curr%grid.A11_inv_scgc.slice(iAlt)+grad_y_curr%grid.A21_inv_scgc.slice(iAlt); - grad_lat.slice(iAlt) = grad_x_curr%grid.A12_inv_scgc.slice(iAlt)+grad_y_curr%grid.A22_inv_scgc.slice(iAlt); + grad_lon.slice(iAlt) = grad_x_curr % grid.A11_inv_scgc.slice( + iAlt) + grad_y_curr % grid.A21_inv_scgc.slice(iAlt); + grad_lat.slice(iAlt) = grad_x_curr % grid.A12_inv_scgc.slice( + iAlt) + grad_y_curr % grid.A22_inv_scgc.slice(iAlt); } // Not initializing with array like procedure in case I get bugs diff --git a/src/tools.cpp b/src/tools.cpp index f8643e17..92b1fb38 100644 --- a/src/tools.cpp +++ b/src/tools.cpp @@ -119,8 +119,10 @@ bool compare(precision_t value1, precision_t value2) { std::string add_cmember(std::string inString) { std::string outString = inString; std::size_t found = outString.rfind("."); - if (found!=std::string::npos) + + if (found != std::string::npos) outString.replace(found, 1, "_" + cMember + "."); + return outString; } diff --git a/src/transform.cpp b/src/transform.cpp index 6c5973dd..33a23bbe 100644 --- a/src/transform.cpp +++ b/src/transform.cpp @@ -13,8 +13,10 @@ std::string mklower(std::string inString) { std::string outString = inString; int64_t nChars = outString.length(); + for (int64_t iChar = 0; iChar < nChars; iChar++) outString[iChar] = tolower(outString[iChar]); + return outString; }