compute_background_PE.cpp 6.67 KB
 Benjamin Storer committed Jul 05, 2018 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 ``````#include "../Science.hpp" #include "math.h" #include "../Par_util.hpp" #include "stdio.h" #include "../Split_reader.hpp" #include "../T_util.hpp" #include "../Parformer.hpp" #include "../Sorter.hpp" #include using blitz::Array; using blitz::cos; using namespace TArrayn; using namespace NSIntegrator; using namespace Transformer; bool compare_pairs( pair a, pair b ) { return a.first < b.first; } // Compute Background Potential Energy (BPE) void compute_Background_PE(double & BPE_tot, TArrayn::DTArray & rho, TArrayn::DTArray & quad3, int Nx, int Ny, int Nz, double Lx, double Ly, double Lz, double g, double rho_0, int iter, bool dimensional_rho, bool mapped, Array hill) { // Tensor variables for indexing blitz::firstIndex ii; blitz::secondIndex jj; blitz::thirdIndex kk; // Set mpi parameters int myrank, numprocs; MPI_Comm_rank(MPI_COMM_WORLD, &myrank); MPI_Comm_size(MPI_COMM_WORLD, &numprocs); // Arrays for sorting static Array sort_rho, sort_quad; static vector sort_hill(Nx), sort_dx(Nx); static vector Lx_partsum(Nx), hillvol_partsum(Nx); double *local_vols = new double[numprocs]; double *local_vols_r = new double[numprocs]; static double hill_max; static double hill_vol; static vector < pair > height_width(Nx); // Stuff to do once at the beginning if (iter == 0) { // adjust if mapped if ( mapped ) { // information about the hill hill_vol = pssum(sum(hill*Ly*(*get_quad_x()))); hill_max = psmax(max(hill)); /* copy and sort the hill */ double *hill_tmp = new double[Nx]; double *hill_tmp2 = new double[Nx]; double *dx_tmp = new double[Nx]; double *dx_tmp2 = new double[Nx]; // create temporary empty arrays for holding the grid and hill for (int II = 0; II < Nx; II++) { if ( (II >= hill.lbound(firstDim)) and (II <= hill.ubound(firstDim)) ) { // populate these arrays with their processor's information hill_tmp2[II] = hill(II); dx_tmp2[II] = (*get_quad_x())(II); } else { // else set to zero hill_tmp2[II] = 0.; dx_tmp2[II] = 0.; } } // share across processors MPI_Allreduce(hill_tmp2, hill_tmp, Nx, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); MPI_Allreduce(dx_tmp2, dx_tmp, Nx, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); // put into a double pair for (int II = 0; II < Nx; II++) { height_width[II].first = hill_tmp[II]; height_width[II].second = dx_tmp[II]; } // sort the height_width pair according to the height sort(height_width.begin(), height_width.end(), compare_pairs); // put sorted pairs into separate vectors for (int II = 0; II < Nx; II++) { sort_hill[II] = height_width[II].first; sort_dx[II] = height_width[II].second; } // Compute cumulative sums of hill volume (reusing Lx_partsum) for (int II = 0; II < Nx; II++) { Lx_partsum[II] = sort_hill[II]*sort_dx[II]; } std::partial_sum(Lx_partsum.begin(), Lx_partsum.end(), hillvol_partsum.begin()); // Compute cumulative sums of dxs std::partial_sum(sort_dx.begin(), sort_dx.end(), Lx_partsum.begin()); // clean-up delete[] dx_tmp; delete[] dx_tmp2; delete[] hill_tmp; delete[] hill_tmp2; } } // Compute sorted rho sortarray(rho, quad3, sort_rho, sort_quad); // Volume of memory-local portion of sorted array double local_vol = sum(sort_quad); // Share local volume information with other processors // ie. processor 0 has lightest fluid for (int II = 0; II < numprocs; II++) { local_vols[II] = (II <= myrank) ? local_vol : 0; } MPI_Allreduce(local_vols, local_vols_r, numprocs, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); // The total volume of fluid with greater or equal density // to the lightest element on the local domain double base_vol = local_vols_r[myrank]; // Find the depth at which base_vol will fill to double tmpH, Area_star; if ( !mapped ) { tmpH = base_vol/Lx/Ly; } else { // check if base_vol will fill above the max hill height if ((base_vol + hill_vol)/Lx/Ly >= hill_max) { tmpH = (base_vol + hill_vol)/Lx/Ly; } else { // if it doesn't, find the depth where the fluid will fill to int II = 0; while ( (base_vol > Lx_partsum[II]*sort_hill[II] - hillvol_partsum[II]) and (II0 && "Something bad happened leading up to the tmpH calculation."); // now subtract off the bit we went over by (draw yourself a picture, it works) tmpH = sort_hill[II] - (sort_hill[II]*Lx_partsum[II] - base_vol - hillvol_partsum[II])/Lx_partsum[II-1]; } } // Now compute BPE from the sorted rho field double dH = 0; double BPE = 0; int LL = Nx-1; for (int II = sort_rho.lbound(firstDim); II <= sort_rho.ubound(firstDim); II++) { for (int KK = sort_rho.lbound(thirdDim); KK <= sort_rho.ubound(thirdDim); KK++) { for (int JJ = sort_rho.lbound(secondDim); JJ <= sort_rho.ubound(secondDim); JJ++) { // Find the surface area of the water at depth tmpH if ( mapped && (tmpH < hill_max) ) { while ( (sort_hill[LL] > tmpH) && (LL > 0) ) { LL--; } Area_star = Lx_partsum[LL]*Ly; } else { Area_star = Lx*Ly; } // spread volume over domain and compute BPE for that cell dH = sort_quad(II,JJ,KK)/Area_star; if (dimensional_rho) { // rho is dimensional density BPE += g*sort_rho(II,JJ,KK)*(tmpH - 0.5*dH)*dH*Area_star; } else { // rho is density anomaly BPE += g*rho_0*(1+sort_rho(II,JJ,KK))*(tmpH - 0.5*dH)*dH*Area_star; } tmpH -= dH; } } } // Share BPE with other processors MPI_Allreduce(&BPE, &BPE_tot, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); // delete variables delete[] local_vols; delete[] local_vols_r; }``````