diff --git a/edu/examples/Advection/cubesphere_equal_angle.cpp b/edu/examples/Advection/cubesphere_equal_angle.cpp new file mode 100644 index 00000000..3c3687e4 --- /dev/null +++ b/edu/examples/Advection/cubesphere_equal_angle.cpp @@ -0,0 +1,1757 @@ +/* + This is an example of a second order 2D solver for the Euler equations. + + to compile: + g++ -I/usr/local/include -I/Users/ridley/Software/Json/json/include -o cubesphere2d cubesphere2d.cpp + +*/ + +#include +#include + +using precision_t = double; + +/// Armadillo type vector (single column) with compile-time precision. +using arma_vec = arma::Col; +/// Armadillo type matrix (two dimension) with compile-time precision. +using arma_mat = arma::Mat; +/// Armadillo type cube (three dimension) with compile-time precision. +using arma_cube = arma::Cube; + +precision_t cPI = 3.141592653589793; +precision_t cTWOPI = 2.0 * cPI; +precision_t cRtoD = 180.0 / cPI; +precision_t cPIdiv2 = cPI / 2; +precision_t cGamma = 5.0 / 3.0; // Specific ratio of heat +precision_t cKb = 1.38e-23; +precision_t mmm = 16.0 * 1.67e-27; + +// --------------------------------------------------------- +// A couple of global variables +// --------------------------------------------------------- + +int64_t verbose = 1; + +struct projection_struct { + arma_mat gradLR; + arma_mat gradDU; + arma_mat R; + arma_mat L; + arma_mat U; + arma_mat D; +}; + +struct grid_struct { + + // sizes: + int64_t nXt, nYt, nGCs; + int64_t iXfirst_, iXlast_; + int64_t iYfirst_, iYlast_; + + // Positions: + arma_mat lon; + arma_mat lat; + + // These are for Ronchi et al., JCP 124, 93-114, 1996 + arma_mat X, Y, Z, C, D, d; + arma_mat dlx, dln, dS; + // xi is the LR direction + // nu is the UD direction + arma_mat xi, nu; + arma_mat x, y, r; + arma_mat Apn, Apx, Atn, Atx; + arma_mat Axt, Axp, Ant, Anp; + precision_t dxi, dnu, R; + arma_mat alpha; + arma_mat sinAlpha; + + arma_mat nXiLon; + arma_mat nXiLat; + arma_mat nNuLon; + arma_mat nNuLat; + // These are eq28 of Nair (g lower ij): + arma_mat gl11, gl12, gl21, gl22; + // These are eq29 of Nair (g upper ij): + arma_mat sqrtg; + arma_mat gu11, gu12, gu21, gu22; + // These are eq32 of Nair (sphere-to-cube): + arma_mat s2c11, s2c12, s2c21, s2c22; + arma_mat c2s11, c2s12, c2s21, c2s22; +}; + +// --------------------------------------------------------- +// +// --------------------------------------------------------- + +precision_t calc_dt(arma_mat dx, + arma_mat dy, + arma_mat &wsLR, + arma_mat &wsDU, + int64_t nGCs) { + + if (verbose > 2) + std::cout << " --> calc_dt\n"; + + int64_t nX = wsLR.n_rows; + int64_t nY = wsLR.n_cols; + + precision_t wsX, wsY, dtX, dtY, dt; + + dt = 1e32; + + for (int64_t j = nGCs; j < nY - nGCs; j++) { + for (int64_t i = nGCs; i < nX - nGCs; i++) { + wsX = (wsLR(i + 1, j) + wsLR(i, j)) / 2; + dtX = dx(i, j) / wsX; + wsY = (wsDU(i, j + 1) + wsDU(i, j)) / 2; + dtY = dy(i, j) / wsY; + + if (dtX < dt) + dt = dtX; + + if (dtY < dt) + dt = dtY; + } + } + + return dt; +} + +/** + * Output function + * + * @param values Values + * @param filename FileName + * @param DoAppend + */ +void output(arma_mat &values, + std::string filename, + bool DoAppend) { + + std::ofstream outfile; + + if (DoAppend) + outfile.open(filename, std::ios_base::app); + else { + outfile.open(filename); + int64_t nX = values.n_rows; + int64_t nY = values.n_cols; + outfile << nX << " " << nY << "\n"; + } + + outfile << values; + outfile.close(); +} + +/** + * Transform spherical coordinates to 3D Cartesian + * + * doi: 10.1016/j.jcp.2007.07.022 + * Section 3, Eqn (23) + * + * @return dh Great Circle Distance between two points + */ +arma_vec sph2cart(precision_t lon, + precision_t lat, + precision_t r) { + arma_vec xyz(3); + xyz(0) = r * std::cos(lat) * std::cos(lon); + xyz(1) = r * std::cos(lat) * std::sin(lon); + xyz(2) = r * std::sin(lat); + return xyz; +} + +grid_struct init_grid_equidistant(int iFace, + int64_t nX, + int64_t nY, + int64_t nGCs, + precision_t R, + precision_t xOff, + precision_t yOff) { + + double a = R / std::sqrt(3); + + grid_struct grid; + int64_t nXt = nX + 2 * nGCs; + int64_t nYt = nY + 2 * nGCs; + + grid.nXt = nXt; + grid.nYt = nYt; + grid.nGCs = nGCs; + grid.iXfirst_ = nGCs; + grid.iYfirst_ = nGCs; + grid.iXlast_ = nX + nGCs; + grid.iYlast_ = nY + nGCs; + + // Positions: + grid.lon.resize(nXt, nYt); + grid.lat.resize(nXt, nYt); + + grid.x.resize(nXt, nYt); + grid.y.resize(nXt, nYt); + grid.r.resize(nXt, nYt); + grid.xi.resize(nXt, nYt); + grid.nu.resize(nXt, nYt); + grid.X.resize(nXt, nYt); + grid.Y.resize(nXt, nYt); + grid.Z.resize(nXt, nYt); + grid.C.resize(nXt, nYt); + grid.D.resize(nXt, nYt); + grid.d.resize(nXt, nYt); + grid.dlx.resize(nXt, nYt); + grid.dln.resize(nXt, nYt); + grid.dS.resize(nXt, nYt); + + grid.s2c11.resize(nXt, nYt); + grid.s2c12.resize(nXt, nYt); + grid.s2c21.resize(nXt, nYt); + grid.s2c22.resize(nXt, nYt); + grid.c2s11.resize(nXt, nYt); + grid.c2s12.resize(nXt, nYt); + grid.c2s21.resize(nXt, nYt); + grid.c2s22.resize(nXt, nYt); + + grid.gl11.resize(nXt, nYt); + grid.gl12.resize(nXt, nYt); + grid.gl21.resize(nXt, nYt); + grid.gl22.resize(nXt, nYt); + grid.gu11.resize(nXt, nYt); + grid.gu12.resize(nXt, nYt); + grid.gu21.resize(nXt, nYt); + grid.gu22.resize(nXt, nYt); + grid.sqrtg.resize(nXt, nYt); + + double iD, iL, x, y, theta, phi, r; + + double dx = 2 * a / nX; + double dy = 2 * a / nY; + double r3, r4; + double R2 = R * R; + + // Loop through each point and derive the coordinate + // DU is y-direction (down-up) + // LR is x-direction (left-right) + for (int iDU = 0; iDU < nYt; iDU++) { + for (int iLR = 0; iLR < nXt; iLR++) { + + // the offsets are so we can find cell centers, edges, and corners + // Centers assume Off = 0.5, which edges assume Off = 0 + double iD = iDU - nGCs + yOff; + double iL = iLR - nGCs + xOff; + + x = -a + iL * dx; + y = -a + iD * dy; + phi = std::atan(x / a); + // y = a * tan(theta) * sec(phi) => y * cos(phi) = a * tan(theta) + theta = std::atan(y * std::cos(phi) / a); + //std::cout << "Grid creation : " + // << iDU << " " + // << iLR << " " + // << y << " " + // << x << " " + // << theta * cRtoD << " " + // << phi * cRtoD << "\n"; + grid.lon(iLR, iDU) = phi; + grid.lat(iLR, iDU) = theta; + grid.X(iLR, iDU) = R * std::cos(theta) * std::cos(phi); + grid.Y(iLR, iDU) = R * std::cos(theta) * std::sin(phi); + grid.Z(iLR, iDU) = R * std::sin(theta); + grid.r(iLR, iDU) = std::sqrt(a * a + x * x + y * y); + + // Equation 28 of Nair: + r4 = r * r * r * r; + grid.gl11(iLR, iDU) = R2 / r4 * (a * a + y * y); + grid.gl12(iLR, iDU) = - R2 / r4 * (x * y); + grid.gl21(iLR, iDU) = - R2 / r4 * (x * y); + grid.gl22(iLR, iDU) = R2 / r4 * (a * a + x * x); + // Equation 29 of Nair: + r3 = r * r * r; + grid.sqrtg(iLR, iDU) = R2 * a / r3; + grid.gu11(iLR, iDU) = grid.gl22(iLR, iDU) / grid.sqrtg(iLR, iDU); + grid.gu12(iLR, iDU) = - grid.gl12(iLR, iDU) / grid.sqrtg(iLR, iDU); + grid.gu21(iLR, iDU) = - grid.gl21(iLR, iDU) / grid.sqrtg(iLR, iDU); + grid.gu22(iLR, iDU) = grid.gl11(iLR, iDU) / grid.sqrtg(iLR, iDU); + + grid.s2c11(iLR, iDU) = a / (R * std::cos(theta) * std::cos(phi)) * + (1 / std::cos(theta)); + grid.s2c12(iLR, iDU) = 0.0; + grid.s2c21(iLR, iDU) = a / (R * std::cos(theta) * std::cos(phi)) * + (std::tan(theta) * std::tan(phi)); + grid.s2c22(iLR, iDU) = a / (R * std::cos(theta) * std::cos(phi)) * + (1 / std::cos(phi)); + grid.c2s11(iLR, iDU) = (R * std::cos(theta) * std::cos(phi)) / a * std::cos( + theta); + grid.c2s12(iLR, iDU) = 0.0; + grid.c2s21(iLR, iDU) = -(R * std::cos(theta) * std::cos(phi)) / a * + std::sin(theta) * std::sin(phi); + grid.c2s22(iLR, iDU) = (R * std::cos(theta) * std::cos(phi)) / a * std::cos( + phi); + + } + } + + //xxx + return grid; +} + + +grid_struct init_grid(int iFace, + int64_t nX, int64_t nY, int64_t nGCs, + precision_t R, + precision_t xOff, precision_t yOff) { + + grid_struct grid; + int64_t nXt = nX + 2 * nGCs; + int64_t nYt = nY + 2 * nGCs; + + grid.nXt = nXt; + grid.nYt = nYt; + grid.nGCs = nGCs; + grid.iXfirst_ = nGCs; + grid.iYfirst_ = nGCs; + grid.iXlast_ = nX + nGCs; + grid.iYlast_ = nY + nGCs; + + // Positions: + grid.lon.resize(nXt, nYt); + grid.lat.resize(nXt, nYt); + + grid.xi.resize(nXt, nYt); + grid.nu.resize(nXt, nYt); + grid.X.resize(nXt, nYt); + grid.Y.resize(nXt, nYt); + grid.C.resize(nXt, nYt); + grid.D.resize(nXt, nYt); + grid.d.resize(nXt, nYt); + grid.dlx.resize(nXt, nYt); + grid.dln.resize(nXt, nYt); + grid.dS.resize(nXt, nYt); + + grid.Axt.resize(nXt, nYt); + grid.Axp.resize(nXt, nYt); + grid.Ant.resize(nXt, nYt); + grid.Anp.resize(nXt, nYt); + + grid.Apn.resize(nXt, nYt); + grid.Apx.resize(nXt, nYt); + grid.Atn.resize(nXt, nYt); + grid.Atx.resize(nXt, nYt); + + precision_t fortyfive = cPI / 4.0; + // Xi is LR (x), Nu is UD (y) + precision_t dxi = 2.0 * fortyfive / (nX - 1 + xOff * 2); + precision_t dnu = 2.0 * fortyfive / (nY - 1 + yOff * 2); + + grid.dxi = dxi; + grid.dnu = dnu; + grid.R = R; + + precision_t latp, lonp; + + // Loop through each point and derive the coordinate + + precision_t total_area = 0.0, det, dmo; + + for (int iDU = 0; iDU < nYt; iDU++) { + for (int iLR = 0; iLR < nXt; iLR++) { + + // the offsets are so we can find cell centers, edges, and corners + double iD = iDU - nGCs + yOff; + double iL = iLR - nGCs + xOff; + + // Define local coordinates: + // Xi is LR (x), Nu is UD (y) + grid.nu(iLR, iDU) = (-fortyfive + dnu * iD); + grid.xi(iLR, iDU) = (-fortyfive + dxi * iL); + + grid.X(iLR, iDU) = tan(grid.xi(iLR, iDU)); + grid.Y(iLR, iDU) = tan(grid.nu(iLR, iDU)); + + // Transformation from 3D Cartesian to LatLong + // lonp = std::atan2(y_cart, x_cart) + cPI/2.0; + if (iFace == 0) { + lonp = std::atan(grid.X(iLR, iDU)); + // Theta in Ronchi is from the north pole, so lat is 90 - theta + latp = std::atan(1.0 / grid.Y(iLR, iDU) / std::cos(lonp)); + } + + if (iFace == 1) { + lonp = std::atan(-1.0 / grid.X(iLR, iDU)); + + if (lonp < 0) + lonp = cPI + lonp; + + // Theta in Ronchi is from the north pole, so lat is 90 - theta + latp = std::atan(1.0 / grid.Y(iLR, iDU) / std::sin(lonp)); + } + + if (iFace == 2) { + lonp = std::atan(grid.X(iLR, iDU)) + cPI; + // Theta in Ronchi is from the north pole, so lat is 90 - theta + latp = std::atan(-1.0 / grid.Y(iLR, iDU) / std::cos(lonp)); + } + + if (iFace == 3) { + lonp = std::atan(-1.0 / grid.X(iLR, iDU)); + + if (lonp > 0) + lonp = lonp + cPI; + else + lonp = 2 * cPI + lonp; + + // Theta in Ronchi is from the north pole, so lat is 90 - theta + latp = std::atan(-1.0 / grid.Y(iLR, iDU) / std::sin(lonp)); + } + + if (iFace == 4) { + lonp = std::atan2(grid.X(iLR, iDU), grid.Y(iLR, iDU)); + latp = std::atan2(-grid.Y(iLR, iDU), cos(lonp) ); + } + + if (iFace == 5) { + lonp = std::atan2(-grid.X(iLR, iDU), grid.Y(iLR, iDU)); + latp = -std::atan2(-grid.Y(iLR, iDU), cos(lonp) ); + } + + if (latp > 0) + latp = cPI / 2 - latp; + else + latp = -(cPI / 2 + latp); + + // Fill Computed coords + grid.lat(iLR, iDU) = latp; + grid.lon(iLR, iDU) = lonp; + + grid.d(iLR, iDU) = + 1 + + grid.X(iLR, iDU) * grid.X(iLR, iDU) + + grid.Y(iLR, iDU) * grid.Y(iLR, iDU); + grid.C(iLR, iDU) = sqrt(1 + + grid.X(iLR, iDU) * grid.X(iLR, iDU)); + grid.D(iLR, iDU) = sqrt(1 + + grid.Y(iLR, iDU) * grid.Y(iLR, iDU)); + + if (iFace < 4) { + grid.Axt(iLR, iDU) = 0.0; + grid.Axp(iLR, iDU) = grid.C(iLR, iDU) * grid.D(iLR, iDU) / + sqrt(grid.d(iLR, iDU)); + grid.Ant(iLR, iDU) = -1.0; + grid.Anp(iLR, iDU) = grid.X(iLR, iDU) * grid.Y(iLR, iDU) / + sqrt(grid.d(iLR, iDU)); + } else { + dmo = 1.0 / std::sqrt(grid.d(iLR, iDU) - 1); + + //if (dmo > 100.0) + // dmo = 100.0; + + if (iFace == 4) { + grid.Axt(iLR, iDU) = - dmo * grid.D(iLR, iDU) * grid.X(iLR, iDU); + grid.Axp(iLR, iDU) = dmo * grid.D(iLR, iDU) * grid.Y(iLR, iDU) / + sqrt(grid.d(iLR, iDU)); + grid.Ant(iLR, iDU) = - dmo * grid.C(iLR, iDU) * grid.Y(iLR, iDU); + grid.Anp(iLR, iDU) = - dmo * grid.C(iLR, iDU) * grid.X(iLR, iDU) / + sqrt(grid.d(iLR, iDU)); + + } else { + // iFace == 5 + grid.Axt(iLR, iDU) = dmo * grid.D(iLR, iDU) * grid.X(iLR, iDU); + grid.Axp(iLR, iDU) = - dmo * grid.D(iLR, iDU) * grid.Y(iLR, iDU) / + sqrt(grid.d(iLR, iDU)); + grid.Ant(iLR, iDU) = dmo * grid.C(iLR, iDU) * grid.Y(iLR, iDU); + grid.Anp(iLR, iDU) = dmo * grid.C(iLR, iDU) * grid.X(iLR, iDU) / + sqrt(grid.d(iLR, iDU)); + } + } + + // Calculate inverse of matrix for calculating Ax and An from At and Ap: + det = 1.0 / (grid.Axt(iLR, iDU) * grid.Anp(iLR, iDU) - + grid.Axp(iLR, iDU) * grid.Ant(iLR, iDU)); + + grid.Atx(iLR, iDU) = det * grid.Anp(iLR, iDU); + grid.Atn(iLR, iDU) = - det * grid.Axp(iLR, iDU); + grid.Apx(iLR, iDU) = - det * grid.Ant(iLR, iDU); + //grid.Apx(iLR, iDU) = grid.d(iLR, iDU) / (grid.C(iLR, iDU) * grid.D(iLR, iDU)); + grid.Apn(iLR, iDU) = det * grid.Axt(iLR, iDU); + + grid.dlx(iLR, iDU) = + R * grid.D(iLR, iDU) * dxi / + grid.d(iLR, iDU) / + (cos(grid.xi(iLR, iDU)) * cos(grid.xi(iLR, iDU))); + grid.dln(iLR, iDU) = + R * grid.C(iLR, iDU) * dnu / + grid.d(iLR, iDU) / + (cos(grid.nu(iLR, iDU)) * cos(grid.nu(iLR, iDU))); + + //grid.dS(iLR, iDU) = R * R * dxi * dnu / + // (sqrt(grid.d(iLR, iDU) * grid.d(iLR, iDU) * grid.d(iLR, iDU))) * + // grid.C(iLR, iDU) * grid.C(iLR, iDU) * + // grid.D(iLR, iDU) * grid.D(iLR, iDU); + grid.dS(iLR, iDU) = R * R * dxi * dnu / + (sqrt(grid.d(iLR, iDU) * grid.d(iLR, iDU) * grid.d(iLR, iDU)) * + cos(grid.xi(iLR, iDU)) * cos(grid.xi(iLR, iDU)) * + cos(grid.nu(iLR, iDU)) * cos(grid.nu(iLR, iDU))); + //grid.dS(iLR, iDU) = R * R * dxi * dnu * + // grid.C(iLR, iDU) * grid.D(iLR, iDU) / + // (grid.d(iLR, iDU) * grid.d(iLR, iDU) * + // cos(grid.xi(iLR, iDU)) * cos(grid.xi(iLR, iDU)) * + // cos(grid.nu(iLR, iDU)) * cos(grid.nu(iLR, iDU))); + + if (iLR > 2 && iLR < nXt - 3 && + iDU > 2 && iDU < nYt - 3) + total_area = total_area + grid.dS(iLR, iDU); + } + } + + std::cout << "Total Area : " << total_area << "; expected : " << 4 * cPI * R * + R / 6.0 << "\n"; + + return grid; +} + +/** + * Calculate Great Circle Distance + * + * doi: 10.1016/j.jcp.2007.07.022 + * Section 3, Eqn (23) + * + * @return dh Great Circle Distance between two points + */ +precision_t calc_great_circle(precision_t lon1, + precision_t lon2, + precision_t lat1, + precision_t lat2) { + + precision_t dlon_2 = (lon2 - lon1) / 2.0; + precision_t dlat_2 = (lat2 - lat1) / 2.0; + + precision_t dh = 2.0 * std::asin(std::sqrt(std::sin(dlat_2) * std::sin(dlat_2) + + std::sin(dlon_2) * std::sin(dlon_2) * std::cos(lat1) * std::cos(lat2))); + + return dh; +} + +// --------------------------------------------------------- +// Angle between three points on a sphere +// --------------------------------------------------------- +precision_t calc_angle_given_three_lon_lat(precision_t p1_lon, + precision_t p1_lat, + precision_t p2_lon, + precision_t p2_lat, + precision_t p3_lon, + precision_t p3_lat, + precision_t R) { + + arma_vec p1 = sph2cart(p1_lon, p1_lat, R); + arma_vec p2 = sph2cart(p2_lon, p2_lat, R); + arma_vec p3 = sph2cart(p3_lon, p3_lat, R); + arma_vec d1 = p1 - p2; + arma_vec d2 = p3 - p2; + precision_t n1 = std::sqrt(d1(0) * d1(0) + + d1(1) * d1(1) + + d1(2) * d1(2)); + precision_t n2 = std::sqrt(d2(0) * d2(0) + + d2(1) * d2(1) + + d2(2) * d2(2)); + d1 = d1 / n1; + d2 = d2 / n2; + precision_t angle = std::acos(d1(0) * d2(0) + + d1(1) * d2(1) + + d1(2) * d2(2)); + + return angle; +} + +// --------------------------------------------------------- +// bin edges +// --------------------------------------------------------- + +arma_vec calc_bin_edges(arma_vec centers) { + + int64_t nPts = centers.n_elem; + arma_vec edges(nPts + 1); + + precision_t dc = centers(1) - centers(0); + + edges(0) = centers(0) - dc / 2.0; + edges(1) = centers(0) + dc / 2.0; + + for (int64_t i = 2; i < nPts + 1; i++) + edges(i) = 2 * centers(i - 1) - edges(i - 1); + + return edges; +} + +// --------------------------------------------------------- +// bin edges +// --------------------------------------------------------- + +arma_mat calc_bin_edges(arma_mat centers, bool DoX) { + + // X is first dimension (row), Y is second dimension (col) + + int64_t nX = centers.n_rows; + int64_t nY = centers.n_cols; + arma_mat edges; + arma_vec centers1d; + + if (DoX) { + if (verbose > 2) + std::cout << " --> x\n"; + + edges.resize(nX + 1, nY); + + for (int64_t j = 0; j < nY; j++) { + centers1d = centers.col(j); + edges.col(j) = calc_bin_edges(centers1d); + } + } else { + if (verbose > 2) + std::cout << " --> y\n"; + + edges.resize(nX, nY + 1); + + for (int64_t i = 0; i < nX; i++) { + centers1d = centers.row(i).as_col(); + edges.row(i) = calc_bin_edges(centers1d).as_row(); + } + } + + return edges; +} + +// --------------------------------------------------------- +// bin widths +// --------------------------------------------------------- + +arma_vec calc_bin_widths(arma_vec edges) { + + int64_t nPts = edges.n_elem - 1; + arma_vec widths(nPts); + + for (int64_t i = 0; i < nPts; i++) + widths(i) = edges(i + 1) - edges(i); + + return widths; +} + +// --------------------------------------------------------- +// bin widths 2d +// --------------------------------------------------------- + +arma_mat calc_bin_widths(arma_mat edges, bool DoX) { + + int64_t nX = edges.n_rows; + int64_t nY = edges.n_cols; + + arma_mat widths; + arma_vec edges1d; + + if (DoX) { + if (verbose > 2) + std::cout << " --> x\n"; + + nX--; + widths.resize(nX, nY); + + for (int64_t j = 0; j < nY; j++) { + edges1d = edges.col(j); + widths.col(j) = calc_bin_widths(edges1d); + } + } else { + if (verbose > 2) + std::cout << " --> y\n"; + + nY--; + widths.resize(nX, nY); + + for (int64_t i = 0; i < nX; i++) { + edges1d = edges.row(i).as_col(); + widths.row(i) = calc_bin_widths(edges1d).as_row(); + } + } + + return widths; +} + +/**SOME PROJECTION AND GRADIENT CODE **/ +// --------------------------------------------------------- +// +// --------------------------------------------------------- + +arma_vec limiter_mc(arma_vec &left, + arma_vec &right, + int64_t nPts, + int64_t nGCs) { + + precision_t beta = 0.8; + + arma_vec s = left % right; + arma_vec combined = (left + right) * 0.5; + + left = left * beta; + right = right * beta; + arma_vec limited = left; + + for (int64_t i = 1; i < nPts + 2 * nGCs - 1; i++) { + if (s(i) < 0) { + // Sign < 0 means opposite signed left and right: + limited(i) = 0.0; + } else { + if (left(i) > 0 && right(i) > 0) { + if (right(i) < limited(i)) + limited(i) = right(i); + + if (combined(i) < limited(i)) + limited(i) = combined(i); + } else { + if (right(i) > limited(i)) + limited(i) = right(i); + + if (combined(i) > limited(i)) + limited(i) = combined(i); + } + } + } + + return limited; +} + +void print(arma_vec values) { + int64_t nP = values.n_elem; + + for (int64_t i = 0; i < nP; i++) + std::cout << values(i) << " "; + + std::cout << "\n"; +} + +// --------------------------------------------------------- +// calc gradients at centers +// - values and x defined at centers +// --------------------------------------------------------- + +arma_vec calc_grad_1d(arma_vec &values, + arma_vec &x, + int64_t nPts, + int64_t nGCs) { + + arma_vec gradients = values * 0.0; + arma_vec gradL = values * 0.0; + arma_vec gradR = values * 0.0; + + precision_t factor1 = 0.625; + precision_t factor2 = 0.0416667; + precision_t h; + + int64_t i; + arma_vec hv = values * 0.0; + + i = nGCs - 1; + h = 2.0 / (x(i + 1) - x(i)); + gradR(i) = h * (factor1 * (values(i + 1) - values(i)) - + factor2 * (values(i + 2) - values(i - 1))); + gradL(i) = (values(i) - values(i - 1)) / (x(i) - x(i - 1)); + + // This is attempting to vectorize the problem, but it seems to be slower? + // int64_t iS = nGCs; + // int64_t iE = nPts + nGCs - 1; + // hv.rows(iS, iE) = 2.0 / (x.rows(iS, iE) - x.rows(iS-1, iE-1)); + // gradL.rows(iS, iE) = hv.rows(iS,iE) % (factor1 * (values.rows(iS, iE) - + // values.rows(iS-1, iE-1)) - + // factor2 * (values.rows(iS+1, iE+1) - + // values.rows(iS-2, iE-2))); + // hv.rows(iS, iE) = 2.0 / (x.rows(iS+1, iE+1) - x.rows(iS, iE)); + // gradR.rows(iS, iE) = hv.rows(iS,iE) % (factor1 * (values.rows(iS+1, iE+1) - + // values.rows(iS, iE)) - + // factor2 * (values.rows(iS+2, iE+2) - + // values.rows(iS-1, iE-1))); + + for (i = nGCs; i < nPts + nGCs; i++) { + h = 2.0 / (x(i) - x(i - 1)); + gradL(i) = h * (factor1 * (values(i) - values(i - 1)) - + factor2 * (values(i + 1) - values(i - 2))); + h = 2.0 / (x(i + 1) - x(i)); + gradR(i) = h * (factor1 * (values(i + 1) - values(i)) - + factor2 * (values(i + 2) - values(i - 1))); + } + + i = nPts + nGCs; + h = 2.0 / (x(i) - x(i - 1)); + gradL(i) = h * (factor1 * (values(i) - values(i - 1)) - + factor2 * (values(i + 1) - values(i - 2))); + gradR(i) = (values(i + 1) - values(i)) / (x(i + 1) - x(i)); + + gradients = limiter_mc(gradL, gradR, nPts, nGCs); + + return gradients; +} + +// --------------------------------------------------------- +// calc gradients at centers for 2d matrices +// - values and x defined at centers +// --------------------------------------------------------- + +arma_mat calc_grad(arma_mat values, + arma_mat x, + int64_t nGCs, + bool DoX) { + + arma_mat v2d, x2d; + + if (DoX) { + v2d = values; + x2d = x; + } else { + v2d = values.t(); + x2d = x.t(); + } + + int64_t nX = v2d.n_rows; + int64_t nY = v2d.n_cols; + arma_mat grad2d = v2d * 0.0; + + int64_t nPts = nX - 2 * nGCs; + arma_vec values1d(nX); + arma_vec x1d(nX); + + for (int64_t j = 1; j < nY - 1; j++) { + values1d = v2d.col(j); + x1d = x2d.col(j); + grad2d.col(j) = calc_grad_1d(values1d, x1d, nPts, nGCs); + } + + arma_mat gradients; + + if (DoX) + gradients = grad2d; + else + gradients = grad2d.t(); + + return gradients; +} + +// --------------------------------------------------------- +// Project gradients + values to the right face, from the left +// returned values are on the i - 1/2 edges +// (between i-1 and i cell center) +// --------------------------------------------------------- + +arma_mat project_from_left(arma_mat values, + arma_mat gradients, + arma_mat x_centers, + arma_mat x_edges, + int64_t nGCs) { + + int64_t nX = values.n_rows; + int64_t nY = values.n_cols; + + // Define at edges: + arma_mat projected(nX + 1, nY); + projected.zeros(); + + // no gradient in the 0 or iEnd cells + for (int64_t j = 0; j < nY; j++) { + for (int64_t i = 1; i < nX - 1; i++) { + projected(i + 1, j) = values(i, j) + + gradients(i, j) * (x_edges(i + 1, j) - x_centers(i, j)); + } + + projected(1, j) = projected(2, j); + projected(0, j) = projected(1, j); + projected(nX, j) = projected(nX - 1, j); + } + + return projected; +} + +// --------------------------------------------------------- +// Project gradients + values to the left face, from the right +// returned values are on the i - 1 edges +// (between i-1 and i cell center) +// --------------------------------------------------------- + +arma_mat project_from_right(arma_mat values, + arma_mat gradients, + arma_mat x_centers, + arma_mat x_edges, + int64_t nGCs) { + int64_t nX = values.n_rows; + int64_t nY = values.n_cols; + + // Define at edges: + arma_mat projected(nX + 1, nY); + projected.zeros(); + + // no gradient in the 0 or iEnd cells + for (int64_t j = 0; j < nY; j++) { + for (int64_t i = 1; i < nX - 1; i++) { + projected(i, j) = values(i, j) + + gradients(i, j) * (x_edges(i, j) - x_centers(i, j)); + } + + projected(0, j) = projected(1, j); + projected(nX - 1, j) = projected(nX - 2, j); + projected(nX, j) = projected(nX - 1, j); + } + + return projected; +} + +// --------------------------------------------------------- +// Limiter on values +// projected is assumed to be on the edge between the +// i-1 and i cell (i-1/2) +// limited is returned at edges +// --------------------------------------------------------- + +arma_vec limiter_value(arma_vec projected, + arma_vec values, + int64_t nPts, + int64_t nGCs) { + + int64_t iStart = 0; + int64_t iEnd = nPts + 2 * nGCs; + + arma_vec limited = projected; + + precision_t mini, maxi; + + for (int64_t i = iStart + 1; i < iEnd - 1; i++) { + + mini = values(i - 1); + + if (values(i) < mini) + mini = values(i); + + maxi = values(i - 1); + + if (values(i) > maxi) + maxi = values(i); + + if (limited(i) < mini) + limited(i) = mini; + + if (limited(i) > maxi) + limited(i) = maxi; + } + + return limited; +} + +// --------------------------------------------------------- +// take gradients and project to all edges +// --------------------------------------------------------- + +projection_struct project_to_edges(arma_mat &values, + arma_mat &x_centers, arma_mat &x_edges, + arma_mat &y_centers, arma_mat &y_edges, + int64_t nGCs) { + + int64_t nX = values.n_rows; + int64_t nY = values.n_cols; + + projection_struct proj; + + proj.gradLR = calc_grad(values, x_centers, nGCs, true); + proj.gradDU = calc_grad(values.t(), y_centers.t(), nGCs, true).t(); + + proj.R = project_from_left(values, proj.gradLR, + x_centers, x_edges, nGCs); + // Left side of edge from left + proj.L = project_from_right(values, proj.gradLR, + x_centers, x_edges, nGCs); + // Up side of edge from down (left) + proj.U = project_from_left(values.t(), proj.gradDU.t(), + y_centers.t(), y_edges.t(), nGCs) + .t(); + // Down side of edge from up (right) + proj.D = project_from_right(values.t(), proj.gradDU.t(), + y_centers.t(), y_edges.t(), nGCs) + .t(); + + return proj; +} + +/**** SOME INITIALIZATION FUNCTION, NOT CORE ****/ +// --------------------------------------------------------- +// initial rho: initialize the whole field to be 2.0 +// --------------------------------------------------------- +arma_mat init_rho(arma_mat &x, + arma_mat &y) { + + int64_t nX = x.n_rows; + int64_t nY = x.n_cols; + + arma_mat rho(nX, nY); + arma_mat r; + + r = sqrt((x - 0.0) % (x - 0.0) + (y - 0.0) % (y - 0.0)); + rho.fill(1.0); + rho.elem( find( r < 0.25)).fill(1.2); + + return rho; +} + +// --------------------------------------------------------- +// initial velocity: initialize zero velocity +// --------------------------------------------------------- + +arma_mat init_vel(arma_mat &x, + arma_mat &y) { + int64_t nX = x.n_rows; + int64_t nY = x.n_cols; + arma_mat vel(nX, nY); + // all cells positive to right: + vel.zeros(); + vel.fill(0.5); + return vel; +} + + +// --------------------------------------------------------- +// initial values +// --------------------------------------------------------- + +arma_mat init_value(arma_mat &x, + arma_mat &y, + precision_t inVal) { + int64_t nX = x.n_rows; + int64_t nY = x.n_cols; + arma_mat val(nX, nY); + val.fill(inVal); + return val; +} + +arma_mat init_vel2(arma_mat &x, + arma_mat &y) { + int64_t nX = x.n_rows; + int64_t nY = x.n_cols; + arma_mat vel(nX, nY); + // all cells positive to right: + vel.zeros(); + //vel.fill(1.0); + return vel; +} + +// --------------------------------------------------------- +// initial temp (E): constant total energy +// THIS IS NOT e but E, the total energy +// --------------------------------------------------------- + +arma_mat init_temp(arma_mat &x, + arma_mat &y) { + int64_t nX = x.n_rows; + int64_t nY = x.n_cols; + + arma_mat temp(nX, nY); + temp.fill(100.0); + return temp; +} + +// --------------------------------------------------------- +// Calculate the max speed in the x and y directions +// --------------------------------------------------------- + +void calc_cmax(arma_mat &xVel, + arma_mat &yVel, + arma_mat &temp, + arma_mat &xMax, + arma_mat &yMax) { + + if (verbose > 2) + std::cout << " --> calc_max\n"; + + arma_mat xVel2, yVel2; + + xVel2 = xVel % xVel; + yVel2 = yVel % yVel; + xMax = sqrt(xVel2) + sqrt(cKb / mmm * temp); + yMax = sqrt(yVel2) + sqrt(cKb / mmm * temp); + + return; +} + + +// --------------------------------------------------------- +// Set Boundary Conditions +// --------------------------------------------------------- + +void set_bcs(arma_mat &rho, + arma_mat &xVel, + arma_mat &yVel, + arma_mat &temp, + grid_struct gridC) { + + if (verbose > 2) + std::cout << " --> set_bcs\n"; + + // ------------------------------------------------ + // Exchange messages (set BCs, really): + for (int64_t i = gridC.iXfirst_; i < gridC.iXlast_; i++) { + for (int64_t j = 0; j < gridC.nGCs; j++) { + // bottom bc: + rho(i, gridC.iYfirst_ - 1 - j) = rho(i, gridC.iYlast_ - 1 - j); + // top bc: + rho(i, gridC.iYlast_ + j) = rho(i, gridC.iYfirst_ + j); + } + } + + for (int64_t j = gridC.iYfirst_; j < gridC.iYlast_; j++) { + for (int64_t i = 0; i < gridC.nGCs; i++) { + // left bc: + rho(gridC.iXfirst_ - 1 - i, j) = rho(gridC.iXlast_ - 1 - i, j); + // right bc: + rho(gridC.iXlast_ + i, j) = rho(gridC.iXfirst_ + i, j); + } + } + + return; +} + +// --------------------------------------------------------- +// Convert vector from Alat, Alon to Axi, Anu +// -> Using equation (7) of Ronchi et al: +// --------------------------------------------------------- + +void convert_vector_ll_to_xn(arma_mat aLon, + arma_mat aLat, + arma_mat &aXi, + arma_mat &aNu, + grid_struct grid) { + + // Ronchi defines aPhi = aLon, aTheta = -aLat + aXi = -grid.Axt % aLat + grid.Axp % aLon; + aNu = -grid.Ant % aLat + grid.Anp % aLon; + return; +} + +// --------------------------------------------------------- +// Convert vector from sphere to cube +// -> Using equation (32) of Nair: +// --------------------------------------------------------- + +void convert_vector_sphere_to_cube(arma_mat u, + arma_mat v, + arma_mat &u1, + arma_mat &u2, + grid_struct grid) { + + // Ronchi defines aPhi = aLon, aTheta = -aLat + u1 = grid.s2c11 % u + grid.s2c12 % v; + u2 = grid.s2c21 % u + grid.s2c22 % v; + return; +} + +// --------------------------------------------------------- +// Convert vector cube to sphere +// -> Using equation (32) of Nair: +// --------------------------------------------------------- + +void convert_vector_cube_to_sphere(arma_mat u1, + arma_mat u2, + arma_mat &u, + arma_mat &v, + grid_struct grid) { + + // Ronchi defines aPhi = aLon, aTheta = -aLat + u = grid.c2s11 % u1 + grid.c2s12 % u2; + v = grid.c2s21 % u1 + grid.c2s22 % u2; + return; +} + +// --------------------------------------------------------- +// Convert vector from Alat, Alon to Axi, Anu +// -> Using equation (7) of Ronchi et al: +// --------------------------------------------------------- + +void convert_vector_xn_to_ll(arma_mat aXi, + arma_mat aNu, + arma_mat &aLon, + arma_mat &aLat, + grid_struct grid) { + + // Ronchi defines aPhi = aLon, aTheta = -aLat + aLat = -(grid.Atx % aXi + grid.Atn % aNu); + aLon = grid.Apx % aXi + grid.Apn % aNu; + + return; +} + +// --------------------------------------------------------- +// Convert vector from Alat, Alon to Axi, Anu +// -> Using equation (7) of Ronchi et al: +// --------------------------------------------------------- + +arma_mat calc_angle_between_coords(grid_struct grid) { + + arma_mat e1Lat, e1Lon, e2Lat, e2Lon, m, one, zero; + arma_mat e1dote2, alpha; + + m.resize(grid.nXt, grid.nYt); + one.resize(grid.nXt, grid.nYt); + one.fill(1.0); + zero.resize(grid.nXt, grid.nYt); + zero.fill(0.0); + + // define e1 as the LR (xi) direction: + e1Lat.resize(grid.nXt, grid.nYt); + e1Lon.resize(grid.nXt, grid.nYt); + convert_vector_xn_to_ll(one, zero, e1Lon, e1Lat, grid); + m = sqrt(e1Lon % e1Lon + e1Lat % e1Lat); + e1Lon = e1Lon / m; + e1Lat = e1Lat / m; + + // define e2 as the DU (nu) direction: + e2Lat.resize(grid.nXt, grid.nYt); + e2Lon.resize(grid.nXt, grid.nYt); + convert_vector_xn_to_ll(zero, one, e2Lon, e2Lat, grid); + m = sqrt(e2Lon % e2Lon + e2Lat % e2Lat); + e2Lon = e2Lon / m; + e2Lat = e2Lat / m; + + alpha = acos(e1Lat % e2Lat + e1Lon % e2Lon); + + return alpha; +} + +// --------------------------------------------------------- +// Convert vector from Alat, Alon to Axi, Anu +// -> Using equation (7) of Ronchi et al: +// --------------------------------------------------------- + +void calc_norms(grid_struct &grid) { + + arma_mat e1Lat, e1Lon, e2Lat, e2Lon, m, one, zero; + + grid.nXiLon.resize(grid.nXt, grid.nYt); + grid.nXiLat.resize(grid.nXt, grid.nYt); + grid.nNuLon.resize(grid.nXt, grid.nYt); + grid.nNuLat.resize(grid.nXt, grid.nYt); + + m.resize(grid.nXt, grid.nYt); + one.resize(grid.nXt, grid.nYt); + one.fill(1.0); + zero.resize(grid.nXt, grid.nYt); + zero.fill(0.0); + + // define e1 as the LR (xi) direction: + e1Lat.resize(grid.nXt, grid.nYt); + e1Lon.resize(grid.nXt, grid.nYt); + convert_vector_xn_to_ll(one, zero, e1Lon, e1Lat, grid); + m = sqrt(e1Lon % e1Lon + e1Lat % e1Lat); + + // Rotate by 90 deg (CCW) to get the norm: + grid.nNuLon = -e1Lat / m; + grid.nNuLat = e1Lon / m; + + // define e2 as the DU (nu) direction: + e2Lat.resize(grid.nXt, grid.nYt); + e2Lon.resize(grid.nXt, grid.nYt); + convert_vector_xn_to_ll(zero, one, e2Lon, e2Lat, grid); + m = sqrt(e2Lon % e2Lon + e2Lat % e2Lat); + // Rotate by 90 deg (CW) to get the norm: + grid.nXiLon = e2Lat / m; + grid.nXiLat = -e2Lon / m; + + return; +} + + +// --------------------------------------------------------- +// Update States +// --------------------------------------------------------- + +void update_states(arma_mat rho, + arma_mat &xVel, + arma_mat &yVel, + arma_mat &temp, + arma_mat &drhodt, + arma_mat &dlonVeldt, + arma_mat &dlatVeldt, + arma_mat &dtempdt, + grid_struct gridC, + grid_struct gridL, + grid_struct gridD, + precision_t dt) { + + arma_mat xMomentum, yMomentum; + arma_mat rhoE, energy, vel2; + + precision_t cv = 1500.0; + + if (verbose > 2) + std::cout << " --> update_states\n"; + + // Derived variables: + xMomentum = rho % xVel; // x1momentum, pure scalar field + yMomentum = rho % yVel; // y1momentum, pure scalar field + rhoE = rho % temp; + + vel2 = xVel % xVel + yVel % yVel; + //energy = rho % (0.5 * vel2 + cv * temp); + energy = cv * rho % temp; + + /** Initialize projection constructs */ + static projection_struct rhoP; + static projection_struct xMomentumP, xVelP; + static projection_struct yMomentumP, yVelP; + static projection_struct energyP; + static projection_struct tempP; + + // They are all pure scalar fields without sqrt(g) + static arma_mat totaleL, totaleR, totaleD, totaleU; + static arma_mat velL2, velR2, velD2, velU2; + static arma_mat pressureL, pressureR, pressureD, pressureU; + + arma_mat dxVeldt = xVel * 0.0; + arma_mat dyVeldt = yVel * 0.0; + + dlonVeldt = dxVeldt * 0.0 + 1; + dlatVeldt = dyVeldt * 0.0 + 1; + + static arma_mat velNormL, velNormR, velNormU, velNormD; + + /** Initialize Flux and Wave Speed Storages */ + static arma_mat eq1FluxLR, eq1FluxDU; + static arma_mat eq1FluxL, eq1FluxR, eq1FluxD, eq1FluxU; + static arma_mat eq2FluxLR, eq2FluxDU; + static arma_mat eq2FluxL, eq2FluxR, eq2FluxD, eq2FluxU; + static arma_mat eq3FluxLR, eq3FluxDU; + static arma_mat eq3FluxL, eq3FluxR, eq3FluxD, eq3FluxU; + static arma_mat eq4FluxLR, eq4FluxDU; + static arma_mat eq4FluxL, eq4FluxR, eq4FluxD, eq4FluxU; + + arma_mat wsL, wsR, wsD, wsU, wsLR, wsDU; + + arma_mat diff; // for Riemann Solver + + if (verbose > 3) + std::cout << " ---> Projecting\n"; + + rhoP = project_to_edges(rho, gridC.xi, gridL.xi, gridC.nu, gridD.nu, + gridC.nGCs); + // project the lon / lat velocities to the edges: + xVelP = project_to_edges(xVel, gridC.xi, gridL.xi, gridC.nu, gridD.nu, + gridC.nGCs); + yVelP = project_to_edges(yVel, gridC.xi, gridL.xi, gridC.nu, gridD.nu, + gridC.nGCs); + xMomentumP = project_to_edges(xMomentum, gridC.xi, gridL.xi, gridC.nu, gridD.nu, + gridC.nGCs); + yMomentumP = project_to_edges(yMomentum, gridC.xi, gridL.xi, gridC.nu, gridD.nu, + gridC.nGCs); + energyP = project_to_edges(energy, gridC.xi, gridL.xi, gridC.nu, gridD.nu, + gridC.nGCs); + tempP = project_to_edges(temp, gridC.xi, gridL.xi, gridC.nu, gridD.nu, + gridC.nGCs); + + if (verbose > 3) + std::cout << " ---> Derived values\n"; + + velL2 = (xVelP.L % xVelP.L + yVelP.L % yVelP.L); + velR2 = (xVelP.R % xVelP.R + yVelP.R % yVelP.R); + velD2 = (xVelP.D % xVelP.D + yVelP.D % yVelP.D); + velU2 = (xVelP.U % xVelP.U + yVelP.U % yVelP.U); + + precision_t k = 1.38e-23; + // let's be Oxygen: + precision_t mass = 16.0 * 1.67e-27; + pressureL = k / mass * (rhoP.L % tempP.L); + pressureR = k / mass * (rhoP.R % tempP.R); + pressureD = k / mass * (rhoP.D % tempP.D); + pressureU = k / mass * (rhoP.U % tempP.U); + + arma_mat pressureLR = (pressureL + pressureR) / 2; + arma_mat pressureDU = (pressureD + pressureU) / 2; + + if (verbose > 3) + std::cout << " ---> Normal Velocities\n"; + + // Calculate the normal velocity at the boundaries: + velNormL = xVelP.L % gridL.nXiLon + yVelP.L % gridL.nXiLat; + velNormR = xVelP.R % gridL.nXiLon + yVelP.R % gridL.nXiLat; + velNormU = xVelP.U % gridD.nNuLon + yVelP.U % gridD.nNuLat; + velNormD = xVelP.D % gridD.nNuLon + yVelP.D % gridD.nNuLat; + + if (verbose > 3) + std::cout << " ---> Fluxes eq 1\n"; + + // Flux calculated from the left of the edge + eq1FluxL = rhoP.L % velNormL; + // Flux calculated from the right of the edge + eq1FluxR = rhoP.R % velNormR; + // Flux calculated from the down of the edge + eq1FluxD = rhoP.D % velNormD; + // Flux calculated from the up of the edge + eq1FluxU = rhoP.U % velNormU; + + if (verbose > 3) + std::cout << " ---> Fluxes eq 2\n"; + + eq2FluxL = (xMomentumP.L % velNormL); + eq2FluxR = (xMomentumP.R % velNormR); + eq2FluxD = (xMomentumP.D % velNormD); + eq2FluxU = (xMomentumP.U % velNormU); + + if (verbose > 3) + std::cout << " ---> Fluxes eq 3\n"; + + eq3FluxL = (yMomentumP.L % velNormL); + eq3FluxR = (yMomentumP.R % velNormR); + eq3FluxD = (yMomentumP.D % velNormD); + eq3FluxU = (yMomentumP.U % velNormU); + + eq4FluxL = energyP.L % velNormL; + eq4FluxR = energyP.R % velNormR; + eq4FluxD = energyP.D % velNormD; + eq4FluxU = energyP.U % velNormU; + + // ------------------------------------------------ + // Calculate the wave speed for the diffusive flux: + // In Reference velocities + if (verbose > 3) + std::cout << " ---> Diffusive Fluxes\n"; + + wsL = sqrt(velL2) + sqrt(cGamma * (cGamma - 1) * tempP.L); + wsR = sqrt(velR2) + sqrt(cGamma * (cGamma - 1) * tempP.R); + wsD = sqrt(velD2) + sqrt(cGamma * (cGamma - 1) * tempP.D); + wsU = sqrt(velU2) + sqrt(cGamma * (cGamma - 1) * tempP.U); + + wsLR = wsR; + + for (int64_t i = 0; i < gridC.nXt + 1; i++) { + for (int64_t j = 0; j < gridC.nYt; j++) { + if (wsL(i, j) > wsLR(i, j)) + wsLR(i, j) = wsL(i, j); + } + } + + wsDU = wsD; + + for (int64_t i = 0; i < gridC.nXt; i++) { + for (int64_t j = 0; j < gridC.nYt + 1; j++) { + if (wsU(i, j) > wsDU(i, j)) + wsDU(i, j) = wsU(i, j); + } + } + + // ------------------------------------------------ + // Calculate average flux at the edges (Rusanov Flux): + + if (verbose > 3) + std::cout << " ---> Averaging fluxes at edges\n"; + + diff = (rhoP.R - rhoP.L); + eq1FluxLR = (eq1FluxL + eq1FluxR) / 2 + 0.5 * wsLR % diff; + diff = (rhoP.U - rhoP.D); + eq1FluxDU = (eq1FluxD + eq1FluxU) / 2 + 0.5 * wsDU % diff; + + diff = (xMomentumP.R - xMomentumP.L); + eq2FluxLR = (eq2FluxL + eq2FluxR) / 2 + 0.5 * wsLR % diff; + diff = (xMomentumP.U - xMomentumP.D); + eq2FluxDU = (eq2FluxD + eq2FluxU) / 2 + 0.5 * wsDU % diff; + + diff = (yMomentumP.R - yMomentumP.L); + eq3FluxLR = (eq3FluxL + eq3FluxR) / 2 + 0.5 * wsLR % diff; + diff = (yMomentumP.U - yMomentumP.D); + eq3FluxDU = (eq3FluxD + eq3FluxU) / 2 + 0.5 * wsDU % diff; + + diff = (energyP.R - energyP.L); + eq4FluxLR = (eq4FluxL + eq4FluxR) / 2 + 0.5 * wsLR % diff; + diff = (energyP.U - energyP.D); + eq4FluxDU = (eq4FluxD + eq4FluxU) / 2 + 0.5 * wsDU % diff; + + // ------------------------------------------------ + // Update values: + if (verbose > 3) + std::cout << " ---> Updating equations of state\n"; + + precision_t dpdx, dpdn, pp, pm; + + arma_mat ax, an; + + ax = xVel * 0.0; + an = yVel * 0.0; + arma_mat dedt = xVel * 0.0; + + arma_mat rhoNew = rho; + + // Only deal with inner cell + for (int64_t j = gridC.iYfirst_; j < gridC.iYlast_; j++) { + for (int64_t i = gridC.iXfirst_; i < gridC.iXlast_; i++) { + precision_t rhoResidual_ij = (gridL.dln(i + 1, j) * eq1FluxLR(i + 1, j) - + gridL.dln(i, j) * eq1FluxLR(i, j) + + gridD.dlx(i, j + 1) * eq1FluxDU(i, j + 1) - + gridD.dlx(i, j) * eq1FluxDU(i, j)); + drhodt(i, j) = rhoResidual_ij / gridC.dS(i, j); + + rhoNew(i, j) = rho(i, j) + dt * drhodt(i, j); + + precision_t xMomentumResidual_ij = (gridL.dln(i + 1, j) * eq2FluxLR(i + 1, j) - + gridL.dln(i, j) * eq2FluxLR(i, j) + + gridD.dlx(i, j + 1) * eq2FluxDU(i, j + 1) - + gridD.dlx(i, j) * eq2FluxDU(i, j)); + dxVeldt(i, j) = xMomentumResidual_ij / gridC.dS(i, j) / rhoNew(i, j); + + precision_t yMomentumResidual_ij = (gridL.dln(i + 1, j) * eq3FluxLR(i + 1, j) - + gridL.dln(i, j) * eq3FluxLR(i, j) + + gridD.dlx(i, j + 1) * eq3FluxDU(i, j + 1) - + gridD.dlx(i, j) * eq3FluxDU(i, j)); + dyVeldt(i, j) = yMomentumResidual_ij / gridC.dS(i, j) / rhoNew(i, j); + + // Calculate the gradient in the potential in the cubesphere + // coordinate system: + dpdx = 1 / gridC.R * gridC.D(i, j) * + (pressureLR(i + 1, j) - pressureLR(i, j)) / gridC.dxi; + dpdn = 1 / gridC.R * gridC.X(i, j) * gridC.Y(i, j) / + gridC.D(i, j) * + (pressureDU(i, j + 1) - pressureDU(i, j)) / gridC.dnu; + ax(i, j) = (dpdx + dpdn) / rhoNew(i, j); + + dpdx = 1 / gridC.R * gridC.X(i, j) * gridC.Y(i, j) / + gridC.C(i, j) * (pressureLR(i + 1, j) - pressureLR(i, j)) / gridC.dxi; + dpdn = 1 / gridC.R * gridC.C(i, j) * + (pressureDU(i, j + 1) - pressureDU(i, j)) / gridC.dnu; + an(i, j) = (dpdx + dpdn) / rhoNew(i, j); + + precision_t energyResidual_ij = (gridL.dln(i + 1, j) * eq4FluxLR(i + 1, j) - + gridL.dln(i, j) * eq4FluxLR(i, j) + + gridD.dlx(i, j + 1) * eq4FluxDU(i, j + 1) - + gridD.dlx(i, j) * eq4FluxDU(i, j)); + dedt(i, j) = energyResidual_ij / gridC.dS(i, j); + + } + } + + dlatVeldt = dyVeldt - (ax % gridC.Atx + an % gridC.Atn); + dlonVeldt = dxVeldt + ax % gridC.Apx + an % gridC.Apn; + dtempdt = dedt / rhoNew / cv; + + return; +} + + +// --------------------------------------------------------- +// Main Code! +// --------------------------------------------------------- + +int main() { + precision_t dt; // Time Step + precision_t current_time = 0.0; // Initial Time 0 + precision_t total_time = 100.0; // Total simulation time + precision_t cfl = 0.75; + int iFace = 5; + + precision_t dtOut = total_time / 50.0; // Output Interval + precision_t dtReport = total_time / 200.0; // Output Interval + + int64_t iStep; // Iterator of Time Step + + int64_t nX = 50; // Number of x grid cells + int64_t nY = 50; // Number of y grid cells + int64_t nGCs = 2; // Number of ghost cells + + // Radius of Sphere + precision_t R = 10000.; + + if (verbose > 0) + std::cout << "> generating cubesphere cell center and metrics\n"; + + grid_struct gridC_eq = init_grid_equidistant(iFace, nX, nY, nGCs, R, 0.5, 0.5); + + output(gridC_eq.lat, "eq_lat.txt", false); + output(gridC_eq.lon, "eq_lon.txt", false); + output(gridC_eq.r, "eq_r.txt", false); + + grid_struct gridC = init_grid(iFace, nX, nY, nGCs, R, 0.5, 0.5); + grid_struct gridL = init_grid(iFace, nX + 1, nY, nGCs, R, 0.0, 0.5); + grid_struct gridD = init_grid(iFace, nX, nY + 1, nGCs, R, 0.5, 0.0); + + gridC.alpha = calc_angle_between_coords(gridC); + gridL.alpha = calc_angle_between_coords(gridL); + gridD.alpha = calc_angle_between_coords(gridD); + + gridC.sinAlpha = sin(gridC.alpha); + gridL.sinAlpha = sin(gridL.alpha); + gridD.sinAlpha = sin(gridD.alpha); + + calc_norms(gridC); + calc_norms(gridL); + calc_norms(gridD); + + /** State Initialization */ + /// Initialize Density + if (verbose > 0) + std::cout << "> initializing rho\n"; + + arma_mat rho = init_rho(gridC.xi, gridC.nu); // rho, pure scalar field + + /// Initialize Velocity and Momentum + if (verbose > 0) + std::cout << "> initializing vel\n"; + + // Initialize spherical velocity + // Supposed to be Longitudinal Velocity: + arma_mat vLon = init_value(gridC.xi, gridC.nu, 0.0); + // Supposed to be Latitudinal Velocity: + arma_mat vLat = init_value(gridC.xi, gridC.nu, 0.0); + + arma_mat vXi, vNu; + convert_vector_ll_to_xn(vLon, vLat, vXi, vNu, gridC); + + + /// Initialize total energy + if (verbose > 0) + std::cout << "> initializing energy\n"; + + arma_mat temp = init_temp(gridC.xi, gridC.nu); + + /** Output some pre-simulation results */ + if (verbose > 0) + std::cout << "-> outputting\n"; + + output(gridC.xi, "xi.txt", false); + output(gridC.dS, "dS.txt", false); + output(gridC.d, "d.txt", false); + output(gridC.nu, "nu.txt", false); + output(gridC.lat, "lat.txt", false); + output(gridC.lon, "lon.txt", false); + output(gridC.alpha, "alpha.txt", false); + output(rho, "rho.txt", false); + output(vLon, "vLon.txt", false); + output(vLat, "vLat.txt", false); + output(vXi, "vXi.txt", false); + output(vNu, "vNu.txt", false); + output(temp, "temp.txt", false); + iStep = 0; + + arma_mat xMax, yMax; + + arma_mat drhodt; + arma_mat dxVeldt, dyVeldt; + arma_mat dtempdt; + + drhodt.resize(gridC.nXt, gridC.nYt); + dxVeldt.resize(gridC.nXt, gridC.nYt); + dyVeldt.resize(gridC.nXt, gridC.nYt); + dtempdt.resize(gridC.nXt, gridC.nYt); + + arma_mat k1rho, k2rho, k3rho, k4rho, rhoInterK1, rhoInterK2, rhoInterK3; + k1rho.resize(gridC.nXt, gridC.nYt); + k2rho.resize(gridC.nXt, gridC.nYt); + k3rho.resize(gridC.nXt, gridC.nYt); + k4rho.resize(gridC.nXt, gridC.nYt); + rhoInterK1.resize(gridC.nXt, gridC.nYt); + rhoInterK2.resize(gridC.nXt, gridC.nYt); + rhoInterK3.resize(gridC.nXt, gridC.nYt); + + arma_mat k1vLon, k2vLon, k3vLon, k4vLon, vLonInterK1, vLonInterK2, vLonInterK3; + k1vLon.resize(gridC.nXt, gridC.nYt); + k2vLon.resize(gridC.nXt, gridC.nYt); + k3vLon.resize(gridC.nXt, gridC.nYt); + k4vLon.resize(gridC.nXt, gridC.nYt); + vLonInterK1.resize(gridC.nXt, gridC.nYt); + vLonInterK2.resize(gridC.nXt, gridC.nYt); + vLonInterK3.resize(gridC.nXt, gridC.nYt); + + arma_mat k1vLat, k2vLat, k3vLat, k4vLat, vLatInterK1, vLatInterK2, vLatInterK3; + k1vLat.resize(gridC.nXt, gridC.nYt); + k2vLat.resize(gridC.nXt, gridC.nYt); + k3vLat.resize(gridC.nXt, gridC.nYt); + k4vLat.resize(gridC.nXt, gridC.nYt); + vLatInterK1.resize(gridC.nXt, gridC.nYt); + vLatInterK2.resize(gridC.nXt, gridC.nYt); + vLatInterK3.resize(gridC.nXt, gridC.nYt); + + arma_mat k1temp, k2temp, k3temp, k4temp, tempInterK1, tempInterK2, tempInterK3; + k1temp.resize(gridC.nXt, gridC.nYt); + k2temp.resize(gridC.nXt, gridC.nYt); + k3temp.resize(gridC.nXt, gridC.nYt); + k4temp.resize(gridC.nXt, gridC.nYt); + tempInterK1.resize(gridC.nXt, gridC.nYt); + tempInterK2.resize(gridC.nXt, gridC.nYt); + tempInterK3.resize(gridC.nXt, gridC.nYt); + + while (current_time < total_time) { + + if (int((current_time - dt) / dtReport) != int((current_time ) / dtReport)) { + std::cout << "step : " << iStep << "; time : " << current_time << "\n"; + arma_vec amin, amax; + precision_t mini, maxi; + amin = rho.min(); + mini = amin.min(); + amax = max(rho, 1); + maxi = arma::max(amax); + std::cout << " -> min/max (rho) : " << mini << " / " << maxi << "\n"; + amin = temp.min(); + mini = amin.min(); + amax = max(temp, 1); + maxi = arma::max(amax); + std::cout << " -> min/max (temp) : " << mini << " / " << maxi << "\n"; + } + + calc_cmax(vLon, vLat, temp, xMax, yMax); + dt = calc_dt(gridC.dlx, gridC.dln, xMax, yMax, gridC.nGCs); + dt = cfl * dt; + std::cout << " dt: " << dt << "\n"; + + // k1 - start at t0, go to t+1/2 to figure out slope at t0 (k1) + update_states(rho, vLon, vLat, temp, + k1rho, k1vLon, k1vLat, k1temp, + gridC, gridL, gridD, dt / 2); + // Take 1/2 step to figure out t+1/2 values, using k1: + rhoInterK1 = rho + k1rho * dt / 2; + vLonInterK1 = vLon + k1vLon * dt / 2; + vLatInterK1 = vLat + k1vLat * dt / 2; + tempInterK1 = temp + k1temp * dt / 2; + set_bcs(rhoInterK1, vLonInterK1, vLatInterK1, tempInterK1, gridC); + + // k2 - start at t+1/2, go to t+1 to figure out slope at t+1/2 (k2) + update_states(rhoInterK1, vLonInterK1, vLatInterK1, tempInterK1, + k2rho, k2vLon, k2vLat, k2temp, + gridC, gridL, gridD, dt / 2); + // Take 1/2 step to figure out t+1/2 values, using k2: + rhoInterK2 = rho + k2rho * dt / 2; + vLonInterK2 = vLon + k2vLon * dt / 2; + vLatInterK2 = vLat + k2vLat * dt / 2; + tempInterK2 = temp + k2temp * dt / 2; + set_bcs(rhoInterK2, vLonInterK2, vLatInterK2, tempInterK2, gridC); + + // k3 - start at t+1/2, go to t+1 to figure out slope at t+1/2 (k3) + update_states(rhoInterK2, vLonInterK2, vLatInterK2, tempInterK2, + k3rho, k3vLon, k3vLat, k3temp, + gridC, gridL, gridD, dt / 2); + // Take full step to figure out k4, using k3 slope: + rhoInterK3 = rho + k3rho * dt; + vLonInterK3 = vLon + k3vLon * dt; + vLatInterK3 = vLat + k3vLat * dt; + tempInterK3 = temp + k3rho * dt; + set_bcs(rhoInterK3, vLonInterK3, vLatInterK3, tempInterK3, gridC); + + // k4 - start at t+1, go to t+2 to figure out slope at t+1 (k4) + update_states(rhoInterK3, vLonInterK3, vLatInterK3, tempInterK3, + k4rho, k4vLon, k4vLat, k4temp, + gridC, gridL, gridD, dt); + + rho = rho - dt / 6 * (k1rho + 2 * k2rho + 2 * k3rho + k4rho); + vLon = vLon - dt / 6 * (k1vLon + 2 * k2vLon + 2 * k3vLon + k4vLon); + vLat = vLat - dt / 6 * (k1vLat + 2 * k2vLat + 2 * k3vLat + k4vLat); + temp = temp - dt / 6 * (k1temp + 2 * k2temp + 2 * k3temp + k4temp); + set_bcs(rho, vLon, vLat, temp, gridC); + + iStep++; + current_time += dt; + + if (verbose > 3) + std::cout << " ---> Outputing\n"; + + if (int((current_time - dt) / dtOut) != int((current_time ) / dtOut)) { + std::cout << "> Outputing at time : " << current_time << "\n"; + output(rho, "rho.txt", true); + output(vLon, "vLon.txt", true); + output(vLat, "vLat.txt", true); + output(temp, "temp.txt", true); + } + } + + return 0; +} diff --git a/edu/examples/Advection/cubesphere_test.cpp b/edu/examples/Advection/cubesphere_test.cpp index f201e00b..f23176aa 100644 --- a/edu/examples/Advection/cubesphere_test.cpp +++ b/edu/examples/Advection/cubesphere_test.cpp @@ -97,7 +97,6 @@ void output(arma_mat &values, outfile << nX << " " << nY << "\n"; } outfile << values; - outfile << "----"; outfile.close(); } @@ -801,7 +800,7 @@ arma_mat init_rho(arma_mat &x, r = sqrt((x - 0.0) % (x - 0.0) + (y - 0.0) % (y - 0.0)); rho.fill(1.0); - // rho.elem( find( r < 0.25)).fill(2.2); + rho.elem( find( r < 0.25)).fill(2.2); // rho.elem( find( r < 0.25)) = 2.25 - r.elem( find( r < 0.25)); return rho; @@ -857,7 +856,7 @@ int main() { precision_t dt = 0.0001; // Time Step precision_t current_time = 0.0; // Initial Time 0 - precision_t total_time = 2.0; // Total simulation time + precision_t total_time = 0.1; // Total simulation time precision_t cfl = 0.1; // CFL Number precision_t gamma = 5.0 / 3.0; // Specific ratio of heat @@ -1274,4 +1273,4 @@ int main() output(xVelSph_output, "xvel_sph.txt", true); output(yVelSph_output, "yvel_sph.txt", true); return 0; -} \ No newline at end of file +} diff --git a/include/cubesphere.h b/include/cubesphere.h index e99036a6..ecca2aed 100644 --- a/include/cubesphere.h +++ b/include/cubesphere.h @@ -15,32 +15,32 @@ namespace CubeSphere { /// The normalized origins of each face of the cube (i.e. corner) static const arma_mat ORIGINS = { - {-1.0, -1.0, -1.0}, - { 1.0, -1.0, -1.0}, - { 1.0, 1.0, -1.0}, - {-1.0, 1.0, -1.0}, - {-1.0, -1.0, -1.0}, - { 1.0, -1.0, 1.0} + {-1.0, -1.0, -1.0}, + { 1.0, -1.0, -1.0}, + { 1.0, 1.0, -1.0}, + {-1.0, 1.0, -1.0}, + {-1.0, -1.0, -1.0}, + { 1.0, -1.0, 1.0} }; /// Normalized right steps in cube static const arma_mat RIGHTS = { - { 2.0, 0.0, 0.0}, - { 0.0, 2.0, 0.0}, - {-2.0, 0.0, 0.0}, - { 0.0, -2.0, 0.0}, - { 0.0, 2.0, 0.0}, - { 0.0, 2.0, 0.0} + { 2.0, 0.0, 0.0}, + { 0.0, 2.0, 0.0}, + {-2.0, 0.0, 0.0}, + { 0.0, -2.0, 0.0}, + { 0.0, 2.0, 0.0}, + { 0.0, 2.0, 0.0} }; /// Normalized right steps in cube static const arma_mat UPS = { - { 0.0, 0.0, 2.0}, - { 0.0, 0.0, 2.0}, - { 0.0, 0.0, 2.0}, - { 0.0, 0.0, 2.0}, - { 2.0, 0.0, 0.0}, - {-2.0, 0.0, 0.0} + { 0.0, 0.0, 2.0}, + { 0.0, 0.0, 2.0}, + { 0.0, 0.0, 2.0}, + { 0.0, 0.0, 2.0}, + { 2.0, 0.0, 0.0}, + {-2.0, 0.0, 0.0} }; } // CubeSphere:: diff --git a/include/grid.h b/include/grid.h index 9bb98b27..daae63d9 100644 --- a/include/grid.h +++ b/include/grid.h @@ -7,6 +7,40 @@ #include #include "mpi.h" +// ---------------------------------------------------------------------------- +// This structure needs to be defined outside of the grid, since we can just +// pass this stuff to the solver. +// ---------------------------------------------------------------------------- + +struct cubesphere_chars { + // For convenience, store the grid size: + int64_t nXt, nYt, nGCs; + int64_t iXfirst_, iXlast_; + int64_t iYfirst_, iYlast_; + + // These are for Ronchi et al., JCP 124, 93-114, 1996 + arma_mat X, Y, Z, C, D, d; + // These are the only things that depend on altitude: + arma_cube dlx, dln, dS; + // In theory, the radius is just a 1D vector: + arma_vec R; + // xi is the LR direction + // nu is the UD direction + arma_mat xi, nu; + // for the equal-angle grid, we can just use these: + precision_t dxi, dnu; + arma_mat Apn, Apx, Atn, Atx; + arma_mat Axt, Axp, Ant, Anp; + + // These are for computing normals to the cell edges (horizontal) + arma_mat nXiLon; + arma_mat nXiLat; + arma_mat nNuLon; + arma_mat nNuLat; + arma_mat lat, lon; +}; + + // ---------------------------------------------------------------------------- // Grid class // ---------------------------------------------------------------------------- @@ -65,6 +99,10 @@ class Grid { arma_cube g11_upper_Down, g12_upper_Down, g21_upper_Down, g22_upper_Down; arma_cube sqrt_g_Down; + cubesphere_chars cubeC, cubeL, cubeD; + + // The magnetic latitude and altitude need to be defined better. This should be the angle between + // magnetic equator and the point, but sometimes it is invariant latitude. // These define the magnetic grid (only defined for a dipole grid): // The magnetic latitude is the angle between the magnetic equator and the point. arma_cube magLat_scgc, magY_scgc; @@ -320,6 +358,29 @@ class Grid { void create_sphere_grid(Quadtree quadtree); void create_cubesphere_connection(Quadtree quadtree); void create_cubesphere_grid(Quadtree quadtree); + + // These two go together, since one builds the angles and the + // other scales by the radius: + void init_cubesphere_grid(Quadtree quadtree, + arma_vec dr, + arma_vec du, + arma_vec ll, + precision_t left_off, + precision_t down_off, + cubesphere_chars &cubeX); + void scale_cube_by_radius(cubesphere_chars &cubeX); + + void convert_vector_xn_to_ll(arma_mat aXi, + arma_mat aNu, + arma_mat &aLon, + arma_mat &aLat, + cubesphere_chars grid); + void convert_vector_ll_to_xn(arma_mat aLon, + arma_mat aLat, + arma_mat &aXi, + arma_mat &aNu, + cubesphere_chars grid); + void create_altitudes(Planets planet); void fill_grid_bfield(Planets planet); bool read_restart(std::string dir); @@ -334,6 +395,27 @@ class Grid { // Support functions: void calc_dipole_grid_spacing(Planets planet); + void calc_alt_dipole_grid_spacing(); + void calc_lat_dipole_grid_spacing(); + void calc_long_dipole_grid_spacing(); + void fill_field_lines(arma_vec baseLats, precision_t min_altRe, + precision_t Gamma, Planets planet, + bool isCorner); + void dipole_alt_edges(Planets planet, precision_t min_altRe); + // get the latitude spacing given the quadtree start & size, and the latitude limits + // extent: quadtree up + // origin: quadtree origin + // upper_lim: upper latitude limit (input) + // lower_lim: lower latitude limit (from min_apex) + // nLats: number of latitudes (nY) + // spacing_factor: (not supported yet), so always 1.0. Will adjust baselat spacing, eventually. + arma_vec baselat_spacing(precision_t extent, + precision_t origin, + precision_t upper_lim, + precision_t lower_lim, + // int16_t nLats, + precision_t spacing_factor); + // Update ghost cells with values from other processors void exchange(arma_cube &data, const bool pole_inverse); @@ -425,8 +507,8 @@ class Grid { bool set_interpolation_coefs(const std::vector &Lons, const std::vector &Lats, const std::vector &Alts, - bool areLocsGeo=true, - bool areLocsIJK=true); + bool areLocsGeo = true, + bool areLocsIJK = true); /** * \brief Set the interpolation coefficients for the dipole grid diff --git a/include/inputs.h b/include/inputs.h index ceae2a5c..ce273f76 100644 --- a/include/inputs.h +++ b/include/inputs.h @@ -19,7 +19,7 @@ class Inputs { -public: + public: int iVerbose; int iVerboseProc; @@ -45,7 +45,7 @@ class Inputs { // - "cubesphere", sets grid.iGridShape_ = iCubesphere_ // - "dipole", sets grid.iGridShape_ = iDipole_ std::string shape; - + // Minimum altitude to simulate: precision_t alt_min; // Some grids allow the specification of the maximum altitude: @@ -61,7 +61,7 @@ class Inputs { bool IsUniformAlt; // Only needed for Mag Field grid: - // min_apex (not used) and LatStretch is used + // min_apex (not used) and LatStretch is used // as lat = min_lat + dlat where dlat = acos(cos(lat^stretch))^(1/stretch) precision_t min_apex, LatStretch, FieldLineStretch, max_blat; @@ -75,41 +75,41 @@ class Inputs { /********************************************************************** - \brief - \param + \brief + \param **/ Inputs() {} - + /********************************************************************** - \brief - \param + \brief + \param **/ Inputs(Times &time); /********************************************************************** - \brief - \param + \brief + \param **/ int read(Times &time); /********************************************************************** - \brief - \param + \brief + \param **/ bool read_inputs_json(Times &time); - + /********************************************************************** - \brief - \param + \brief + \param **/ bool set_verbose(json in); - + // -------------------------------------------------------------------- // get functions: // - These functions offer access to specific parts of the settings json. // - They call the general functions that check whether the key(s) exists. // - If the key does not exist, an error flag (in report) is set. - // - + // - // -------------------------------------------------------------------- // --------------------- @@ -121,13 +121,13 @@ class Inputs { \param none **/ int get_verbose(); - + /********************************************************************** \brief returns settings["Debug"]["iProc"] \param none **/ int get_verbose_proc(); - + /********************************************************************** \brief returns settings["Debug"]["dt"] \param none @@ -144,13 +144,13 @@ class Inputs { \param none **/ precision_t get_n_outputs(); - + /********************************************************************** \brief returns settings["Outputs"]["dt"][iOutput] \param iOutput int specifying which output file type to report on **/ precision_t get_dt_output(int iOutput); - + /********************************************************************** \brief returns settings["Outputs"]["type"][iOutput] \param iOutput int specifying which output file type to report on @@ -162,7 +162,7 @@ class Inputs { \param none **/ precision_t get_dt_euv(); - + /********************************************************************** \brief returns settings["Euv"]["IncludePhotoElectrons"] \param none @@ -175,263 +175,263 @@ class Inputs { \param none **/ std::string get_diffuse_auroral_model(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ std::string get_potential_model(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ std::string get_electrodynamics_dir(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ std::string get_electrodynamics_file(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ std::string get_electrodynamics_north_file(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ std::string get_electrodynamics_south_file(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ precision_t get_euv_heating_eff_neutrals(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ std::string get_euv_model(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ std::string get_euv_file(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ bool get_euv_douse(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ std::string get_aurora_file(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ std::string get_chemistry_file(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ std::string get_indices_lookup_file(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ std::vector get_omniweb_files(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ int get_number_of_omniweb_files(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ std::string get_f107_file(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ std::string get_planet(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ std::string get_planetary_file(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ std::string get_planet_species_file(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ std::string get_collision_file(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ bool get_do_calc_bulk_ion_temp(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ precision_t get_eddy_coef(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ precision_t get_eddy_bottom(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ precision_t get_eddy_top(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ bool get_use_eddy_momentum(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ bool get_use_eddy_energy(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ std::string get_bfield_type(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ bool get_do_restart(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ std::string get_restartout_dir(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ std::string get_restartin_dir(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ precision_t get_dt_write_restarts(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ int get_original_seed(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ int get_updated_seed(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ void set_seed(int seed); - + /********************************************************************** \brief returns settings[" - \param + \param **/ bool write_restart(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ - json get_perturb_values(); - + json get_perturb_values(); + /********************************************************************** \brief returns settings[" - \param + \param **/ bool get_do_lat_dependent_radius(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ bool get_do_J2(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ bool get_check_for_nans(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ bool get_nan_test(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ std::string get_nan_test_variable(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ bool get_is_cubesphere(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ bool get_NO_cooling(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ bool get_O_cooling(); @@ -474,53 +474,54 @@ class Inputs { /********************************************************************** \brief returns settings[" - \param + \param **/ bool get_use_centripetal(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ bool get_use_coriolis(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ bool get_cent_acc(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ std::string get_student_name(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ bool get_is_student(); - - + + /********************************************************************** \brief returns settings[" - \param + \param **/ json get_initial_condition_types(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ json get_boundary_condition_types(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ std::string get_advection_neutrals_vertical(); + std::string get_advection_neutrals_horizontal(); bool get_advection_neutrals_bulkwinds(); bool get_advection_neutrals_implicitfriction(); @@ -528,81 +529,81 @@ class Inputs { /********************************************************************** \brief returns settings[" - \param + \param **/ int get_nLons(std::string gridtype); - + /********************************************************************** \brief returns settings[" - \param + \param **/ int get_nLats(std::string gridtype); - + /********************************************************************** \brief returns settings[" - \param + \param **/ int get_nAlts(std::string gridtype); /********************************************************************** \brief returns settings[gridtype, "shape"] - \param + \param **/ std::string get_grid_shape(std::string gridtype); - + /********************************************************************** \brief returns settings[" - \param + \param **/ int get_nMembers(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ std::string get_logfile(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ std::string get_logfile(int64_t iLog); - + /********************************************************************** \brief returns settings[" - \param + \param **/ std::vector get_species_vector(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ bool get_logfile_append(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ precision_t get_logfile_dt(); // Satellites - + /********************************************************************** \brief returns settings[" - \param + \param **/ std::vector get_satellite_files(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ std::vector get_satellite_names(); - + /********************************************************************** \brief returns settings[" - \param + \param **/ std::vector get_satellite_dts(); @@ -613,148 +614,148 @@ class Inputs { json get_tests(); // General get_setting functions with error checks: - + /********************************************************************** - \brief - \param + \brief + \param **/ std::string get_setting_str(std::string key1); - + /********************************************************************** - \brief - \param + \brief + \param **/ std::string get_setting_str(std::string key1, std::string key2); - + /********************************************************************** - \brief - \param + \brief + \param **/ std::string get_setting_str(std::string key1, std::string key2, std::string key3); - + /********************************************************************** - \brief - \param + \brief + \param **/ json get_setting_json(std::string key1); - + /********************************************************************** - \brief - \param + \brief + \param **/ json get_setting_json(std::string key1, std::string key2); - + /********************************************************************** - \brief - \param + \brief + \param **/ bool get_setting_bool(std::string key1); - + /********************************************************************** - \brief - \param + \brief + \param **/ bool get_setting_bool(std::string key1, std::string key2); - + /********************************************************************** - \brief - \param + \brief + \param **/ bool get_setting_bool(std::string key1, std::string key2, std::string key3); - + /********************************************************************** - \brief - \param + \brief + \param **/ precision_t get_setting_float(std::string key1); - + /********************************************************************** - \brief - \param + \brief + \param **/ precision_t get_setting_float(std::string key1, std::string key2); - + /********************************************************************** - \brief - \param + \brief + \param **/ int64_t get_setting_int(std::string key1); - + /********************************************************************** - \brief - \param + \brief + \param **/ int64_t get_setting_int(std::string key1, std::string key2); - + /********************************************************************** - \brief - \param + \brief + \param **/ std::vector get_setting_intarr(std::string key1); - + /********************************************************************** - \brief - \param + \brief + \param **/ std::vector get_setting_intarr(std::string key1, std::string key2); /********************************************************************** - \brief - \param + \brief + \param **/ std::vector get_setting_timearr(std::string key1); // Check settings functions: - + /********************************************************************** - \brief - \param + \brief + \param **/ bool check_settings(std::string key1); - + /********************************************************************** - \brief - \param + \brief + \param **/ bool check_settings(std::string key1, std::string key2); - + /********************************************************************** - \brief - \param + \brief + \param **/ std::string check_settings_str(std::string key1); - + /********************************************************************** - \brief - \param + \brief + \param **/ std::string check_settings_str(std::string key1, std::string key2); - + /********************************************************************** - \brief - \param + \brief + \param **/ precision_t check_settings_pt(std::string key1, std::string key2); - + /********************************************************************** \brief Check to see if internal state of class is ok **/ bool is_ok(); - -private: + + private: // This is the main variable that contains all of the settings in Aether: json settings; - + // These are a bunch of misc strings that should go away: std::string euv_file = "UA/inputs/euv.csv"; std::string aurora_file = "UA/inputs/aurora_earth.csv"; @@ -784,7 +785,7 @@ class Inputs { std::string restart_in_directory = "UA/restartIn"; bool DoRestart; - + precision_t dt_euv; precision_t dt_report; @@ -793,7 +794,7 @@ class Inputs { int nAltsGeo; int updated_seed; - + /// An internal variable to hold the state of the class bool isOk; diff --git a/include/ions.h b/include/ions.h index 05ab7bbe..b43a5f89 100644 --- a/include/ions.h +++ b/include/ions.h @@ -11,13 +11,13 @@ * \class Ions * * \brief Defines the ion states - * + * * The Ion class defines the ion states as well as a bunch - * of derived states and source/loss terms. + * of derived states and source/loss terms. * * \author Aaron Ridley * - * \date 2021/03/28 + * \date 2021/03/28 * **************************************************************/ @@ -29,7 +29,7 @@ class Ions { // species of ion. We will then have a vector of these species. int64_t nSpecies = 8; - + struct species_chars { /// Name of the species @@ -68,7 +68,7 @@ class Ions { /// Ion - Electron collision frequencies: std::vector nu_ion_electron; - + // Sources and Losses: /// Number density of species (/m3) @@ -187,12 +187,25 @@ class Ions { /// Average energy of diffuse electron aurora (keV, tbc): arma_mat avee; + // Some variables that we are going to use in calc_ion_v: + std::vector gravity_vcgc; + std::vector wind_acc; + std::vector total_acc; + std::vector efield_acc; + std::vector a_par; + std::vector a_perp; + std::vector a_x_b; + std::vector grad_Pi_plus_Pe; + arma_cube rho, nuin, nuin_sum, Nie, sum_rho; + arma_cube top, bottom; + + /// Number of species to advect: int nSpeciesAdvect; - + /// IDs of species to advect: std::vector species_to_advect; - + // names and units const std::string density_name = "Neutral Bulk Density"; const std::string density_unit = "/m3"; @@ -207,7 +220,7 @@ class Ions { const std::string potential_name = "Potential"; const std::string potential_unit = "Volts"; - + // -------------------------------------------------------------------- // Functions: @@ -225,11 +238,11 @@ class Ions { species_chars create_species(Grid grid); /********************************************************************** - \brief + \brief \param planet contains information about the species to simulate **/ int read_planet_file(Planets planet); - + /********************************************************************** \brief Initialize the ion temperature (to the neutral temperature) \param neutrals the neutral class to grab the temperature from @@ -325,16 +338,16 @@ class Ions { \param dt the delta-t for the current time **/ void calc_ion_drift(Neutrals neutrals, - Grid grid, - precision_t dt); - + Grid grid, + precision_t dt); + /********************************************************************** \brief Calculate the ion + electron pressure gradient \param iIon which ion to act upon \param grid this is the grid to solve the equation on **/ std::vector calc_ion_electron_pressure_gradient(int64_t iIon, - Grid grid); + Grid grid); /********************************************************************** \brief Calculates the ion temperature(s) on the given grid @@ -356,7 +369,7 @@ class Ions { /// @brief Calculate epsilon /// @details intermediate variable used in photoelectron & ionization heating /// From (Smithro & Solomon, 2008). - /// @param neutrals + /// @param neutrals /// @return epsilon **/ arma_cube calc_epsilon(Neutrals &neutrals); @@ -366,17 +379,17 @@ class Ions { \details Based on (Swartz & Nisbet, 1972) & (Smithro & Solomon, 2008) Uses equations 9-12 from (Zhu & Ridley, 2016) https://doi.org/10.1016/j.jastp.2016.01.005 - \param epsilon - \return Qphe + \param epsilon + \return Qphe **/ arma_cube calc_photoelectron_heating(arma_cube epsilon); /********************************************************************** \brief Calculates auroral heating - \details NOTE: in GITM this is solved separately for ion precipitation & auroral + \details NOTE: in GITM this is solved separately for ion precipitation & auroral ionization. In Aether these are both in ions.species[iIon].ionization_scgc... - \param epsilon - \return Qaurora + \param epsilon + \return Qaurora **/ arma_cube calc_ionization_heating(arma_cube epsilon); @@ -394,26 +407,28 @@ class Ions { /********************************************************************** \brief Calculates electron-neutral elastic collisional heating \details From Schunk and Nagy 2009 - \param neutrals + \param neutrals \return vector **/ - std::vector calc_electron_neutral_elastic_collisions(Neutrals &neutrals); + std::vector calc_electron_neutral_elastic_collisions( + Neutrals &neutrals); /********************************************************************** \brief Calculates the electron-neutral inelastic collisional heating \details From Schunk and Nagy 2009 pages 277, 282. This includes N2, O2 rotation, fine structure, O(1D) exitation & vibration, N2 vibration. See equation 15 from (Zhu, Ridley, Deng, 2016) https://doi.org/10.1016/j.jastp.2016.01.005 - \param neutrals + \param neutrals \return vector **/ - std::vector calc_electron_neutral_inelastic_collisions(Neutrals &neutrals); + std::vector calc_electron_neutral_inelastic_collisions( + Neutrals &neutrals); /********************************************************************** \brief Calculate the thermoelectric current (same at all altitudes) \details Use eq. 6 of https://doi.org/10.1016/j.jastp.2016.01.005 - Since we do not know e- parallel velocity, the dipole needs to do it this way too. - \param grid + \param grid \return arma_mat JParaAlt **/ arma_mat calc_thermoelectric_current(Grid &grid); @@ -422,7 +437,7 @@ class Ions { \brief Check all of the variables for nonfinites, such as nans \param none **/ - bool check_for_nonfinites(); + bool check_for_nonfinites(std::string location); /********************************************************************** \brief Run through a test of an arma_cube to see if it contains nans @@ -445,7 +460,7 @@ class Ions { bool exchange_old(Grid &grid); /********************************************************************** - \brief Vertical advection solver - Rusanov + \brief Vertical advection solver - Rusanov \param grid The grid to define the neutrals on \param time contains information about the current time **/ diff --git a/include/neutrals.h b/include/neutrals.h index d7ac5321..34f67e34 100644 --- a/include/neutrals.h +++ b/include/neutrals.h @@ -11,7 +11,7 @@ * \class Neutrals * * \brief Defines the neutral states - * + * * The Neutrals class defines the neutrals states as well as a bunch * of derived states and source/loss terms. The initial temperature * structure as well as the lower boundary densities can be set @@ -19,7 +19,7 @@ * * \author Aaron Ridley * - * \date 2021/03/28 + * \date 2021/03/28 * **************************************************************/ @@ -52,7 +52,7 @@ class Neutrals { /// Number density of species (/m3) arma_cube density_scgc; arma_cube newDensity_scgc; - + /// Velocity of each species (m/s). For all below: /// Index 0 = longitudinal component of velocity /// Index 1 = latitudinal @@ -65,14 +65,14 @@ class Neutrals { /// Coefficient for the friction term (sum of friction coefs with others) arma_cube neutral_friction_coef; - + /// Acceleration of each species based on Eddy contribution. /// Only in vertical direction. arma_cube acc_eddy; - + /// Acceleration of each species due to ion drag. std::vector acc_ion_drag; - + /// concentration (density of species / total density) arma_cube concentration_scgc; // mass concentration (mass * density of species / rho) @@ -104,7 +104,7 @@ class Neutrals { std::vector iEuvPeiId_; /// Which ion species results from the ionization? std::vector iEuvPeiSpecies_; - + int nAuroraIonSpecies; std::vector iAuroraIonSpecies_; float Aurora_Coef; @@ -127,7 +127,7 @@ class Neutrals { /// Chemistry source rate (/m3/s) arma_cube sources_scgc; - + /// Chemistry loss rate (/m3/s) arma_cube losses_scgc; @@ -145,7 +145,7 @@ class Neutrals { /// sound speed + abs(bulk velocity (m/s)) std::vector cMax_vcgc; - + /// bunk temperature (K) arma_cube temperature_scgc; arma_cube newTemperature_scgc; @@ -177,7 +177,7 @@ class Neutrals { /// Viscosity arma_cube viscosity_scgc; - /// O cooling + /// O cooling arma_cube O_cool_scgc; /// NO cooling @@ -225,10 +225,10 @@ class Neutrals { std::vector initial_altitudes; std::vector initial_temperatures; int64_t nInitial_temps = 0; - + /// Number of species to advect: int nSpeciesAdvect; - + /// IDs of species to advect: std::vector species_to_advect; @@ -257,9 +257,9 @@ class Neutrals { \param indices used to help set initial conditions **/ Neutrals(Grid grid, - Planets planet, - Times time, - Indices indices); + Planets planet, + Times time, + Indices indices); /********************************************************************** \brief Creates the variables within the species_chars structure @@ -270,7 +270,7 @@ class Neutrals { /********************************************************************** \brief Read in the planet-specific file - This file specifies the species to model, their masses, + This file specifies the species to model, their masses, diffusion coefficients and all of the other things needed for specifying the neutrals. @@ -285,8 +285,8 @@ class Neutrals { \param indices used to help set initial conditions **/ bool initial_conditions(Grid grid, - Times time, - Indices indices); + Times time, + Indices indices); /********************************************************************** \brief temporary function to set neutral densities with in the model @@ -301,14 +301,14 @@ class Neutrals { \param grid The grid to define the neutrals on **/ void fill_with_hydrostatic(int64_t iStart, - int64_t iEnd, - Grid grid); + int64_t iEnd, + Grid grid); void fill_with_hydrostatic(int64_t iSpecies, - int64_t iStart, - int64_t iEnd, - Grid grid); - + int64_t iStart, + int64_t iEnd, + Grid grid); + /********************************************************************** \brief Limit the density to a floor and a ceiling **/ @@ -324,38 +324,38 @@ class Neutrals { \param grid The grid to define the neutrals on **/ void calc_scale_height(Grid grid); - + /********************************************************************** \brief Calculate the viscosity coefficient **/ void calc_viscosity(); - + /********************************************************************** \brief Calculate the eddy diffusion coefficient in valid pressure **/ void calc_kappa_eddy(); - + /********************************************************************** \brief Calculate the concentration for each species (species ndensity / total ndensity) **/ void calc_concentration(); /********************************************************************** - \brief Calculate the density of each species from the mass concentration + \brief Calculate the density of each species from the mass concentration for each species and rho (ndensity = con * rho / mass) **/ void calc_density_from_mass_concentration(); - + /********************************************************************** \brief Calculate the bulk mean major mass **/ void calc_mean_major_mass(); - + /********************************************************************** \brief Calculate the mean pressure **/ void calc_pressure(); - + /********************************************************************** \brief Calculate bulk velocity **/ @@ -421,8 +421,8 @@ class Neutrals { \param indices used to help set initial conditions **/ bool set_bcs(Grid grid, - Times time, - Indices indices); + Times time, + Indices indices); /********************************************************************** \brief Set lower boundary conditions for the neutrals @@ -431,8 +431,8 @@ class Neutrals { \param indices used to help set initial conditions **/ bool set_lower_bcs(Grid grid, - Times time, - Indices indices); + Times time, + Indices indices); /********************************************************************** \brief Set upper boundary conditions for the neutrals @@ -448,7 +448,7 @@ class Neutrals { \param grid The grid to define the neutrals on **/ bool set_horizontal_bcs(int64_t iDir, Grid grid); - + /********************************************************************** \brief Get the species ID number (int) given the species name (string) \param name string holding the species name (e.g., "O+") @@ -470,7 +470,7 @@ class Neutrals { \param dir directory to write restart files \param DoRead read the restart files if true, write if false **/ - bool restart_file(std::string dir, std::string cGridtype, bool DoRead); + bool restart_file(std::string dir, std::string cGridtype, bool DoRead); /********************************************************************** \brief Exchange messages between processors @@ -482,9 +482,9 @@ class Neutrals { /********************************************************************** \brief add eddy contributions to vertical acceleration \param grid The grid to define the neutrals on - **/ + **/ void vertical_momentum_eddy(Grid &grid); - + /********************************************************************** \brief Exchange one face for the NEUTRALS @@ -503,29 +503,29 @@ class Neutrals { **/ bool exchange_one_face(int iReceiver, int iSender, - precision_t *buffer, - int64_t iTotalSize, - int nG, int iDir); + precision_t *buffer, + int64_t iTotalSize, + int nG, int iDir); bool pack_one_face(int iReceiver, - precision_t *buffer, - int nG, int iDir, - bool IsPole); + precision_t *buffer, + int nG, int iDir, + bool IsPole); bool unpack_one_face(int iSender, - precision_t *buffer, - int nG, int iDir, - bool DoReverseX, - bool DoReverseY, - bool XbecomesY); + precision_t *buffer, + int nG, int iDir, + bool DoReverseX, + bool DoReverseY, + bool XbecomesY); /********************************************************************** - \brief Vertical advection solver - Rusanov + \brief Vertical advection solver - Rusanov \param grid The grid to define the neutrals on \param time contains information about the current time **/ void solver_vertical_rusanov(Grid grid, - Times time); - + Times time); + /********************************************************************** \brief Call the correct vertical advection scheme \param grid The grid to define the neutrals on @@ -542,18 +542,144 @@ class Neutrals { \param vels updated velocity, which acts as a source term for the implicit solve **/ arma_vec calc_friction_one_cell(int64_t iLong, int64_t iLat, int64_t iAlt, - precision_t dt, arma_vec &vels); + precision_t dt, arma_vec &vels); /********************************************************************** \brief Calculate the neutral friction in all cells (calls one_cell above) \param dt time step **/ - void calc_neutral_friction_implicit(precision_t dt); + void calc_neutral_friction_implicit(precision_t dt); /********************************************************************** \brief Calculate the neutral friction coefficients for semi-implicit solver **/ - void calc_neutral_friction_coefs(); + void calc_neutral_friction_coefs(); + + /********************************************************************** + \brief Residuals for **fluid motion** horizontally with Rusanov + \brief It actually updates the weighted residuals (-1/Area*R) for efficiency + + \param grid + \param time + \param states + **/ + std::vector residual_horizontal_rusanov(std::vector& states, + Grid& grid, Times& time, int64_t iAlt); + + /********************************************************************** + \brief Solves for **fluid motion** horizontally with RK4 + + \param grid + \param time + \param report + **/ + void solver_horizontal_RK4(Grid& grid, Times& time); + + /********************************************************************** + \brief Solves for **fluid motion** horizontally with RK1 + + \param grid + \param time + \param report + **/ + void solver_horizontal_RK1(Grid& grid, Times& time); + void solver_horizontal_RK1_rochi(Grid& grid, Times& time); + + /********************************************************************** + \brief Call the correct horizontal advection scheme with CE eqn + \param grid The grid to define the neutrals on + \param time contains information about the current time + **/ + bool advect_horizontal(Grid& grid, Times& time); + + /********************************************************************** + \brief Solves for fluid motion (pure advect) horizontally with Rusanov + + \param grid + \param time + **/ + void solver_horizontal_rusanov_advection(Grid& grid, Times& time); + void advect_sphere(Grid &grid, Times &time); + + /********************************************************************** + \brief Solves for fluid motion (pure advect) horizontally with RK1 + + \param grid + \param time + **/ + void solver_horizontal_RK1_advection(Grid& grid, Times& time); + + /********************************************************************** + \brief Solves for fluid motion (pure advect) horizontally with RK2 + + \param grid + \param time + **/ + void solver_horizontal_RK2_advection(Grid& grid, Times& time); + + /********************************************************************** + \brief Solves for fluid motion (pure advect) horizontally with RK4 + + \param grid + \param time + **/ + void solver_horizontal_RK4_advection(Grid& grid, Times& time); + + /********************************************************************** + \brief Residuals for fluid motion (pure advect) horizontally with HLLE + \brief It actually updates the weighted residuals (-1/Area*R) for efficiency + + \param grid + \param time + \param states + **/ + std::vector residual_horizontal_hlle_advection( + std::vector& states, Grid& grid, Times& time); + + /********************************************************************** + \brief Residuals for fluid motion (pure advect) horizontally with Rusanov + \brief It actually updates the weighted residuals (-1/Area*R) for efficiency + + \param grid + \param time + \param states + **/ + std::vector residual_horizontal_rusanov_advection( + std::vector& states, Grid& grid, Times& time); + + /********************************************************************** + \brief Call the horizontal advection scheme with only advection + \param grid The grid to define the neutrals on + \param time contains information about the current time + **/ + bool advect_horizontal_advection(Grid& grid, Times& time); + + /********************************************************************** + \brief Setup initial condition for the cosine bell test + \brief For advection test + \param grid The grid to define the neutrals on + \param time contains information about the current time + \param indices used to help set initial conditions + \param planet planet data for extracting the radius + **/ + bool cosine_bell_ic(Grid grid, + Times time, + Indices indices, + Planets planet); + + /********************************************************************** + \brief Setup initial condition for the blob test + \brief For Actual Cubesphere fluid solver + \param grid The grid to define the neutrals on + \param time contains information about the current time + \param indices used to help set initial conditions + \param planet planet data for extracting the radius + **/ + bool blob_ic(Grid grid, + Times time, + Indices indices, + Planets planet); + }; #endif // INCLUDE_NEUTRALS_H_ diff --git a/include/solvers.h b/include/solvers.h index 531ea23c..d7d2412a 100644 --- a/include/solvers.h +++ b/include/solvers.h @@ -21,32 +21,68 @@ struct projection_struct { arma_mat grad_edge_DU; }; -arma_vec limiter_mc(arma_vec &left, arma_vec &right, int64_t nPts, int64_t nGCs); +arma_vec limiter_mc(arma_vec &left, arma_vec &right, int64_t nPts, + int64_t nGCs); arma_vec calc_grad_1d(arma_vec &values, - arma_vec &x, - int64_t nPts, - int64_t nGCs); + arma_vec &x, + int64_t nPts, + int64_t nGCs); arma_mat calc_grad(arma_mat values, arma_mat x, int64_t nGCs, bool DoX); + +projection_struct project_to_edges(arma_mat &values, + arma_mat &x_centers, arma_mat &x_edges, + arma_mat &y_centers, arma_mat &y_edges, + int64_t nGCs); + +namespace Cubesphere_tools { +/* +struct projection_struct { + arma_mat gradLR; + arma_mat gradDU; + arma_mat R; + arma_mat L; + arma_mat U; + arma_mat D; +}; +*/ +arma_vec limiter_mc(arma_vec &left, arma_vec &right, int64_t nPts, + int64_t nGCs); +void print(arma_vec values); +arma_vec calc_grad_1d(arma_vec &values, arma_vec &x, int64_t nPts, + int64_t nGCs); +arma_mat calc_grad(arma_mat values, arma_mat x, int64_t nGCs, bool DoX); + +arma_mat project_from_left(arma_mat values, arma_mat gradients, + arma_mat x_centers, arma_mat x_edges, int64_t nGCs); +arma_mat project_from_right(arma_mat values, arma_mat gradients, + arma_mat x_centers, arma_mat x_edges, int64_t nGCs); +arma_vec limiter_value(arma_vec projected, arma_vec values, int64_t nPts, + int64_t nGCs); +//projection_struct project_to_edges(arma_mat &values, arma_mat &x_centers, +// arma_mat &x_edges, arma_mat &y_centers, arma_mat &y_edges, int64_t nGCs); +} + + void advect(Grid &grid, Times &time, Neutrals &neutrals); arma_vec solver_conduction( - arma_vec value, - arma_vec lambda, - arma_vec front, - arma_vec source, - arma_vec dx, - precision_t dt, - int64_t nGCs, - bool return_diff = false, - arma_vec source2 = arma_vec()); + arma_vec value, + arma_vec lambda, + arma_vec front, + arma_vec source, + arma_vec dx, + precision_t dt, + int64_t nGCs, + bool return_diff = false, + arma_vec source2 = arma_vec()); arma_cube solver_chemistry(arma_cube density, - arma_cube source, - arma_cube loss, - precision_t dt); + arma_cube source, + arma_cube loss, + precision_t dt); std::vector coriolis(std::vector velocity, precision_t rotation_rate, @@ -56,32 +92,32 @@ std::vector coriolis(std::vector velocity, /// or an interpolated value should be used. const int iPrevious_ = 1; const int iNext_ = 2; -const int iClosest_ = 3; +const int iClosest_ = 3; const int iInterp_ = 4; double interpolate_1d(double outX, - std::vector inXs, - std::vector inValues); + std::vector inXs, + std::vector inValues); double interpolate_1d_get_index_doubles(double intime, - std::vector times); + std::vector times); // Overloading the interpolation function: double interpolate_1d_w_index(std::vector values, - double interpolation_index, - int interpolation_type); + double interpolation_index, + int interpolation_type); double interpolate_1d_w_index(std::vector values, - double interpolation_index, - int interpolation_type); + double interpolation_index, + int interpolation_type); double interpolate_1d_w_index(std::vector values, - float interpolation_index, - int interpolation_type); + float interpolation_index, + int interpolation_type); double interpolate_1d_w_index(arma_vec values, - double interpolation_index, - int interpolation_type); + double interpolation_index, + int interpolation_type); fmat interpolate_1d_w_index(std::vector values, - double interpolation_index, - int interpolation_type); + double interpolation_index, + int interpolation_type); arma_cube calc_gradient_lon(arma_cube value, Grid grid); arma_cube calc_gradient_lat(arma_cube value, Grid grid); @@ -114,18 +150,18 @@ precision_t interpolate_unit_cube(const arma_cube &data, const precision_t zRatio); precision_t limiter_mc(precision_t dUp, - precision_t dDown, - precision_t beta); - - - /********************************************************************** - \brief Calculate dt (cell size / cMax) in each direction, and take min - \param dt returns the neutral time-step - \param grid The grid to define the neutrals on - **/ - precision_t calc_dt(Grid grid, std::vector cMax_vcgc); - precision_t calc_dt_sphere(Grid grid, std::vector cMax_vcgc); - precision_t calc_dt_cubesphere(Grid grid, std::vector cMax_vcgc); - precision_t calc_dt_vertical(Grid grid, std::vector cMax_vcgc); + precision_t dDown, + precision_t beta); + + +/********************************************************************** + \brief Calculate dt (cell size / cMax) in each direction, and take min + \param dt returns the neutral time-step + \param grid The grid to define the neutrals on + **/ +precision_t calc_dt(Grid grid, std::vector cMax_vcgc); +precision_t calc_dt_sphere(Grid grid, std::vector cMax_vcgc); +precision_t calc_dt_cubesphere(Grid grid, std::vector cMax_vcgc); +precision_t calc_dt_vertical(Grid grid, std::vector cMax_vcgc); #endif // INCLUDE_SOLVERS_H_ diff --git a/share/run/UA/inputs/defaults.json b/share/run/UA/inputs/defaults.json index fd7bfc84..91ce86a8 100644 --- a/share/run/UA/inputs/defaults.json +++ b/share/run/UA/inputs/defaults.json @@ -24,7 +24,7 @@ "Advection" : { "Neutrals" : { "Vertical" : "rusanov", - "Horizontal" : "default", + "Horizontal" : "fv", "useBulkWinds" : true, "useImplicitFriction" : true}, "Ions" : { diff --git a/share/run/aether.json b/share/run/aether.json index cc3e330d..b192b5b5 100644 --- a/share/run/aether.json +++ b/share/run/aether.json @@ -10,21 +10,21 @@ "dt" : 10.0, "check_for_nans" : true}, - "EndTime" : [2011, 3, 20, 0, 10, 0], + "EndTime" : [2011, 3, 20, 0, 1, 0], "neuGrid" : { "Shape": "sphere4", - "nLonsPerBlock" : 24, - "nLatsPerBlock" : 22, + "nLonsPerBlock" : 20, + "nLatsPerBlock" : 18, "nAlts" : 40, - "dAltScale" : 0.25, + "dAltScale" : 0.3, "IsUniformAlt" : false}, "ionGrid": { "Shape": "dipole4", - "nLonsPerBlock": 36, + "nLonsPerBlock": 24, "nLatsPerBlock": 18, - "nAlts": 100, + "nAlts": 50, "LatRange": [10, 80], "AltRange": [80.0, 1000], "LonRange": [0.0, 360.0]}, diff --git a/src/advance.cpp b/src/advance.cpp index 134a0365..96230431 100644 --- a/src/advance.cpp +++ b/src/advance.cpp @@ -74,7 +74,6 @@ bool advance(Planets &planet, didWork = neutralsMag.check_for_nonfinites("Ion Grid: After extras"); - ions.fill_electrons(); ions.calc_sound_speed(); ions.calc_cMax(); @@ -111,7 +110,6 @@ bool advance(Planets &planet, didWork = neutralsMag.check_for_nonfinites("Ion Grid: set bcs"); - // advect in the 3rd dimension (vertical), but only if we have it: if (gGrid.get_nAlts(false) > 1) { neutrals.advect_vertical(gGrid, time); @@ -119,7 +117,11 @@ bool advance(Planets &planet, if (didWork & input.get_check_for_nans()) didWork = neutrals.check_for_nonfinites("After Vertical Neutral Advection"); - ions.advect_vertical(gGrid, time); + // ajr - ions.advect_vertical(gGrid, time); + + if (didWork & input.get_check_for_nans()) + didWork = ions.check_for_nonfinites("After Vertical Ion Advection"); + } // advect in the 1st and 2nd dimensions (horizontal), but only if @@ -127,13 +129,19 @@ bool advance(Planets &planet, if (gGrid.get_HasXdim() || gGrid.get_HasYdim()) { neutrals.exchange_old(gGrid); ions.exchange_old(gGrid); + neutrals.advect_horizontal(gGrid, time); ionsMag.exchange_old(mGrid); - advect(gGrid, time, neutrals); + //advect(gGrid, time, neutrals); } - if (didWork & input.get_check_for_nans()) { + if (input.get_check_for_nans()) { didWork = neutrals.check_for_nonfinites("Geo Grid: After Horizontal Advection"); didWork = neutralsMag.check_for_nonfinites("Ion Grid: After Horizontal Advection"); + + if (!didWork) { + report.exit(function); + return didWork; + } } // ------------------------------------ @@ -200,10 +208,10 @@ bool advance(Planets &planet, didWork = neutralsMag.check_for_nonfinites("Ion Grid: After Add Sources"); } - ions.calc_ion_temperature(neutrals, gGrid, time); - // ions.calc_electron_temperature(neutrals, gGrid, time); + //ions.calc_ion_temperature(neutrals, gGrid, time); + //ions.calc_electron_temperature(neutrals, gGrid, time); //ionsMag.calc_ion_temperature(neutralsMag, mGrid, time); - ionsMag.calc_electron_temperature(neutralsMag, mGrid, time); + //ionsMag.calc_electron_temperature(neutralsMag, mGrid, time); if (didWork & input.get_check_for_nans()) didWork = neutrals.check_for_nonfinites("After Vertical Advection"); diff --git a/src/calc_ion_drift.cpp b/src/calc_ion_drift.cpp index ec64aec6..0cc9821a 100644 --- a/src/calc_ion_drift.cpp +++ b/src/calc_ion_drift.cpp @@ -13,7 +13,8 @@ void Ions::calc_efield(Grid grid) { efield_vcgc = calc_gradient_vector(-1.0 * potential_scgc, grid); // Remove component along b-field (should be zero, anyways!) - arma_cube edotb = dot_product(efield_vcgc, grid.bfield_unit_vcgc); + arma_cube edotb; + edotb = dot_product(efield_vcgc, grid.bfield_unit_vcgc); for (int64_t iComp = 0; iComp < 3; iComp++) efield_vcgc[iComp] = @@ -25,9 +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); @@ -96,21 +99,10 @@ void Ions::calc_ion_drift(Neutrals neutrals, report.print(5, "going into calc_exb_drift"); calc_exb_drift(grid); - std::vector gravity_vcgc = make_cube_vector(nX, nY, nZ, 3); - std::vector wind_acc = make_cube_vector(nX, nY, nZ, 3); - std::vector total_acc = make_cube_vector(nX, nY, nZ, 3); - std::vector efield_acc = make_cube_vector(nX, nY, nZ, 3); + int64_t iIon, iNeutral, iDim; + int64_t iComp; - int64_t iIon, iNeutral, iComp; - - std::vector grad_Pi_plus_Pe; - arma_cube rho, nuin, nuin_sum, Nie, sum_rho; - arma_cube top, bottom; - - nuin_sum.set_size(nX, nY, nZ); nuin_sum.zeros(); - - sum_rho.set_size(nX, nY, nZ); sum_rho.zeros(); fill_electrons(); @@ -118,10 +110,6 @@ void Ions::calc_ion_drift(Neutrals neutrals, for (iComp = 0; iComp < 3; iComp++) velocity_vcgc[iComp].zeros(); - std::vector a_par = make_cube_vector(nX, nY, nZ, 3); - std::vector a_perp = make_cube_vector(nX, nY, nZ, 3); - std::vector a_x_b; - for (iIon = 0; iIon < nSpecies; iIon++) { for (iComp = 0; iComp < 3; iComp++) diff --git a/src/cubesphere_tools.cpp b/src/cubesphere_tools.cpp new file mode 100644 index 00000000..d94011bb --- /dev/null +++ b/src/cubesphere_tools.cpp @@ -0,0 +1,300 @@ +// Copyright 2024, the Aether Development Team (see doc/dev_team.md for members) +// Full license can be found in License.md + +// Initial version: F. Cheng, Feb 2024 + +#include "aether.h" + +arma_vec Cubesphere_tools::limiter_mc(arma_vec &left, + arma_vec &right, + int64_t nPts, + int64_t nGCs) { + + precision_t beta = 0.8; + + arma_vec s = left % right; + arma_vec combined = (left + right) * 0.5; + + left = left * beta; + right = right * beta; + arma_vec limited = left; + + for (int64_t i = 1; i < nPts + 2 * nGCs - 1; i++) { + if (s(i) < 0) { + // Sign < 0 means opposite signed left and right: + limited(i) = 0.0; + } else { + if (left(i) > 0 && right(i) > 0) { + if (right(i) < limited(i)) + limited(i) = right(i); + + if (combined(i) < limited(i)) + limited(i) = combined(i); + } else { + if (right(i) > limited(i)) + limited(i) = right(i); + + if (combined(i) > limited(i)) + limited(i) = combined(i); + } + } + } + + return limited; +} + +void Cubesphere_tools::print(arma_vec values) { + int64_t nP = values.n_elem; + + for (int64_t i = 0; i < nP; i++) + std::cout << values(i) << " "; + + std::cout << "\n"; +} + +// --------------------------------------------------------- +// calc gradients at centers +// - values and x defined at centers +// --------------------------------------------------------- + +arma_vec Cubesphere_tools::calc_grad_1d(arma_vec &values, + arma_vec &x, + int64_t nPts, + int64_t nGCs) { + + arma_vec gradients = values * 0.0; + arma_vec gradL = values * 0.0; + arma_vec gradR = values * 0.0; + + precision_t factor1 = 0.625; + precision_t factor2 = 0.0416667; + precision_t h; + + int64_t i; + arma_vec hv = values * 0.0; + + i = nGCs - 1; + h = 2.0 / (x(i + 1) - x(i)); + gradR(i) = h * (factor1 * (values(i + 1) - values(i)) - + factor2 * (values(i + 2) - values(i - 1))); + gradL(i) = (values(i) - values(i - 1)) / (x(i) - x(i - 1)); + + // This is attempting to vectorize the problem, but it seems to be slower? + // int64_t iS = nGCs; + // int64_t iE = nPts + nGCs - 1; + // hv.rows(iS, iE) = 2.0 / (x.rows(iS, iE) - x.rows(iS-1, iE-1)); + // gradL.rows(iS, iE) = hv.rows(iS,iE) % (factor1 * (values.rows(iS, iE) - + // values.rows(iS-1, iE-1)) - + // factor2 * (values.rows(iS+1, iE+1) - + // values.rows(iS-2, iE-2))); + // hv.rows(iS, iE) = 2.0 / (x.rows(iS+1, iE+1) - x.rows(iS, iE)); + // gradR.rows(iS, iE) = hv.rows(iS,iE) % (factor1 * (values.rows(iS+1, iE+1) - + // values.rows(iS, iE)) - + // factor2 * (values.rows(iS+2, iE+2) - + // values.rows(iS-1, iE-1))); + + for (i = nGCs; i < nPts + nGCs; i++) { + h = 2.0 / (x(i) - x(i - 1)); + gradL(i) = h * (factor1 * (values(i) - values(i - 1)) - + factor2 * (values(i + 1) - values(i - 2))); + h = 2.0 / (x(i + 1) - x(i)); + gradR(i) = h * (factor1 * (values(i + 1) - values(i)) - + factor2 * (values(i + 2) - values(i - 1))); + } + + i = nPts + nGCs; + h = 2.0 / (x(i) - x(i - 1)); + gradL(i) = h * (factor1 * (values(i) - values(i - 1)) - + factor2 * (values(i + 1) - values(i - 2))); + gradR(i) = (values(i + 1) - values(i)) / (x(i + 1) - x(i)); + + gradients = Cubesphere_tools::limiter_mc(gradL, gradR, nPts, nGCs); + + return gradients; +} + +// --------------------------------------------------------- +// calc gradients at centers for 2d matrices +// - values and x defined at centers +// --------------------------------------------------------- + +arma_mat Cubesphere_tools::calc_grad(arma_mat values, + arma_mat x, + int64_t nGCs, + bool DoX) { + + arma_mat v2d, x2d; + + if (DoX) { + v2d = values; + x2d = x; + } else { + v2d = values.t(); + x2d = x.t(); + } + + int64_t nX = v2d.n_rows; + int64_t nY = v2d.n_cols; + arma_mat grad2d = v2d * 0.0; + + int64_t nPts = nX - 2 * nGCs; + arma_vec values1d(nX); + arma_vec x1d(nX); + + for (int64_t j = 1; j < nY - 1; j++) { + values1d = v2d.col(j); + x1d = x2d.col(j); + grad2d.col(j) = calc_grad_1d(values1d, x1d, nPts, nGCs); + } + + arma_mat gradients; + + if (DoX) + gradients = grad2d; + else + gradients = grad2d.t(); + + return gradients; +} + +// --------------------------------------------------------- +// Project gradients + values to the right face, from the left +// returned values are on the i - 1/2 edges +// (between i-1 and i cell center) +// --------------------------------------------------------- + +arma_mat Cubesphere_tools::project_from_left(arma_mat values, + arma_mat gradients, + arma_mat x_centers, + arma_mat x_edges, + int64_t nGCs) { + + int64_t nX = values.n_rows; + int64_t nY = values.n_cols; + + // Define at edges: + arma_mat projected(nX + 1, nY); + projected.zeros(); + + // no gradient in the 0 or iEnd cells + for (int64_t j = 0; j < nY; j++) { + for (int64_t i = 1; i < nX - 1; i++) { + projected(i + 1, j) = values(i, j) + + gradients(i, j) * (x_edges(i + 1, j) - x_centers(i, j)); + } + + projected(1, j) = projected(2, j); + projected(0, j) = projected(1, j); + projected(nX, j) = projected(nX - 1, j); + } + + return projected; +} + +// --------------------------------------------------------- +// Project gradients + values to the left face, from the right +// returned values are on the i - 1 edges +// (between i-1 and i cell center) +// --------------------------------------------------------- + +arma_mat Cubesphere_tools::project_from_right(arma_mat values, + arma_mat gradients, + arma_mat x_centers, + arma_mat x_edges, + int64_t nGCs) { + int64_t nX = values.n_rows; + int64_t nY = values.n_cols; + + // Define at edges: + arma_mat projected(nX + 1, nY); + projected.zeros(); + + // no gradient in the 0 or iEnd cells + for (int64_t j = 0; j < nY; j++) { + for (int64_t i = 1; i < nX - 1; i++) { + projected(i, j) = values(i, j) + + gradients(i, j) * (x_edges(i, j) - x_centers(i, j)); + } + + projected(0, j) = projected(1, j); + projected(nX - 1, j) = projected(nX - 2, j); + projected(nX, j) = projected(nX - 1, j); + } + + return projected; +} + +// --------------------------------------------------------- +// Limiter on values +// projected is assumed to be on the edge between the +// i-1 and i cell (i-1/2) +// limited is returned at edges +// --------------------------------------------------------- + +arma_vec Cubesphere_tools::limiter_value(arma_vec projected, + arma_vec values, + int64_t nPts, + int64_t nGCs) { + + int64_t iStart = 0; + int64_t iEnd = nPts + 2 * nGCs; + + arma_vec limited = projected; + + precision_t mini, maxi; + + for (int64_t i = iStart + 1; i < iEnd - 1; i++) { + + mini = values(i - 1); + + if (values(i) < mini) + mini = values(i); + + maxi = values(i - 1); + + if (values(i) > maxi) + maxi = values(i); + + if (limited(i) < mini) + limited(i) = mini; + + if (limited(i) > maxi) + limited(i) = maxi; + } + + return limited; +} + +// // --------------------------------------------------------- +// // take gradients and project to all edges +// // --------------------------------------------------------- + +// projection_struct Cubesphere_tools::project_to_edges(arma_mat &values, +// arma_mat &x_centers, arma_mat &x_edges, +// arma_mat &y_centers, arma_mat &y_edges, +// int64_t nGCs) { + +// int64_t nX = values.n_rows; +// int64_t nY = values.n_cols; + +// projection_struct proj; + +// proj.gradLR = calc_grad(values, x_centers, nGCs, true); +// proj.gradDU = calc_grad(values.t(), y_centers.t(), nGCs, true).t(); + +// proj.R = project_from_left(values, proj.gradLR, +// x_centers, x_edges, nGCs); +// // Left side of edge from left +// proj.L = project_from_right(values, proj.gradLR, +// x_centers, x_edges, nGCs); +// // Up side of edge from down (left) +// proj.U = project_from_left(values.t(), proj.gradDU.t(), +// y_centers.t(), y_edges.t(), nGCs) +// .t(); +// // Down side of edge from up (right) +// proj.D = project_from_right(values.t(), proj.gradDU.t(), +// y_centers.t(), y_edges.t(), nGCs) +// .t(); + +// return proj; +// } diff --git a/src/exchange_messages.cpp b/src/exchange_messages.cpp index 895bb854..09f5ed1e 100644 --- a/src/exchange_messages.cpp +++ b/src/exchange_messages.cpp @@ -1109,8 +1109,6 @@ bool exchange_one_var(Grid &grid, if (!grid.isExchangeInitialized) { DidWork = exchange_sides_init(grid, nVarsToPass); grid.isExchangeInitialized = true; - std::cout << "initializing : " << grid.get_gridtype() << " " << nZ << " " << - iProc << "\n"; } var_scgc.set_size(nX, nY, nZ); diff --git a/src/grid_cubesphere.cpp b/src/grid_cubesphere.cpp index d5c6742d..4be8c25c 100644 --- a/src/grid_cubesphere.cpp +++ b/src/grid_cubesphere.cpp @@ -72,6 +72,404 @@ void Grid::create_cubesphere_connection(Quadtree quadtree) { return; } + +// ---------------------------------------------------------------------- +// This function takes the normalized coordinates and makes latitude +// and longitude arrays from them. It can do this for the corners or +// edges, depending on the offset. +// ---------------------------------------------------------------------- + +void Grid::init_cubesphere_grid(Quadtree quadtree, + arma_vec dr, + arma_vec du, + arma_vec ll, + precision_t left_off, + precision_t down_off, + cubesphere_chars &cubeX) { + + std::string function = "Grid::init_cubesphere_grid"; + static int iFunction = -1; + report.enter(function, iFunction); + + precision_t dnu, dxi, nu0, xi0; + + // The du, dr, and ll were meant to be used on the cube + // and not really on the equal-angle grid. So, we probably + // want to rethink these... + if (quadtree.iSide == 0) { + dnu = du[2]; + dxi = dr[0]; + nu0 = ll[0]; + xi0 = ll[2]; + } + + if (quadtree.iSide == 1) { + dnu = du[2]; + dxi = dr[1]; + nu0 = ll[1]; + xi0 = ll[2]; + } + + if (quadtree.iSide == 2) { + dnu = du[2]; + dxi = -dr[0]; + nu0 = -ll[0]; + xi0 = ll[2]; + } + + if (quadtree.iSide == 3) { + dnu = du[2]; + dxi = -dr[1]; + nu0 = -ll[1]; + xi0 = ll[2]; + } + + if (quadtree.iSide == 4) { + dnu = du[0]; + dxi = dr[1]; + nu0 = ll[1]; + xi0 = ll[0]; + } + + if (quadtree.iSide == 5) { + dnu = -du[0]; + dxi = dr[1]; + nu0 = ll[1]; + xi0 = -ll[0]; + } + + // Normalized from -1 to 1 -> -pi/4 to pi/4 + dnu = dnu * cPI / 4.0; + dxi = dxi * cPI / 4.0; + nu0 = nu0 * cPI / 4.0; + xi0 = xi0 * cPI / 4.0; + + cubeX.dnu = dnu; + cubeX.dxi = dxi; + + int64_t iDU, iLR; + precision_t iD, iL; + int64_t nXp = nX, nYp = nY; + + // If we are shifting the grid over and doing edges, we + // need to increase the number of points by 1 in that + // direction: + if (left_off < cSmall) + nXp++; + + if (down_off < cSmall) + nYp++; + + // These are convenient for the solver: + cubeX.nXt = nXp; + cubeX.nYt = nYp; + cubeX.nGCs = nGCs; + cubeX.iXfirst_ = nGCs; + cubeX.iXlast_ = nXp - nGCs; + cubeX.iYfirst_ = nGCs; + cubeX.iYlast_ = nYp - nGCs; + + // these are coordinates: + cubeX.lat.resize(nXp, nYp); + cubeX.lon.resize(nXp, nYp); + cubeX.nu.resize(nXp, nYp); + cubeX.xi.resize(nXp, nYp); + + cubeX.X.resize(nXp, nYp); + cubeX.Y.resize(nXp, nYp); + cubeX.Z.resize(nXp, nYp); + cubeX.C.resize(nXp, nYp); + cubeX.D.resize(nXp, nYp); + cubeX.d.resize(nXp, nYp); + + // These are dependent on radius, + // but that is not included at this time: + cubeX.dlx.resize(nXp, nYp, nZ); + cubeX.dln.resize(nXp, nYp, nZ); + cubeX.dS.resize(nXp, nYp, nZ); + cubeX.R.resize(nZ); + + // These are matricies for rotating vectors: + cubeX.Apn.resize(nXp, nYp); + cubeX.Apx.resize(nXp, nYp); + cubeX.Atn.resize(nXp, nYp); + cubeX.Atx.resize(nXp, nYp); + cubeX.Axt.resize(nXp, nYp); + cubeX.Axp.resize(nXp, nYp); + cubeX.Ant.resize(nXp, nYp); + cubeX.Anp.resize(nXp, nYp); + + // These are for computing normals to the cell edges (horizontal) + cubeX.nXiLon.resize(nXp, nYp); + cubeX.nXiLat.resize(nXp, nYp); + cubeX.nNuLon.resize(nXp, nYp); + cubeX.nNuLat.resize(nXp, nYp); + + precision_t det, dmo, latp, lonp; + + // Loop through each point and derive the coordinate + for (iDU = 0; iDU < nY; iDU++) { + for (iLR = 0; iLR < nX; iLR++) { + + // the offsets are so we can find cell centers, edges, and corners + iD = iDU - nGCs + down_off; + iL = iLR - nGCs + left_off; + + // Define local coordinates: + // Xi is LR (x), Nu is UD (y) + cubeX.nu(iLR, iDU) = (nu0 + dnu * iD); + cubeX.xi(iLR, iDU) = (xi0 + dxi * iL); + + cubeX.X(iLR, iDU) = tan(cubeX.xi(iLR, iDU)); + cubeX.Y(iLR, iDU) = tan(cubeX.nu(iLR, iDU)); + + // Transformation from 3D Cartesian to LatLong + // lonp = std::atan2(y_cart, x_cart) + cPI/2.0; + if (quadtree.iSide == 0) { + lonp = std::atan(cubeX.X(iLR, iDU)); + // Theta in Ronchi is from the north pole, so lat is 90 - theta + latp = std::atan(1.0 / cubeX.Y(iLR, iDU) / std::cos(lonp)); + } + + if (quadtree.iSide == 1) { + lonp = std::atan(-1.0 / cubeX.X(iLR, iDU)); + + if (lonp < 0) + lonp = cPI + lonp; + + // Theta in Ronchi is from the north pole, so lat is 90 - theta + latp = std::atan(1.0 / cubeX.Y(iLR, iDU) / std::sin(lonp)); + } + + if (quadtree.iSide == 2) { + lonp = std::atan(cubeX.X(iLR, iDU)) + cPI; + // Theta in Ronchi is from the north pole, so lat is 90 - theta + latp = std::atan(-1.0 / cubeX.Y(iLR, iDU) / std::cos(lonp)); + } + + if (quadtree.iSide == 3) { + lonp = std::atan(-1.0 / cubeX.X(iLR, iDU)); + + if (lonp > 0) + lonp = lonp + cPI; + else + lonp = 2 * cPI + lonp; + + // Theta in Ronchi is from the north pole, so lat is 90 - theta + latp = std::atan(-1.0 / cubeX.Y(iLR, iDU) / std::sin(lonp)); + } + + if (quadtree.iSide == 4) { + lonp = std::atan2(cubeX.X(iLR, iDU), cubeX.Y(iLR, iDU)); + latp = std::atan2(-cubeX.Y(iLR, iDU), cos(lonp) ); + } + + if (quadtree.iSide == 5) { + lonp = std::atan2(-cubeX.X(iLR, iDU), cubeX.Y(iLR, iDU)); + latp = -std::atan2(-cubeX.Y(iLR, iDU), cos(lonp) ); + } + + if (latp > 0) + latp = cPI / 2 - latp; + else + latp = -(cPI / 2 + latp); + + if (lonp > cTWOPI) + lonp = lonp - cTWOPI; + + if (lonp < 0.0) + lonp = lonp + cTWOPI; + + // Fill Computed coords + cubeX.lat(iLR, iDU) = latp; + cubeX.lon(iLR, iDU) = lonp; + + cubeX.d(iLR, iDU) = + 1 + + cubeX.X(iLR, iDU) * cubeX.X(iLR, iDU) + + cubeX.Y(iLR, iDU) * cubeX.Y(iLR, iDU); + + cubeX.C(iLR, iDU) = + sqrt(1 + cubeX.X(iLR, iDU) * cubeX.X(iLR, iDU)); + cubeX.D(iLR, iDU) = + sqrt(1 + cubeX.Y(iLR, iDU) * cubeX.Y(iLR, iDU)); + + if (quadtree.iSide < 4) { + cubeX.Axt(iLR, iDU) = 0.0; + cubeX.Axp(iLR, iDU) = + cubeX.C(iLR, iDU) * cubeX.D(iLR, iDU) / + sqrt(cubeX.d(iLR, iDU)); + cubeX.Ant(iLR, iDU) = -1.0; + cubeX.Anp(iLR, iDU) = + cubeX.X(iLR, iDU) * cubeX.Y(iLR, iDU) / + sqrt(cubeX.d(iLR, iDU)); + } else { + if (cubeX.d(iLR, iDU) < 1.0001) + cubeX.d(iLR, iDU) = 1.0001; + + dmo = 1.0 / std::sqrt(cubeX.d(iLR, iDU) - 1); + + if (quadtree.iSide == 4) { + cubeX.Axt(iLR, iDU) = + - dmo * cubeX.D(iLR, iDU) * cubeX.X(iLR, iDU); + cubeX.Axp(iLR, iDU) = + dmo * cubeX.D(iLR, iDU) * cubeX.Y(iLR, iDU) / + sqrt(cubeX.d(iLR, iDU)); + cubeX.Ant(iLR, iDU) = + - dmo * cubeX.C(iLR, iDU) * cubeX.Y(iLR, iDU); + cubeX.Anp(iLR, iDU) = + - dmo * cubeX.C(iLR, iDU) * cubeX.X(iLR, iDU) / + sqrt(cubeX.d(iLR, iDU)); + + } else { + // iFace == 5 + cubeX.Axt(iLR, iDU) = + dmo * cubeX.D(iLR, iDU) * cubeX.X(iLR, iDU); + cubeX.Axp(iLR, iDU) = + - dmo * cubeX.D(iLR, iDU) * + cubeX.Y(iLR, iDU) / + sqrt(cubeX.d(iLR, iDU)); + cubeX.Ant(iLR, iDU) = + dmo * cubeX.C(iLR, iDU) * cubeX.Y(iLR, iDU); + cubeX.Anp(iLR, iDU) = + dmo * cubeX.C(iLR, iDU) * + cubeX.X(iLR, iDU) / + sqrt(cubeX.d(iLR, iDU)); + } + } + + // Calculate inverse of matrix for calculating Ax and An from At and Ap: + det = 1.0 / (cubeX.Axt(iLR, iDU) * cubeX.Anp(iLR, iDU) - + cubeX.Axp(iLR, iDU) * cubeX.Ant(iLR, iDU)); + + cubeX.Atx(iLR, iDU) = det * cubeX.Anp(iLR, iDU); + cubeX.Atn(iLR, iDU) = - det * cubeX.Axp(iLR, iDU); + cubeX.Apx(iLR, iDU) = - det * cubeX.Ant(iLR, iDU); + cubeX.Apn(iLR, iDU) = det * cubeX.Axt(iLR, iDU); + + // These (dlx and dln) need to be multiplied by radius + cubeX.dlx(iLR, iDU, 0) = + cubeX.D(iLR, iDU) * dxi / + cubeX.d(iLR, iDU) / + (cos(cubeX.xi(iLR, iDU)) * cos(cubeX.xi(iLR, iDU))); + cubeX.dln(iLR, iDU, 0) = + cubeX.C(iLR, iDU) * dnu / + cubeX.d(iLR, iDU) / + (cos(cubeX.nu(iLR, iDU)) * cos(cubeX.nu(iLR, iDU))); + + // Need to multiply dS * radius ^ 2 + cubeX.dS(iLR, iDU, 0) = + dxi * dnu / + (sqrt(cubeX.d(iLR, iDU) * cubeX.d(iLR, iDU) * cubeX.d(iLR, iDU)) * + cos(cubeX.xi(iLR, iDU)) * cos(cubeX.xi(iLR, iDU)) * + cos(cubeX.nu(iLR, iDU)) * cos(cubeX.nu(iLR, iDU))); + + } + } + + // Calculate norms given the values above: + arma_mat e1Lat, e1Lon, e2Lat, e2Lon, m, one, zero; + m.resize(nXp, nYp); + one.resize(nXp, nYp); + one.fill(1.0); + zero.resize(nXp, nYp); + zero.fill(0.0); + + // define e1 as the LR (xi) direction: + e1Lat.resize(nXp, nYp); + e1Lon.resize(nXp, nYp); + convert_vector_xn_to_ll(one, zero, e1Lon, e1Lat, cubeX); + m = sqrt(e1Lon % e1Lon + e1Lat % e1Lat); + + // Rotate by 90 deg (CCW) to get the norm: + cubeX.nNuLon = -e1Lat / m; + cubeX.nNuLat = e1Lon / m; + + // define e2 as the DU (nu) direction: + e2Lat.resize(nXp, nYp); + e2Lon.resize(nXp, nYp); + convert_vector_xn_to_ll(zero, one, e2Lon, e2Lat, cubeX); + m = sqrt(e2Lon % e2Lon + e2Lat % e2Lat); + // Rotate by 90 deg (CW) to get the norm: + cubeX.nXiLon = e2Lat / m; + cubeX.nXiLat = -e2Lon / m; + + report.exit(function); + return; +} + +// --------------------------------------------------------- +// Convert vector from Alat, Alon to Axi, Anu +// -> Using equation (7) of Ronchi et al: +// --------------------------------------------------------- + +void Grid::convert_vector_xn_to_ll(arma_mat aXi, + arma_mat aNu, + arma_mat &aLon, + arma_mat &aLat, + cubesphere_chars grid) { + + // Ronchi defines aPhi = aLon, aTheta = -aLat + aLat = -(grid.Atx % aXi + grid.Atn % aNu); + aLon = grid.Apx % aXi + grid.Apn % aNu; + + return; +} + +// --------------------------------------------------------- +// Convert vector from Alat, Alon to Axi, Anu +// -> Using equation (7) of Ronchi et al: +// --------------------------------------------------------- + +void Grid::convert_vector_ll_to_xn(arma_mat aLon, + arma_mat aLat, + arma_mat &aXi, + arma_mat &aNu, + cubesphere_chars grid) { + + // Ronchi defines aPhi = aLon, aTheta = -aLat + aXi = -grid.Axt % aLat + grid.Axp % aLon; + aNu = -grid.Ant % aLat + grid.Anp % aLon; + return; +} + + + +// ---------------------------------------------------------------------- +// This function scales the deltas in the grid by the radius +// - This assumes that radius is not dependent on lat / lon!!! +// ---------------------------------------------------------------------- + +void Grid::scale_cube_by_radius(cubesphere_chars &cubeX) { + + int64_t iZ; + + for (iZ = 1; iZ < nZ; iZ++) { + cubeX.R(iZ) = radius_scgc(nGCs, nGCs, iZ); + // These are distances: + cubeX.dlx.slice(iZ) = + cubeX.dlx.slice(0) * cubeX.R(iZ); + cubeX.dln.slice(iZ) = + cubeX.dln.slice(0) * cubeX.R(iZ); + // This is an area: + cubeX.dS.slice(iZ) = + cubeX.dS.slice(0) * cubeX.R(iZ) * cubeX.R(iZ); + } + + // Lastly, scale the 0th slice + iZ = 0; + cubeX.R(iZ) = radius_scgc(nGCs, nGCs, iZ); + cubeX.dlx.slice(iZ) = + cubeX.dlx.slice(0) * cubeX.R(iZ); + cubeX.dln.slice(iZ) = + cubeX.dln.slice(0) * cubeX.R(iZ); + // This is an area: + cubeX.dS.slice(iZ) = + cubeX.dS.slice(0) * cubeX.R(iZ) * cubeX.R(iZ); + + return; +} + // ---------------------------------------------------------------------- // This function takes the normalized coordinates and makes latitude // and longitude arrays from them. It can do this for the corners or @@ -85,10 +483,10 @@ void fill_cubesphere_lat_lon_from_norms(Quadtree quadtree, int64_t nGCs, precision_t left_off, precision_t down_off, - arma_mat &lat2d, - arma_mat &lon2d, - arma_mat &refx, - arma_mat &refy) { + arma_mat & lat2d, + arma_mat & lon2d, + arma_mat & refx, + arma_mat & refy) { int64_t nX = lat2d.n_rows; int64_t nY = lat2d.n_cols; @@ -172,25 +570,25 @@ void fill_cubesphere_lat_lon_from_norms(Quadtree quadtree, // generate transformation and metric tensors // ---------------------------------------------------------------------- void transformation_metrics(Quadtree quadtree, - arma_mat &lat2d, - arma_mat &lon2d, - arma_mat &refx, - arma_mat &refy, - arma_mat &A11, - arma_mat &A12, - arma_mat &A21, - arma_mat &A22, - arma_mat &A11_inv, - arma_mat &A12_inv, - arma_mat &A21_inv, - arma_mat &A22_inv, - arma_mat &g11_upper, - arma_mat &g12_upper, - arma_mat &g21_upper, - arma_mat &g22_upper, - arma_mat &sqrt_g, - arma_mat &refx_angle, - arma_mat &refy_angle) { + arma_mat & lat2d, + arma_mat & lon2d, + arma_mat & refx, + arma_mat & refy, + arma_mat & A11, + arma_mat & A12, + arma_mat & A21, + arma_mat & A22, + arma_mat & A11_inv, + arma_mat & A12_inv, + arma_mat & A21_inv, + arma_mat & A22_inv, + arma_mat & g11_upper, + arma_mat & g12_upper, + arma_mat & g21_upper, + arma_mat & g22_upper, + arma_mat & sqrt_g, + arma_mat & refx_angle, + arma_mat & refy_angle) { int64_t nX = lat2d.n_rows; int64_t nY = lat2d.n_cols; @@ -343,6 +741,12 @@ void Grid::create_cubesphere_grid(Quadtree quadtree) { du = size_up_norm / (nLats - 2 * nGCs); ll = lower_left_norm; + // This function builds the equal-angle grid, but doesn't + // scale them with altitude, since that has not been created, yet: + init_cubesphere_grid(quadtree, dr, du, ll, 0.5, 0.5, cubeC); + init_cubesphere_grid(quadtree, dr, du, ll, 0.0, 0.5, cubeL); + init_cubesphere_grid(quadtree, dr, du, ll, 0.5, 0.0, cubeD); + int64_t iAlt, iLon, iLat; // --------------------------------------------- diff --git a/src/init_geo_grid.cpp b/src/init_geo_grid.cpp index 95138d5b..d61fb623 100644 --- a/src/init_geo_grid.cpp +++ b/src/init_geo_grid.cpp @@ -206,12 +206,10 @@ bool Grid::init_geo_grid(Quadtree quadtree, // report.print(1, "Restarting! Reading grid files!"); // DidWork = read_restart(input.get_restartin_dir()); //} else { - if (iGridShape_ == iCubesphere_) { - //if (input.get_do_restart()) - // report.print(0, "Not restarting the grid - it is too complicated!"); - + if (iGridShape_ == iCubesphere_) create_cubesphere_grid(quadtree); - } else + + else create_sphere_grid(quadtree); //MPI_Barrier(aether_comm); @@ -227,8 +225,15 @@ bool Grid::init_geo_grid(Quadtree quadtree, // Correct the reference grid with correct length scale: // (with R = actual radius) - if (iGridShape_ == iCubesphere_) + if (iGridShape_ == iCubesphere_) { correct_xy_grid(planet); + // New functions for equal-angular grid (center, left, down): + report.print(3, "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"); + } if (IsMagGrid) { report.print(0, "--> Grid is Magnetic, so rotating"); diff --git a/src/inputs.cpp b/src/inputs.cpp index b2bae93d..d4303538 100644 --- a/src/inputs.cpp +++ b/src/inputs.cpp @@ -778,7 +778,8 @@ bool Inputs::get_do_ionization_heating() { // ----------------------------------------------------------------------- bool Inputs::get_do_electron_ion_collisional_heating() { - return get_setting_bool("Sources", "Ions", "IncludeElectronIonCollisionalHeating"); + return get_setting_bool("Sources", "Ions", + "IncludeElectronIonCollisionalHeating"); } // ----------------------------------------------------------------------- @@ -786,7 +787,8 @@ bool Inputs::get_do_electron_ion_collisional_heating() { // ----------------------------------------------------------------------- bool Inputs::get_do_electron_neutral_elastic_collisional_heating() { - return get_setting_bool("Sources", "Ions", "IncludeElectronNeutralElasticCollisionalHeating"); + return get_setting_bool("Sources", "Ions", + "IncludeElectronNeutralElasticCollisionalHeating"); } // ----------------------------------------------------------------------- @@ -794,7 +796,8 @@ bool Inputs::get_do_electron_neutral_elastic_collisional_heating() { // ----------------------------------------------------------------------- bool Inputs::get_do_electron_neutral_inelastic_collisional_heating() { - return get_setting_bool("Sources", "Ions", "IncludeElectronNeutralInelasticCollisionalHeating"); + return get_setting_bool("Sources", "Ions", + "IncludeElectronNeutralInelasticCollisionalHeating"); } // ----------------------------------------------------------------------- @@ -1199,6 +1202,10 @@ std::string Inputs::get_advection_neutrals_vertical() { return get_setting_str("Advection", "Neutrals", "Vertical"); } +std::string Inputs::get_advection_neutrals_horizontal() { + return get_setting_str("Advection", "Neutrals", "Horizontal"); +} + std::string Inputs::get_advection_ions_along() { return get_setting_str("Advection", "Ions", "Along"); } diff --git a/src/ions.cpp b/src/ions.cpp index 075c8df0..bb372a33 100644 --- a/src/ions.cpp +++ b/src/ions.cpp @@ -107,6 +107,25 @@ Ions::Ions(Grid grid, Planets planet) { velocity_vcgc = make_cube_vector(nLons, nLats, nAlts, 3); cMax_vcgc = make_cube_vector(nLons, nLats, nAlts, 3); + // Some variables that will be used in calc_ion_v: + gravity_vcgc = make_cube_vector(nLons, nLats, nAlts, 3); + wind_acc = make_cube_vector(nLons, nLats, nAlts, 3); + total_acc = make_cube_vector(nLons, nLats, nAlts, 3); + efield_acc = make_cube_vector(nLons, nLats, nAlts, 3); + a_par = make_cube_vector(nLons, nLats, nAlts, 3); + a_perp = make_cube_vector(nLons, nLats, nAlts, 3); + a_x_b = make_cube_vector(nLons, nLats, nAlts, 3); + grad_Pi_plus_Pe = make_cube_vector(nLons, nLats, nAlts, 3); + rho.set_size(nLons, nLats, nAlts); + nuin.set_size(nLons, nLats, nAlts); + nuin_sum.set_size(nLons, nLats, nAlts); + Nie.set_size(nLons, nLats, nAlts); + sum_rho.set_size(nLons, nLats, nAlts); + top.set_size(nLons, nLats, nAlts); + bottom.set_size(nLons, nLats, nAlts); + + + Cv_scgc.set_size(nLons, nLats, nAlts); Cv_scgc.zeros(); lambda.set_size(nLons, nLats, nAlts); @@ -260,7 +279,7 @@ void Ions::nan_test(std::string variable) { // Checks for nans and +/- infinities in density, temp, and velocity //---------------------------------------------------------------------- -bool Ions::check_for_nonfinites() { +bool Ions::check_for_nonfinites(std::string location) { bool didWork = true; if (!all_finite(density_scgc, "density_scgc") || @@ -269,7 +288,7 @@ bool Ions::check_for_nonfinites() { didWork = false; if (!didWork) - throw std::string("Check for nonfinites failed!!!\n"); + report.error("ions are nan from location : " + location); return didWork; } @@ -495,7 +514,7 @@ void Ions::fill_electrons() { // Will return nSpecies for electrons //---------------------------------------------------------------------- -int Ions::get_species_id(const std::string &name) const{ +int Ions::get_species_id(const std::string &name) const { std::string function = "Ions::get_species_id"; static int iFunction = -1; diff --git a/src/main/main.cpp b/src/main/main.cpp index f5e6957a..fcb62ce0 100644 --- a/src/main/main.cpp +++ b/src/main/main.cpp @@ -9,7 +9,7 @@ // ----------------------------------------------------------------------------- int main() { - + int iErr = 0; std::string sError; bool didWork = true, testsPassing=true; @@ -24,6 +24,7 @@ int main() { try { // Create inputs (reading the input file): input = Inputs(time); + if (!input.is_ok()) throw std::string("input initialization failed!"); @@ -39,27 +40,32 @@ int main() { // cubesphere (6 root) Quadtree quadtree(input.get_grid_shape("neuGrid")); Quadtree quadtree_ion(input.get_grid_shape("ionGrid")); + if (!quadtree.is_ok()) throw std::string("quadtree initialization failed!"); // Initialize MPI and parallel aspects of the code: didWork = init_parallel(quadtree, quadtree_ion); + if (!didWork) throw std::string("init_parallel failed!"); // Everything should be set for the inputs now, so write a restart file: didWork = input.write_restart(); + if (!didWork) throw std::string("input.write_restart failed!"); // Initialize the EUV system: Euv euv; + if (!euv.is_ok()) throw std::string("EUV initialization failed!"); // Initialize the planet: Planets planet; MPI_Barrier(aether_comm); + if (!planet.is_ok()) throw std::string("planet initialization failed!"); @@ -67,6 +73,7 @@ int main() { Indices indices; didWork = read_and_store_indices(indices); MPI_Barrier(aether_comm); + if (!didWork) throw std::string("read_and_store_indices failed!"); @@ -77,6 +84,7 @@ int main() { Grid gGrid("neuGrid"); didWork = gGrid.init_geo_grid(quadtree, planet); MPI_Barrier(aether_comm); + if (!didWork) throw std::string("init_geo_grid failed!"); @@ -95,6 +103,8 @@ int main() { if (mGrid.iGridShape_ == mGrid.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"; mGrid.set_IsDipole(false); @@ -136,9 +146,12 @@ int main() { if (input.get_check_for_nans()) { didWork = neutrals.check_for_nonfinites("After Inputs"); + if (!didWork) throw std::string("NaNs found in Neutrals in Initialize!\n"); - didWork = ions.check_for_nonfinites(); + + didWork = ions.check_for_nonfinites("NaNs found in Ions in Initialize!\n"); + if (!didWork) throw std::string("NaNs found in Ions in Initialize!\n"); } @@ -164,9 +177,12 @@ int main() { // Initialize electrodynamics and check if electrodynamics times // works with input time Electrodynamics electrodynamics(time); + if (!electrodynamics.is_ok()) throw std::string("electrodynamics on geo grid initialization failed!"); + Electrodynamics electrodynamicsMag(time); + if (!electrodynamicsMag.is_ok()) throw std::string("electrodynamics on mag grid initialization failed!"); @@ -184,6 +200,7 @@ int main() { didWork = output(neutrals, ions, gGrid, time, planet); didWork = output(neutralsMag, ionsMag, mGrid, time, planet); } + if (!didWork) throw std::string("Initial output failed!"); @@ -244,14 +261,18 @@ int main() { if (!time.check_time_gate(input.get_dt_write_restarts())) { report.print(3, "Writing restart files"); - didWork = neutrals.restart_file(input.get_restartout_dir(), gGrid.get_gridtype(), DoWrite); - didWork = neutralsMag.restart_file(input.get_restartout_dir(), mGrid.get_gridtype(), DoWrite); + didWork = neutrals.restart_file(input.get_restartout_dir(), + gGrid.get_gridtype(), DoWrite); + didWork = neutralsMag.restart_file(input.get_restartout_dir(), + mGrid.get_gridtype(), DoWrite); if (!didWork) throw std::string("Writing Restart for Neutrals Failed!!!\n"); - didWork = ions.restart_file(input.get_restartout_dir(), gGrid.get_gridtype(), DoWrite); - didWork = ionsMag.restart_file(input.get_restartout_dir(), mGrid.get_gridtype(), DoWrite); + didWork = ions.restart_file(input.get_restartout_dir(), gGrid.get_gridtype(), + DoWrite); + didWork = ionsMag.restart_file(input.get_restartout_dir(), mGrid.get_gridtype(), + DoWrite); if (!didWork) throw std::string("Writing Restart for Ions Failed!!!\n"); diff --git a/src/neutrals_advect.cpp b/src/neutrals_advect.cpp index a80dc274..f2b11652 100644 --- a/src/neutrals_advect.cpp +++ b/src/neutrals_advect.cpp @@ -33,3 +33,43 @@ bool Neutrals::advect_vertical(Grid grid, Times time) { return didWork; } +bool Neutrals::advect_horizontal(Grid & grid, Times & time) { + bool didWork = true; + + std::string function = "Neutrals::advance_horizontal"; + static int iFunction = -1; + report.enter(function, iFunction); + + if (grid.iGridShape_ == grid.iCubesphere_) { + if (input.get_advection_neutrals_horizontal() == "advect_test") + solver_horizontal_RK4_advection(grid, time); + else if (input.get_advection_neutrals_horizontal() == "fv") + solver_horizontal_RK1_rochi(grid, time); + else { + std::cout << "Horizontal solver not found!\n"; + std::cout << " ==> Requested : " + << input.get_advection_neutrals_horizontal() + << "\n"; + didWork = false; + } + } else + advect_sphere(grid, time); + + + report.exit(function); + return didWork; +} + +bool Neutrals::advect_horizontal_advection(Grid & grid, Times & time) { + bool didWork = true; + + std::string function = "Neutrals::advance_horizontal_advection"; + static int iFunction = -1; + report.enter(function, iFunction); + + //solver_horizontal_rusanov_advection(grid, time); + solver_horizontal_RK4_advection(grid, time); + + report.exit(function); + return didWork; +} diff --git a/src/neutrals_ics.cpp b/src/neutrals_ics.cpp index 8da58901..7cefe289 100644 --- a/src/neutrals_ics.cpp +++ b/src/neutrals_ics.cpp @@ -35,6 +35,8 @@ bool Neutrals::initial_conditions(Grid grid, precision_t alt, r; int64_t nAlts = grid.get_nZ(true); int64_t nGCs = grid.get_nGCs(); + int64_t nLons = grid.get_nLons(); + int64_t nLats = grid.get_nLats(); report.print(3, "Creating Neutrals initial_condition"); @@ -120,10 +122,6 @@ bool Neutrals::initial_conditions(Grid grid, // temperature profile in the planet.in file. // --------------------------------------------------------------------- - int64_t nLons = grid.get_nLons(); - int64_t nLats = grid.get_nLats(); - int64_t nAlts = grid.get_nAlts(); - // Let's assume that the altitudes are not dependent on lat/lon: arma_vec alt1d(nAlts); @@ -194,6 +192,41 @@ bool Neutrals::initial_conditions(Grid grid, } // type = planet } + /* + This section is for putting an initial blob into the simulation + to test the advection solver. + precision_t lon_0 = 0.0; + precision_t lat_0 = 0.0; + precision_t r_0 = 150.0 * 1000.0 * 10.0; + + for (iAlt = 0; iAlt < nAlts; iAlt++) { + + for (int64_t iLat = 0; iLat < nLats; iLat++) { + for (int64_t iLon = 0; iLon < nLons; iLon++) { + precision_t curr_lat = grid.geoLat_scgc(iLon, iLat, iAlt); + precision_t curr_lon = grid.geoLon_scgc(iLon, iLat, iAlt); + precision_t R = grid.radius_scgc(iLon, iLat, iAlt); + + // Calculate great circle distance + precision_t dlon_2 = (curr_lon - lon_0) / 2.0; + precision_t dlat_2 = (curr_lat - lat_0) / 2.0; + + precision_t r_d = 2.0 * R * asin(sqrt(sin(dlat_2) * sin(dlat_2) + sin( + dlon_2) * sin(dlon_2) * cos(curr_lat) * cos(lat_0))); + + if (r_d < r_0) { + for (int iSpecies = 0; iSpecies < nSpecies; iSpecies++) { + species[iSpecies].density_scgc(iLon, iLat, + iAlt) = species[iSpecies].density_scgc(iLon, iLat, iAlt) * 10.; + std::cout << "increasing density!\n"; + } + } + } + } + } + */ + + // ensure that the densities are all within bounds: clamp_density(); @@ -204,5 +237,221 @@ bool Neutrals::initial_conditions(Grid grid, return didWork; } +bool Neutrals::cosine_bell_ic(Grid grid, + Times time, + Indices indices, + Planets planet) { + std::string function = "Neutrals::cosine_bell_ic"; + static int iFunction = -1; + report.enter(function, iFunction); + + // Planet.get_radius() takes in latitude + // but at current stage is unimplemented + // Anyway, we use equator radius as assumption for CubeSphere + // CubeSphere must be a perfect sphere!! + precision_t planet_R = planet.get_radius(0); + + // radius of planet + altitude + // just pick alt at (0,0) loction + arma_vec R_Alts = grid.geoAlt_scgc.tube(0, 0) + planet_R; + + + /** Get a bunch of constants for setting up the ic **/ + precision_t R = R_Alts(2); // select R in the middle + //7.37128e+06 meters Earth radius + 1000km height + //std::cout << R << std::endl; + + // Determine flow direction + // 0 - Equatorial + // pi/2 - Meridional + // pi/4 - NE direction + precision_t alpha_0 = 0.; + + // Scaling factor for physical velocity + // 12 day period in miliseconds + precision_t u_0 = cTWOPI * R / (12.*24.*60.*60.); + + // Radius of the cosine bell + precision_t r_0 = R / 3.; + + // Center of the cosine bell + precision_t lon_0 = 7.*cPI / 4. + cPI / 8; + precision_t lat_0 = 0.; + + // Maximum height for the cosine bell + precision_t h_0 = 1000.; + + // Some grid dimensions and coordinates + int64_t nLats = grid.get_nLats(); + int64_t nLons = grid.get_nLons(); + int64_t nAlts = grid.get_nAlts(); + arma_mat lat_grid = grid.geoLat_scgc.slice(2); + arma_mat lon_grid = grid.geoLon_scgc.slice(2); + + // Calculate for physical velocity for every altitude + // First we prepare velocities for one slice + arma_mat slice_u = velocity_vcgc[0].slice(2); + arma_mat slice_v = velocity_vcgc[1].slice(2); + + // Fill velocities in one slice + for (int64_t iLat = 0; iLat < nLats; iLat++) { + for (int64_t iLon = 0; iLon < nLons; iLon++) { + precision_t curr_lat = lat_grid(iLon, iLat); + precision_t curr_lon = lon_grid(iLon, iLat); + slice_u(iLon, iLat) = u_0 * (cos(alpha_0) * cos(curr_lat) + sin(alpha_0) * cos( + curr_lon) * sin(curr_lat)); + slice_v(iLon, iLat) = -u_0 * sin(alpha_0) * sin(curr_lon); + } + } + + // Update this slice of velocity to all slices (for completeness) + for (int64_t iAlt = 0; iAlt < nAlts; iAlt++) { + velocity_vcgc[0].slice(iAlt) = slice_u; + velocity_vcgc[1].slice(iAlt) = slice_v; + } + + // Calculate the cosine bell or rho_scgc + // First, again take a slice + arma_mat slice_rho = rho_scgc.slice(2); + + // Fill rho in one slice + for (int64_t iLat = 0; iLat < nLats; iLat++) { + for (int64_t iLon = 0; iLon < nLons; iLon++) { + precision_t curr_lat = lat_grid(iLon, iLat); + precision_t curr_lon = lon_grid(iLon, iLat); + + // Calculate great circle distance + precision_t dlon_2 = (curr_lon - lon_0) / 2.0; + precision_t dlat_2 = (curr_lat - lat_0) / 2.0; + + precision_t r_d = 2.0 * R * asin(sqrt(sin(dlat_2) * sin(dlat_2) + sin( + dlon_2) * sin(dlon_2) * cos(curr_lat) * cos(lat_0))); + + if (r_d < r_0) + slice_rho(iLon, iLat) = (h_0 / 2) * (1 + cos(cPI * r_d / r_0)); + + else + slice_rho(iLon, iLat) = 0.; + } + } + + // Update this slice of rho to all slices (for completeness) + for (int64_t iAlt = 0; iAlt < nAlts; iAlt++) { + rho_scgc.slice(iAlt) = slice_rho; + + // Do zero concentration conversion + for (int64_t iSpec = 0; iSpec < nSpecies; iSpec++) + species[iSpec].density_scgc.slice(iAlt) = slice_rho / species[iSpec].mass; + } + + // Add some velocity pertubation + //std::cout << velocity_vcgc[0].slice(2) << std::endl; + //std::cout << rho_scgc.slice(2) << std::endl; + + return 1; +} + +bool Neutrals::blob_ic(Grid grid, + Times time, + Indices indices, + Planets planet) { + std::string function = "Neutrals::blob_ic"; + static int iFunction = -1; + report.enter(function, iFunction); + + // Planet.get_radius() takes in latitude + // but at current stage is unimplemented + // Anyway, we use equator radius as assumption for CubeSphere + // CubeSphere must be a perfect sphere!! + precision_t planet_R = planet.get_radius(0); + + // radius of planet + altitude + // just pick alt at (0,0) loction + arma_vec R_Alts = grid.geoAlt_scgc.tube(0, 0) + planet_R; + + + /** Get a bunch of constants for setting up the ic **/ + precision_t R = R_Alts(2); // select R in the middle + //7.37128e+06 meters Earth radius + 1000km height + //std::cout << R << std::endl; + + // Radius of the blob + // Hardcoded + precision_t r_0 = 111321 * 10; + + // Center of the blob + precision_t lon_0 = 7.*cPI / 4. - cPI / 8.; + precision_t lat_0 = 0.; + + // Some grid dimensions and coordinates + int64_t nLats = grid.get_nLats(); + int64_t nLons = grid.get_nLons(); + int64_t nAlts = grid.get_nAlts(); + arma_mat lat_grid = grid.geoLat_scgc.slice(2); + arma_mat lon_grid = grid.geoLon_scgc.slice(2); + + // Calculate for physical velocity for every altitude + // First we prepare velocities for one slice + arma_mat slice_u = velocity_vcgc[0].slice(2); + arma_mat slice_v = velocity_vcgc[1].slice(2); + + // Fill velocities in one slice + for (int64_t iLat = 0; iLat < nLats; iLat++) { + for (int64_t iLon = 0; iLon < nLons; iLon++) { + slice_u(iLon, iLat) = 0.; + slice_v(iLon, iLat) = 0.; + } + } + + // Update this slice of velocity to all slices (for completeness) + for (int64_t iAlt = 0; iAlt < nAlts; iAlt++) { + velocity_vcgc[0].slice(iAlt) = slice_u; + velocity_vcgc[1].slice(iAlt) = slice_v; + } + + // Calculate the cosine bell or rho_scgc + // First, again take a slice + arma_mat slice_rho = rho_scgc.slice(2); + + // Fill rho in one slice + for (int64_t iLat = 0; iLat < nLats; iLat++) { + for (int64_t iLon = 0; iLon < nLons; iLon++) { + precision_t curr_lat = lat_grid(iLon, iLat); + precision_t curr_lon = lon_grid(iLon, iLat); + + // Calculate great circle distance + precision_t dlon_2 = (curr_lon - lon_0) / 2.0; + precision_t dlat_2 = (curr_lat - lat_0) / 2.0; + + precision_t r_d = 2.0 * R * asin(sqrt(sin(dlat_2) * sin(dlat_2) + sin( + dlon_2) * sin(dlon_2) * cos(curr_lat) * cos(lat_0))); + + if (r_d < r_0) + slice_rho(iLon, iLat) = 5e-12; + + else + slice_rho(iLon, iLat) = 1e-12; + } + } + + // Update this slice of rho to all slices (for completeness) + for (int64_t iAlt = 0; iAlt < nAlts; iAlt++) { + rho_scgc.slice(iAlt) = slice_rho; + + // Do zero concentration conversion + for (int64_t iSpec = 0; iSpec < nSpecies; iSpec++) + species[iSpec].density_scgc.slice(iAlt) = slice_rho / species[iSpec].mass; + } + + // Temperature setup + for (int64_t iAlt = 0; iAlt < nAlts; iAlt++) + temperature_scgc.slice(iAlt) = 600.*arma_mat(nLons, nLats, fill::ones); + + // Add some velocity pertubation + //std::cout << velocity_vcgc[0].slice(2) << std::endl; + //std::cout << rho_scgc.slice(2) << std::endl; + + return 1; +} diff --git a/src/solver_advection.cpp b/src/solver_advection.cpp index 1e69c7bc..48e88688 100644 --- a/src/solver_advection.cpp +++ b/src/solver_advection.cpp @@ -282,9 +282,7 @@ precision_t calc_dt(arma_mat &xWidth, // // --------------------------------------------------------- -void advect(Grid &grid, - Times &time, - Neutrals &neutrals) { +void Neutrals::advect_sphere(Grid &grid, Times &time) { std::string function = "advect"; static int iFunction = -1; @@ -337,9 +335,9 @@ void advect(Grid &grid, arma_mat gamma2d; // These are all needed by the solver: - neutrals.calc_mass_density(); - neutrals.calc_mean_major_mass(); - neutrals.calc_specific_heat(); + calc_mass_density(); + calc_mean_major_mass(); + calc_specific_heat(); arma_mat t_to_e; @@ -350,14 +348,14 @@ void advect(Grid &grid, if (report.test_verbose(3)) std::cout << "Advection: Working with iAlt: " << iAlt << "\n"; - xVel = neutrals.velocity_vcgc[0].slice(iAlt); - yVel = neutrals.velocity_vcgc[1].slice(iAlt); - rho = neutrals.rho_scgc.slice(iAlt); + xVel = velocity_vcgc[0].slice(iAlt); + yVel = velocity_vcgc[1].slice(iAlt); + rho = rho_scgc.slice(iAlt); // this is "e", or temperature expressed as an energy - gamma2d = neutrals.gamma_scgc.slice(iAlt); - t_to_e = 1.0 / (gamma2d - 1.0) * cKB / neutrals.mean_major_mass_scgc.slice( + gamma2d = gamma_scgc.slice(iAlt); + t_to_e = 1.0 / (gamma2d - 1.0) * cKB / mean_major_mass_scgc.slice( iAlt); - temp = t_to_e % neutrals.temperature_scgc.slice(iAlt); + temp = t_to_e % temperature_scgc.slice(iAlt); // ------------------------------------------------ // Calculate derived equations (at cell centers - these will be updated) @@ -542,8 +540,8 @@ void advect(Grid &grid, xVel = xMomentum / rho; yVel = yMomentum / rho; - neutrals.velocity_vcgc[0].slice(iAlt) = xVel; - neutrals.velocity_vcgc[1].slice(iAlt) = yVel; + velocity_vcgc[0].slice(iAlt) = xVel; + velocity_vcgc[1].slice(iAlt) = yVel; temp = (totalE / rho - 0.5 * (xVel % xVel + yVel % yVel)) / t_to_e; temp.clamp(200, 2000); @@ -555,35 +553,35 @@ void advect(Grid &grid, //if (cos(grid.geoLat_scgc(i,j,iAlt)) < 0.2) { // fac = fac * (0.2 - cos(grid.geoLat_scgc(i,j,iAlt))); //} - //dm = (1.0 - fac) * neutrals.temperature_scgc(i,j,iAlt); - //dp = (1.0 + fac) * neutrals.temperature_scgc(i,j,iAlt); + //dm = (1.0 - fac) * temperature_scgc(i,j,iAlt); + //dp = (1.0 + fac) * temperature_scgc(i,j,iAlt); //if (temp(i,j) < dm) temp(i,j) = dm; //if (temp(i,j) > dp) temp(i,j) = dp; - neutrals.temperature_scgc(i, j, iAlt) = temp(i, j); + temperature_scgc(i, j, iAlt) = temp(i, j); - //dm = (1.0 - fac) * neutrals.rho_scgc(i,j,iAlt); - //dp = (1.0 + fac) * neutrals.rho_scgc(i,j,iAlt); + //dm = (1.0 - fac) * rho_scgc(i,j,iAlt); + //dp = (1.0 + fac) * rho_scgc(i,j,iAlt); //if (rho(i,j) < dm) rho(i,j) = dm; //if (rho(i,j) > dp) rho(i,j) = dp; - neutrals.rho_scgc(i, j, iAlt) = rho(i, j); + rho_scgc(i, j, iAlt) = rho(i, j); } } if (report.test_verbose(3) && iAlt == 8) { - std::cout << "end t : " << neutrals.temperature_scgc.slice( - iAlt).min() << " " << neutrals.temperature_scgc.slice(iAlt).max() << "\n"; + std::cout << "end t : " << temperature_scgc.slice( + iAlt).min() << " " << temperature_scgc.slice(iAlt).max() << "\n"; std::cout << "end temp : " << temp.min() << " " << temp.max() << "\n"; std::cout << "end xVel : " << xVel.min() << " " << xVel.max() << "\n"; std::cout << "end yVel : " << yVel.min() << " " << yVel.max() << "\n"; } } - neutrals.calc_density_from_mass_concentration(); + calc_density_from_mass_concentration(); // Assign bulk horizontal velocity to all species: - for (int64_t iSpecies = 0; iSpecies < neutrals.nSpecies; iSpecies++) + for (int64_t iSpecies = 0; iSpecies < nSpecies; iSpecies++) for (int64_t iDir = 0; iDir < 2; iDir++) - neutrals.species[iSpecies].velocity_vcgc[iDir] = neutrals.velocity_vcgc[iDir]; + species[iSpecies].velocity_vcgc[iDir] = velocity_vcgc[iDir]; report.exit(function); return; diff --git a/src/solver_horizontal_cubesphere.cpp b/src/solver_horizontal_cubesphere.cpp index ebce8b2c..98868b1a 100644 --- a/src/solver_horizontal_cubesphere.cpp +++ b/src/solver_horizontal_cubesphere.cpp @@ -2,5 +2,1215 @@ // Full license can be found in License.md // Initial version: F. Cheng, July 2023 +// Moved to new solver: August 2025 -#include "../include/aether.h" \ No newline at end of file +#include "aether.h" + +// --------------------------------------------------------- +// Update States +// --------------------------------------------------------- + +void update_states_cubesphere(arma_mat rho, + arma_mat &xVel, + arma_mat &yVel, + arma_mat &temp, + arma_mat &drhodt, + arma_mat &dlonVeldt, + arma_mat &dlatVeldt, + arma_mat &dtempdt, + cubesphere_chars gridC, + cubesphere_chars gridL, + cubesphere_chars gridD, + precision_t dt, + int64_t iZ) { + + arma_mat xMomentum, yMomentum; + arma_mat rhoE, energy, vel2; + + precision_t cv = 1500.0; + + if (report.test_verbose(2)) + std::cout << " --> update_states\n"; + + // Derived variables: + xMomentum = rho % xVel; // x1momentum, pure scalar field + yMomentum = rho % yVel; // y1momentum, pure scalar field + rhoE = rho % temp; + + vel2 = xVel % xVel + yVel % yVel; + //energy = rho % (0.5 * vel2 + cv * temp); + energy = cv * rho % temp; + + /** Initialize projection constructs */ + static projection_struct rhoP; + static projection_struct xMomentumP, xVelP; + static projection_struct yMomentumP, yVelP; + static projection_struct energyP; + static projection_struct tempP; + + // They are all pure scalar fields without sqrt(g) + static arma_mat totaleL, totaleR, totaleD, totaleU; + static arma_mat velL2, velR2, velD2, velU2; + static arma_mat pressureL, pressureR, pressureD, pressureU; + + arma_mat dxVeldt = xVel * 0.0; + arma_mat dyVeldt = yVel * 0.0; + + dlonVeldt = dxVeldt * 0.0 + 1; + dlatVeldt = dyVeldt * 0.0 + 1; + + static arma_mat velNormL, velNormR, velNormU, velNormD; + + /** Initialize Flux and Wave Speed Storages */ + static arma_mat eq1FluxLR, eq1FluxDU; + static arma_mat eq1FluxL, eq1FluxR, eq1FluxD, eq1FluxU; + static arma_mat eq2FluxLR, eq2FluxDU; + static arma_mat eq2FluxL, eq2FluxR, eq2FluxD, eq2FluxU; + static arma_mat eq3FluxLR, eq3FluxDU; + static arma_mat eq3FluxL, eq3FluxR, eq3FluxD, eq3FluxU; + static arma_mat eq4FluxLR, eq4FluxDU; + static arma_mat eq4FluxL, eq4FluxR, eq4FluxD, eq4FluxU; + + arma_mat wsL, wsR, wsD, wsU, wsLR, wsDU; + + arma_mat diff; // for Riemann Solver + + if (report.test_verbose(3)) + std::cout << " ---> Projecting\n"; + + rhoP = project_to_edges(rho, gridC.xi, gridL.xi, gridC.nu, gridD.nu, + gridC.nGCs); + // project the lon / lat velocities to the edges: + xVelP = project_to_edges(xVel, gridC.xi, gridL.xi, gridC.nu, gridD.nu, + gridC.nGCs); + yVelP = project_to_edges(yVel, gridC.xi, gridL.xi, gridC.nu, gridD.nu, + gridC.nGCs); + xMomentumP = project_to_edges(xMomentum, gridC.xi, gridL.xi, gridC.nu, gridD.nu, + gridC.nGCs); + yMomentumP = project_to_edges(yMomentum, gridC.xi, gridL.xi, gridC.nu, gridD.nu, + gridC.nGCs); + energyP = project_to_edges(energy, gridC.xi, gridL.xi, gridC.nu, gridD.nu, + gridC.nGCs); + tempP = project_to_edges(temp, gridC.xi, gridL.xi, gridC.nu, gridD.nu, + gridC.nGCs); + + if (report.test_verbose(3)) + std::cout << " ---> Derived values\n"; + + velL2 = (xVelP.L % xVelP.L + yVelP.L % yVelP.L); + velR2 = (xVelP.R % xVelP.R + yVelP.R % yVelP.R); + velD2 = (xVelP.D % xVelP.D + yVelP.D % yVelP.D); + velU2 = (xVelP.U % xVelP.U + yVelP.U % yVelP.U); + + precision_t k = 1.38e-23; + // let's be Oxygen: + precision_t mass = 16.0 * 1.67e-27; + pressureL = k / mass * (rhoP.L % tempP.L); + pressureR = k / mass * (rhoP.R % tempP.R); + pressureD = k / mass * (rhoP.D % tempP.D); + pressureU = k / mass * (rhoP.U % tempP.U); + + arma_mat pressureLR = (pressureL + pressureR) / 2; + arma_mat pressureDU = (pressureD + pressureU) / 2; + + if (report.test_verbose(3)) + std::cout << " ---> Normal Velocities\n"; + + // Calculate the normal velocity at the boundaries: + velNormL = xVelP.L % gridL.nXiLon + yVelP.L % gridL.nXiLat; + velNormR = xVelP.R % gridL.nXiLon + yVelP.R % gridL.nXiLat; + velNormU = xVelP.U % gridD.nNuLon + yVelP.U % gridD.nNuLat; + velNormD = xVelP.D % gridD.nNuLon + yVelP.D % gridD.nNuLat; + + if (report.test_verbose(3)) + std::cout << " ---> Fluxes eq 1\n"; + + // Flux calculated from the left of the edge + eq1FluxL = rhoP.L % velNormL; + // Flux calculated from the right of the edge + eq1FluxR = rhoP.R % velNormR; + // Flux calculated from the down of the edge + eq1FluxD = rhoP.D % velNormD; + // Flux calculated from the up of the edge + eq1FluxU = rhoP.U % velNormU; + + if (report.test_verbose(3)) + std::cout << " ---> Fluxes eq 2\n"; + + eq2FluxL = (xMomentumP.L % velNormL); + eq2FluxR = (xMomentumP.R % velNormR); + eq2FluxD = (xMomentumP.D % velNormD); + eq2FluxU = (xMomentumP.U % velNormU); + + if (report.test_verbose(3)) + std::cout << " ---> Fluxes eq 3\n"; + + eq3FluxL = (yMomentumP.L % velNormL); + eq3FluxR = (yMomentumP.R % velNormR); + eq3FluxD = (yMomentumP.D % velNormD); + eq3FluxU = (yMomentumP.U % velNormU); + + eq4FluxL = energyP.L % velNormL; + eq4FluxR = energyP.R % velNormR; + eq4FluxD = energyP.D % velNormD; + eq4FluxU = energyP.U % velNormU; + + // ------------------------------------------------ + // Calculate the wave speed for the diffusive flux: + // In Reference velocities + if (report.test_verbose(3)) + std::cout << " ---> Diffusive Fluxes\n"; + + precision_t cGamma = 5.0 / 3.0; + + wsL.resize(gridC.nXt + 1, gridC.nYt); + wsR.resize(gridC.nXt + 1, gridC.nYt); + wsD.resize(gridC.nXt, gridC.nYt + 1); + wsU.resize(gridC.nXt, gridC.nYt + 1); + + wsL.zeros(); + wsR.zeros(); + wsU.zeros(); + wsD.zeros(); + + for (int64_t i = 0; i < gridC.nXt; i++) { + for (int64_t j = 0; j < gridC.nYt; j++) { + wsL(i, j) = sqrt(velL2(i, j)) + sqrt(cGamma * (cGamma - 1) * tempP.L(i, j)); + wsR(i, j) = sqrt(velR2(i, j)) + sqrt(cGamma * (cGamma - 1) * tempP.R(i, j)); + wsD(i, j) = sqrt(velD2(i, j)) + sqrt(cGamma * (cGamma - 1) * tempP.D(i, j)); + wsU(i, j) = sqrt(velU2(i, j)) + sqrt(cGamma * (cGamma - 1) * tempP.U(i, j)); + } + } + + wsLR = wsR; + + for (int64_t i = 0; i < gridC.nXt; i++) { + for (int64_t j = 0; j < gridC.nYt; j++) { + if (wsL(i, j) > wsLR(i, j)) + wsLR(i, j) = wsL(i, j); + } + } + + wsDU = wsD; + + for (int64_t i = 0; i < gridC.nXt; i++) { + for (int64_t j = 0; j < gridC.nYt; j++) { + if (wsU(i, j) > wsDU(i, j)) + wsDU(i, j) = wsU(i, j); + } + } + + // ------------------------------------------------ + // Calculate average flux at the edges (Rusanov Flux): + + if (report.test_verbose(3)) + std::cout << " ---> Averaging fluxes at edges\n"; + + diff = (rhoP.R - rhoP.L); + eq1FluxLR = (eq1FluxL + eq1FluxR) / 2 + 0.5 * wsLR % diff; + diff = (rhoP.U - rhoP.D); + eq1FluxDU = (eq1FluxD + eq1FluxU) / 2 + 0.5 * wsDU % diff; + + diff = (xMomentumP.R - xMomentumP.L); + eq2FluxLR = (eq2FluxL + eq2FluxR) / 2 + 0.5 * wsLR % diff; + diff = (xMomentumP.U - xMomentumP.D); + eq2FluxDU = (eq2FluxD + eq2FluxU) / 2 + 0.5 * wsDU % diff; + + diff = (yMomentumP.R - yMomentumP.L); + eq3FluxLR = (eq3FluxL + eq3FluxR) / 2 + 0.5 * wsLR % diff; + diff = (yMomentumP.U - yMomentumP.D); + eq3FluxDU = (eq3FluxD + eq3FluxU) / 2 + 0.5 * wsDU % diff; + + diff = (energyP.R - energyP.L); + eq4FluxLR = (eq4FluxL + eq4FluxR) / 2 + 0.5 * wsLR % diff; + diff = (energyP.U - energyP.D); + eq4FluxDU = (eq4FluxD + eq4FluxU) / 2 + 0.5 * wsDU % diff; + + // ------------------------------------------------ + // Update values: + if (report.test_verbose(3)) + std::cout << " ---> Updating equations of state\n"; + + precision_t dpdx, dpdn, pp, pm; + + arma_mat ax(gridC.nXt, gridC.nYt), an(gridC.nXt, gridC.nYt); + + ax.zeros(); + an.zeros(); + arma_mat dedt(gridC.nXt, gridC.nYt); + dedt.zeros(); + + arma_mat rhoNew = rho; + + // Only deal with inner cell + for (int64_t j = gridC.iYfirst_; j < gridC.iYlast_; j++) { + for (int64_t i = gridC.iXfirst_; i < gridC.iXlast_; i++) { + precision_t rhoResidual_ij = (gridL.dln(i + 1, j, iZ) * eq1FluxLR(i + 1, j) - + gridL.dln(i, j, iZ) * eq1FluxLR(i, j) + + gridD.dlx(i, j + 1, iZ) * eq1FluxDU(i, j + 1) - + gridD.dlx(i, j, iZ) * eq1FluxDU(i, j)); + drhodt(i, j) = rhoResidual_ij / gridC.dS(i, j, iZ); + + rhoNew(i, j) = rho(i, j) + dt * drhodt(i, j); + + precision_t xMomentumResidual_ij = (gridL.dln(i + 1, j, iZ) * eq2FluxLR(i + 1, + j) - + gridL.dln(i, j, iZ) * eq2FluxLR(i, j) + + gridD.dlx(i, j + 1, iZ) * eq2FluxDU(i, j + 1) - + gridD.dlx(i, j, iZ) * eq2FluxDU(i, j)); + dxVeldt(i, j) = xMomentumResidual_ij / gridC.dS(i, j, iZ) / rhoNew(i, j); + + precision_t yMomentumResidual_ij = (gridL.dln(i + 1, j, iZ) * eq3FluxLR(i + 1, + j) - + gridL.dln(i, j, iZ) * eq3FluxLR(i, j) + + gridD.dlx(i, j + 1, iZ) * eq3FluxDU(i, j + 1) - + gridD.dlx(i, j, iZ) * eq3FluxDU(i, j)); + dyVeldt(i, j) = yMomentumResidual_ij / gridC.dS(i, j, iZ) / rhoNew(i, j); + + // Calculate the gradient in the potential in the cubesphere + // coordinate system: + dpdx = 1 / gridC.R(iZ) * gridC.D(i, j) * + (pressureLR(i + 1, j) - pressureLR(i, j)) / gridC.dxi; + dpdn = 1 / gridC.R(iZ) * gridC.X(i, j) * gridC.Y(i, j) / + gridC.D(i, j) * + (pressureDU(i, j + 1) - pressureDU(i, j)) / gridC.dnu; + ax(i, j) = (dpdx + dpdn) / rhoNew(i, j); + + dpdx = 1 / gridC.R(iZ) * gridC.X(i, j) * gridC.Y(i, j) / + gridC.C(i, j) * (pressureLR(i + 1, j) - pressureLR(i, j)) / gridC.dxi; + dpdn = 1 / gridC.R(iZ) * gridC.C(i, j) * + (pressureDU(i, j + 1) - pressureDU(i, j)) / gridC.dnu; + an(i, j) = (dpdx + dpdn) / rhoNew(i, j); + + precision_t energyResidual_ij = (gridL.dln(i + 1, j, iZ) * eq4FluxLR(i + 1, j) - + gridL.dln(i, j, iZ) * eq4FluxLR(i, j) + + gridD.dlx(i, j + 1, iZ) * eq4FluxDU(i, j + 1) - + gridD.dlx(i, j, iZ) * eq4FluxDU(i, j)); + dedt(i, j) = energyResidual_ij / gridC.dS(i, j, iZ); + + } + } + + // lat is negative because of the Rochi definition of theta: + dlatVeldt = dyVeldt - (ax % gridC.Atx + an % gridC.Atn); + dlonVeldt = dxVeldt + ax % gridC.Apx + an % gridC.Apn; + dtempdt = dedt / rhoNew / cv; + + return; +} + + +// using namespace Cubesphere_tools; + +std::vector Neutrals::residual_horizontal_rusanov( + std::vector& states, + Grid & grid, + Times & time, + int64_t iAlt) { + + // Dimensions of Spatial Discretization + int64_t nXs = grid.get_nX(); + int64_t nYs = grid.get_nY(); + int64_t nGCs = grid.get_nGCs(); + int64_t nAlts = grid.get_nAlts(); + + /** Extract Grid Features **/ + arma_mat x = grid.refx_scgc.slice(iAlt); + arma_mat xEdges = grid.refx_Left.slice(iAlt); + arma_mat y = grid.refy_scgc.slice(iAlt); + arma_mat yEdges = grid.refy_Down.slice(iAlt); + + // Get reference grid dimensions (Assume dx = dy and equidistant) + arma_vec x_vec = x.col(0); + precision_t dx = x_vec(1) - x_vec(0); + precision_t area = dx * dx; + arma_mat jacobian = grid.sqrt_g_scgc.slice(iAlt); + + /** States preprocessing **/ + /* MASS DENSITY */ + arma_mat rho = states[0]; + + /* VELOCITY */ + // Get contravariant velocity + //arma_mat xVel = states[1]; // u^1 + //arma_mat yVel = states[2]; // u^2 + + // Generate contravriant momentum + arma_mat xMomentum = states[1]; // x1momentum + arma_mat yMomentum = states[2]; // x2momentum + + // Resolve to contravariant velocity + arma_mat xVel = xMomentum / rho; // u^1 + arma_mat yVel = yMomentum / rho; // u^2 + + // Generate velocity magnitude squared + arma_mat vel2 = xVel % xVel + yVel % yVel; + + /* TEMP and ENERGY */ + // Generate total energy (rhoE) (TODO: Verify) + arma_mat rhoE = states[3]; + + /** Advancing **/ + /* Initialize projection constructs and storages */ + projection_struct rhoP; + projection_struct xMomentumP; + projection_struct yMomentumP; + projection_struct rhoEP; + projection_struct gammaP; + projection_struct tempP; + projection_struct numberDensityP; + + // They are all pure scalar fields without sqrt(g) + arma_mat rhoL, rhoR, rhoD, rhoU; + arma_mat xVelL, xVelR, xVelD, xVelU; + arma_mat yVelL, yVelR, yVelD, yVelU; + arma_mat totalEL, totalER, totalED, totalEU; + + arma_mat velL2, velR2, velD2, velU2; + arma_mat internaleL, internaleR, internaleD, internaleU; + arma_mat pressureL, pressureR, pressureD, pressureU; + + /** Initialize Flux and Wave Speed Storages */ + arma_mat eq1FluxLR, eq1FluxDU; + arma_mat eq1FluxL, eq1FluxR, eq1FluxD, eq1FluxU; + + arma_mat eq2FluxLR, eq2FluxDU; + arma_mat eq2FluxL, eq2FluxR, eq2FluxD, eq2FluxU; + + arma_mat eq3FluxLR, eq3FluxDU; + arma_mat eq3FluxL, eq3FluxR, eq3FluxD, eq3FluxU; + + arma_mat eq4FluxLR, eq4FluxDU; + arma_mat eq4FluxL, eq4FluxR, eq4FluxD, eq4FluxU; + + arma_mat wsL, wsR, wsD, wsU, wsLR, wsDU; + + arma_mat diff; // for Riemann Solver + + /* Projection */ + rhoP = project_to_edges(rho, x, xEdges, y, yEdges, nGCs); + xMomentumP = project_to_edges(xMomentum, x, xEdges, y, yEdges, nGCs); + yMomentumP = project_to_edges(yMomentum, x, xEdges, y, yEdges, nGCs); + rhoEP = project_to_edges(rhoE, x, xEdges, y, yEdges, nGCs); + // Also need to project gamma and temp - these should be passed, since + // they need to be updated for the RK4 scheme: + gammaP = project_to_edges(gamma_scgc.slice(iAlt), x, xEdges, y, yEdges, nGCs); + tempP = project_to_edges(temperature_scgc.slice(iAlt), x, xEdges, y, yEdges, + nGCs); + numberDensityP = project_to_edges(density_scgc.slice(iAlt), x, xEdges, y, + yEdges, + nGCs); + + // Resolve Scalar Fields into rho, xVel, yVel, and totalE (without rho) + rhoL = rhoP.L; + rhoR = rhoP.R; + rhoD = rhoP.D; + rhoU = rhoP.U; + + xVelL = xMomentumP.L / rhoL; + xVelR = xMomentumP.R / rhoR; + xVelD = xMomentumP.D / rhoD; + xVelU = xMomentumP.U / rhoU; + + yVelL = yMomentumP.L / rhoL; + yVelR = yMomentumP.R / rhoR; + yVelD = yMomentumP.D / rhoD; + yVelU = yMomentumP.U / rhoU; + + totalEL = rhoEP.L / rhoL; + totalER = rhoEP.R / rhoR; + totalED = rhoEP.D / rhoD; + totalEU = rhoEP.U / rhoU; + + velL2 = xVelL % xVelL + yVelL % yVelL; + velR2 = xVelR % xVelR + yVelR % yVelR; + velD2 = xVelD % xVelD + yVelD % yVelD; + velU2 = xVelU % xVelU + yVelU % yVelU; + + internaleL = totalEL - 0.5 * velL2; + internaleR = totalER - 0.5 * velR2; + internaleD = totalED - 0.5 * velD2; + internaleU = totalEU - 0.5 * velU2; + + //pressureL = (gammaP.L - 1) % (rhoP.L % internaleL); + //pressureR = (gammaP.R - 1) % (rhoP.R % internaleR); + //pressureD = (gammaP.D - 1) % (rhoP.D % internaleD); + //pressureU = (gammaP.U - 1) % (rhoP.U % internaleU); + + pressureL = cKB * (numberDensityP.L % tempP.L); + pressureR = cKB * (numberDensityP.R % tempP.R); + pressureD = cKB * (numberDensityP.D % tempP.D); + pressureU = cKB * (numberDensityP.U % tempP.U); + + /* Calculate Edge Fluxes */ + // Note that dot product between normal vector at edge and flux vector + // resolves into a pure one component flux or either hat{x} or hat{y} + // Flux calculated from the left of the edge + eq1FluxL = rhoL % xVelL % grid.sqrt_g_Left.slice(iAlt); + // Flux calculated from the right of the edge + eq1FluxR = rhoR % xVelR % grid.sqrt_g_Left.slice(iAlt); + // Flux calculated from the down of the edge + eq1FluxD = rhoD % yVelD % grid.sqrt_g_Down.slice(iAlt); + // Flux calculated from the up of the edge + eq1FluxU = rhoU % yVelU % grid.sqrt_g_Down.slice(iAlt); + /* + eq2FluxL = (rhoL % xVelL % xVelL + + pressureL % grid.g11_upper_Left.slice(iAlt)) % + grid.sqrt_g_Left.slice(iAlt); + eq2FluxR = (rhoR % xVelR % xVelR + + pressureR % grid.g11_upper_Left.slice(iAlt)) % + grid.sqrt_g_Left.slice(iAlt); + eq2FluxD = (rhoD % yVelD % xVelD + + pressureD % grid.g12_upper_Down.slice(iAlt)) % + grid.sqrt_g_Down.slice(iAlt); + eq2FluxU = (rhoU % yVelU % xVelU + + pressureU % grid.g12_upper_Down.slice(iAlt)) % + grid.sqrt_g_Down.slice(iAlt); + */ + eq2FluxL = (rhoL % xVelL % xVelL + + pressureL) % + grid.sqrt_g_Left.slice(iAlt); + eq2FluxR = (rhoR % xVelR % xVelR + + pressureR) % + grid.sqrt_g_Left.slice(iAlt); + eq2FluxD = (rhoD % yVelD % xVelD + + pressureD) % + grid.sqrt_g_Down.slice(iAlt); + eq2FluxU = (rhoU % yVelU % xVelU + + pressureU) % + grid.sqrt_g_Down.slice(iAlt); + /* + eq3FluxL = (rhoL % xVelL % yVelL + + pressureL % grid.g21_upper_Left.slice(iAlt)) % + grid.sqrt_g_Left.slice(iAlt); + eq3FluxR = (rhoR % xVelR % yVelR + + pressureR % grid.g21_upper_Left.slice(iAlt)) % + grid.sqrt_g_Left.slice(iAlt); + eq3FluxD = (rhoD % yVelD % yVelD + + pressureD % grid.g22_upper_Down.slice(iAlt)) % + grid.sqrt_g_Down.slice(iAlt); + eq3FluxU = (rhoU % yVelU % yVelU + + pressureU % grid.g22_upper_Down.slice(iAlt)) % + grid.sqrt_g_Down.slice(iAlt); + */ + eq3FluxL = (rhoL % xVelL % yVelL + + pressureL) % + grid.sqrt_g_Left.slice(iAlt); + eq3FluxR = (rhoR % xVelR % yVelR + + pressureR) % + grid.sqrt_g_Left.slice(iAlt); + eq3FluxD = (rhoD % yVelD % yVelD + + pressureD) % + grid.sqrt_g_Down.slice(iAlt); + eq3FluxU = (rhoU % yVelU % yVelU + + pressureU) % + grid.sqrt_g_Down.slice(iAlt); + + eq4FluxL = (rhoEP.L + pressureL) % xVelL % grid.sqrt_g_Left.slice(iAlt); + eq4FluxR = (rhoEP.R + pressureR) % xVelR % grid.sqrt_g_Left.slice(iAlt); + eq4FluxD = (rhoEP.D + pressureD) % yVelD % grid.sqrt_g_Down.slice(iAlt); + eq4FluxU = (rhoEP.U + pressureU) % yVelU % grid.sqrt_g_Down.slice(iAlt); + + /* Wave Speed Calculation */ + wsL = sqrt(velL2) + sqrt(gammaP.L % (gammaP.L - 1.) % tempP.L); + wsR = sqrt(velR2) + sqrt(gammaP.R % (gammaP.R - 1.) % tempP.R); + wsD = sqrt(velD2) + sqrt(gammaP.D % (gammaP.D - 1.) % tempP.D); + wsU = sqrt(velU2) + sqrt(gammaP.U % (gammaP.U - 1.) % tempP.U); + + //wsL = abs(xVelL) + sqrt(gammaP.L % (gammaP.L - 1.) % internaleL); + //wsR = abs(xVelR) + sqrt(gammaP.R % (gammaP.R - 1.) % internaleR); + //wsD = abs(yVelD) + sqrt(gammaP.D % (gammaP.D - 1.) % internaleD); + //wsU = abs(yVelU) + sqrt(gammaP.U % (gammaP.U - 1.) % internaleU); + + // Find the maximum wave speed + wsLR = wsR; + + for (int i = 0; i < nXs + 1; i++) { + for (int j = 0; j < nYs; j++) { + if (wsL(i, j) > wsLR(i, j)) + wsLR(i, j) = wsL(i, j); + } + } + + wsDU = wsD; + + for (int i = 0; i < nXs; i++) { + for (int j = 0; j < nYs + 1; j++) { + if (wsU(i, j) > wsDU(i, j)) + wsDU(i, j) = wsU(i, j); + } + } + + /* Calculate average flux at the edges (Rusanov Flux) */ + /* Why is it + instead of - for the state difference? + * Because the projection actually works backwards + * Left states are actually right + * Right states are actually left + * Due to the convention in the past codes + * We keep it this way for consistency + */ + + // State difference, need to add sqrt(g) + diff = (rhoR - rhoL) % grid.sqrt_g_Left.slice(iAlt); + eq1FluxLR = (eq1FluxL + eq1FluxR) / 2 + 0.5 * wsLR % diff; + diff = (rhoU - rhoD) % grid.sqrt_g_Down.slice(iAlt); + eq1FluxDU = (eq1FluxD + eq1FluxU) / 2 + 0.5 * wsDU % diff; + + if (iAlt == -1) { + std::cout << "in solver: " << iProc << " " << + wsDU(13, 23) << " " << + wsU(13, 23) << " " << + wsD(13, 23) << " " << + gammaP.D(13, 23) << " " << + internaleD(13, 23) << " " << + gammaP.U(14, 23) << " " << + internaleU(13, 23) << " " << + diff(13, 22) << "\n"; + } + + diff = (rhoR % xVelR - rhoL % xVelL) % grid.sqrt_g_Left.slice(iAlt); + eq2FluxLR = (eq2FluxL + eq2FluxR) / 2 + 0.5 * wsLR % diff; + diff = (rhoU % xVelU - rhoD % xVelD) % grid.sqrt_g_Down.slice(iAlt); + eq2FluxDU = (eq2FluxD + eq2FluxU) / 2 + 0.5 * wsDU % diff; + + diff = (rhoR % yVelR - rhoL % yVelL) % grid.sqrt_g_Left.slice(iAlt); + eq3FluxLR = (eq3FluxL + eq3FluxR) / 2 + 0.5 * wsLR % diff; + diff = (rhoU % yVelU - rhoD % yVelD) % grid.sqrt_g_Down.slice(iAlt); + eq3FluxDU = (eq3FluxD + eq3FluxU) / 2 + 0.5 * wsDU % diff; + + diff = (rhoR % totalER - rhoL % totalEL) % grid.sqrt_g_Left.slice(iAlt); + eq4FluxLR = (eq4FluxL + eq4FluxR) / 2 + 0.5 * wsLR % diff; + diff = (rhoU % totalEU - rhoD % totalED) % grid.sqrt_g_Down.slice(iAlt); + eq4FluxDU = (eq4FluxD + eq4FluxU) / 2 + 0.5 * wsDU % diff; + + // Setup residual storage for return + arma_mat eq1_residual(nXs, nYs, fill::zeros); + arma_mat eq2_residual(nXs, nYs, fill::zeros); + arma_mat eq3_residual(nXs, nYs, fill::zeros); + arma_mat eq4_residual(nXs, nYs, fill::zeros); + + // State Update + // Note the ghost cells WILL NOT BE UPDATED + for (int j = nGCs; j < nYs - nGCs; j++) { + for (int i = nGCs; i < nXs - nGCs; i++) { + precision_t rhoResidual_ij = dx * eq1FluxLR(i + 1, j) - + dx * eq1FluxLR(i, j) + + dx * eq1FluxDU(i, j + 1) - + dx * eq1FluxDU(i, j); + eq1_residual(i, j) = -1 / area * rhoResidual_ij; + precision_t xMomentumResidual_ij = dx * eq2FluxLR(i + 1, j) - + dx * eq2FluxLR(i, j) + + dx * eq2FluxDU(i, j + 1) - + dx * eq2FluxDU(i, j); + eq2_residual(i, j) = -1 / area * xMomentumResidual_ij; + precision_t yMomentumResidual_ij = dx * eq3FluxLR(i + 1, j) - + dx * eq3FluxLR(i, j) + + dx * eq3FluxDU(i, j + 1) - + dx * eq3FluxDU(i, j); + eq3_residual(i, j) = -1 / area * yMomentumResidual_ij; + precision_t rhoEResidual_ij = dx * eq4FluxLR(i + 1, j) - + dx * eq4FluxLR(i, j) + + dx * eq4FluxDU(i, j + 1) - + dx * eq4FluxDU(i, j); + eq4_residual(i, j) = -1 / area * rhoEResidual_ij; + } + } + + if (iAlt == -1) { + std::cout << "in solver2: " << iProc << " " << + eq1_residual(13, 22) << " " << + eq2_residual(13, 22) << " " << + eq3_residual(13, 22) << " " << + eq4_residual(13, 22) << " " << + area << "\n"; + } + + + // Setup return vector + std::vector return_vector; + return_vector.push_back(eq1_residual); + return_vector.push_back(eq2_residual); + return_vector.push_back(eq3_residual); + return_vector.push_back(eq4_residual); + + return return_vector; +} + + +//-------------------------------------------------------------------- +// New solver using Rochi +//-------------------------------------------------------------------- + +void Neutrals::solver_horizontal_RK1_rochi(Grid & grid, Times & time) { + + std::string function = "Neutrals::solver_horizontal_RK1_rochi"; + static int iFunction = -1; + report.enter(function, iFunction); + + precision_t dt = time.get_dt(); + + int64_t nXs = grid.get_nX(); + int64_t nYs = grid.get_nY(); + int64_t nGCs = grid.get_nGCs(); + int64_t nAlts = grid.get_nAlts(); + + calc_concentration(); + + arma_mat temp(nXs, nYs), rho(nXs, nYs), vLon(nXs, nYs), vLat(nXs, nYs); + + arma_mat k1rho(nXs, nYs); + arma_mat k1vLon(nXs, nYs), k1vLat(nXs, nYs); + arma_mat k1temp(nXs, nYs); + + int64_t iAlt; + + for (iAlt = nGCs; iAlt < nAlts - nGCs; iAlt++) { + + /** States preprocessing **/ + /* MASS DENSITY */ + rho = rho_scgc.slice(iAlt); + vLon = velocity_vcgc[0].slice(iAlt); + vLat = velocity_vcgc[1].slice(iAlt); + temp = temperature_scgc.slice(iAlt); + + // k1 - start at t0, go to t+1/2 to figure out slope at t0 (k1) + update_states_cubesphere( + rho, vLon, vLat, temp, + k1rho, k1vLon, k1vLat, k1temp, + grid.cubeC, grid.cubeL, grid.cubeD, dt, iAlt); + // Take full step using k1: + rho_scgc.slice(iAlt) = rho - k1rho * dt; + velocity_vcgc[0].slice(iAlt) = vLon - k1vLon * dt; + velocity_vcgc[1].slice(iAlt) = vLat - k1vLat * dt; + temperature_scgc.slice(iAlt) = temp - k1temp * dt; + + } + + + calc_density_from_mass_concentration(); + + report.exit(function); + return; + +} + + +void Neutrals::solver_horizontal_RK1(Grid & grid, Times & time) { + // Function Reporting + std::string function = "Neutrals::solver_horizontal_RK1"; + static int iFunction = -1; + report.enter(function, iFunction); + + // Dimensions of Spatial Discretization + int64_t nXs = grid.get_nX(); + int64_t nYs = grid.get_nY(); + int64_t nGCs = grid.get_nGCs(); + int64_t nAlts = grid.get_nAlts(); + int iAlt, iSpec; + + // Time Discretization (TODO: change dt calculation method) + precision_t dt = time.get_dt(); + + arma_mat x(nXs, nYs), y(nXs, nYs); + arma_mat jacobian(nXs, nYs), rho(nXs, nYs), rhoE(nXs, nYs), vel2(nXs, nYs); + arma_mat uVel(nXs, nYs), vVel(nXs, nYs), xVel(nXs, nYs), yVel(nXs, nYs); + arma_mat xMomentum(nXs, nYs), yMomentum(nXs, nYs); + arma_mat xMomentum_0(nXs, nYs), yMomentum_0(nXs, nYs); + arma_mat rho_0(nXs, nYs), rhoE_0(nXs, nYs); + arma_mat f_0_eq1(nXs, nYs), f_0_eq2(nXs, nYs); + arma_mat f_0_eq3(nXs, nYs), f_0_eq4(nXs, nYs); + + calc_concentration(); + + // Advance for bulk calculation first, calculate for every altitude + + for (iAlt = nGCs; iAlt < nAlts - nGCs; iAlt++) { + /** Extract Grid Features **/ + x = grid.refx_scgc.slice(iAlt); + arma_mat xEdges = grid.refx_Left.slice(iAlt); + y = grid.refy_scgc.slice(iAlt); + arma_mat yEdges = grid.refy_Down.slice(iAlt); + + // Get reference grid dimensions (Assume dx = dy and equidistant) + arma_vec x_vec = x.col(0); + precision_t dx = x_vec(1) - x_vec(0); + precision_t area = dx * dx; + jacobian = grid.sqrt_g_scgc.slice(iAlt); + + /** States preprocessing **/ + /* MASS DENSITY */ + rho = rho_scgc.slice(iAlt); + + /* VELOCITY */ + // Get spherical velocity + uVel = velocity_vcgc[0].slice(iAlt); + vVel = velocity_vcgc[1].slice(iAlt); + // Convert to contravariant (reference) velocity + xVel = uVel % grid.A11_inv_scgc.slice(iAlt) + vVel % + grid.A12_inv_scgc.slice(iAlt); // u^1 + yVel = uVel % grid.A21_inv_scgc.slice(iAlt) + vVel % + grid.A22_inv_scgc.slice(iAlt); // u^2 + vel2 = xVel % xVel + yVel % yVel; + // Generate contravriant momentum (no sqrt(g)) + xMomentum = rho % xVel; // x1momentum + yMomentum = rho % yVel; // x2momentum + + /* TEMP and ENERGY */ + // Generate total energy (rhoE (no sqrt(g))) + // (TODO: Verify units) + rhoE = rho % (temperature_scgc.slice(iAlt) % Cv_scgc.slice( + iAlt) + 0.5 * vel2); + + + if (iAlt == -1) { + std::cout << "before solve: " << iProc << " " << temperature_scgc(13, 22, + 2) << " " << + rhoE(13, 22) << " " << xVel(13, 22) << " " << yVel(13, 22) << " " << + rho_scgc(13, 22, 2) << "\n"; + } + + if (iAlt == -2) { + std::cout << "before solve: " << + temperature_scgc(13, 22, 2) << " " << + xVel(13, 22) << " " << yVel(13, 22) << " " << + rho_scgc(13, 22, 2) << "\n"; + } + + + /** Advancing with RK4 **/ + // Setup Containers + rho_0 = rho; + xMomentum_0 = xMomentum; + yMomentum_0 = yMomentum; + rhoE_0 = rhoE; + + // FIRST (1) STEP, Compute F_0-> State_1 + // Pass in state vector + std::vector state_0; + state_0.push_back(rho_0); + state_0.push_back(xMomentum_0); + state_0.push_back(yMomentum_0); + state_0.push_back(rhoE_0); + std::vector f_0_vec = residual_horizontal_rusanov(state_0, grid, time, + iAlt); + // Extract Gradients + f_0_eq1 = f_0_vec[0]; + f_0_eq2 = f_0_vec[1]; + f_0_eq3 = f_0_vec[2]; + f_0_eq4 = f_0_vec[3]; + + /* Update Bulk Scalars and Contravariant velocity */ + // Euler State Update + for (int j = nGCs; j < nYs - nGCs; j++) { + for (int i = nGCs; i < nXs - nGCs; i++) { + rho(i, j) = rho_0(i, j) + dt * f_0_eq1(i, j) / jacobian(i, j); + xMomentum(i, j) = xMomentum_0(i, j) - dt * f_0_eq2(i, j) / jacobian(i, j); + yMomentum(i, j) = yMomentum_0(i, j) - dt * f_0_eq3(i, j) / jacobian(i, j); + rhoE(i, j) = rhoE_0(i, j) + dt * f_0_eq4(i, j) / jacobian(i, j); + } + } + + if (iAlt == -1) { + std::cout << "after solve: " << iProc << " " << + dt << " " << + f_0_eq1(13, 22) << " " << + f_0_eq2(13, 22) << " " << + f_0_eq3(13, 22) << " " << + f_0_eq4(13, 22) << " " << + rhoE(13, 22) << " " << xMomentum(13, 22) << " " << yMomentum(13, 22) << " " << + rho_scgc(13, 22, 2) << "\n"; + } + + + /* Re-derive Spherical Velocity and Bulk States */ + // Density + rho_scgc.slice(iAlt) = rho; + + // Bulk Velocity + xVel = xMomentum / rho; // u^1 + yVel = yMomentum / rho; // u^2 + vel2 = xVel % xVel + yVel % yVel; // Squared Magnitude of Contravariant + + if (iProc == 0 && iAlt == -2) { + std::cout << "a11 : " << grid.A11_scgc(13, 20, 2) << " " + << grid.A12_scgc(13, 20, 2) << "\n"; + std::cout << "inv : " << grid.A11_inv_scgc(13, 20, 2) << " " + << grid.A12_inv_scgc(13, 20, 2) << "\n"; + } + + velocity_vcgc[0].slice(iAlt) = xVel % grid.A11_scgc.slice( + iAlt) + yVel % grid.A12_scgc.slice(iAlt); + velocity_vcgc[1].slice(iAlt) = + xVel % grid.A21_scgc.slice(iAlt) + + yVel % grid.A22_scgc.slice(iAlt); + + /* Update temperature */ + temperature_scgc.slice(iAlt) = (rhoE / rho - 0.5 * vel2) / Cv_scgc.slice(iAlt); + + + //if (iAlt == 10) { + std::cout << "after solve: " << iAlt << " " << + temperature_scgc(13, 22, 10) << " " << + velocity_vcgc[0](13, 22, 10) << " " << velocity_vcgc[1](13, 22, 10) << " " << + rho_scgc(13, 22, 10) << "\n"; + //} + + calc_density_from_mass_concentration(); + //assign_bulk_velocity(); + + } + + report.exit(function); + return; +} + +/* +void Neutrals::solver_horizontal_RK1(Grid& grid, Times& time) { + // Function Reporting + std::string function = "Neutrals::solver_horizontal_RK1"; + static int iFunction = -1; + report.enter(function, iFunction); + + // Dimensions of Spatial Discretization + int64_t nXs = grid.get_nX(); + int64_t nYs = grid.get_nY(); + int64_t nGCs = grid.get_nGCs(); + int64_t nAlts = grid.get_nAlts(); + int iAlt, iSpec; + + // Time Discretization (TODO: change dt calculation method) + precision_t dt = time.get_dt(); + + iAlt = 0; + + arma_mat x(nXs, nYs), y(nXs, nYs); + arma_mat jacobian(nXs, nYs), rho(nXs, nYs), rhoE(nXs, nYs), vel2(nXs, nYs); + arma_mat uVel(nXs, nYs), vVel(nXs, nYs), xVel(nXs, nYs), yVel(nXs, nYs); + arma_mat xMomentum(nXs, nYs), yMomentum(nXs, nYs); + arma_mat xMomentum_0(nXs, nYs), yMomentum_0(nXs, nYs); + arma_mat rho_0(nXs, nYs), rhoE_0(nXs, nYs); + arma_mat f_0_eq1(nXs, nYs), f_0_eq2(nXs, nYs); + arma_mat f_0_eq3(nXs, nYs), f_0_eq4(nXs, nYs); + +// Advance for bulk calculation first, calculate for every altitude +for (iAlt = nGCs; iAlt < nAlts - nGCs; iAlt++) { + // Extract Grid Features +arma_mat x = grid.refx_scgc.slice(iAlt); +arma_mat xEdges = grid.refx_Left.slice(iAlt); +arma_mat y = grid.refy_scgc.slice(iAlt); +arma_mat yEdges = grid.refy_Down.slice(iAlt); + +// Get reference grid dimensions (Assume dx = dy and equidistant) +arma_vec x_vec = x.col(0); +precision_t dx = x_vec(1) - x_vec(0); +precision_t area = dx * dx; +arma_mat jacobian = grid.sqrt_g_scgc.slice(iAlt); + +// States preprocessing +// MASS DENSITY +arma_mat rho = rho_scgc.slice(iAlt); + +// VELOCITY +// Get spherical velocity +arma_mat uVel = velocity_vcgc[0].slice(iAlt); +arma_mat vVel = velocity_vcgc[1].slice(iAlt); +// Convert to contravariant (reference) velocity +arma_mat xVel = uVel % grid.A11_inv_scgc.slice(iAlt) + vVel % + grid.A12_inv_scgc.slice(iAlt); // u^1 +arma_mat yVel = uVel % grid.A21_inv_scgc.slice(iAlt) + vVel % + grid.A22_inv_scgc.slice(iAlt); // u^2 +arma_mat vel2 = xVel % xVel + yVel % yVel; +// Generate contravriant momentum (no sqrt(g)) +arma_mat xMomentum = rho % xVel; // x1momentum +arma_mat yMomentum = rho % yVel; // x2momentum + +// TEMP and ENERGY +// Generate total energy (rhoE (no sqrt(g))) +// (TODO: Verify units) +arma_mat rhoE = rho % (temperature_scgc.slice(iAlt) % + Cv_scgc.slice(iAlt) + 0.5 * vel2); + +// Advancing with RK4 +// Setup Containers +arma_mat rho_0 = rho; +arma_mat xMomentum_0 = xMomentum; +arma_mat yMomentum_0 = yMomentum; +arma_mat rhoE_0 = rhoE; + +// FIRST (1) STEP, Compute F_0-> State_1 +// Pass in state vector +std::vector state_0; +state_0.push_back(rho_0); +state_0.push_back(xMomentum_0); +state_0.push_back(yMomentum_0); +state_0.push_back(rhoE_0); +std::vector f_0_vec = residual_horizontal_rusanov(state_0, grid, time, + iAlt); +// Extract Gradients +arma_mat f_0_eq1 = f_0_vec[0]; +arma_mat f_0_eq2 = f_0_vec[1]; +arma_mat f_0_eq3 = f_0_vec[2]; +arma_mat f_0_eq4 = f_0_vec[3]; + +// Update Bulk Scalars and Contravariant velocity +// Euler State Update +for (int j = nGCs; j < nYs - nGCs; j++) { + for (int i = nGCs; i < nXs - nGCs; i++) { + rho(i, j) = rho_0(i, j) + dt * f_0_eq1(i, j) / jacobian(i, j); + xMomentum(i, j) = xMomentum_0(i, j) - dt * f_0_eq2(i, j) / jacobian(i, j); + yMomentum(i, j) = yMomentum_0(i, j) - dt * f_0_eq3(i, j) / jacobian(i, j); + rhoE(i, j) = rhoE_0(i, j) + dt * f_0_eq4(i, j) / jacobian(i, j); + } +} + +// Re-derive Spherical Velocity and Bulk States +// Density +rho_scgc.slice(iAlt) = rho; + +// Bulk Velocity +xVel = xMomentum / rho; // u^1 +yVel = yMomentum / rho; // u^2 +vel2 = xVel % xVel + yVel % yVel; // Squared Magnitude of Contravariant +velocity_vcgc[0].slice(iAlt) = xVel % grid.A11_scgc.slice( + iAlt) + yVel % grid.A12_scgc.slice(iAlt); +velocity_vcgc[1].slice(iAlt) = xVel % grid.A21_scgc.slice( + iAlt) + yVel % grid.A22_scgc.slice(iAlt); + +// Update specie number density and velocity +for (iSpec = 0; iSpec < nSpecies; iSpec++) { + //species[iSpec].density_scgc.slice(iAlt) = rho % species[iSpec].mass_concentration_scgc.slice(iAlt); + species[iSpec].density_scgc.slice(iAlt) = rho / species[iSpec].mass; + species[iSpec].velocity_vcgc[0].slice(iAlt) = velocity_vcgc[0].slice(iAlt); + species[iSpec].velocity_vcgc[1].slice(iAlt) = velocity_vcgc[1].slice(iAlt); +} + +// Update temperature +temperature_scgc.slice(iAlt) = (rhoE / rho - 0.5 * vel2) / Cv_scgc.slice(iAlt); + +report.exit(function); +return; +} + +*/ + +void Neutrals::solver_horizontal_RK4(Grid & grid, Times & time) { + // Function Reporting + std::string function = "Neutrals::solver_horizontal_RK4"; + static int iFunction = -1; + report.enter(function, iFunction); + + // Dimensions of Spatial Discretization + int64_t nXs = grid.get_nX(); + int64_t nYs = grid.get_nY(); + int64_t nGCs = grid.get_nGCs(); + int64_t nAlts = grid.get_nAlts(); + int iAlt, iSpec; + + // Time Discretization (TODO: change dt calculation method) + precision_t dt = time.get_dt() / 10; + + // Advance for bulk calculation first, calculate for every altitude + for (iAlt = nGCs; iAlt < nAlts - nGCs; iAlt++) { + /** Extract Grid Features **/ + arma_mat x = grid.refx_scgc.slice(iAlt); + arma_mat xEdges = grid.refx_Left.slice(iAlt); + arma_mat y = grid.refy_scgc.slice(iAlt); + arma_mat yEdges = grid.refy_Down.slice(iAlt); + + // Get reference grid dimensions (Assume dx = dy and equidistant) + arma_vec x_vec = x.col(0); + precision_t dx = x_vec(1) - x_vec(0); + precision_t area = dx * dx; + arma_mat jacobian = grid.sqrt_g_scgc.slice(iAlt); + + /** States preprocessing **/ + /* MASS DENSITY */ + arma_mat rho = rho_scgc.slice(iAlt); + + /* VELOCITY */ + // Get spherical velocity + arma_mat uVel = velocity_vcgc[0].slice(iAlt); + arma_mat vVel = velocity_vcgc[1].slice(iAlt); + // Convert to contravariant (reference) velocity + arma_mat xVel = uVel % grid.A11_inv_scgc.slice(iAlt) + vVel % + grid.A12_inv_scgc.slice(iAlt); // u^1 + arma_mat yVel = uVel % grid.A21_inv_scgc.slice(iAlt) + vVel % + grid.A22_inv_scgc.slice(iAlt); // u^2 + arma_mat vel2 = xVel % xVel + yVel % yVel; + // Generate contravriant momentum (no sqrt(g)) + arma_mat xMomentum = rho % xVel; // x1momentum + arma_mat yMomentum = rho % yVel; // x2momentum + + /* TEMP and ENERGY */ + // Generate total energy (rhoE (no sqrt(g))) + // (TODO: Verify units) + //arma_mat rhoE = rho % (temperature_scgc.slice(iAlt) % Cv_scgc.slice( + // iAlt) + 0.5 * vel2); + arma_mat rhoE = rho % (temperature_scgc.slice(iAlt) % + Cv_scgc.slice(iAlt) + + 0.5 * vel2); + + /** Advancing with RK4 **/ + // Setup Containers + arma_mat rho_0 = rho; + arma_mat rho_1(nXs, nYs, fill::zeros); // corresponding f_1 + arma_mat rho_2(nXs, nYs, fill::zeros); // corresponding f_2 + arma_mat rho_3(nXs, nYs, fill::zeros); // corresponding f_3 + + arma_mat xMomentum_0 = xMomentum; + arma_mat xMomentum_1(nXs, nYs, fill::zeros); // corresponding f_1 + arma_mat xMomentum_2(nXs, nYs, fill::zeros); // corresponding f_2 + arma_mat xMomentum_3(nXs, nYs, fill::zeros); // corresponding f_3 + + arma_mat yMomentum_0 = yMomentum; + arma_mat yMomentum_1(nXs, nYs, fill::zeros); // corresponding f_1 + arma_mat yMomentum_2(nXs, nYs, fill::zeros); // corresponding f_2 + arma_mat yMomentum_3(nXs, nYs, fill::zeros); // corresponding f_3 + + arma_mat rhoE_0 = rhoE; + arma_mat rhoE_1(nXs, nYs, fill::zeros); // corresponding f_1 + arma_mat rhoE_2(nXs, nYs, fill::zeros); // corresponding f_2 + arma_mat rhoE_3(nXs, nYs, fill::zeros); // corresponding f_3 + + // FIRST (1) STEP, Compute F_0-> State_1 + // Pass in state vector + std::vector state_0; + state_0.push_back(rho_0); + state_0.push_back(xMomentum_0); + state_0.push_back(yMomentum_0); + state_0.push_back(rhoE_0); + std::vector f_0_vec = residual_horizontal_rusanov(state_0, grid, time, + iAlt); + // Extract Gradients + arma_mat f_0_eq1 = f_0_vec[0]; + arma_mat f_0_eq2 = f_0_vec[1]; + arma_mat f_0_eq3 = f_0_vec[2]; + arma_mat f_0_eq4 = f_0_vec[3]; + + /* Update Bulk Scalars and Contravariant velocity */ + // Euler State Update + for (int j = nGCs; j < nYs - nGCs; j++) { + for (int i = nGCs; i < nXs - nGCs; i++) { + rho_1(i, j) = rho_0(i, j) + 0.5 * dt * f_0_eq1(i, j) / jacobian(i, j); + xMomentum_1(i, j) = xMomentum_0(i, j) + 0.5 * dt * f_0_eq2(i, j) / jacobian(i, + j); + yMomentum_1(i, j) = yMomentum_0(i, j) + 0.5 * dt * f_0_eq3(i, j) / jacobian(i, + j); + rhoE_1(i, j) = rhoE_0(i, j) + 0.5 * dt * f_0_eq4(i, j) / jacobian(i, j); + } + } + + // SECOND (2) STEP, Compute F_1-> State_2 + // Pass in state vector + std::vector state_1; + state_1.push_back(rho_1); + state_1.push_back(xMomentum_1); + state_1.push_back(yMomentum_1); + state_1.push_back(rhoE_1); + std::vector f_1_vec = residual_horizontal_rusanov(state_1, grid, time, + iAlt); + // Extract Gradients + arma_mat f_1_eq1 = f_1_vec[0]; + arma_mat f_1_eq2 = f_1_vec[1]; + arma_mat f_1_eq3 = f_1_vec[2]; + arma_mat f_1_eq4 = f_1_vec[3]; + + /* Update Bulk Scalars and Contravariant velocity */ + // Euler State Update + for (int j = nGCs; j < nYs - nGCs; j++) { + for (int i = nGCs; i < nXs - nGCs; i++) { + rho_2(i, j) = rho_0(i, j) + 0.5 * dt * f_1_eq1(i, j) / jacobian(i, j); + xMomentum_2(i, j) = xMomentum_0(i, j) + 0.5 * dt * f_1_eq2(i, j) / jacobian(i, + j); + yMomentum_2(i, j) = yMomentum_0(i, j) + 0.5 * dt * f_1_eq3(i, j) / jacobian(i, + j); + rhoE_2(i, j) = rhoE_0(i, j) + 0.5 * dt * f_1_eq4(i, j) / jacobian(i, j); + } + } + + // THIRD (3) STEP, Compute F_2-> State_3 + // Pass in state vector + std::vector state_2; + state_2.push_back(rho_2); + state_2.push_back(xMomentum_2); + state_2.push_back(yMomentum_2); + state_2.push_back(rhoE_2); + std::vector f_2_vec = residual_horizontal_rusanov(state_2, grid, time, + iAlt); + // Extract Gradients + arma_mat f_2_eq1 = f_2_vec[0]; + arma_mat f_2_eq2 = f_2_vec[1]; + arma_mat f_2_eq3 = f_2_vec[2]; + arma_mat f_2_eq4 = f_2_vec[3]; + + /* Update Bulk Scalars and Contravariant velocity */ + // Euler State Update + for (int j = nGCs; j < nYs - nGCs; j++) { + for (int i = nGCs; i < nXs - nGCs; i++) { + rho_3(i, j) = rho_0(i, j) + dt * f_2_eq1(i, j) / jacobian(i, j); + xMomentum_3(i, j) = xMomentum_0(i, j) + dt * f_2_eq2(i, j) / jacobian(i, j); + yMomentum_3(i, j) = yMomentum_0(i, j) + dt * f_2_eq3(i, j) / jacobian(i, j); + rhoE_3(i, j) = rhoE_0(i, j) + dt * f_2_eq4(i, j) / jacobian(i, j); + } + } + + // FOURTH (4) STEP, Compute F_3 + // Pass in state vector + std::vector state_3; + state_3.push_back(rho_3); + state_3.push_back(xMomentum_3); + state_3.push_back(yMomentum_3); + state_3.push_back(rhoE_3); + std::vector f_3_vec = residual_horizontal_rusanov(state_3, grid, time, + iAlt); + // Extract Gradients + arma_mat f_3_eq1 = f_3_vec[0]; + arma_mat f_3_eq2 = f_3_vec[1]; + arma_mat f_3_eq3 = f_3_vec[2]; + arma_mat f_3_eq4 = f_3_vec[3]; + + // Summing all steps for final update + arma_mat f_sum_eq1 = f_0_eq1 + 2 * f_1_eq1 + 2 * f_2_eq1 + f_3_eq1; + arma_mat f_sum_eq2 = f_0_eq2 + 2 * f_1_eq2 + 2 * f_2_eq2 + f_3_eq2; + arma_mat f_sum_eq3 = f_0_eq3 + 2 * f_1_eq3 + 2 * f_2_eq3 + f_3_eq3; + arma_mat f_sum_eq4 = f_0_eq4 + 2 * f_1_eq4 + 2 * f_2_eq4 + f_3_eq4; + + /* Update Bulk Scalars and Contravariant velocity */ + // Euler State Update + for (int j = nGCs; j < nYs - nGCs; j++) { + for (int i = nGCs; i < nXs - nGCs; i++) { + rho(i, j) = rho(i, j) + dt / 6 * f_sum_eq1(i, j) / jacobian(i, j); + xMomentum(i, j) = xMomentum(i, j) + dt / 6 * f_sum_eq2(i, j) / jacobian(i, j); + yMomentum(i, j) = yMomentum(i, j) + dt / 6 * f_sum_eq3(i, j) / jacobian(i, j); + rhoE(i, j) = rhoE(i, j) + dt / 6 * f_sum_eq4(i, j) / jacobian(i, j); + } + } + + /* Re-derive Spherical Velocity and Bulk States */ + // Density + rho_scgc.slice(iAlt) = rho; + + // Bulk Velocity + xVel = xMomentum / rho; // u^1 + yVel = yMomentum / rho; // u^2 + vel2 = xVel % xVel + yVel % yVel; // Squared Magnitude of Contravariant + velocity_vcgc[0].slice(iAlt) = xVel % grid.A11_scgc.slice( + iAlt) + yVel % grid.A12_scgc.slice(iAlt); + velocity_vcgc[1].slice(iAlt) = xVel % grid.A21_scgc.slice( + iAlt) + yVel % grid.A22_scgc.slice(iAlt); + + /* Update specie number density and velocity */ + for (iSpec = 0; iSpec < nSpecies; iSpec++) { + //species[iSpec].density_scgc.slice(iAlt) = rho % species[iSpec].concentration_scgc.slice(iAlt); + species[iSpec].density_scgc.slice(iAlt) = rho / species[iSpec].mass; + species[iSpec].velocity_vcgc[0].slice(iAlt) = velocity_vcgc[0].slice(iAlt); + species[iSpec].velocity_vcgc[1].slice(iAlt) = velocity_vcgc[1].slice(iAlt); + } + + /* Update temperature */ + temperature_scgc.slice(iAlt) = (rhoE / rho - 0.5 * vel2) / Cv_scgc.slice(iAlt); + + report.exit(function); + return; + } +} diff --git a/src/solver_horizontal_cubesphere_advection.cpp b/src/solver_horizontal_cubesphere_advection.cpp new file mode 100644 index 00000000..e96b1e88 --- /dev/null +++ b/src/solver_horizontal_cubesphere_advection.cpp @@ -0,0 +1,993 @@ +// Copyright 2023, the Aether Development Team (see doc/dev_team.md for members) +// Full license can be found in License.md + +// Initial version: F. Cheng, July 2023 + +#include "aether.h" + +using namespace Cubesphere_tools; + +// DOES NOT WORK WELL +std::vector Neutrals::residual_horizontal_hlle_advection( + std::vector& states, Grid& grid, Times& time) { + // Dimensions of Spatial Discretization + int64_t nXs = grid.get_nX(); + int64_t nYs = grid.get_nY(); + int64_t nGCs = grid.get_nGCs(); + int64_t nAlts = grid.get_nAlts(); + int iAlt, iSpec; + + /** Extract Grid Features **/ + arma_mat x = grid.refx_scgc.slice(iAlt); + arma_mat xEdges = grid.refx_Left.slice(iAlt); + arma_mat y = grid.refy_scgc.slice(iAlt); + arma_mat yEdges = grid.refy_Down.slice(iAlt); + + // Get reference grid dimensions (Assume dx = dy and equidistant) + arma_vec x_vec = x.col(0); + precision_t dx = x_vec(1) - x_vec(0); + precision_t area = dx * dx; + arma_mat jacobian = grid.sqrt_g_scgc.slice(iAlt); + + /** State/Velocity extraction **/ + /* MASS DENSITY */ + arma_mat rho = states[0]; + + /* VELOCITY */ + // Convert to contravariant (reference) velocity + arma_mat xVel = states[1]; // u^1 + arma_mat yVel = states[2]; // u^2 + + // Generate velocity magnitude squared + arma_mat vel2 = xVel % xVel + yVel % yVel; + + /** Advancing **/ + /* Initialize projection constructs and storages */ + projection_struct rhoP; + projection_struct xVelP; + projection_struct yVelP; + + // They are all pure scalar fields without sqrt(g) + arma_mat rhoL, rhoR, rhoD, rhoU; + arma_mat xVelL, xVelR, xVelD, xVelU; + arma_mat yVelL, yVelR, yVelD, yVelU; + + arma_mat velL2, velR2, velD2, velU2; + + /** Initialize Flux and Wave Speed Storages */ + arma_mat eq1FluxLR_left, eq1FluxDU_down; + arma_mat eq1FluxLR_right, eq1FluxDU_upper; + arma_mat eq1FluxL, eq1FluxR, eq1FluxD, eq1FluxU; + + arma_mat wsL, wsR, wsD, wsU; + arma_mat wsL_min, wsL_max, wsR_min, wsR_max; + arma_mat wsD_min, wsD_max, wsU_min, wsU_max; + arma_mat wsLR_max, wsDU_max, wsLR_min, wsDU_min; + + arma_mat diff; // for Riemann Solver + + /* Projection */ + rhoP = project_to_edges(rho, x, xEdges, y, yEdges, nGCs); + xVelP = project_to_edges(xVel, x, xEdges, y, yEdges, nGCs); + yVelP = project_to_edges(yVel, x, xEdges, y, yEdges, nGCs); + + // Resolve Scalar Fields into rho, xVel, yVel, and totalE (without rho) + rhoL = rhoP.L; + rhoR = rhoP.R; + rhoD = rhoP.D; + rhoU = rhoP.U; + + xVelL = xVelP.L; + xVelR = xVelP.R; + xVelD = xVelP.D; + xVelU = xVelP.U; + + yVelL = yVelP.L; + yVelR = yVelP.R; + yVelD = yVelP.D; + yVelU = yVelP.U; + + //velL2 = xVelL % xVelL + yVelL % yVelL; + //velR2 = xVelR % xVelR + yVelR % yVelR; + //velD2 = xVelD % xVelD + yVelD % yVelD; + //velU2 = xVelU % xVelU + yVelU % yVelU; + + /* Calculate Edge Fluxes */ + // Note that dot product between normal vector at edge and flux vector + // resolves into a pure one component flux or either hat{x} or hat{y} + + // Flux calculated from the left of the edge + eq1FluxL = rhoL % xVelL % grid.sqrt_g_Left.slice(iAlt); + // Flux calculated from the right of the edge + eq1FluxR = rhoR % xVelR % grid.sqrt_g_Left.slice(iAlt); + // Flux calculated from the down of the edge + eq1FluxD = rhoD % yVelD % grid.sqrt_g_Down.slice(iAlt); + // Flux calculated from the up of the edge + eq1FluxU = rhoU % yVelU % grid.sqrt_g_Down.slice(iAlt); + + /* Wave Speed Calculation (Left/Down) */ + wsL = xVelL; + wsR = xVelR; + wsD = yVelD; + wsU = yVelU; + + wsL_max = wsL; + wsL_min = wsL; + wsR_max = wsR; + wsR_min = wsR; + wsD_max = wsD; + wsD_min = wsD; + wsU_max = wsU; + wsU_min = wsU; + + // Process wave speeds from each direction first + for (int i = 0; i < nXs + 1; i++) { + for (int j = 0; j < nYs; j++) { + + if (wsL(i, j) > 0.) + wsL_min(i, j) = 0.; + + else + wsL_max(i, j) = 0.; + + if (wsR(i, j) > 0.) + wsR_min(i, j) = 0.; + + else + wsR_max(i, j) = 0.; + } + } + + for (int i = 0; i < nXs; i++) { + for (int j = 0; j < nYs + 1; j++) { + if (wsD(i, j) > 0.) + wsD_min(i, j) = 0.; + + else + wsD_max(i, j) = 0.; + + if (wsU(i, j) > 0.) + wsU_min(i, j) = 0.; + + else + wsU_max(i, j) = 0.; + } + } + + // Process edge wave speeds + wsLR_max = wsR_max; + + for (int i = 0; i < nXs + 1; i++) { + for (int j = 0; j < nYs; j++) { + if (wsL_max(i, j) > wsLR_max(i, j)) + wsLR_max(i, j) = wsL_max(i, j); + } + } + + wsDU_max = wsD_max; + + for (int i = 0; i < nXs; i++) { + for (int j = 0; j < nYs + 1; j++) { + if (wsU_max(i, j) > wsDU_max(i, j)) + wsDU_max(i, j) = wsU_max(i, j); + } + } + + wsLR_min = wsR_min; + + for (int i = 0; i < nXs + 1; i++) { + for (int j = 0; j < nYs; j++) { + if (wsL_min(i, j) < wsLR_min(i, j)) + wsLR_min(i, j) = wsL_min(i, j); + } + } + + wsDU_min = wsD_min; + + for (int i = 0; i < nXs; i++) { + for (int j = 0; j < nYs + 1; j++) { + if (wsU_min(i, j) < wsDU_min(i, j)) + wsDU_min(i, j) = wsU_min(i, j); + } + } + + /* Calculate average flux at the edges (HLLE Flux) */ + arma_mat wsLR_sum = wsLR_max + wsLR_min; + arma_mat wsLR_diff = wsLR_max - wsLR_min; + diff = (rhoR - rhoL) % grid.sqrt_g_Left.slice( + iAlt); // State difference, need to add sqrt(g) + eq1FluxLR_left = 0.5 * (eq1FluxL + eq1FluxR) + 0.5 * (wsLR_sum / wsLR_diff) % + (eq1FluxR - eq1FluxL) - (wsLR_max % wsLR_min) / wsLR_diff % diff; + + arma_mat wsDU_sum = wsDU_max + wsDU_min; + arma_mat wsDU_diff = wsDU_max - wsDU_min; + diff = (rhoU - rhoD) % grid.sqrt_g_Down.slice(iAlt); + eq1FluxDU_down = 0.5 * (eq1FluxU + eq1FluxD) + 0.5 * (wsDU_sum / wsDU_diff) % + (eq1FluxD - eq1FluxU) - (wsDU_max % wsDU_min) / wsDU_diff % diff; + + /* Wave Speed Calculation (Right/Up) */ + wsL = -xVelL; + wsR = -xVelR; + wsD = -yVelD; + wsU = -yVelU; + + wsL_max = wsL; + wsL_min = wsL; + wsR_max = wsR; + wsR_min = wsR; + wsD_max = wsD; + wsD_min = wsD; + wsU_max = wsU; + wsU_min = wsU; + + // Process wave speeds from each direction first + for (int i = 0; i < nXs + 1; i++) { + for (int j = 0; j < nYs; j++) { + + if (wsL(i, j) > 0.) + wsL_min(i, j) = 0.; + + else + wsL_max(i, j) = 0.; + + if (wsR(i, j) > 0.) + wsR_min(i, j) = 0.; + + else + wsR_max(i, j) = 0.; + } + } + + for (int i = 0; i < nXs; i++) { + for (int j = 0; j < nYs + 1; j++) { + if (wsD(i, j) > 0.) + wsD_min(i, j) = 0.; + + else + wsD_max(i, j) = 0.; + + if (wsU(i, j) > 0.) + wsU_min(i, j) = 0.; + + else + wsU_max(i, j) = 0.; + } + } + + // Process edge wave speeds + wsLR_max = wsR_max; + + for (int i = 0; i < nXs + 1; i++) { + for (int j = 0; j < nYs; j++) { + if (wsL_max(i, j) > wsLR_max(i, j)) + wsLR_max(i, j) = wsL_max(i, j); + } + } + + wsDU_max = wsD_max; + + for (int i = 0; i < nXs; i++) { + for (int j = 0; j < nYs + 1; j++) { + if (wsU_max(i, j) > wsDU_max(i, j)) + wsDU_max(i, j) = wsU_max(i, j); + } + } + + wsLR_min = wsR_min; + + for (int i = 0; i < nXs + 1; i++) { + for (int j = 0; j < nYs; j++) { + if (wsL_min(i, j) < wsLR_min(i, j)) + wsLR_min(i, j) = wsL_min(i, j); + } + } + + wsDU_min = wsD_min; + + for (int i = 0; i < nXs; i++) { + for (int j = 0; j < nYs + 1; j++) { + if (wsU_min(i, j) < wsDU_min(i, j)) + wsDU_min(i, j) = wsU_min(i, j); + } + } + + /* Calculate average flux at the edges (HLLE Flux) */ + wsLR_sum = wsLR_max + wsLR_min; + wsLR_diff = wsLR_max - wsLR_min; + diff = (rhoR - rhoL) % grid.sqrt_g_Left.slice( + iAlt); // State difference, need to add sqrt(g) + eq1FluxLR_right = 0.5 * (eq1FluxL + eq1FluxR) + 0.5 * (wsLR_sum / wsLR_diff) % + (eq1FluxR - eq1FluxL) - (wsLR_max % wsLR_min) / wsLR_diff % diff; + + wsDU_sum = wsDU_max + wsDU_min; + wsDU_diff = wsDU_max - wsDU_min; + diff = (rhoU - rhoD) % grid.sqrt_g_Down.slice(iAlt); + eq1FluxDU_upper = 0.5 * (eq1FluxU + eq1FluxD) + 0.5 * (wsDU_sum / wsDU_diff) % + (eq1FluxD - eq1FluxU) - (wsDU_max % wsDU_min) / wsDU_diff % diff; + + + // Setup residual storage for return + arma_mat eq1_residual(nXs, nYs, fill::zeros); + + // State Update + // Note the ghost cells WILL NOT BE UPDATED + for (int j = nGCs; j < nYs - nGCs; j++) { + for (int i = nGCs; i < nXs - nGCs; i++) { + precision_t rhoResidual_ij = dx * eq1FluxLR_right(i + 1, j) - + dx * eq1FluxLR_left(i, j) + + dx * eq1FluxDU_upper(i, j + 1) - + dx * eq1FluxDU_down(i, j); + eq1_residual(i, j) = -1 / area * rhoResidual_ij; + } + } + + // Setup return vector + std::vector return_vector; + return_vector.push_back(eq1_residual); + + return return_vector; +} + +// WORKS, but diffusive +std::vector Neutrals::residual_horizontal_rusanov_advection( + std::vector& states, Grid& grid, Times& time) { + // Dimensions of Spatial Discretization + int64_t nXs = grid.get_nX(); + int64_t nYs = grid.get_nY(); + int64_t nGCs = grid.get_nGCs(); + int64_t nAlts = grid.get_nAlts(); + int iAlt, iSpec; + + /** Extract Grid Features **/ + arma_mat x = grid.refx_scgc.slice(iAlt); + arma_mat xEdges = grid.refx_Left.slice(iAlt); + arma_mat y = grid.refy_scgc.slice(iAlt); + arma_mat yEdges = grid.refy_Down.slice(iAlt); + + // Get reference grid dimensions (Assume dx = dy and equidistant) + arma_vec x_vec = x.col(0); + precision_t dx = x_vec(1) - x_vec(0); + precision_t area = dx * dx; + arma_mat jacobian = grid.sqrt_g_scgc.slice(iAlt); + + /** State/Velocity extraction **/ + /* MASS DENSITY */ + arma_mat rho = states[0]; + + /* VELOCITY */ + // Convert to contravariant (reference) velocity + arma_mat xVel = states[1]; // u^1 + arma_mat yVel = states[2]; // u^2 + + // Generate velocity magnitude squared + arma_mat vel2 = xVel % xVel + yVel % yVel; + + /** Advancing **/ + /* Initialize projection constructs and storages */ + projection_struct rhoP; + projection_struct xVelP; + projection_struct yVelP; + + // They are all pure scalar fields without sqrt(g) + arma_mat rhoL, rhoR, rhoD, rhoU; + arma_mat xVelL, xVelR, xVelD, xVelU; + arma_mat yVelL, yVelR, yVelD, yVelU; + + arma_mat velL2, velR2, velD2, velU2; + + /** Initialize Flux and Wave Speed Storages */ + arma_mat eq1FluxLR, eq1FluxDU; + arma_mat eq1FluxL, eq1FluxR, eq1FluxD, eq1FluxU; + + arma_mat wsL, wsR, wsD, wsU, wsLR, wsDU; + + arma_mat diff; // for Riemann Solver + + /* Projection */ + rhoP = project_to_edges(rho, x, xEdges, y, yEdges, nGCs); + xVelP = project_to_edges(xVel, x, xEdges, y, yEdges, nGCs); + yVelP = project_to_edges(yVel, x, xEdges, y, yEdges, nGCs); + + // Resolve Scalar Fields into rho, xVel, yVel, and totalE (without rho) + rhoL = rhoP.L; + rhoR = rhoP.R; + rhoD = rhoP.D; + rhoU = rhoP.U; + + xVelL = xVelP.L; + xVelR = xVelP.R; + xVelD = xVelP.D; + xVelU = xVelP.U; + + yVelL = yVelP.L; + yVelR = yVelP.R; + yVelD = yVelP.D; + yVelU = yVelP.U; + + velL2 = xVelL % xVelL + yVelL % yVelL; + velR2 = xVelR % xVelR + yVelR % yVelR; + velD2 = xVelD % xVelD + yVelD % yVelD; + velU2 = xVelU % xVelU + yVelU % yVelU; + + /* Calculate Edge Fluxes */ + // Note that dot product between normal vector at edge and flux vector + // resolves into a pure one component flux or either hat{x} or hat{y} + + // Flux calculated from the left of the edge + eq1FluxL = rhoL % xVelL % grid.sqrt_g_Left.slice(iAlt); + // Flux calculated from the right of the edge + eq1FluxR = rhoR % xVelR % grid.sqrt_g_Left.slice(iAlt); + // Flux calculated from the down of the edge + eq1FluxD = rhoD % yVelD % grid.sqrt_g_Down.slice(iAlt); + // Flux calculated from the up of the edge + eq1FluxU = rhoU % yVelU % grid.sqrt_g_Down.slice(iAlt); + + /* Wave Speed Calculation */ + wsL = sqrt(velL2); + wsR = sqrt(velR2); + wsD = sqrt(velD2); + wsU = sqrt(velU2); + + wsLR = wsR; + + for (int i = 0; i < nXs + 1; i++) { + for (int j = 0; j < nYs; j++) { + if (wsL(i, j) > wsLR(i, j)) + wsLR(i, j) = wsL(i, j); + } + } + + wsDU = wsD; + + for (int i = 0; i < nXs; i++) { + for (int j = 0; j < nYs + 1; j++) { + if (wsU(i, j) > wsDU(i, j)) + wsDU(i, j) = wsU(i, j); + } + } + + /* Calculate average flux at the edges (Rusanov Flux) */ + diff = (rhoR - rhoL) % grid.sqrt_g_Left.slice( + iAlt); // State difference, need to add sqrt(g) + eq1FluxLR = (eq1FluxL + eq1FluxR) / 2 + 0.5 * wsLR % diff; + diff = (rhoU - rhoD) % grid.sqrt_g_Down.slice(iAlt); + eq1FluxDU = (eq1FluxD + eq1FluxU) / 2 + 0.5 * wsDU % diff; + + // Setup residual storage for return + arma_mat eq1_residual(nXs, nYs, fill::zeros); + + // State Update + // Note the ghost cells WILL NOT BE UPDATED + for (int j = nGCs; j < nYs - nGCs; j++) { + for (int i = nGCs; i < nXs - nGCs; i++) { + precision_t rhoResidual_ij = dx * eq1FluxLR(i + 1, j) - + dx * eq1FluxLR(i, j) + + dx * eq1FluxDU(i, j + 1) - + dx * eq1FluxDU(i, j); + eq1_residual(i, j) = -1 / area * rhoResidual_ij; + } + } + + // Setup return vector + std::vector return_vector; + return_vector.push_back(eq1_residual); + + return return_vector; +} + +void Neutrals::solver_horizontal_rusanov_advection(Grid& grid, Times& time) { + // Function Reporting + std::string function = "Neutrals::solver_horizontal_rusanov_advection"; + static int iFunction = -1; + report.enter(function, iFunction); + + // Dimensions of Spatial Discretization + int64_t nXs = grid.get_nX(); + int64_t nYs = grid.get_nY(); + int64_t nGCs = grid.get_nGCs(); + int64_t nAlts = grid.get_nAlts(); + int iAlt, iSpec; + + // Time Discretization (TODO: change dt calculation method) + precision_t dt = time.get_dt(); + + + // Advance for bulk calculation first, calculate for every altitude + for (iAlt = nGCs; iAlt < nAlts - nGCs; iAlt++) { + /** Extract Grid Features **/ + arma_mat x = grid.refx_scgc.slice(iAlt); + arma_mat xEdges = grid.refx_Left.slice(iAlt); + arma_mat y = grid.refy_scgc.slice(iAlt); + arma_mat yEdges = grid.refy_Down.slice(iAlt); + + // Get reference grid dimensions (Assume dx = dy and equidistant) + arma_vec x_vec = x.col(0); + precision_t dx = x_vec(1) - x_vec(0); + precision_t area = dx * dx; + arma_mat jacobian = grid.sqrt_g_scgc.slice(iAlt); + + /** States preprocessing **/ + /* MASS DENSITY */ + arma_mat rho = rho_scgc.slice(iAlt); + + /* VELOCITY */ + // Get spherical velocity + arma_mat uVel = velocity_vcgc[0].slice(iAlt); + arma_mat vVel = velocity_vcgc[1].slice(iAlt); + // Convert to contravariant (reference) velocity + arma_mat xVel = uVel % grid.A11_inv_scgc.slice(iAlt) + vVel % + grid.A12_inv_scgc.slice(iAlt); // u^1 + arma_mat yVel = uVel % grid.A21_inv_scgc.slice(iAlt) + vVel % + grid.A22_inv_scgc.slice(iAlt); // u^2 + + // Generate velocity magnitude squared + arma_mat vel2 = xVel % xVel + yVel % yVel; + + /** Advancing **/ + /* Initialize projection constructs and storages */ + projection_struct rhoP; + projection_struct xVelP; + projection_struct yVelP; + + // They are all pure scalar fields without sqrt(g) + arma_mat rhoL, rhoR, rhoD, rhoU; + arma_mat xVelL, xVelR, xVelD, xVelU; + arma_mat yVelL, yVelR, yVelD, yVelU; + + arma_mat velL2, velR2, velD2, velU2; + + /** Initialize Flux and Wave Speed Storages */ + arma_mat eq1FluxLR, eq1FluxDU; + arma_mat eq1FluxL, eq1FluxR, eq1FluxD, eq1FluxU; + + arma_mat wsL, wsR, wsD, wsU, wsLR, wsDU; + + arma_mat diff; // for Riemann Solver + + /* Projection */ + rhoP = project_to_edges(rho, x, xEdges, y, yEdges, nGCs); + xVelP = project_to_edges(xVel, x, xEdges, y, yEdges, nGCs); + yVelP = project_to_edges(yVel, x, xEdges, y, yEdges, nGCs); + + // Resolve Scalar Fields into rho, xVel, yVel, and totalE (without rho) + rhoL = rhoP.L; + rhoR = rhoP.R; + rhoD = rhoP.D; + rhoU = rhoP.U; + + xVelL = xVelP.L; + xVelR = xVelP.R; + xVelD = xVelP.D; + xVelU = xVelP.U; + + yVelL = yVelP.L; + yVelR = yVelP.R; + yVelD = yVelP.D; + yVelU = yVelP.U; + + velL2 = xVelL % xVelL + yVelL % yVelL; + velR2 = xVelR % xVelR + yVelR % yVelR; + velD2 = xVelD % xVelD + yVelD % yVelD; + velU2 = xVelU % xVelU + yVelU % yVelU; + + + /* Calculate Edge Fluxes */ + // Note that dot product between normal vector at edge and flux vector + // resolves into a pure one component flux or either hat{x} or hat{y} + + // Flux calculated from the left of the edge + eq1FluxL = rhoL % xVelL % grid.sqrt_g_Left.slice(iAlt); + // Flux calculated from the right of the edge + eq1FluxR = rhoR % xVelR % grid.sqrt_g_Left.slice(iAlt); + // Flux calculated from the down of the edge + eq1FluxD = rhoD % yVelD % grid.sqrt_g_Down.slice(iAlt); + // Flux calculated from the up of the edge + eq1FluxU = rhoU % yVelU % grid.sqrt_g_Down.slice(iAlt); + + /* Wave Speed Calculation */ + wsL = sqrt(velL2); + wsR = sqrt(velR2); + wsD = sqrt(velD2); + wsU = sqrt(velU2); + + wsLR = wsR; + + for (int i = 0; i < nXs + 1; i++) { + for (int j = 0; j < nYs; j++) { + if (wsL(i, j) > wsLR(i, j)) + wsLR(i, j) = wsL(i, j); + } + } + + wsDU = wsD; + + for (int i = 0; i < nXs; i++) { + for (int j = 0; j < nYs + 1; j++) { + if (wsU(i, j) > wsDU(i, j)) + wsDU(i, j) = wsU(i, j); + } + } + + /* Calculate average flux at the edges (Rusanov Flux) */ + diff = (rhoR - rhoL) % grid.sqrt_g_Left.slice( + iAlt); // State difference, need to add sqrt(g) + eq1FluxLR = (eq1FluxL + eq1FluxR) / 2 + 0.5 * wsLR % diff; + diff = (rhoU - rhoD) % grid.sqrt_g_Down.slice(iAlt); + eq1FluxDU = (eq1FluxD + eq1FluxU) / 2 + 0.5 * wsDU % diff; + + /* Update Bulk Scalars and Contravariant velocity */ + // Euler State Update + for (int j = nGCs; j < nYs - nGCs; j++) { + for (int i = nGCs; i < nXs - nGCs; i++) { + precision_t rhoResidual_ij = dx * eq1FluxLR(i + 1, j) - + dx * eq1FluxLR(i, j) + + dx * eq1FluxDU(i, j + 1) - + dx * eq1FluxDU(i, j); + rho(i, j) = rho(i, j) - dt / area / jacobian(i, j) * rhoResidual_ij; + } + } + + /* Re-derive Spherical Velocity and Bulk States */ + // Density + rho_scgc.slice(iAlt) = rho; + + // Bulk Velocity + //vel2 = xVel % xVel + yVel % yVel; // Squared Magnitude of Contravariant + //velocity_vcgc[0].slice(iAlt) = xVel%grid.A11_scgc.slice(iAlt) + yVel%grid.A12_scgc.slice(iAlt); + //velocity_vcgc[1].slice(iAlt) = xVel%grid.A21_scgc.slice(iAlt) + yVel%grid.A22_scgc.slice(iAlt); + + /* Update specie density */ + for (iSpec = 0; iSpec < nSpecies; iSpec++) { + //species[iSpec].density_scgc.slice(iAlt) = rho % species[iSpec].concentration_scgc.slice(iAlt); + species[iSpec].density_scgc.slice(iAlt) = rho / species[iSpec].mass; + } + + + report.exit(function); + return; + } +} + +void Neutrals::solver_horizontal_RK1_advection(Grid& grid, Times& time) { + // Function Reporting + std::string function = "Neutrals::solver_horizontal_RK1_advection"; + static int iFunction = -1; + report.enter(function, iFunction); + + // Dimensions of Spatial Discretization + int64_t nXs = grid.get_nX(); + int64_t nYs = grid.get_nY(); + int64_t nGCs = grid.get_nGCs(); + int64_t nAlts = grid.get_nAlts(); + int iAlt, iSpec; + + // Time Discretization (TODO: change dt calculation method) + precision_t dt = time.get_dt(); + + // Advance for bulk calculation first, calculate for every altitude + for (iAlt = nGCs; iAlt < nAlts - nGCs; iAlt++) { + /** Extract Grid Features **/ + arma_mat x = grid.refx_scgc.slice(iAlt); + arma_mat xEdges = grid.refx_Left.slice(iAlt); + arma_mat y = grid.refy_scgc.slice(iAlt); + arma_mat yEdges = grid.refy_Down.slice(iAlt); + + // Get reference grid dimensions (Assume dx = dy and equidistant) + arma_vec x_vec = x.col(0); + precision_t dx = x_vec(1) - x_vec(0); + precision_t area = dx * dx; + arma_mat jacobian = grid.sqrt_g_scgc.slice(iAlt); + + /** States preprocessing **/ + /* MASS DENSITY */ + arma_mat rho = rho_scgc.slice(iAlt); + + /* VELOCITY */ + // Get spherical velocity + arma_mat uVel = velocity_vcgc[0].slice(iAlt); + arma_mat vVel = velocity_vcgc[1].slice(iAlt); + // Convert to contravariant (reference) velocity + arma_mat xVel = uVel % grid.A11_inv_scgc.slice(iAlt) + vVel % + grid.A12_inv_scgc.slice(iAlt); // u^1 + arma_mat yVel = uVel % grid.A21_inv_scgc.slice(iAlt) + vVel % + grid.A22_inv_scgc.slice(iAlt); // u^2 + + /** Advancing with RK4 **/ + // Setup Containers + arma_mat rho_0 = rho; + + // FIRST (1) STEP, Compute F_0-> State_1 + // Pass in state vector + std::vector state_0; + state_0.push_back(rho_0); + state_0.push_back(xVel); + state_0.push_back(yVel); + std::vector f_0_vec = residual_horizontal_rusanov_advection(state_0, + grid, time); + + // Extract Gradients + arma_mat f_0_eq1 = f_0_vec[0]; + + /* Update Bulk Scalars and Contravariant velocity */ + // Euler State Update + for (int j = nGCs; j < nYs - nGCs; j++) { + for (int i = nGCs; i < nXs - nGCs; i++) + rho(i, j) = rho_0(i, j) + dt * f_0_eq1(i, j) / jacobian(i, j); + } + + /* Re-derive Spherical Velocity and Bulk States */ + // Density + rho_scgc.slice(iAlt) = rho; + + // Bulk Velocity + //vel2 = xVel % xVel + yVel % yVel; // Squared Magnitude of Contravariant + //velocity_vcgc[0].slice(iAlt) = xVel%grid.A11_scgc.slice(iAlt) + yVel%grid.A12_scgc.slice(iAlt); + //velocity_vcgc[1].slice(iAlt) = xVel%grid.A21_scgc.slice(iAlt) + yVel%grid.A22_scgc.slice(iAlt); + + /* Update specie density */ + for (iSpec = 0; iSpec < nSpecies; iSpec++) { + //species[iSpec].density_scgc.slice(iAlt) = rho % species[iSpec].concentration_scgc.slice(iAlt); + species[iSpec].density_scgc.slice(iAlt) = rho / species[iSpec].mass; + } + + + report.exit(function); + return; + } +} + +void Neutrals::solver_horizontal_RK2_advection(Grid& grid, Times& time) { + // Function Reporting + std::string function = "Neutrals::solver_horizontal_RK2_advection"; + static int iFunction = -1; + report.enter(function, iFunction); + + // Dimensions of Spatial Discretization + int64_t nXs = grid.get_nX(); + int64_t nYs = grid.get_nY(); + int64_t nGCs = grid.get_nGCs(); + int64_t nAlts = grid.get_nAlts(); + int iAlt, iSpec; + + // Time Discretization (TODO: change dt calculation method) + precision_t dt = time.get_dt(); + + // Advance for bulk calculation first, calculate for every altitude + for (iAlt = nGCs; iAlt < nAlts - nGCs; iAlt++) { + /** Extract Grid Features **/ + arma_mat x = grid.refx_scgc.slice(iAlt); + arma_mat xEdges = grid.refx_Left.slice(iAlt); + arma_mat y = grid.refy_scgc.slice(iAlt); + arma_mat yEdges = grid.refy_Down.slice(iAlt); + + // Get reference grid dimensions (Assume dx = dy and equidistant) + arma_vec x_vec = x.col(0); + precision_t dx = x_vec(1) - x_vec(0); + precision_t area = dx * dx; + arma_mat jacobian = grid.sqrt_g_scgc.slice(iAlt); + + /** States preprocessing **/ + /* MASS DENSITY */ + arma_mat rho = rho_scgc.slice(iAlt); + + /* VELOCITY */ + // Get spherical velocity + arma_mat uVel = velocity_vcgc[0].slice(iAlt); + arma_mat vVel = velocity_vcgc[1].slice(iAlt); + // Convert to contravariant (reference) velocity + arma_mat xVel = uVel % grid.A11_inv_scgc.slice(iAlt) + vVel % + grid.A12_inv_scgc.slice(iAlt); // u^1 + arma_mat yVel = uVel % grid.A21_inv_scgc.slice(iAlt) + vVel % + grid.A22_inv_scgc.slice(iAlt); // u^2 + + /** Advancing with RK4 **/ + // Setup Containers + arma_mat rho_0 = rho; + arma_mat rho_1(nXs, nYs, fill::zeros); // corresponding f_1 + + // FIRST (1) STEP, Compute F_0-> State_1 + // Pass in state vector + std::vector state_0; + state_0.push_back(rho_0); + state_0.push_back(xVel); + state_0.push_back(yVel); + std::vector f_0_vec = residual_horizontal_hlle_advection(state_0, + grid, time); + // Extract Gradients + arma_mat f_0_eq1 = f_0_vec[0]; + + /* Update Bulk Scalars and Contravariant velocity */ + // Euler State Update + for (int j = nGCs; j < nYs - nGCs; j++) { + for (int i = nGCs; i < nXs - nGCs; i++) + rho_1(i, j) = rho_0(i, j) + dt * f_0_eq1(i, j) / jacobian(i, j); + } + + // SECOND (2) STEP, Compute F_1-> State_2 + // Pass in state vector + std::vector state_1; + state_1.push_back(rho_1); + state_1.push_back(xVel); + state_1.push_back(yVel); + std::vector f_1_vec = residual_horizontal_hlle_advection(state_1, + grid, time); + // Extract Gradients + arma_mat f_1_eq1 = f_1_vec[0]; + + // Summing all steps for final update + arma_mat f_sum_eq1 = f_0_eq1 + f_1_eq1; + + /* Update Bulk Scalars and Contravariant velocity */ + // Euler State Update + for (int j = nGCs; j < nYs - nGCs; j++) { + for (int i = nGCs; i < nXs - nGCs; i++) + rho(i, j) = rho(i, j) + 0.5 * dt * f_sum_eq1(i, j) / jacobian(i, j); + } + + /* Re-derive Spherical Velocity and Bulk States */ + // Density + rho_scgc.slice(iAlt) = rho; + + // Bulk Velocity + //vel2 = xVel % xVel + yVel % yVel; // Squared Magnitude of Contravariant + //velocity_vcgc[0].slice(iAlt) = xVel%grid.A11_scgc.slice(iAlt) + yVel%grid.A12_scgc.slice(iAlt); + //velocity_vcgc[1].slice(iAlt) = xVel%grid.A21_scgc.slice(iAlt) + yVel%grid.A22_scgc.slice(iAlt); + + /* Update specie density */ + for (iSpec = 0; iSpec < nSpecies; iSpec++) { + //species[iSpec].density_scgc.slice(iAlt) = rho % species[iSpec].concentration_scgc.slice(iAlt); + species[iSpec].density_scgc.slice(iAlt) = rho / species[iSpec].mass; + } + + + report.exit(function); + return; + } +} + +void Neutrals::solver_horizontal_RK4_advection(Grid& grid, Times& time) { + // Function Reporting + std::string function = "Neutrals::solver_horizontal_RK4_advection"; + static int iFunction = -1; + report.enter(function, iFunction); + + // Dimensions of Spatial Discretization + int64_t nXs = grid.get_nX(); + int64_t nYs = grid.get_nY(); + int64_t nGCs = grid.get_nGCs(); + int64_t nAlts = grid.get_nAlts(); + int iAlt, iSpec; + + // Time Discretization (TODO: change dt calculation method) + precision_t dt = time.get_dt(); + + // Advance for bulk calculation first, calculate for every altitude + for (iAlt = nGCs; iAlt < nAlts - nGCs; iAlt++) { + /** Extract Grid Features **/ + arma_mat x = grid.refx_scgc.slice(iAlt); + arma_mat xEdges = grid.refx_Left.slice(iAlt); + arma_mat y = grid.refy_scgc.slice(iAlt); + arma_mat yEdges = grid.refy_Down.slice(iAlt); + + // Get reference grid dimensions (Assume dx = dy and equidistant) + arma_vec x_vec = x.col(0); + precision_t dx = x_vec(1) - x_vec(0); + precision_t area = dx * dx; + arma_mat jacobian = grid.sqrt_g_scgc.slice(iAlt); + + /** States preprocessing **/ + /* MASS DENSITY */ + arma_mat rho = rho_scgc.slice(iAlt); + + /* VELOCITY */ + // Get spherical velocity + arma_mat uVel = velocity_vcgc[0].slice(iAlt); + arma_mat vVel = velocity_vcgc[1].slice(iAlt); + // Convert to contravariant (reference) velocity + arma_mat xVel = uVel % grid.A11_inv_scgc.slice(iAlt) + vVel % + grid.A12_inv_scgc.slice(iAlt); // u^1 + arma_mat yVel = uVel % grid.A21_inv_scgc.slice(iAlt) + vVel % + grid.A22_inv_scgc.slice(iAlt); // u^2 + + /** Advancing with RK4 **/ + // Setup Containers + arma_mat rho_0 = rho; + arma_mat rho_1(nXs, nYs, fill::zeros); // corresponding f_1 + arma_mat rho_2(nXs, nYs, fill::zeros); // corresponding f_2 + arma_mat rho_3(nXs, nYs, fill::zeros); // corresponding f_3 + + // FIRST (1) STEP, Compute F_0-> State_1 + // Pass in state vector + std::vector state_0; + state_0.push_back(rho_0); + state_0.push_back(xVel); + state_0.push_back(yVel); + std::vector f_0_vec = residual_horizontal_rusanov_advection(state_0, + grid, time); + // Extract Gradients + arma_mat f_0_eq1 = f_0_vec[0]; + + /* Update Bulk Scalars and Contravariant velocity */ + // Euler State Update + for (int j = nGCs; j < nYs - nGCs; j++) { + for (int i = nGCs; i < nXs - nGCs; i++) + rho_1(i, j) = rho_0(i, j) + 0.5 * dt * f_0_eq1(i, j) / jacobian(i, j); + } + + // SECOND (2) STEP, Compute F_1-> State_2 + // Pass in state vector + std::vector state_1; + state_1.push_back(rho_1); + state_1.push_back(xVel); + state_1.push_back(yVel); + std::vector f_1_vec = residual_horizontal_rusanov_advection(state_1, + grid, time); + // Extract Gradients + arma_mat f_1_eq1 = f_1_vec[0]; + + /* Update Bulk Scalars and Contravariant velocity */ + // Euler State Update + for (int j = nGCs; j < nYs - nGCs; j++) { + for (int i = nGCs; i < nXs - nGCs; i++) + rho_2(i, j) = rho_0(i, j) + 0.5 * dt * f_1_eq1(i, j) / jacobian(i, j); + } + + // THIRD (3) STEP, Compute F_2-> State_3 + // Pass in state vector + std::vector state_2; + state_2.push_back(rho_2); + state_2.push_back(xVel); + state_2.push_back(yVel); + std::vector f_2_vec = residual_horizontal_rusanov_advection(state_2, + grid, time); + // Extract Gradients + arma_mat f_2_eq1 = f_2_vec[0]; + + /* Update Bulk Scalars and Contravariant velocity */ + // Euler State Update + for (int j = nGCs; j < nYs - nGCs; j++) { + for (int i = nGCs; i < nXs - nGCs; i++) + rho_3(i, j) = rho_0(i, j) + dt * f_2_eq1(i, j) / jacobian(i, j); + } + + // FOURTH (4) STEP, Compute F_3 + // Pass in state vector + std::vector state_3; + state_3.push_back(rho_3); + state_3.push_back(xVel); + state_3.push_back(yVel); + std::vector f_3_vec = residual_horizontal_rusanov_advection(state_3, + grid, time); + // Extract Gradients + arma_mat f_3_eq1 = f_3_vec[0]; + + // Summing all steps for final update + arma_mat f_sum_eq1 = f_0_eq1 + 2 * f_1_eq1 + 2 * f_2_eq1 + f_3_eq1; + + /* Update Bulk Scalars and Contravariant velocity */ + // Euler State Update + for (int j = nGCs; j < nYs - nGCs; j++) { + for (int i = nGCs; i < nXs - nGCs; i++) + rho(i, j) = rho(i, j) + dt / 6 * f_sum_eq1(i, j) / jacobian(i, j); + } + + /* Re-derive Spherical Velocity and Bulk States */ + // Density + rho_scgc.slice(iAlt) = rho; + + // Bulk Velocity + //vel2 = xVel % xVel + yVel % yVel; // Squared Magnitude of Contravariant + //velocity_vcgc[0].slice(iAlt) = xVel%grid.A11_scgc.slice(iAlt) + yVel%grid.A12_scgc.slice(iAlt); + //velocity_vcgc[1].slice(iAlt) = xVel%grid.A21_scgc.slice(iAlt) + yVel%grid.A22_scgc.slice(iAlt); + + /* Update specie density */ + for (iSpec = 0; iSpec < nSpecies; iSpec++) { + //species[iSpec].density_scgc.slice(iAlt) = rho % species[iSpec].concentration_scgc.slice(iAlt); + species[iSpec].density_scgc.slice(iAlt) = rho / species[iSpec].mass; + } + + + report.exit(function); + return; + } +} \ No newline at end of file diff --git a/src/solver_vertical_rusanov.cpp b/src/solver_vertical_rusanov.cpp index 50e376d8..d57a6c4b 100644 --- a/src/solver_vertical_rusanov.cpp +++ b/src/solver_vertical_rusanov.cpp @@ -86,23 +86,6 @@ void calc_facevalues_alts_rusanov(Grid &grid, 0.5 * dVarLimited.slice(iZ) % grid.dk_edge_m.slice(iZ); } - /* - if (iProc == 11) - std::cout << "facevalues : " - << inVar(7,19,17) << " " - << inVar(7,19,18) << " " - << inVar(7,19,19) << " " - << inVar(7,19,20) << " " - << dVarLimited(7,19,18) << " " - << grid.dk_edge_m(7,19,17) << " " - << outRight(7, 19, 17) << " " - << outRight(7, 19, 18) << " " - << outLeft(7, 19, 17) << " " - << outLeft(7, 19, 18) << " " - << dVarUp(7, 19) << " " - << dVarDown(7, 19) << "\n"; - */ - return; } @@ -144,16 +127,6 @@ void calc_grad_and_diff_alts_rusanov(Grid &grid, varLeft.slice(iZ) - varRight.slice(iZ)) / grid.dk_center_m_scgc.slice(iZ); - /* - if (iProc == 11) - std::cout << "calc_grad : " - << varLeft(7, 19, 17) << " " - << varLeft(7, 19, 18) << " " - << varRight(7, 19, 17) << " " - << varRight(7, 19, 18) << " " - << grid.dk_edge_m(7, 19, 17) << " " - << outGrad(7, 19, 17) << "\n"; - */ for (iZ = nGCs; iZ < nZs - nGCs + 1; iZ++) { for (iX = nGCs; iX < nXs - nGCs; iX++) for (iY = nGCs; iY < nYs - nGCs; iY++) { @@ -164,11 +137,7 @@ void calc_grad_and_diff_alts_rusanov(Grid &grid, diffFlux(iX, iY, iZ) = 0.5 * cMaxLocal * (varRight(iX, iY, iZ) - varLeft(iX, iY, iZ)); - //if (iZ <= 10 && iX == 4 && iY == 4) - // std::cout << "diff flux : " << diffFlux(iX, iY, iZ) - // << " " << cMaxLocal - // << " " << varRight(iX, iY, iZ) - // << " " << varLeft(iX, iY, iZ) << "\n"; + } } @@ -311,17 +280,7 @@ void Neutrals::solver_vertical_rusanov(Grid grid, species[iSpecies].velocity_vcgc[2] % gradLogN_s[iSpecies]) + dt * diffLogN_s[iSpecies]; species[iSpecies].newDensity_scgc = exp(log_s); - /* - std::cout << iSpecies << " " << log_s(2,2,19) << " " - << dt << " " - << divVertVel_s[iSpecies](2,2,19) << " " - << species[iSpecies].velocity_vcgc[2](2,2,19) << " " - << gradLogN_s[iSpecies](2,2,19) << " " - << species[iSpecies].velocity_vcgc[2](2,2,19) * gradLogN_s[iSpecies](2,2,19) << " " - << diffLogN_s[iSpecies](2,2,19) << " " - << species[iSpecies].density_scgc(2,2,19) << " " - << species[iSpecies].newDensity_scgc(2,2,19) << "\n"; - */ + accTotal = dt * grid.gravity_vcgc[2] - dt * temperature_scgc % gradLogN_s[iSpecies] * cKB / mass @@ -386,36 +345,6 @@ void Neutrals::solver_vertical_rusanov(Grid grid, } } - bool doPrintThis = false; - - if (doPrintThis) { - iX = 2; - iY = 2; - iSpecies = 0; - mass = species[iSpecies].mass; - - for (int iAlt = 19; iAlt < 20; iAlt++) { - std::cout << iAlt << " " - << log(species[iSpecies].density_scgc(iX, iY, iAlt)) << " " - << temperature_scgc(iX, iY, iAlt) << " " - << species[iSpecies].velocity_vcgc[2](iX, iY, iAlt) << " " - << temperature_scgc(iX, iY, iAlt) * gradLogN_s[iSpecies](iX, iY, - iAlt) * cKB / mass << " " - << gradTemp(iX, iY, iAlt) * cKB / mass << " " - << grid.gravity_vcgc[2](iX, iY, iAlt) << "\n"; - } - } - - //calc_neutral_friction(); - /* - for (iSpecies = 0; iSpecies < nSpecies; iSpecies++) { - if (species[iSpecies].DoAdvect) { - species[iSpecies].velocity_vcgc[2] = - species[iSpecies].velocity_vcgc[2] + dt * - species[iSpecies].acc_neutral_friction[2]; - } - } - */ calc_mass_density(); // Calculate bulk vertical winds: velocity_vcgc[2].zeros();