| Line | % of fetches | Source |
|---|---|---|
| 1 | #include "grid.hpp" | |
| 2 | #ifdef DO_OUTPUT | |
| 3 | #include <silo.h> | |
| 4 | #endif | |
| 5 | #include <fstream> | |
| 6 | #include <thread> | |
| 7 | #include <cmath> | |
| 8 | ||
| 9 | #include <hpx/include/lcos.hpp> | |
| 10 | //#define EQ_ONLY | |
| 11 | //#define RHO_ONLY | |
| 12 | ||
| 13 | namespace hpx { | |
| 14 | using mutex = hpx::lcos::local::spinlock; | |
| 15 | } | |
| 16 | ||
| 17 | std::vector<std::vector<real>>& TLS_V(); | |
| 18 | ||
| 19 | #include <unordered_map> | |
| 20 | ||
| 21 | inline bool float_eq(xpoint_type a, xpoint_type b) { | |
| 22 | constexpr static xpoint_type eps = 0.00000011920928955078125; // std::pow(xpoint_type(2), -23); | |
| 23 | // const xpoint_type eps = std::pow(xpoint_type(2), -23); | |
| 24 | return std::abs(a - b) < eps; | |
| 25 | } | |
| 26 | ||
| 27 | bool grid::xpoint_eq(const xpoint& a, const xpoint& b) { | |
| 28 | bool rc = true; | |
| 29 | for (integer d = 0; d != NDIM; ++d) { | |
| 30 | if (!float_eq(a[d], b[d])) { | |
| 31 | rc = false; | |
| 32 | break; | |
| 33 | } | |
| 34 | } | |
| 35 | return rc; | |
| 36 | } | |
| 37 | ||
| 38 | bool grid::node_point::operator==(const node_point& other) const { | |
| 39 | return xpoint_eq(other.pt, pt); | |
| 40 | } | |
| 41 | ||
| 42 | bool grid::node_point::operator<(const node_point& other) const { | |
| 43 | bool rc = false; | |
| 44 | for (integer d = 0; d != NDIM; ++d) { | |
| 45 | if (!float_eq(pt[d], other.pt[d])) { | |
| 46 | rc = (pt[d] < other.pt[d]); | |
| 47 | break; | |
| 48 | } | |
| 49 | } | |
| 50 | return rc; | |
| 51 | } | |
| 52 | ||
| 53 | void grid::merge_output_lists(grid::output_list_type& l1, grid::output_list_type&& l2) { | |
| 54 | ||
| 55 | std::unordered_map<zone_int_type, zone_int_type> index_map; | |
| 56 | ||
| 57 | if (l2.zones.size() > l1.zones.size()) { | |
| 58 | std::swap(l1, l2); | |
| 59 | } | |
| 60 | for (auto i = l2.nodes.begin(); i != l2.nodes.end(); ++i) { | |
| 61 | zone_int_type index, oindex; | |
| 62 | auto this_x = *i; | |
| 63 | oindex = this_x.index; | |
| 64 | auto j = l1.nodes.find(this_x); | |
| 65 | if (j != l1.nodes.end()) { | |
| 66 | index = j->index; | |
| 67 | } else { | |
| 68 | index = l1.nodes.size(); | |
| 69 | this_x.index = index; | |
| 70 | l1.nodes.insert(this_x); | |
| 71 | } | |
| 72 | index_map[oindex] = index; | |
| 73 | } | |
| 74 | integer zzz = l1.zones.size(); | |
| 75 | l1.zones.resize(zzz + l2.zones.size()); | |
| 76 | for (auto i = l2.zones.begin(); i != l2.zones.end(); ++i) { | |
| 77 | l1.zones[zzz] = index_map[*i]; | |
| 78 | ++zzz; | |
| 79 | } | |
| 80 | for (integer field = 0; field < NF + NGF + NPF; ++field) { | |
| 81 | const auto l1sz = l1.data[field].size(); | |
| 82 | l1.data[field].resize(l1sz + l2.data[field].size()); | |
| 83 | std::move(l2.data[field].begin(), l2.data[field].end(), l1.data[field].begin() + l1sz); | |
| 84 | } | |
| 85 | if (l1.analytic.size()) { | |
| 86 | for (integer field = 0; field < NF; ++field) { | |
| 87 | const auto l1sz = l1.analytic[field].size(); | |
| 88 | l1.analytic[field].resize(l1sz + l2.analytic[field].size()); | |
| 89 | std::move(l2.analytic[field].begin(), l2.analytic[field].end(), l1.analytic[field].begin() + l1sz); | |
| 90 | } | |
| 91 | } | |
| 92 | } | |
| 93 | ||
| 94 | grid::output_list_type grid::get_output_list(bool analytic) const { | |
| 95 | auto& V = TLS_V(); | |
| 96 | compute_primitives(); | |
| 97 | output_list_type rc; | |
| 98 | const integer vertex_order[8] = { 0, 1, 3, 2, 4, 5, 7, 6 }; | |
| 99 | ||
| 100 | std::set<node_point>& node_list = rc.nodes; | |
| 101 | std::vector<zone_int_type>& zone_list = rc.zones; | |
| 102 | std::array<std::vector<real>, NF + NGF + NPF> &data = rc.data; | |
| 103 | std::array<std::vector<real>, NF + NGF + NPF> &A = rc.analytic; | |
| 104 | ||
| 105 | for (integer field = 0; field != NF + NGF + NPF; ++field) { | |
| 106 | data[field].reserve(INX * INX * INX); | |
| 107 | if (analytic) { | |
| 108 | A[field].reserve(INX * INX * INX); | |
| 109 | } | |
| 110 | } | |
| 111 | const integer this_bw = H_BW; | |
| 112 | zone_list.reserve(cube(H_NX - 2 * this_bw) * NCHILD); | |
| 113 | for (integer i = this_bw; i != H_NX - this_bw; ++i) { | |
| 114 | for (integer j = this_bw; j != H_NX - this_bw; ++j) { | |
| 115 | for (integer k = this_bw; k != H_NX - this_bw; ++k) { | |
| 116 | const integer iii = hindex(i, j, k); | |
| 117 | const integer iiig = gindex(i - H_BW, j - H_BW, k - H_BW); | |
| 118 | #ifdef EQ_ONLY | |
| 119 | if (!(std::abs(X[ZDIM][iii]) < dx) && !(std::abs(X[YDIM][iii]) < dx)) { | |
| 120 | continue; | |
| 121 | } | |
| 122 | #endif | |
| 123 | for (integer ci = 0; ci != NVERTEX; ++ci) { | |
| 124 | const integer vi = vertex_order[ci]; | |
| 125 | const integer xi = (vi >> 0) & 1; | |
| 126 | const integer yi = (vi >> 1) & 1; | |
| 127 | ||
| 128 | const integer zi = (vi >> 2) & 1; | |
| 129 | node_point this_x; | |
| 130 | this_x.pt[XDIM] = X[XDIM][iii] + (real(xi) - HALF) * dx; | |
| 131 | this_x.pt[YDIM] = X[YDIM][iii] + (real(yi) - HALF) * dx; | |
| 132 | this_x.pt[ZDIM] = X[ZDIM][iii] + (real(zi) - HALF) * dx; | |
| 133 | auto iter = node_list.find(this_x); | |
| 134 | integer index; | |
| 135 | if (iter != node_list.end()) { | |
| 136 | index = iter->index; | |
| 137 | } else { | |
| 138 | index = node_list.size(); | |
| 139 | this_x.index = index; | |
| 140 | node_list.insert(this_x); | |
| 141 | } | |
| 142 | zone_list.push_back(index); | |
| 143 | } | |
| 144 | if (analytic) { | |
| 145 | for (integer field = 0; field != NF; ++field) { | |
| 146 | A[field].push_back(Ua[field][iii]); | |
| 147 | } | |
| 148 | } | |
| 149 | for (integer field = 0; field != NF; ++field) { | |
| 150 | data[field].push_back(U[field][iii]); | |
| 151 | } | |
| 152 | for (integer field = 0; field != NGF; ++field) { | |
| 153 | data[field + NF].push_back(G[iiig][field]); | |
| 154 | } | |
| 155 | data[NGF + NF + 0].push_back(V[vx_i][iii]); | |
| 156 | data[NGF + NF + 1].push_back(V[vy_i][iii]); | |
| 157 | data[NGF + NF + 2].push_back(V[vz_i][iii]); | |
| 158 | if (V[egas_i][iii] * V[rho_i][iii] < de_switch2 * U[egas_i][iii]) { | |
| 159 | data[NGF + NF + 3].push_back(std::pow(V[tau_i][iii], fgamma) / V[rho_i][iii]); | |
| 160 | } else { | |
| 161 | data[NGF + NF + 3].push_back(V[egas_i][iii]); | |
| 162 | } | |
| 163 | data[NGF + NF + 4].push_back(V[zz_i][iii]); | |
| 164 | } | |
| 165 | } | |
| 166 | } | |
| 167 | ||
| 168 | return rc; | |
| 169 | } | |
| 170 | ||
| 171 | void grid::output(const output_list_type& olists, std::string _filename, real _t, int cycle, bool analytic) { | |
| 172 | #ifdef DO_OUTPUT | |
| 173 | ||
| 174 | std::thread( | |
| 175 | [&](const std::string& filename, real t) { | |
| 176 | printf( "t = %e\n", t); | |
| 177 | const std::set<node_point>& node_list = olists.nodes; | |
| 178 | const std::vector<zone_int_type>& zone_list = olists.zones; | |
| 179 | ||
| 180 | const int nzones = zone_list.size() / NVERTEX; | |
| 181 | std::vector<int> zone_nodes; | |
| 182 | zone_nodes = std::move(zone_list); | |
| 183 | ||
| 184 | const int nnodes = node_list.size(); | |
| 185 | std::vector<double> x_coord(nnodes); | |
| 186 | std::vector<double> y_coord(nnodes); | |
| 187 | std::vector<double> z_coord(nnodes); | |
| 188 | std::array<double*, NDIM> node_coords = {x_coord.data(), y_coord.data(), z_coord.data()}; | |
| 189 | for (auto iter = node_list.begin(); iter != node_list.end(); ++iter) { | |
| 190 | const integer i = iter->index; | |
| 191 | x_coord[i] = iter->pt[0]; | |
| 192 | y_coord[i] = iter->pt[1]; | |
| 193 | z_coord[i] = iter->pt[2]; | |
| 194 | } | |
| 195 | ||
| 196 | constexpr | |
| 197 | int nshapes = 1; | |
| 198 | int shapesize[1] = {NVERTEX}; | |
| 199 | int shapetype[1] = {DB_ZONETYPE_HEX}; | |
| 200 | int shapecnt[1] = {nzones}; | |
| 201 | const char* coord_names[NDIM] = {"x", "y", "z"}; | |
| 202 | ||
| 203 | #ifndef __MIC__ | |
| 204 | auto olist = DBMakeOptlist(1); | |
| 205 | double time = double(t); | |
| 206 | int ndim = 3; | |
| 207 | DBAddOption(olist, DBOPT_CYCLE, &cycle); | |
| 208 | DBAddOption(olist, DBOPT_DTIME, &time); | |
| 209 | DBAddOption(olist, DBOPT_NSPACE, &ndim ); | |
| 210 | DBfile *db = DBCreateReal(filename.c_str(), DB_CLOBBER, DB_LOCAL, "Euler Mesh", DB_PDB); | |
| 211 | assert(db); | |
| 212 | DBPutZonelist2(db, "zones", nzones, int(NDIM), zone_nodes.data(), nzones * NVERTEX, 0, 0, 0, shapetype, shapesize, | |
| 213 | shapecnt, nshapes, olist); | |
| 214 | DBPutUcdmesh(db, "mesh", int(NDIM), const_cast<char**>(coord_names), node_coords.data(), nnodes, nzones, "zones", nullptr, DB_DOUBLE, | |
| 215 | olist); | |
| 216 | const char* analytic_names[] = {"rho_a", "egas_a", "sx_a", "sy_a", "sz_a", "tau_a"}; | |
| 217 | DBFreeOptlist(olist); | |
| 218 | for (int field = 0; field != NF + NGF + NPF; ++field) { | |
| 219 | auto olist = DBMakeOptlist(1); | |
| 220 | double time = double(t); | |
| 221 | int istrue = 1; | |
| 222 | int isfalse = 0; | |
| 223 | DBAddOption(olist, DBOPT_CYCLE, &cycle); | |
| 224 | DBAddOption(olist, DBOPT_DTIME, &time); | |
| 225 | DBAddOption(olist, DBOPT_NSPACE, &ndim ); | |
| 226 | if( field == rho_i || field == sx_i || field == sy_i || field == sz_i || field == spc_ac_i || field == spc_ae_i || field == spc_dc_i || field == spc_de_i || field == spc_vac_i ) { | |
| 227 | DBAddOption(olist, DBOPT_CONSERVED, &istrue); | |
| 228 | } else { | |
| 229 | DBAddOption(olist, DBOPT_CONSERVED, &isfalse ); | |
| 230 | } | |
| 231 | if( field < NF ) { | |
| 232 | DBAddOption(olist, DBOPT_EXTENSIVE, &istrue); | |
| 233 | } else { | |
| 234 | DBAddOption(olist, DBOPT_EXTENSIVE, &isfalse); | |
| 235 | } | |
| 236 | DBAddOption(olist, DBOPT_EXTENSIVE, &istrue); | |
| 237 | DBPutUcdvar1(db, field_names[field], "mesh", const_cast<void*>(reinterpret_cast<const void*>(olists.data[field].data())), nzones, nullptr, 0, DB_DOUBLE, DB_ZONECENT, | |
| 238 | olist); | |
| 239 | if( analytic && field < 6) { | |
| 240 | DBPutUcdvar1(db, analytic_names[field], "mesh", const_cast<void*>(reinterpret_cast<const void*>(olists.analytic[field].data())), nzones, nullptr, 0, DB_DOUBLE, DB_ZONECENT, | |
| 241 | olist); | |
| 242 | } | |
| 243 | DBFreeOptlist(olist); | |
| 244 | #ifdef RHO_ONLY | |
| 245 | break; | |
| 246 | #endif | |
| 247 | } | |
| 248 | DBClose(db); | |
| 249 | #endif | |
| 250 | }, _filename, _t).join(); | |
| 251 | #endif | |
| 252 | } | |
| 253 | ||
| 254 | std::size_t grid::load(FILE* fp) { | |
| 255 | std::size_t cnt = 0; | |
| 256 | auto foo = std::fread; | |
| 257 | { | |
| 258 | static hpx::mutex mtx; | |
| 259 | std::lock_guard<hpx::mutex> lock(mtx); | |
| 260 | cnt += foo(&scaling_factor, sizeof(real), 1, fp) * sizeof(real); | |
| 261 | cnt += foo(&max_level, sizeof(integer), 1, fp) * sizeof(integer); | |
| 262 | cnt += foo(&Acons, sizeof(real), 1, fp) * sizeof(real); | |
| 263 | cnt += foo(&Bcons, sizeof(integer), 1, fp) * sizeof(integer); | |
| 264 | } | |
| 265 | cnt += foo(&is_leaf, sizeof(bool), 1, fp) * sizeof(bool); | |
| 266 | cnt += foo(&is_root, sizeof(bool), 1, fp) * sizeof(bool); | |
| 267 | ||
| 268 | allocate(); | |
| 269 | ||
| 270 | for (integer f = 0; f != NF; ++f) { | |
| 271 | for (integer i = H_BW; i < H_NX - H_BW; ++i) { | |
| 272 | for (integer j = H_BW; j < H_NX - H_BW; ++j) { | |
| 273 | const integer iii = hindex(i, j, H_BW); | |
| 274 | cnt += foo(&(U[f][iii]), sizeof(real), INX, fp) * sizeof(real); | |
| 275 | } | |
| 276 | } | |
| 277 | } | |
| 278 | for (integer i = 0; i < G_NX; ++i) { | |
| 279 | for (integer j = 0; j < G_NX; ++j) { | |
| 280 | for (integer k = 0; k < G_NX; ++k) { | |
| 281 | const integer iii = gindex(i, j, k); | |
| 282 | real g[NGF]; | |
| 283 | cnt += foo(g, sizeof(real), NGF, fp) * sizeof(real); | |
| 284 | for( integer f = 0; f != NGF; ++f) { | |
| 285 | G[iii][f] = g[f]; | |
| 286 | } | |
| 287 | } | |
| 288 | } | |
| 289 | } | |
| 290 | cnt += foo(U_out.data(), sizeof(real), U_out.size(), fp) * sizeof(real); | |
| 291 | set_coordinates(); | |
| 292 | return cnt; | |
| 293 | } | |
| 294 | ||
| 295 | std::size_t grid::save(FILE* fp) const { | |
| 296 | std::size_t cnt = 0; | |
| 297 | auto foo = std::fwrite; | |
| 298 | { | |
| 299 | static hpx::mutex mtx; | |
| 300 | std::lock_guard<hpx::mutex> lock(mtx); | |
| 301 | cnt += foo(&scaling_factor, sizeof(real), 1, fp) * sizeof(real); | |
| 302 | cnt += foo(&max_level, sizeof(integer), 1, fp) * sizeof(integer); | |
| 303 | cnt += foo(&Acons, sizeof(real), 1, fp) * sizeof(real); | |
| 304 | cnt += foo(&Bcons, sizeof(integer), 1, fp) * sizeof(integer); | |
| 305 | } | |
| 306 | cnt += foo(&is_leaf, sizeof(bool), 1, fp) * sizeof(bool); | |
| 307 | cnt += foo(&is_root, sizeof(bool), 1, fp) * sizeof(bool); | |
| 308 | for (integer f = 0; f != NF; ++f) { | |
| 309 | for (integer i = H_BW; i < H_NX - H_BW; ++i) { | |
| 310 | for (integer j = H_BW; j < H_NX - H_BW; ++j) { | |
| 311 | const integer iii = hindex(i, j, H_BW); | |
| 312 | cnt += foo(&U[f][iii], sizeof(real), INX, fp) * sizeof(real); | |
| 313 | } | |
| 314 | } | |
| 315 | } | |
| 316 | for (integer i = 0; i < G_NX; ++i) { | |
| 317 | for (integer j = 0; j < G_NX; ++j) { | |
| 318 | for (integer k = 0; k < G_NX; ++k) { | |
| 319 | const integer iii = gindex(i, j, k); | |
| 320 | const auto d = G[iii][0]; | |
| 321 | cnt += foo(&d, sizeof(real), NGF, fp) * sizeof(real); | |
| 322 | } | |
| 323 | } | |
| 324 | } | |
| 325 | cnt += foo(U_out.data(), sizeof(real), U_out.size(), fp) * sizeof(real); | |
| 326 | return cnt; | |
| 327 | } | |
| 328 | ||
| 329 |
Copyright (c) 2006-2012 Rogue Wave Software, Inc. All Rights Reserved.
Patents pending.