Line | % of fetches | Source | |
---|---|---|---|
1 | #include "future.hpp" | ||
2 | #include "grid.hpp" | ||
3 | #include "problem.hpp" | ||
4 | #include "options.hpp" | ||
5 | #include "profiler.hpp" | ||
6 | #include "taylor.hpp" | ||
7 | //#include "helmholtz.hpp" | ||
8 | #include "node_server.hpp" | ||
9 | #include "exact_sod.hpp" | ||
10 | |||
11 | #include <cmath> | ||
12 | #include <cassert> | ||
13 | |||
14 | #include <hpx/include/runtime.hpp> | ||
15 | |||
16 | extern options opts; | ||
17 | |||
18 | char const* grid::field_names[] = { "rho", "egas", "sx", "sy", "sz", "tau", "pot", "zx", "zy", "zz", "primary_core", "primary_envelope", "secondary_core", | ||
19 | "secondary_envelope", "vacuum", "phi", "gx", "gy", "gz", "vx", "vy", "vz", "eint", "zzs" }; | ||
20 | |||
21 | real grid::omega = ZERO; | ||
22 | space_vector grid::pivot(ZERO); | ||
23 | real grid::scaling_factor = 1.0; | ||
24 | |||
25 | integer grid::max_level = 0; | ||
26 | |||
27 | struct tls_data_t { | ||
28 | std::vector<std::vector<real>> v; | ||
29 | std::vector<std::vector<std::vector<real>>>dvdx; | ||
30 | std::vector<std::vector<std::vector<real>>> dudx; | ||
31 | std::vector<std::vector<std::vector<real>>> uf; | ||
32 | std::vector<std::vector<real>> zz; | ||
33 | }; | ||
34 | |||
35 | real grid::Acons = 1.0; | ||
36 | real grid::Bcons = 1.0; | ||
37 | |||
38 | #if !defined(_MSC_VER) | ||
39 | |||
40 | #include <boost/thread/tss.hpp> | ||
41 | |||
42 | class tls_t { | ||
43 | private: | ||
44 | pthread_key_t key; | ||
45 | public: | ||
46 | static void cleanup(void* ptr) { | ||
47 | tls_data_t* _ptr = (tls_data_t*) ptr; | ||
48 | delete _ptr; | ||
49 | } | ||
50 | tls_t() { | ||
51 | pthread_key_create(&key, cleanup); | ||
52 | } | ||
53 | tls_data_t* get_ptr() { | ||
54 | tls_data_t* ptr = (tls_data_t*) pthread_getspecific(key); | ||
55 | if (ptr == nullptr) { | ||
56 | ptr = new tls_data_t; | ||
57 | ptr->v.resize(NF, std::vector < real > (H_N3)); | ||
58 | ptr->zz.resize(NDIM, std::vector < real > (H_N3)); | ||
59 | ptr->dvdx.resize(NDIM, std::vector < std::vector < real >> (NF, std::vector < real > (H_N3))); | ||
60 | ptr->dudx.resize(NDIM, std::vector < std::vector < real >> (NF, std::vector < real > (H_N3))); | ||
61 | ptr->uf.resize(NFACE, std::vector < std::vector < real >> (NF, std::vector < real > (H_N3))); | ||
62 | pthread_setspecific(key, ptr); | ||
63 | } | ||
64 | return ptr; | ||
65 | } | ||
66 | }; | ||
67 | |||
68 | #else | ||
69 | #include <hpx/util/thread_specific_ptr.hpp> | ||
70 | |||
71 | class tls_t | ||
72 | { | ||
73 | private: | ||
74 | struct tls_data_tag {}; | ||
75 | static hpx::util::thread_specific_ptr<tls_data_t, tls_data_tag> data; | ||
76 | |||
77 | public: | ||
78 | // static void cleanup(void* ptr) | ||
79 | // { | ||
80 | // tls_data_t* _ptr = (tls_data_t*) ptr; | ||
81 | // delete _ptr; | ||
82 | // } | ||
83 | |||
84 | tls_data_t* get_ptr() | ||
85 | { | ||
86 | tls_data_t* ptr = data.get(); | ||
87 | if (ptr == nullptr) { | ||
88 | ptr = new tls_data_t; | ||
89 | ptr->v.resize(NF, std::vector < real > (H_N3)); | ||
90 | ptr->zz.resize(NDIM, std::vector < real > (H_N3)); | ||
91 | ptr->dvdx.resize(NDIM, std::vector < std::vector < real >> (NF, std::vector < real > (H_N3))); | ||
92 | ptr->dudx.resize(NDIM, std::vector < std::vector < real >> (NF, std::vector < real > (H_N3))); | ||
93 | ptr->uf.resize(NFACE, std::vector < std::vector < real >> (NF, std::vector < real > (H_N3))); | ||
94 | data.reset(ptr); | ||
95 | } | ||
96 | return ptr; | ||
97 | } | ||
98 | }; | ||
99 | |||
100 | hpx::util::thread_specific_ptr<tls_data_t, tls_t::tls_data_tag> tls_t::data; | ||
101 | |||
102 | #endif | ||
103 | |||
104 | static tls_t tls; | ||
105 | |||
106 | std::vector<std::vector<real>>& TLS_V() { | ||
107 | return tls.get_ptr()->v; | ||
108 | } | ||
109 | |||
110 | static std::vector<std::vector<std::vector<real>>>& TLS_dVdx() { | ||
111 | return tls.get_ptr()->dvdx; | ||
112 | } | ||
113 | |||
114 | static std::vector<std::vector<std::vector<real>>>& TLS_dUdx() { | ||
115 | return tls.get_ptr()->dudx; | ||
116 | } | ||
117 | |||
118 | static std::vector<std::vector<real>>& TLS_zz() { | ||
119 | return tls.get_ptr()->zz; | ||
120 | } | ||
121 | |||
122 | static std::vector<std::vector<std::vector<real>>>& TLS_Uf() { | ||
123 | return tls.get_ptr()->uf; | ||
124 | } | ||
125 | |||
126 | space_vector grid::get_cell_center(integer i, integer j, integer k) { | ||
127 | const integer iii0 = hindex(H_BW,H_BW,H_BW); | ||
128 | space_vector c; | ||
129 | c[XDIM] = X[XDIM][iii0] + (i) * dx; | ||
130 | c[YDIM] = X[XDIM][iii0] + (j) * dx; | ||
131 | c[ZDIM] = X[XDIM][iii0] + (k) * dx; | ||
132 | return c; | ||
133 | } | ||
134 | |||
135 | void grid::set_hydro_boundary(const std::vector<real>& data, const geo::direction& dir, integer width, bool etot_only) { | ||
136 | PROF_BEGIN; | ||
137 | std::array<integer, NDIM> lb, ub; | ||
138 | if (!etot_only) { | ||
139 | get_boundary_size(lb, ub, dir, OUTER, INX, width); | ||
140 | } else { | ||
141 | get_boundary_size(lb, ub, dir, OUTER, INX, width); | ||
142 | } | ||
143 | integer iter = 0; | ||
144 | |||
145 | for (integer field = 0; field != NF; ++field) { | ||
146 | if (!etot_only || (etot_only && field == egas_i)) { | ||
147 | for (integer i = lb[XDIM]; i < ub[XDIM]; ++i) { | ||
148 | for (integer j = lb[YDIM]; j < ub[YDIM]; ++j) { | ||
149 | for (integer k = lb[ZDIM]; k < ub[ZDIM]; ++k) { | ||
150 ![]() | 3.7% | U[field][hindex( i, j, k)] = data[iter]; | |
![]() ![]() ![]() | |||
151 | ++iter; | ||
152 | } | ||
153 | } | ||
154 | } | ||
155 | } | ||
156 | } | ||
157 | PROF_END; | ||
158 | } | ||
159 | |||
160 | std::vector<real> grid::get_hydro_boundary(const geo::direction& dir, integer width, bool etot_only) { | ||
161 | PROF_BEGIN; | ||
162 | std::array<integer, NDIM> lb, ub; | ||
163 | std::vector < real > data; | ||
164 | integer size; | ||
165 | if (!etot_only) { | ||
166 | size = NF * get_boundary_size(lb, ub, dir, INNER, INX, width); | ||
167 | } else { | ||
168 | size = get_boundary_size(lb, ub, dir, INNER, INX, width); | ||
169 | } | ||
170 | data.resize(size); | ||
171 | integer iter = 0; | ||
172 | |||
173 | for (integer field = 0; field != NF; ++field) { | ||
174 | if (!etot_only || (etot_only && field == egas_i)) { | ||
175 | for (integer i = lb[XDIM]; i < ub[XDIM]; ++i) { | ||
176 | for (integer j = lb[YDIM]; j < ub[YDIM]; ++j) { | ||
177 | for (integer k = lb[ZDIM]; k < ub[ZDIM]; ++k) { | ||
178 | data[iter] = U[field][hindex( i, j, k)]; | ||
179 | ++iter; | ||
180 | } | ||
181 | } | ||
182 | } | ||
183 | } | ||
184 | } | ||
185 | PROF_END; | ||
186 | return data; | ||
187 | |||
188 | } | ||
189 | |||
190 | line_of_centers_t grid::line_of_centers(const std::pair<space_vector, space_vector>& line) { | ||
191 | PROF_BEGIN; | ||
192 | line_of_centers_t loc; | ||
193 | for (integer i = H_BW; i != H_NX - H_BW; ++i) { | ||
194 | for (integer j = H_BW; j != H_NX - H_BW; ++j) { | ||
195 | for (integer k = H_BW; k != H_NX - H_BW; ++k) { | ||
196 | const integer iii = hindex(i, j, k); | ||
197 | const integer iiig = gindex(i-H_BW, j-H_BW, k-H_BW); | ||
198 | space_vector a = line.first; | ||
199 | const space_vector& o = line.second; | ||
200 | space_vector b; | ||
201 | real aa = 0.0; | ||
202 | real bb = 0.0; | ||
203 | real ab = 0.0; | ||
204 | for (integer d = 0; d != NDIM; ++d) { | ||
205 | a[d] -= o[d]; | ||
206 | b[d] = X[d][iii] - o[d]; | ||
207 | } | ||
208 | for (integer d = 0; d != NDIM; ++d) { | ||
209 | aa += a[d] * a[d]; | ||
210 | bb += b[d] * b[d]; | ||
211 | ab += a[d] * b[d]; | ||
212 | } | ||
213 | const real d = std::sqrt((aa * bb - ab * ab) / aa); | ||
214 | real p = ab / std::sqrt(aa); | ||
215 | std::vector < real > data(NF + NGF); | ||
216 | if (d < std::sqrt(3.0) * dx / 2.0) { | ||
217 | for (integer ui = 0; ui != NF; ++ui) { | ||
218 | data[ui] = U[ui][iii]; | ||
219 | } | ||
220 | for (integer gi = 0; gi != NGF; ++gi) { | ||
221 | data[NF + gi] = G[iiig][gi]; | ||
222 | } | ||
223 | loc.resize(loc.size() + 1); | ||
224 | loc[loc.size() - 1].first = p; | ||
225 | loc[loc.size() - 1].second = std::move(data); | ||
226 | } | ||
227 | } | ||
228 | } | ||
229 | } | ||
230 | PROF_END; | ||
231 | return loc; | ||
232 | } | ||
233 | |||
234 | std::pair<std::vector<real>, std::vector<real>> grid::diagnostic_error() const { | ||
235 | PROF_BEGIN; | ||
236 | std::pair < std::vector<real>, std::vector < real >> e; | ||
237 | const real dV = dx * dx * dx; | ||
238 | if (opts.problem == SOLID_SPHERE) { | ||
239 | e.first.resize(8, ZERO); | ||
240 | e.second.resize(8, ZERO); | ||
241 | } | ||
242 | for (integer i = 0; i != G_NX; ++i) { | ||
243 | for (integer j = 0; j != G_NX; ++j) { | ||
244 | for (integer k = 0; k != G_NX; ++k) { | ||
245 | const integer iii = gindex(i, j, k); | ||
246 | const integer bdif = H_BW; | ||
247 | const integer iiih = hindex(i + bdif, j + bdif, k + bdif); | ||
248 | const real x = X[XDIM][iiih]; | ||
249 | const real y = X[YDIM][iiih]; | ||
250 | const real z = X[ZDIM][iiih]; | ||
251 | if (opts.problem == SOLID_SPHERE) { | ||
252 | const auto a = solid_sphere_analytic_phi(x, y, z, 0.25); | ||
253 | std::vector < real > n(4); | ||
254 | n[phi_i] = G[iii][phi_i]; | ||
255 | n[gx_i] = G[iii][gx_i]; | ||
256 | n[gy_i] = G[iii][gy_i]; | ||
257 | n[gz_i] = G[iii][gz_i]; | ||
258 | const real rho = U[rho_i][iiih]; | ||
259 | for (integer l = 0; l != 4; ++l) { | ||
260 | e.first[l] += std::abs(a[l] - n[l]) * dV * rho; | ||
261 | e.first[4 + l] += std::abs(a[l]) * dV * rho; | ||
262 | e.second[l] += sqr((a[l] - n[l]) * rho) * dV; | ||
263 | e.second[4 + l] += sqr(a[l] * rho) * dV; | ||
264 | } | ||
265 | } | ||
266 | } | ||
267 | } | ||
268 | } | ||
269 | // printf("%e\n", e[0]); | ||
270 | PROF_END; | ||
271 | return e; | ||
272 | } | ||
273 | |||
274 | real grid::get_A() { | ||
275 | return Acons; | ||
276 | } | ||
277 | |||
278 | real grid::get_B() { | ||
279 | return Bcons; | ||
280 | } | ||
281 | |||
282 | analytic_func_type grid::analytic = nullptr; | ||
283 | |||
284 | void grid::set_analytic_func(const analytic_func_type& func) { | ||
285 | analytic = func; | ||
286 | } | ||
287 | |||
288 | real grid::get_omega() { | ||
289 | return omega; | ||
290 | } | ||
291 | |||
292 | void grid::velocity_inc(const space_vector& dv) { | ||
293 | |||
294 | for (integer iii = 0; iii != H_N3; ++iii) { | ||
295 | const real rho = U[rho_i][iii]; | ||
296 | if (rho != ZERO) { | ||
297 | const real rhoinv = ONE / rho; | ||
298 | real& sx = U[sx_i][iii]; | ||
299 | real& sy = U[sy_i][iii]; | ||
300 | real& sz = U[sz_i][iii]; | ||
301 | real& egas = U[egas_i][iii]; | ||
302 | egas -= HALF * (sx * sx + sy * sy + sz * sz) * rhoinv; | ||
303 | sx += dv[XDIM] * rho; | ||
304 | sy += dv[YDIM] * rho; | ||
305 | sz += dv[ZDIM] * rho; | ||
306 | egas += HALF * (sx * sx + sy * sy + sz * sz) * rhoinv; | ||
307 | } | ||
308 | } | ||
309 | |||
310 | } | ||
311 | |||
312 | void grid::set_pivot(const space_vector& p) { | ||
313 | pivot = p; | ||
314 | } | ||
315 | |||
316 | inline real minmod(real a, real b) { | ||
317 | // return (std::copysign(HALF, a) + std::copysign(HALF, b)) * std::min(std::abs(a), std::abs(b)); | ||
318 | bool a_is_neg = a < 0; | ||
319 | bool b_is_neg = b < 0; | ||
320 | if (a_is_neg != b_is_neg) | ||
321 | return ZERO; | ||
322 | |||
323 | real val = std::min(std::abs(a), std::abs(b)); | ||
324 | return a_is_neg ? -val : val; | ||
325 | } | ||
326 | |||
327 | inline real minmod_theta(real a, real b, real c, real theta) { | ||
328 | return minmod(theta * minmod(a, b), c); | ||
329 | } | ||
330 | |||
331 | inline real minmod_theta(real a, real b, real theta = 1.0) { | ||
332 | return minmod(theta * minmod(a, b), HALF * (a + b)); | ||
333 | } | ||
334 | |||
335 | std::vector<real> grid::get_flux_restrict(const std::array<integer, NDIM>& lb, const std::array<integer, NDIM>& ub, const geo::dimension& dim) const { | ||
336 | PROF_BEGIN; | ||
337 | std::vector < real > data; | ||
338 | integer size = 1; | ||
339 | for (auto& dim : geo::dimension::full_set()) { | ||
340 | size *= (ub[dim] - lb[dim]); | ||
341 | } | ||
342 | size /= (NCHILD / 2); | ||
343 | size *= NF; | ||
344 | data.reserve(size); | ||
345 | const integer stride1 = (dim == XDIM) ? (INX + 1) : (INX + 1) * (INX + 1); | ||
346 | const integer stride2 = (dim == ZDIM) ? (INX + 1) : 1; | ||
347 | for (integer field = 0; field != NF; ++field) { | ||
348 | for (integer i = lb[XDIM]; i < ub[XDIM]; i += 2) { | ||
349 | for (integer j = lb[YDIM]; j < ub[YDIM]; j += 2) { | ||
350 | for (integer k = lb[ZDIM]; k < ub[ZDIM]; k += 2) { | ||
351 | const integer i00 = findex(i, j, k); | ||
352 | const integer i10 = i00 + stride1; | ||
353 | const integer i01 = i00 + stride2; | ||
354 | const integer i11 = i00 + stride1 + stride2; | ||
355 | real value = ZERO; | ||
356 | value += F[dim][field][i00]; | ||
357 | value += F[dim][field][i10]; | ||
358 | value += F[dim][field][i01]; | ||
359 | value += F[dim][field][i11]; | ||
360 | const real f = dx / TWO; | ||
361 | if (field == zx_i) { | ||
362 | if (dim == YDIM) { | ||
363 | value += F[dim][sy_i][i00] * f; | ||
364 | value += F[dim][sy_i][i10] * f; | ||
365 | value -= F[dim][sy_i][i01] * f; | ||
366 | value -= F[dim][sy_i][i11] * f; | ||
367 | } else if (dim == ZDIM) { | ||
368 | value -= F[dim][sz_i][i00] * f; | ||
369 | value -= F[dim][sz_i][i10] * f; | ||
370 | value += F[dim][sz_i][i01] * f; | ||
371 | value += F[dim][sz_i][i11] * f; | ||
372 | } else if (dim == XDIM) { | ||
373 | value += F[dim][sy_i][i00] * f; | ||
374 | value += F[dim][sy_i][i10] * f; | ||
375 | value -= F[dim][sy_i][i01] * f; | ||
376 | value -= F[dim][sy_i][i11] * f; | ||
377 | value -= F[dim][sz_i][i00] * f; | ||
378 | value += F[dim][sz_i][i10] * f; | ||
379 | value -= F[dim][sz_i][i01] * f; | ||
380 | value += F[dim][sz_i][i11] * f; | ||
381 | } | ||
382 | } else if (field == zy_i) { | ||
383 | if (dim == XDIM) { | ||
384 | value -= F[dim][sx_i][i00] * f; | ||
385 | value -= F[dim][sx_i][i10] * f; | ||
386 | value += F[dim][sx_i][i01] * f; | ||
387 | value += F[dim][sx_i][i11] * f; | ||
388 | } else if (dim == ZDIM) { | ||
389 | value += F[dim][sz_i][i00] * f; | ||
390 | value -= F[dim][sz_i][i10] * f; | ||
391 | value += F[dim][sz_i][i01] * f; | ||
392 | value -= F[dim][sz_i][i11] * f; | ||
393 | } else if (dim == YDIM) { | ||
394 | value -= F[dim][sx_i][i00] * f; | ||
395 | value -= F[dim][sx_i][i10] * f; | ||
396 | value += F[dim][sx_i][i01] * f; | ||
397 | value += F[dim][sx_i][i11] * f; | ||
398 | value += F[dim][sz_i][i00] * f; | ||
399 | value -= F[dim][sz_i][i10] * f; | ||
400 | value += F[dim][sz_i][i01] * f; | ||
401 | value -= F[dim][sz_i][i11] * f; | ||
402 | } | ||
403 | } else if (field == zz_i) { | ||
404 | if (dim == XDIM) { | ||
405 | value += F[dim][sx_i][i00] * f; | ||
406 | value -= F[dim][sx_i][i10] * f; | ||
407 | value += F[dim][sx_i][i01] * f; | ||
408 | value -= F[dim][sx_i][i11] * f; | ||
409 | } else if (dim == YDIM) { | ||
410 | value -= F[dim][sy_i][i00] * f; | ||
411 | value += F[dim][sy_i][i10] * f; | ||
412 | value -= F[dim][sy_i][i01] * f; | ||
413 | value += F[dim][sy_i][i11] * f; | ||
414 | } else if (dim == ZDIM) { | ||
415 | value -= F[dim][sy_i][i00] * f; | ||
416 | value += F[dim][sy_i][i10] * f; | ||
417 | value -= F[dim][sy_i][i01] * f; | ||
418 | value += F[dim][sy_i][i11] * f; | ||
419 | value += F[dim][sx_i][i00] * f; | ||
420 | value += F[dim][sx_i][i10] * f; | ||
421 | value -= F[dim][sx_i][i01] * f; | ||
422 | value -= F[dim][sx_i][i11] * f; | ||
423 | } | ||
424 | } | ||
425 | value /= real(4); | ||
426 | data.push_back(value); | ||
427 | } | ||
428 | } | ||
429 | } | ||
430 | } | ||
431 | PROF_END; | ||
432 | return data; | ||
433 | } | ||
434 | |||
435 | void grid::set_flux_restrict(const std::vector<real>& data, const std::array<integer, NDIM>& lb, const std::array<integer, NDIM>& ub, | ||
436 | const geo::dimension& dim) { | ||
437 | PROF_BEGIN; | ||
438 | integer index = 0; | ||
439 | for (integer field = 0; field != NF; ++field) { | ||
440 | for (integer i = lb[XDIM]; i < ub[XDIM]; ++i) { | ||
441 | for (integer j = lb[YDIM]; j < ub[YDIM]; ++j) { | ||
442 | for (integer k = lb[ZDIM]; k < ub[ZDIM]; ++k) { | ||
443 | const integer iii = findex(i, j, k); | ||
444 | F[dim][field][iii] = data[index]; | ||
445 | ++index; | ||
446 | } | ||
447 | } | ||
448 | } | ||
449 | } | ||
450 | PROF_END; | ||
451 | } | ||
452 | |||
453 | void grid::set_outflows(std::vector<real>&& u) { | ||
454 | U_out = std::move(u); | ||
455 | } | ||
456 | |||
457 | void grid::set_prolong(const std::vector<real>& data, std::vector<real>&& outflows) { | ||
458 | PROF_BEGIN; | ||
459 | integer index = 0; | ||
460 | U_out = std::move(outflows); | ||
461 | for (integer field = 0; field != NF; ++field) { | ||
462 | for (integer i = H_BW; i != H_NX - H_BW; ++i) { | ||
463 | for (integer j = H_BW; j != H_NX - H_BW; ++j) { | ||
464 | for (integer k = H_BW; k != H_NX - H_BW; ++k) { | ||
465 | const integer iii = hindex(i, j, k); | ||
466 | auto& value = U[field][iii]; | ||
467 | value = data[index]; | ||
468 | ++index; | ||
469 | } | ||
470 | } | ||
471 | } | ||
472 | } | ||
473 | PROF_END; | ||
474 | } | ||
475 | |||
476 | std::vector<real> grid::get_prolong(const std::array<integer, NDIM>& lb, const std::array<integer, NDIM>& ub, bool etot_only) { | ||
477 | PROF_BEGIN; | ||
478 | auto& dUdx = TLS_dUdx(); | ||
479 | auto& tmpz = TLS_zz(); | ||
480 | std::vector < real > data; | ||
481 | |||
482 | integer size = NF; | ||
483 | for (integer dim = 0; dim != NDIM; ++dim) { | ||
484 | size *= (ub[dim] - lb[dim]); | ||
485 | } | ||
486 | data.reserve(size); | ||
487 | auto lb0 = lb; | ||
488 | auto ub0 = ub; | ||
489 | for (integer d = 0; d != NDIM; ++d) { | ||
490 | lb0[d] /= 2; | ||
491 | ub0[d] = (ub[d] - 1) / 2 + 1; | ||
492 | } | ||
493 | compute_primitives(lb0, ub0, etot_only); | ||
494 | compute_primitive_slopes(1.0, lb0, ub0, etot_only); | ||
495 | compute_conserved_slopes(lb0, ub0, etot_only); | ||
496 | |||
497 | if (!etot_only) { | ||
498 | // #if !defined(HPX_HAVE_DATAPAR) | ||
499 | for (integer i = lb0[XDIM]; i != ub0[XDIM]; ++i) { | ||
500 | for (integer j = lb0[YDIM]; j != ub0[YDIM]; ++j) { | ||
501 | #pragma GCC ivdep | ||
502 | for (integer k = lb0[ZDIM]; k != ub0[ZDIM]; ++k) { | ||
503 | const integer iii = hindex(i,j,k); | ||
504 | tmpz[XDIM][iii] = U[zx_i][iii]; | ||
505 | tmpz[YDIM][iii] = U[zy_i][iii]; | ||
506 | tmpz[ZDIM][iii] = U[zz_i][iii]; | ||
507 | } | ||
508 | } | ||
509 | } | ||
510 | // #else | ||
511 | // #endif | ||
512 | } | ||
513 | |||
514 | for (integer field = 0; field != NF; ++field) { | ||
515 | if (!etot_only || (etot_only && field == egas_i)) { | ||
516 | for (integer i = lb[XDIM]; i != ub[XDIM]; ++i) { | ||
517 | const real xsgn = (i % 2) ? +1 : -1; | ||
518 | for (integer j = lb[YDIM]; j != ub[YDIM]; ++j) { | ||
519 | const real ysgn = (j % 2) ? +1 : -1; | ||
520 | #pragma GCC ivdep | ||
521 | for (integer k = lb[ZDIM]; k != ub[ZDIM]; ++k) { | ||
522 | const integer iii = hindex(i / 2, j / 2, k / 2); | ||
523 | const real zsgn = (k % 2) ? +1 : -1; | ||
524 | real value = U[field][iii]; | ||
525 | value += xsgn * dUdx[XDIM][field][iii] * 0.25; | ||
526 | value += ysgn * dUdx[YDIM][field][iii] * 0.25; | ||
527 | value += zsgn * dUdx[ZDIM][field][iii] * 0.25; | ||
528 | if (field == sx_i) { | ||
529 | U[zy_i][iii] -= 0.25 * zsgn * value * dx / 8.0; | ||
530 | U[zz_i][iii] += 0.25 * ysgn * value * dx / 8.0; | ||
531 | } else if (field == sy_i) { | ||
532 | U[zx_i][iii] += 0.25 * zsgn * value * dx / 8.0; | ||
533 | U[zz_i][iii] -= 0.25 * xsgn * value * dx / 8.0; | ||
534 | } else if (field == sz_i) { | ||
535 | U[zx_i][iii] -= 0.25 * ysgn * value * dx / 8.0; | ||
536 | U[zy_i][iii] += 0.25 * xsgn * value * dx / 8.0; | ||
537 | } | ||
538 | data.push_back(value); | ||
539 | } | ||
540 | } | ||
541 | } | ||
542 | } | ||
543 | } | ||
544 | |||
545 | if (!etot_only) { | ||
546 | for (integer i = lb0[XDIM]; i != ub0[XDIM]; ++i) { | ||
547 | for (integer j = lb0[YDIM]; j != ub0[YDIM]; ++j) { | ||
548 | #pragma GCC ivdep | ||
549 | for (integer k = lb0[ZDIM]; k != ub0[ZDIM]; ++k) { | ||
550 | const integer iii = hindex(i,j,k); | ||
551 | U[zx_i][iii] = tmpz[XDIM][iii]; | ||
552 | U[zy_i][iii] = tmpz[YDIM][iii]; | ||
553 | U[zz_i][iii] = tmpz[ZDIM][iii]; | ||
554 | } | ||
555 | } | ||
556 | } | ||
557 | } | ||
558 | |||
559 | PROF_END; | ||
560 | return data; | ||
561 | } | ||
562 | |||
563 | std::vector<real> grid::get_restrict() const { | ||
564 | PROF_BEGIN; | ||
565 | constexpr | ||
566 | integer Size = NF * INX * INX * INX / NCHILD + NF; | ||
567 | std::vector < real > data; | ||
568 | data.reserve(Size); | ||
569 | for (integer field = 0; field != NF; ++field) { | ||
570 | for (integer i = H_BW; i < H_NX - H_BW; i += 2) { | ||
571 | for (integer j = H_BW; j < H_NX - H_BW; j += 2) { | ||
572 | for (integer k = H_BW; k < H_NX - H_BW; k += 2) { | ||
573 | const integer iii = hindex(i, j, k); | ||
574 | real pt = ZERO; | ||
575 | for (integer x = 0; x != 2; ++x) { | ||
576 | for (integer y = 0; y != 2; ++y) { | ||
577 | for (integer z = 0; z != 2; ++z) { | ||
578 | const integer jjj = iii + x * H_DNX + y * H_DNY + z * H_DNZ; | ||
579 | pt += U[field][jjj]; | ||
580 | if (field == zx_i) { | ||
581 | pt += X[YDIM][jjj] * U[sz_i][jjj]; | ||
582 | pt -= X[ZDIM][jjj] * U[sy_i][jjj]; | ||
583 | } else if (field == zy_i) { | ||
584 | pt -= X[XDIM][jjj] * U[sz_i][jjj]; | ||
585 | pt += X[ZDIM][jjj] * U[sx_i][jjj]; | ||
586 | } else if (field == zz_i) { | ||
587 | pt += X[XDIM][jjj] * U[sy_i][jjj]; | ||
588 | pt -= X[YDIM][jjj] * U[sx_i][jjj]; | ||
589 | } | ||
590 | } | ||
591 | } | ||
592 | } | ||
593 | pt /= real(NCHILD); | ||
594 | data.push_back(pt); | ||
595 | } | ||
596 | } | ||
597 | } | ||
598 | } | ||
599 | for (integer field = 0; field != NF; ++field) { | ||
600 | data.push_back(U_out[field]); | ||
601 | } | ||
602 | PROF_END; | ||
603 | return data; | ||
604 | } | ||
605 | |||
606 | void grid::set_restrict(const std::vector<real>& data, const geo::octant& octant) { | ||
607 | PROF_BEGIN; | ||
608 | integer index = 0; | ||
609 | const integer i0 = octant.get_side(XDIM) * (INX / 2); | ||
610 | const integer j0 = octant.get_side(YDIM) * (INX / 2); | ||
611 | const integer k0 = octant.get_side(ZDIM) * (INX / 2); | ||
612 | for (integer field = 0; field != NF; ++field) { | ||
613 | for (integer i = H_BW; i != H_NX / 2; ++i) { | ||
614 | for (integer j = H_BW; j != H_NX / 2; ++j) { | ||
615 | for (integer k = H_BW; k != H_NX / 2; ++k) { | ||
616 | const integer iii = (i + i0) * H_DNX + (j + j0) * H_DNY + (k + k0) * H_DNZ; | ||
617 | auto& v = U[field][iii]; | ||
618 | v = data[index]; | ||
619 | if (field == zx_i) { | ||
620 | v -= X[YDIM][iii] * U[sz_i][iii]; | ||
621 | v += X[ZDIM][iii] * U[sy_i][iii]; | ||
622 | } else if (field == zy_i) { | ||
623 | v += X[XDIM][iii] * U[sz_i][iii]; | ||
624 | v -= X[ZDIM][iii] * U[sx_i][iii]; | ||
625 | } else if (field == zz_i) { | ||
626 | v -= X[XDIM][iii] * U[sy_i][iii]; | ||
627 | v += X[YDIM][iii] * U[sx_i][iii]; | ||
628 | } | ||
629 | ++index; | ||
630 | } | ||
631 | } | ||
632 | } | ||
633 | } | ||
634 | PROF_END; | ||
635 | } | ||
636 | |||
637 | std::pair<std::vector<real>, std::vector<real> > grid::field_range() const { | ||
638 | PROF_BEGIN; | ||
639 | std::pair < std::vector<real>, std::vector<real> > minmax; | ||
640 | minmax.first.resize(NF); | ||
641 | minmax.second.resize(NF); | ||
642 | for (integer field = 0; field != NF; ++field) { | ||
643 | minmax.first[field] = +std::numeric_limits < real > ::max(); | ||
644 | minmax.second[field] = -std::numeric_limits < real > ::max(); | ||
645 | } | ||
646 | for (integer i = H_BW; i != H_NX - H_BW; ++i) { | ||
647 | for (integer j = H_BW; j != H_NX - H_BW; ++j) { | ||
648 | for (integer k = H_BW; k != H_NX - H_BW; ++k) { | ||
649 | const integer iii = hindex(i, j, k); | ||
650 | for (integer field = 0; field != NF; ++field) { | ||
651 | minmax.first[field] = std::min(minmax.first[field], U[field][iii]); | ||
652 | minmax.second[field] = std::max(minmax.second[field], U[field][iii]); | ||
653 | } | ||
654 | } | ||
655 | } | ||
656 | } | ||
657 | PROF_END; | ||
658 | return minmax; | ||
659 | } | ||
660 | |||
661 | HPX_PLAIN_ACTION(grid::set_AB, set_AB_action); | ||
662 | |||
663 | void grid::set_AB(real a, real b) { | ||
664 | |||
665 | // FIXME: use proper broadcasting... | ||
666 | |||
667 | if (hpx::get_locality_id() == 0) { | ||
668 | std::vector<hpx::future<void>> futs; | ||
669 | auto remotes = hpx::find_remote_localities(); | ||
670 | futs.reserve(remotes.size()); | ||
671 | for (auto& l : remotes) { | ||
672 | futs.push_back(hpx::async < set_AB_action > (l, a, b)); | ||
673 | } | ||
674 | |||
675 | wait_all_and_propagate_exceptions(futs); | ||
676 | } | ||
677 | grid::Acons = a; | ||
678 | grid::Bcons = b; | ||
679 | } | ||
680 | |||
681 | HPX_PLAIN_ACTION(grid::set_omega, set_omega_action); | ||
682 | |||
683 | void grid::set_omega(real omega) { | ||
684 | |||
685 | // FIXME: use proper broadcasting... | ||
686 | |||
687 | if (hpx::get_locality_id() == 0) { | ||
688 | std::vector<hpx::future<void>> futs; | ||
689 | auto remotes = hpx::find_remote_localities(); | ||
690 | futs.reserve(remotes.size()); | ||
691 | for (auto& l : remotes) { | ||
692 | futs.push_back(hpx::async < set_omega_action > (l, omega)); | ||
693 | } | ||
694 | |||
695 | wait_all_and_propagate_exceptions(futs); | ||
696 | } | ||
697 | grid::omega = omega; | ||
698 | } | ||
699 | |||
700 | real grid::roche_volume(const std::pair<space_vector, space_vector>& axis, const std::pair<real, real>& l1, real cx, bool donor) const { | ||
701 | PROF_BEGIN; | ||
702 | const real dV = dx * dx * dx; | ||
703 | real V = 0.0; | ||
704 | for (integer i = H_BW; i != H_NX - H_BW; ++i) { | ||
705 | for (integer j = H_BW; j != H_NX - H_BW; ++j) { | ||
706 | for (integer k = H_BW; k != H_NX - H_BW; ++k) { | ||
707 | const integer D = 0 - H_BW; | ||
708 | const integer iii = hindex(i, j, k); | ||
709 | const integer iiig = gindex(i + D, j + D, k + D); | ||
710 | real x0 = X[XDIM][iii]; | ||
711 | real x = x0 - cx; | ||
712 | real y = X[YDIM][iii]; | ||
713 | real z = X[ZDIM][iii]; | ||
714 | const real R = std::sqrt(x0 * x0 + y * y); | ||
715 | real phi_eff = G[iiig][phi_i] - 0.5 * sqr(omega * R); | ||
716 | // real factor = axis.first[0] == l1.first ? 0.5 : 1.0; | ||
717 | if ((x0 <= l1.first && !donor) || (x0 >= l1.first && donor)) { | ||
718 | if (phi_eff <= l1.second) { | ||
719 | const real fx = G[iiig][gx_i] + x0 * sqr(omega); | ||
720 | const real fy = G[iiig][gy_i] + y * sqr(omega); | ||
721 | const real fz = G[iiig][gz_i]; | ||
722 | real g = x * fx + y * fy + z * fz; | ||
723 | if (g <= 0.0) { | ||
724 | V += dV; | ||
725 | } | ||
726 | } | ||
727 | } | ||
728 | } | ||
729 | } | ||
730 | } | ||
731 | PROF_END; | ||
732 | return V; | ||
733 | } | ||
734 | |||
735 | std::vector<real> grid::frac_volumes() const { | ||
736 | PROF_BEGIN; | ||
737 | std::vector < real > V(NSPECIES, 0.0); | ||
738 | const real dV = dx * dx * dx; | ||
739 | for (integer i = H_BW; i != H_NX - H_BW; ++i) { | ||
740 | for (integer j = H_BW; j != H_NX - H_BW; ++j) { | ||
741 | for (integer k = H_BW; k != H_NX - H_BW; ++k) { | ||
742 | const integer iii = hindex(i, j, k); | ||
743 | for (integer si = 0; si != NSPECIES; ++si) { | ||
744 | if (U[spc_i + si][iii] > 1.0e-5) { | ||
745 | V[si] += (U[spc_i + si][iii] / U[rho_i][iii]) * dV; | ||
746 | } | ||
747 | } | ||
748 | } | ||
749 | } | ||
750 | } | ||
751 | // printf( "%e", V[0]); | ||
752 | PROF_END; | ||
753 | return V; | ||
754 | } | ||
755 | |||
756 | bool grid::is_in_star(const std::pair<space_vector, space_vector>& axis, const std::pair<real, real>& l1, integer frac, integer iii) const { | ||
757 | bool use = false; | ||
758 | if (frac == 0) { | ||
759 | use = true; | ||
760 | } else { | ||
761 | space_vector a = axis.first; | ||
762 | const space_vector& o = axis.second; | ||
763 | space_vector b; | ||
764 | real aa = 0.0; | ||
765 | real ab = 0.0; | ||
766 | for (integer d = 0; d != NDIM; ++d) { | ||
767 | a[d] -= o[d]; | ||
768 | b[d] = X[d][iii] - o[d]; | ||
769 | } | ||
770 | for (integer d = 0; d != NDIM; ++d) { | ||
771 | aa += a[d] * a[d]; | ||
772 | ab += a[d] * b[d]; | ||
773 | } | ||
774 | real p = ab / std::sqrt(aa); | ||
775 | // printf( "%e\n", l1.first); | ||
776 | if (p < l1.first && frac == +1) { | ||
777 | use = true; | ||
778 | } else if (p >= l1.first && frac == -1) { | ||
779 | use = true; | ||
780 | } | ||
781 | } | ||
782 | return use; | ||
783 | } | ||
784 | |||
785 | real grid::z_moments(const std::pair<space_vector, space_vector>& axis, const std::pair<real, real>& l1, integer frac) const { | ||
786 | PROF_BEGIN; | ||
787 | real mom = 0.0; | ||
788 | const real dV = dx * dx * dx; | ||
789 | for (integer i = H_BW; i != H_NX - H_BW; ++i) { | ||
790 | for (integer j = H_BW; j != H_NX - H_BW; ++j) { | ||
791 | for (integer k = H_BW; k != H_NX - H_BW; ++k) { | ||
792 | const integer iii = hindex(i, j, k); | ||
793 | if (is_in_star(axis, l1, frac, iii)) { | ||
794 | mom += (sqr(X[XDIM][iii]) + sqr(dx) / 6.0) * U[rho_i][iii] * dV; | ||
795 | mom += (sqr(X[YDIM][iii]) + sqr(dx) / 6.0) * U[rho_i][iii] * dV; | ||
796 | } | ||
797 | } | ||
798 | } | ||
799 | } | ||
800 | PROF_END; | ||
801 | return mom; | ||
802 | } | ||
803 | |||
804 | std::vector<real> grid::conserved_sums(space_vector& com, space_vector& com_dot, const std::pair<space_vector, space_vector>& axis, | ||
805 | const std::pair<real, real>& l1, integer frac) const { | ||
806 | PROF_BEGIN; | ||
807 | std::vector < real > sum(NF, ZERO); | ||
808 | com[0] = com[1] = com[2] = 0.0; | ||
809 | com_dot[0] = com_dot[1] = com_dot[2] = 0.0; | ||
810 | const real dV = dx * dx * dx; | ||
811 | for (integer i = H_BW; i != H_NX - H_BW; ++i) { | ||
812 | for (integer j = H_BW; j != H_NX - H_BW; ++j) { | ||
813 | for (integer k = H_BW; k != H_NX - H_BW; ++k) { | ||
814 | const integer iii = hindex(i, j, k); | ||
815 | if (is_in_star(axis, l1, frac, iii)) { | ||
816 | com[0] += X[XDIM][iii] * U[rho_i][iii] * dV; | ||
817 | com[1] += X[YDIM][iii] * U[rho_i][iii] * dV; | ||
818 | com[2] += X[ZDIM][iii] * U[rho_i][iii] * dV; | ||
819 | com_dot[0] += U[sx_i][iii] * dV; | ||
820 | com_dot[1] += U[sy_i][iii] * dV; | ||
821 | com_dot[2] += U[sz_i][iii] * dV; | ||
822 | for (integer field = 0; field != NF; ++field) { | ||
823 | sum[field] += U[field][iii] * dV; | ||
824 | } | ||
825 | if (node_server::is_gravity_on()) { | ||
826 | sum[egas_i] += U[pot_i][iii] * HALF * dV; | ||
827 | } | ||
828 | sum[zx_i] += X[YDIM][iii] * U[sz_i][iii] * dV; | ||
829 | sum[zx_i] -= X[ZDIM][iii] * U[sy_i][iii] * dV; | ||
830 | sum[zy_i] -= X[XDIM][iii] * U[sz_i][iii] * dV; | ||
831 | sum[zy_i] += X[ZDIM][iii] * U[sx_i][iii] * dV; | ||
832 | sum[zz_i] += X[XDIM][iii] * U[sy_i][iii] * dV; | ||
833 | sum[zz_i] -= X[YDIM][iii] * U[sx_i][iii] * dV; | ||
834 | } | ||
835 | } | ||
836 | } | ||
837 | } | ||
838 | if (sum[rho_i] > 0.0) { | ||
839 | for (integer d = 0; d != NDIM; ++d) { | ||
840 | com[d] /= sum[rho_i]; | ||
841 | com_dot[d] /= sum[rho_i]; | ||
842 | } | ||
843 | } | ||
844 | PROF_END; | ||
845 | return sum; | ||
846 | } | ||
847 | |||
848 | std::vector<real> grid::gforce_sum(bool torque) const { | ||
849 | PROF_BEGIN; | ||
850 | std::vector < real > sum(NDIM, ZERO); | ||
851 | const real dV = dx * dx * dx; | ||
852 | for (integer i = H_BW; i != H_NX - H_BW; ++i) { | ||
853 | for (integer j = H_BW; j != H_NX - H_BW; ++j) { | ||
854 | for (integer k = H_BW; k != H_NX - H_BW; ++k) { | ||
855 | const auto D = 0 - H_BW; | ||
856 | const integer iii = hindex(i, j, k); | ||
857 | const integer iiig = gindex(i + D, j + D, k + D); | ||
858 | const real& rho = U[rho_i][iii]; | ||
859 | const real x = X[XDIM][iii]; | ||
860 | const real y = X[YDIM][iii]; | ||
861 | const real z = X[ZDIM][iii]; | ||
862 | const real fx = rho * G[iiig][gx_i] * dV; | ||
863 | const real fy = rho * G[iiig][gy_i] * dV; | ||
864 | const real fz = rho * G[iiig][gz_i] * dV; | ||
865 | if (!torque) { | ||
866 | sum[XDIM] += fx; | ||
867 | sum[YDIM] += fy; | ||
868 | sum[ZDIM] += fz; | ||
869 | } else { | ||
870 | sum[XDIM] -= z * fy - y * fz; | ||
871 | sum[YDIM] += z * fx - x * fz; | ||
872 | sum[ZDIM] -= y * fx - x * fy; | ||
873 | } | ||
874 | } | ||
875 | } | ||
876 | } | ||
877 | PROF_END; | ||
878 | return sum; | ||
879 | } | ||
880 | |||
881 | std::vector<real> grid::l_sums() const { | ||
882 | PROF_BEGIN; | ||
883 | std::vector < real > sum(NDIM); | ||
884 | const real dV = dx * dx * dx; | ||
885 | std::fill(sum.begin(), sum.end(), ZERO); | ||
886 | for (integer i = H_BW; i != H_NX - H_BW; ++i) { | ||
887 | for (integer j = H_BW; j != H_NX - H_BW; ++j) { | ||
888 | for (integer k = H_BW; k != H_NX - H_BW; ++k) { | ||
889 | const integer iii = hindex(i, j, k); | ||
890 | sum[XDIM] += X[YDIM][iii] * U[sz_i][iii] * dV; | ||
891 | sum[XDIM] -= X[ZDIM][iii] * U[sy_i][iii] * dV; | ||
892 | |||
893 | sum[YDIM] -= X[XDIM][iii] * U[sz_i][iii] * dV; | ||
894 | sum[YDIM] += X[ZDIM][iii] * U[sx_i][iii] * dV; | ||
895 | |||
896 | sum[ZDIM] += X[XDIM][iii] * U[sy_i][iii] * dV; | ||
897 | sum[ZDIM] -= X[YDIM][iii] * U[sx_i][iii] * dV; | ||
898 | |||
899 | } | ||
900 | } | ||
901 | } | ||
902 | PROF_END; | ||
903 | return sum; | ||
904 | } | ||
905 | |||
906 | bool grid::refine_me(integer lev) const { | ||
907 | PROF_BEGIN; | ||
908 | auto test = get_refine_test(); | ||
909 | if (lev < 2) { | ||
910 | PROF_END; | ||
911 | return true; | ||
912 | } | ||
913 | bool rc = false; | ||
914 | std::vector < real > state(NF); | ||
915 | std::array<std::vector<real>, NDIM> dud; | ||
916 | std::vector<real>& dudx = dud[0]; | ||
917 | std::vector<real>& dudy = dud[1]; | ||
918 | std::vector<real>& dudz = dud[2]; | ||
919 | dudx.resize(NF); | ||
920 | dudy.resize(NF); | ||
921 | dudz.resize(NF); | ||
922 | for (integer i = H_BW - R_BW; i != H_NX - H_BW + R_BW; ++i) { | ||
923 | for (integer j = H_BW - R_BW; j != H_NX - H_BW + R_BW; ++j) { | ||
924 | for (integer k = H_BW - R_BW; k != H_NX - H_BW + R_BW; ++k) { | ||
925 | int cnt = 0; | ||
926 | if (i < H_BW || i >= H_NX - H_BW) { | ||
927 | ++cnt; | ||
928 | } | ||
929 | if (j < H_BW || j >= H_NX - H_BW) { | ||
930 | ++cnt; | ||
931 | } | ||
932 | if (k < H_BW || k >= H_NX - H_BW) { | ||
933 | ++cnt; | ||
934 | } | ||
935 | if (cnt > 1) { | ||
936 | continue; | ||
937 | } | ||
938 | const integer iii = hindex(i, j, k); | ||
939 | for (integer i = 0; i != NF; ++i) { | ||
940 | state[i] = U[i][iii]; | ||
941 | dudx[i] = (U[i][iii + H_DNX] - U[i][iii - H_DNX]) / 2.0; | ||
942 | dudy[i] = (U[i][iii + H_DNY] - U[i][iii - H_DNY]) / 2.0; | ||
943 | dudz[i] = (U[i][iii + H_DNZ] - U[i][iii - H_DNZ]) / 2.0; | ||
944 | } | ||
945 | if (test(lev, max_level, X[XDIM][iii], X[YDIM][iii], X[ZDIM][iii], state, dud)) { | ||
946 | rc = true; | ||
947 | break; | ||
948 | } | ||
949 | } | ||
950 | if (rc) { | ||
951 | break; | ||
952 | } | ||
953 | } | ||
954 | if (rc) { | ||
955 | break; | ||
956 | } | ||
957 | } | ||
958 | PROF_END; | ||
959 | return rc; | ||
960 | } | ||
961 | |||
962 | grid::~grid() { | ||
963 | |||
964 | } | ||
965 | |||
966 | void grid::rho_mult(real f0, real f1) { | ||
967 | for (integer i = 0; i != H_NX; ++i) { | ||
968 | for (integer j = 0; j != H_NX; ++j) { | ||
969 | for (integer k = 0; k != H_NX; ++k) { | ||
970 | U[spc_ac_i][hindex(i,j,k)] *= f0; | ||
971 | U[spc_dc_i][hindex(i,j,k)] *= f1; | ||
972 | U[spc_ae_i][hindex(i,j,k)] *= f0; | ||
973 | U[spc_de_i][hindex(i,j,k)] *= f1; | ||
974 | U[rho_i][hindex(i,j,k)] = 0.0; | ||
975 | for (integer si = 0; si != NSPECIES; ++si) { | ||
976 | U[rho_i][hindex(i,j,k)] += U[spc_i + si][hindex(i, j, k)]; | ||
977 | } | ||
978 | } | ||
979 | } | ||
980 | } | ||
981 | |||
982 | } | ||
983 | |||
984 | void grid::rho_move(real x) { | ||
985 | real w = x / dx; | ||
986 | const real rho_floor = 1.0e-15; | ||
987 | |||
988 | w = std::max(-0.5, std::min(0.5, w)); | ||
989 | for (integer i = 1; i != H_NX - 1; ++i) { | ||
990 | for (integer j = 1; j != H_NX - 1; ++j) { | ||
991 | for (integer k = 1; k != H_NX - 1; ++k) { | ||
992 | for (integer si = spc_i; si != NSPECIES + spc_i; ++si) { | ||
993 | U[si][hindex(i,j,k)] += w * U[si][hindex(i+1,j,k)]; | ||
994 | U[si][hindex(i,j,k)] -= w * U[si][hindex(i-1,j,k)]; | ||
995 | U[si][hindex(i,j,k)] = std::max(U[si][hindex(i,j,k)], 0.0); | ||
996 | } | ||
997 | U[rho_i][hindex(i,j,k)] = 0.0; | ||
998 | for (integer si = 0; si != NSPECIES; ++si) { | ||
999 | U[rho_i][hindex(i,j,k)] += U[spc_i + si][hindex(i, j, k)]; | ||
1000 | } | ||
1001 | U[rho_i][hindex(i,j,k)] = std::max(U[rho_i][hindex(i,j,k)], rho_floor); | ||
1002 | } | ||
1003 | } | ||
1004 | } | ||
1005 | } | ||
1006 | /* | ||
1007 | space_vector& grid::center_of_mass_value(integer i, integer j, integer k) { | ||
1008 | return com[0][gindex(i, j, k)]; | ||
1009 | } | ||
1010 | |||
1011 | const space_vector& grid::center_of_mass_value(integer i, integer j, integer k) const { | ||
1012 | return com[0][gindex(i, j, k)]; | ||
1013 | }*/ | ||
1014 | |||
1015 | space_vector grid::center_of_mass() const { | ||
1016 | auto& M = *M_ptr; | ||
1017 | auto& mon = *mon_ptr; | ||
1018 | PROF_BEGIN; | ||
1019 | space_vector this_com; | ||
1020 | this_com[0] = this_com[1] = this_com[2] = ZERO; | ||
1021 | real m = ZERO; | ||
1022 | auto& com0 = *(com_ptr)[0]; | ||
1023 | for (integer i = 0; i != INX + 0; ++i) { | ||
1024 | for (integer j = 0; j != INX + 0; ++j) { | ||
1025 | for (integer k = 0; k != INX + 0; ++k) { | ||
1026 | const integer iii = gindex(i, j, k); | ||
1027 | const real this_m = is_leaf ? mon[iii] : M[iii](); | ||
1028 | for (auto& dim : geo::dimension::full_set()) { | ||
1029 | this_com[dim] += this_m * com0[iii][dim]; | ||
1030 | } | ||
1031 | m += this_m; | ||
1032 | } | ||
1033 | } | ||
1034 | } | ||
1035 | if (m != ZERO) { | ||
1036 | for (auto& dim : geo::dimension::full_set()) { | ||
1037 | this_com[dim] /= m; | ||
1038 | } | ||
1039 | } | ||
1040 | PROF_END; | ||
1041 | return this_com; | ||
1042 | } | ||
1043 | |||
1044 | grid::grid(real _dx, std::array<real, NDIM> _xmin) : | ||
1045 | U(NF), U0(NF), dUdt(NF), F(NDIM), X(NDIM), G(NGF), is_root(false), is_leaf(true) { | ||
1046 | dx = _dx; | ||
1047 | xmin = _xmin; | ||
1048 | allocate(); | ||
1049 | } | ||
1050 | |||
1051 | void grid::compute_primitives(const std::array<integer, NDIM> lb, const std::array<integer, NDIM> ub, bool etot_only) const { | ||
1052 | PROF_BEGIN; | ||
1053 | auto& V = TLS_V(); | ||
1054 | if (!etot_only) { | ||
1055 | for (integer i = lb[XDIM] - 1; i != ub[XDIM] + 1; ++i) { | ||
1056 | for (integer j = lb[YDIM] - 1; j != ub[YDIM] + 1; ++j) { | ||
1057 | #pragma GCC ivdep | ||
1058 | for (integer k = lb[ZDIM] - 1; k != ub[ZDIM] + 1; ++k) { | ||
1059 | const integer iii = hindex(i, j, k); | ||
1060 | V[rho_i][iii] = U[rho_i][iii]; | ||
1061 | V[tau_i][iii] = U[tau_i][iii]; | ||
1062 | const real rhoinv = 1.0 / V[rho_i][iii]; | ||
1063 | |||
1064 | V[egas_i][iii] = (U[egas_i][iii] | ||
1065 | #ifdef WD_EOS | ||
1066 | - ztwd_energy(U[rho_i][iii]) | ||
1067 | #endif | ||
1068 | ) * rhoinv; | ||
1069 | for (integer si = 0; si != NSPECIES; ++si) { | ||
1070 ![]() | 1.2% | V[spc_i + si][iii] = U[spc_i + si][iii] * rhoinv; | |
![]() | |||
1071 | } | ||
1072 | if (node_server::is_gravity_on()) { | ||
1073 | V[pot_i][iii] = U[pot_i][iii] * rhoinv; | ||
1074 | } | ||
1075 | for (integer d = 0; d != NDIM; ++d) { | ||
1076 | auto& v = V[sx_i + d][iii]; | ||
1077 | v = U[sx_i + d][iii] * rhoinv; | ||
1078 | V[egas_i][iii] -= 0.5 * v * v; | ||
1079 | V[zx_i + d][iii] = U[zx_i + d][iii] * rhoinv; | ||
1080 | } | ||
1081 | |||
1082 | V[sx_i][iii] += X[YDIM][iii] * omega; | ||
1083 | V[sy_i][iii] -= X[XDIM][iii] * omega; | ||
1084 | V[zz_i][iii] -= sqr(dx) * omega / 6.0; | ||
1085 | } | ||
1086 | } | ||
1087 | } | ||
1088 | } else { | ||
1089 | for (integer i = lb[XDIM] - 1; i != ub[XDIM] + 1; ++i) { | ||
1090 | for (integer j = lb[YDIM] - 1; j != ub[YDIM] + 1; ++j) { | ||
1091 | #pragma GCC ivdep | ||
1092 | for (integer k = lb[ZDIM] - 1; k != ub[ZDIM] + 1; ++k) { | ||
1093 | const integer iii = hindex(i, j, k); | ||
1094 | V[rho_i][iii] = U[rho_i][iii]; | ||
1095 | const real rhoinv = 1.0 / V[rho_i][iii]; | ||
1096 | V[egas_i][iii] = (U[egas_i][iii] | ||
1097 | #ifdef WD_EOS | ||
1098 | - ztwd_energy(U[rho_i][iii]) | ||
1099 | #endif | ||
1100 | ) * rhoinv; | ||
1101 | for (integer d = 0; d != NDIM; ++d) { | ||
1102 | auto& v = V[sx_i + d][iii]; | ||
1103 | v = U[sx_i + d][iii] * rhoinv; | ||
1104 | V[egas_i][iii] -= 0.5 * v * v; | ||
1105 | V[zx_i + d][iii] = U[zx_i + d][iii] * rhoinv; | ||
1106 | } | ||
1107 | V[sx_i][iii] += X[YDIM][iii] * omega; | ||
1108 | V[sy_i][iii] -= X[XDIM][iii] * omega; | ||
1109 | V[zz_i][iii] -= sqr(dx) * omega / 6.0; | ||
1110 | } | ||
1111 | } | ||
1112 | } | ||
1113 | } | ||
1114 | PROF_END; | ||
1115 | } | ||
1116 | |||
1117 | void grid::compute_primitive_slopes(real theta, const std::array<integer, NDIM> lb, const std::array<integer, NDIM> ub, bool etot_only) { | ||
1118 | PROF_BEGIN; | ||
1119 | auto& dVdx = TLS_dVdx(); | ||
1120 | auto& V = TLS_V(); | ||
1121 | for (integer f = 0; f != NF; ++f) { | ||
1122 | if (etot_only && (f == tau_i || f == pot_i || (f >= spc_i && f < spc_i + NSPECIES))) { | ||
1123 | continue; | ||
1124 | } | ||
1125 | const auto& v = V[f]; | ||
1126 | for (integer i = lb[XDIM]; i != ub[XDIM]; ++i) { | ||
1127 | for (integer j = lb[YDIM]; j != ub[YDIM]; ++j) { | ||
1128 | #pragma GCC ivdep | ||
1129 | for (integer k = lb[ZDIM]; k != ub[ZDIM]; ++k) { | ||
1130 | const integer iii = hindex(i,j,k); | ||
1131 | const auto v0 = v[iii]; | ||
1132 | dVdx[XDIM][f][iii] = minmod_theta(v[iii + H_DNX] - v0, v0 - v[iii - H_DNX], theta); | ||
1133 | dVdx[YDIM][f][iii] = minmod_theta(v[iii + H_DNY] - v0, v0 - v[iii - H_DNY], theta); | ||
1134 | dVdx[ZDIM][f][iii] = minmod_theta(v[iii + H_DNZ] - v0, v0 - v[iii - H_DNZ], theta); | ||
1135 | } | ||
1136 | } | ||
1137 | } | ||
1138 | } | ||
1139 | for (integer i = lb[XDIM]; i != ub[XDIM]; ++i) { | ||
1140 | for (integer j = lb[YDIM]; j != ub[YDIM]; ++j) { | ||
1141 | #pragma GCC ivdep | ||
1142 | for (integer k = lb[ZDIM]; k != ub[ZDIM]; ++k) { | ||
1143 | const integer iii = hindex(i,j,k); | ||
1144 | real dV_sym[3][3]; | ||
1145 | real dV_ant[3][3]; | ||
1146 | for (integer d0 = 0; d0 != NDIM; ++d0) { | ||
1147 | for (integer d1 = 0; d1 != NDIM; ++d1) { | ||
1148 | dV_sym[d1][d0] = (dVdx[d0][sx_i + d1][iii] + dVdx[d1][sx_i + d0][iii]) / 2.0; | ||
1149 | dV_ant[d1][d0] = 0.0; | ||
1150 | } | ||
1151 | } | ||
1152 | dV_ant[XDIM][YDIM] = +6.0 * V[zz_i][iii] / dx; | ||
1153 | dV_ant[XDIM][ZDIM] = -6.0 * V[zy_i][iii] / dx; | ||
1154 | dV_ant[YDIM][ZDIM] = +6.0 * V[zx_i][iii] / dx; | ||
1155 | dV_ant[YDIM][XDIM] = -dV_ant[XDIM][YDIM]; | ||
1156 | dV_ant[ZDIM][XDIM] = -dV_ant[XDIM][ZDIM]; | ||
1157 | dV_ant[ZDIM][YDIM] = -dV_ant[YDIM][ZDIM]; | ||
1158 | for (integer d0 = 0; d0 != NDIM; ++d0) { | ||
1159 | for (integer d1 = 0; d1 != NDIM; ++d1) { | ||
1160 | const real tmp = dV_sym[d0][d1] + dV_ant[d0][d1]; | ||
1161 | dVdx[d0][sx_i + d1][iii] = minmod(tmp, 2.0 / theta * dVdx[d0][sx_i + d1][iii]); | ||
1162 | } | ||
1163 | } | ||
1164 | } | ||
1165 | } | ||
1166 | } | ||
1167 | PROF_END; | ||
1168 | } | ||
1169 | |||
1170 | void grid::compute_conserved_slopes(const std::array<integer, NDIM> lb, const std::array<integer, NDIM> ub, bool etot_only) { | ||
1171 | PROF_BEGIN; | ||
1172 | auto& dVdx = TLS_dVdx(); | ||
1173 | auto& dUdx = TLS_dUdx(); | ||
1174 | auto& V = TLS_V(); | ||
1175 | const real theta = 1.0; | ||
1176 | if (!etot_only) { | ||
1177 | for (integer i = lb[XDIM]; i != ub[XDIM]; ++i) { | ||
1178 | for (integer j = lb[YDIM]; j != ub[YDIM]; ++j) { | ||
1179 | #pragma GCC ivdep | ||
1180 | for (integer k = lb[ZDIM]; k != ub[ZDIM]; ++k) { | ||
1181 | const integer iii = hindex(i,j,k); | ||
1182 | V[sx_i][iii] -= X[YDIM][iii] * omega; | ||
1183 | V[sy_i][iii] += X[XDIM][iii] * omega; | ||
1184 | V[zz_i][iii] += sqr(dx) * omega / 6.0; | ||
1185 | dVdx[YDIM][sx_i][iii] -= dx * omega; | ||
1186 | dVdx[XDIM][sy_i][iii] += dx * omega; | ||
1187 | } | ||
1188 | } | ||
1189 | } | ||
1190 | for (integer d0 = 0; d0 != NDIM; ++d0) { | ||
1191 | auto& dV = dVdx[d0]; | ||
1192 | auto& dU = dUdx[d0]; | ||
1193 | for (integer i = lb[XDIM]; i != ub[XDIM]; ++i) { | ||
1194 | for (integer j = lb[YDIM]; j != ub[YDIM]; ++j) { | ||
1195 | #pragma GCC ivdep | ||
1196 | for (integer k = lb[ZDIM]; k != ub[ZDIM]; ++k) { | ||
1197 | const integer iii = hindex(i,j,k); | ||
1198 | dU[rho_i][iii] = dV[rho_i][iii]; | ||
1199 | for (integer si = 0; si != NSPECIES; ++si) { | ||
1200 | dU[spc_i + si][iii] = V[spc_i + si][iii] * dV[rho_i][iii] + dV[spc_i + si][iii] * V[rho_i][iii]; | ||
1201 | } | ||
1202 | if (node_server::is_gravity_on()) { | ||
1203 | dU[pot_i][iii] = V[pot_i][iii] * dV[rho_i][iii] + dV[pot_i][iii] * V[rho_i][iii]; | ||
1204 | } | ||
1205 | dU[egas_i][iii] = V[egas_i][iii] * dV[rho_i][iii] + dV[egas_i][iii] * V[rho_i][iii]; | ||
1206 | for (integer d1 = 0; d1 != NDIM; ++d1) { | ||
1207 | dU[sx_i + d1][iii] = V[sx_i + d1][iii] * dV[rho_i][iii] + dV[sx_i + d1][iii] * V[rho_i][iii]; | ||
1208 | dU[egas_i][iii] += V[rho_i][iii] * (V[sx_i + d1][iii] * dV[sx_i + d1][iii]); | ||
1209 | dU[egas_i][iii] += dV[rho_i][iii] * 0.5 * sqr(V[sx_i + d1][iii]); | ||
1210 | dU[zx_i + d1][iii] = V[zx_i + d1][iii] * dV[rho_i][iii]; // + dV[zx_i + d1][iii] * V[rho_i][iii]; | ||
1211 | } | ||
1212 | #ifdef WD_EOS | ||
1213 | V[egas_i][iii] += ztwd_enthalpy(V[rho_i][iii]) * dV[rho_i][iii]; | ||
1214 | #endif | ||
1215 | dU[tau_i][iii] = dV[tau_i][iii]; | ||
1216 | } | ||
1217 | } | ||
1218 | } | ||
1219 | } | ||
1220 | } else { | ||
1221 | for (integer d0 = 0; d0 != NDIM; ++d0) { | ||
1222 | auto& dV = dVdx[d0]; | ||
1223 | auto& dU = dUdx[d0]; | ||
1224 | for (integer i = lb[XDIM]; i != ub[XDIM]; ++i) { | ||
1225 | for (integer j = lb[YDIM]; j != ub[YDIM]; ++j) { | ||
1226 | #pragma GCC ivdep | ||
1227 | for (integer k = lb[ZDIM]; k != ub[ZDIM]; ++k) { | ||
1228 | const integer iii = hindex(i,j,k); | ||
1229 | dU[egas_i][iii] = V[egas_i][iii] * dV[rho_i][iii] + dV[egas_i][iii] * V[rho_i][iii]; | ||
1230 | for (integer d1 = 0; d1 != NDIM; ++d1) { | ||
1231 | dU[egas_i][iii] += V[rho_i][iii] * (V[sx_i + d1][iii] * dV[sx_i + d1][iii]); | ||
1232 | dU[egas_i][iii] += dV[rho_i][iii] * 0.5 * sqr(V[sx_i + d1][iii]); | ||
1233 | } | ||
1234 | } | ||
1235 | } | ||
1236 | } | ||
1237 | } | ||
1238 | } | ||
1239 | |||
1240 | PROF_END; | ||
1241 | } | ||
1242 | |||
1243 | void grid::set_root(bool flag) { | ||
1244 | is_root = flag; | ||
1245 | } | ||
1246 | |||
1247 | void grid::set_leaf(bool flag) { | ||
1248 | if (is_leaf != flag) { | ||
1249 | is_leaf = flag; | ||
1250 | } | ||
1251 | } | ||
1252 | |||
1253 | void grid::set_fgamma(real fg) { | ||
1254 | fgamma = fg; | ||
1255 | } | ||
1256 | |||
1257 | real grid::get_fgamma() { | ||
1258 | return fgamma; | ||
1259 | } | ||
1260 | |||
1261 | real grid::fgamma = 5.0 / 3.0; | ||
1262 | |||
1263 | void grid::set_scaling_factor(real f) { | ||
1264 | scaling_factor = f; | ||
1265 | } | ||
1266 | |||
1267 | real grid::get_scaling_factor() { | ||
1268 | return scaling_factor; | ||
1269 | } | ||
1270 | |||
1271 | bool grid::get_leaf() const { | ||
1272 | return is_leaf; | ||
1273 | } | ||
1274 | |||
1275 | space_vector grid::get_pivot() { | ||
1276 | return pivot; | ||
1277 | } | ||
1278 | |||
1279 | real grid::get_source(integer i, integer j, integer k) const { | ||
1280 | return U[rho_i][hindex(i + H_BW, j + H_BW, k + H_BW)] * dx * dx * dx; | ||
1281 | } | ||
1282 | |||
1283 | std::vector<real> grid::get_outflows() { | ||
1284 | return U_out; | ||
1285 | } | ||
1286 | |||
1287 | void grid::set_coordinates() { | ||
1288 | PROF_BEGIN; | ||
1289 | for (integer i = 0; i != H_NX; ++i) { | ||
1290 | for (integer j = 0; j != H_NX; ++j) { | ||
1291 | for (integer k = 0; k != H_NX; ++k) { | ||
1292 | const integer iii = hindex(i, j, k); | ||
1293 | X[XDIM][iii] = (real(i - H_BW) + HALF) * dx + xmin[XDIM] - pivot[XDIM]; | ||
1294 | X[YDIM][iii] = (real(j - H_BW) + HALF) * dx + xmin[YDIM] - pivot[YDIM]; | ||
1295 | X[ZDIM][iii] = (real(k - H_BW) + HALF) * dx + xmin[ZDIM] - pivot[ZDIM]; | ||
1296 | } | ||
1297 | } | ||
1298 | } | ||
1299 | PROF_END; | ||
1300 | } | ||
1301 | |||
1302 | analytic_t grid::compute_analytic(real t) { | ||
1303 | analytic_t a; | ||
1304 | const real dv = dx * dx * dx; | ||
1305 | if (analytic != nullptr) { | ||
1306 | for (integer i = H_BW; i != H_NX - H_BW; ++i) | ||
1307 | for (integer j = H_BW; j != H_NX - H_BW; ++j) | ||
1308 | for (integer k = H_BW; k != H_NX - H_BW; ++k) { | ||
1309 | const integer iii = hindex(i,j,k); | ||
1310 | const auto A = analytic(X[XDIM][iii], X[YDIM][iii], X[ZDIM][iii], t); | ||
1311 | for (integer field = 0; field != NF; ++field) { | ||
1312 | Ua[field][iii] = A[field]; | ||
1313 | real dif = std::abs(Ua[field][iii] - U[field][iii]); | ||
1314 | a.l1[field] += dif * dv; | ||
1315 | a.l2[field] += dif * dif * dv; | ||
1316 | a.l1a[field] += std::abs(Ua[field][iii]) * dv; | ||
1317 | a.l2a[field] += Ua[field][iii] * Ua[field][iii] * dv; | ||
1318 | } | ||
1319 | } | ||
1320 | } | ||
1321 | return a; | ||
1322 | } | ||
1323 | |||
1324 | void grid::allocate() { | ||
1325 | PROF_BEGIN; | ||
1326 | U_out0 = std::vector < real > (NF, ZERO); | ||
1327 | U_out = std::vector < real > (NF, ZERO); | ||
1328 | dphi_dt = std::vector < real > (INX * INX * INX); | ||
1329 | G.resize(G_N3); | ||
1330 | for (integer dim = 0; dim != NDIM; ++dim) { | ||
1331 | X[dim].resize(H_N3); | ||
1332 | } | ||
1333 | for (integer field = 0; field != NF; ++field) { | ||
1334 | U0[field].resize(INX * INX * INX); | ||
1335 | U[field].resize(H_N3, 0.0); | ||
1336 | dUdt[field].resize(INX * INX * INX); | ||
1337 | for (integer dim = 0; dim != NDIM; ++dim) { | ||
1338 | F[dim][field].resize(F_N3); | ||
1339 | } | ||
1340 | } | ||
1341 | Ua = U; | ||
1342 | L.resize(G_N3); | ||
1343 | L_c.resize(G_N3); | ||
1344 | integer nlevel = 0; | ||
1345 | com_ptr.resize(2); | ||
1346 | |||
1347 | set_coordinates(); | ||
1348 | |||
1349 | L_mtx.reset(new hpx::lcos::local::spinlock); | ||
1350 | |||
1351 | PROF_END; | ||
1352 | } | ||
1353 | |||
1354 | grid::grid() : | ||
1355 | U(NF), U0(NF), dUdt(NF), F(NDIM), X(NDIM), G(NGF), is_root(false), is_leaf(true), U_out(NF, ZERO), U_out0(NF, ZERO), dphi_dt(H_N3) { | ||
1356 | // allocate(); | ||
1357 | } | ||
1358 | |||
1359 | grid::grid(const init_func_type& init_func, real _dx, std::array<real, NDIM> _xmin) : | ||
1360 | U(NF), U0(NF), dUdt(NF), F(NDIM), X(NDIM), G(NGF), is_root(false), is_leaf(true), U_out(NF, ZERO), U_out0(NF, ZERO), dphi_dt(H_N3) { | ||
1361 | PROF_BEGIN; | ||
1362 | dx = _dx; | ||
1363 | xmin = _xmin; | ||
1364 | allocate(); | ||
1365 | for (integer i = H_BW; i != H_NX - H_BW; ++i) { | ||
1366 | for (integer j = H_BW; j != H_NX - H_BW; ++j) { | ||
1367 | for (integer k = H_BW; k != H_NX - H_BW; ++k) { | ||
1368 | const integer iii = hindex(i, j, k); | ||
1369 | std::vector < real > this_u = init_func(X[XDIM][iii], X[YDIM][iii], X[ZDIM][iii], dx); | ||
1370 | for (integer field = 0; field != NF; ++field) { | ||
1371 | U[field][iii] = this_u[field]; | ||
1372 | } | ||
1373 | } | ||
1374 | } | ||
1375 | } | ||
1376 | if (node_server::is_gravity_on()) { | ||
1377 | for (integer i = 0; i != G_N3; ++i) { | ||
1378 | for (integer field = 0; field != NGF; ++field) { | ||
1379 | G[i][field] = 0.0; | ||
1380 | } | ||
1381 | } | ||
1382 | } | ||
1383 | PROF_END; | ||
1384 | } | ||
1385 | |||
1386 | inline real limit_range(real a, real b, real& c) { | ||
1387 | const real max = std::max(a, b); | ||
1388 | const real min = std::min(a, b); | ||
1389 | c = std::min(max, std::max(min, c)); | ||
1390 | } | ||
1391 | ; | ||
1392 | |||
1393 | inline real limit_range_all(real am, real ap, real& bl, real& br) { | ||
1394 | real avg = (br + bl) / 2.0; | ||
1395 | limit_range(am, ap, avg); | ||
1396 | limit_range(am, avg, bl); | ||
1397 | limit_range(ap, avg, br); | ||
1398 | } | ||
1399 | ; | ||
1400 | |||
1401 | inline void limit_slope(real& ql, real q0, real& qr) { | ||
1402 | const real tmp1 = qr - ql; | ||
1403 | const real tmp2 = qr + ql; | ||
1404 | // if ((qr - q0) * (q0 - ql) <= 0.0) { | ||
1405 | // if (qr < q0 || q0 < ql) { | ||
1406 | if (bool(qr < q0) != bool(q0 < ql)) { | ||
1407 | qr = ql = q0; | ||
1408 | } else if (tmp1 * (q0 - 0.5 * tmp2) > sqr(tmp1) / 6.0) { | ||
1409 | ql = 3.0 * q0 - 2.0 * qr; | ||
1410 | } else if (-(sqr(tmp1) / 6.0) > tmp1 * (q0 - 0.5 * tmp2)) { | ||
1411 | qr = 3.0 * q0 - 2.0 * ql; | ||
1412 | } | ||
1413 | } | ||
1414 | ; | ||
1415 | |||
1416 | void grid::reconstruct() { | ||
1417 | |||
1418 | PROF_BEGIN; | ||
1419 | auto& Uf = TLS_Uf(); | ||
1420 | auto& dUdx = TLS_dUdx(); | ||
1421 | auto& dVdx = TLS_dVdx(); | ||
1422 | auto& V = TLS_V(); | ||
1423 | |||
1424 | auto& slpx = dUdx[XDIM]; | ||
1425 | auto& slpy = dUdx[YDIM]; | ||
1426 | auto& slpz = dUdx[ZDIM]; | ||
1427 | |||
1428 | compute_primitives(); | ||
1429 | |||
1430 | for (integer field = 0; field != NF; ++field) { | ||
1431 | if (field >= zx_i && field <= zz_i) { | ||
1432 | continue; | ||
1433 | } | ||
1434 | // printf("%i\n", int(field)); | ||
1435 | const real theta_x = (field == sy_i || field == sz_i) ? 1.0 : 2.0; | ||
1436 | const real theta_y = (field == sx_i || field == sz_i) ? 1.0 : 2.0; | ||
1437 | const real theta_z = (field == sx_i || field == sy_i) ? 1.0 : 2.0; | ||
1438 | std::vector<real> const& Vfield = V[field]; | ||
1439 | #pragma GCC ivdep | ||
1440 | for (integer iii = H_NX * H_NX; iii != H_N3 - H_NX * H_NX; ++iii) { | ||
1441 | if (field == 1) { | ||
1442 | // printf("%i %i %i\n", int(iii / (H_NX * H_NX)), int((iii / H_NX) % H_NX), int(iii % H_NX)); | ||
1443 | // printf("%e %e %e\n", X[XDIM][iii], X[YDIM][iii], X[ZDIM][iii]); | ||
1444 | // printf("%e %e %e\n", V[field][iii + H_DNX], V[field][iii], V[field][iii - H_DNX]); | ||
1445 | // printf("%e %e %e\n", V[field][iii + H_DNY], V[field][iii], V[field][iii - H_DNY]); | ||
1446 | // printf("%e %e %e\n", V[field][iii + H_DNZ], V[field][iii], V[field][iii - H_DNZ]); | ||
1447 | // printf("%e %e %e\n", V[rho_i][iii + H_DNX], V[rho_i][iii], V[rho_i][iii - H_DNX]); | ||
1448 | // printf("%e %e %e\n", V[rho_i][iii + H_DNY], V[rho_i][iii], V[rho_i][iii - H_DNY]); | ||
1449 | // printf("%e %e %e\n", V[rho_i][iii + H_DNZ], V[rho_i][iii], V[rho_i][iii - H_DNZ]); | ||
1450 | } | ||
1451 ![]() | 1.5% | slpx[field][iii] = minmod_theta(Vfield[iii + H_DNX] - Vfield[iii], Vfield[iii] - Vfield[iii - H_DNX], theta_x); | |
![]() | |||
1452 | slpy[field][iii] = minmod_theta(Vfield[iii + H_DNY] - Vfield[iii], Vfield[iii] - Vfield[iii - H_DNY], theta_y); | ||
1453 ![]() | 1.2% | slpz[field][iii] = minmod_theta(Vfield[iii + H_DNZ] - Vfield[iii], Vfield[iii] - Vfield[iii - H_DNZ], theta_z); | |
![]() | |||
1454 | } | ||
1455 | } | ||
1456 | |||
1457 | if (opts.ang_con) { | ||
1458 | //#pragma GCC ivdep | ||
1459 | auto average = [](real& s1, real& s2) { s1 = s2 = 0.5 * (s1 + s2); }; | ||
1460 | auto step1 = [&](real& lhs, real const& rhs) { lhs += 6.0 * rhs / dx; }; | ||
1461 | auto step2 = [&](real& lhs, real const& rhs) { lhs -= 6.0 * rhs / dx; }; | ||
1462 | auto minmod_step = | ||
1463 | [](real& lhs, real const& r1, real const& r2, real const& r3) | ||
1464 | { | ||
1465 | lhs = minmod(lhs, 2.0 * minmod(r1 - r2, r2 - r3)); | ||
1466 | }; | ||
1467 | |||
1468 | for (integer iii = H_NX * H_NX; iii != H_N3 - H_NX * H_NX; ++iii) { | ||
1469 | average(slpx[sy_i][iii], slpy[sx_i][iii]); | ||
1470 | average(slpx[sz_i][iii], slpz[sx_i][iii]); | ||
1471 | average(slpy[sz_i][iii], slpz[sy_i][iii]); | ||
1472 | |||
1473 | step1(slpx[sy_i][iii], V[zz_i][iii]); | ||
1474 | step1(slpy[sz_i][iii], V[zx_i][iii]); | ||
1475 | step1(slpz[sx_i][iii], V[zy_i][iii]); | ||
1476 | |||
1477 | step2(slpy[sx_i][iii], V[zz_i][iii]); | ||
1478 | step2(slpz[sy_i][iii], V[zx_i][iii]); | ||
1479 | step2(slpx[sz_i][iii], V[zy_i][iii]); | ||
1480 | |||
1481 | minmod_step(slpx[sy_i][iii], V[sy_i][iii + H_DNX], V[sy_i][iii], V[sy_i][iii - H_DNX]); | ||
1482 | minmod_step(slpx[sz_i][iii], V[sz_i][iii + H_DNX], V[sz_i][iii], V[sz_i][iii - H_DNX]); | ||
1483 | minmod_step(slpy[sx_i][iii], V[sx_i][iii + H_DNY], V[sx_i][iii], V[sx_i][iii - H_DNY]); | ||
1484 | minmod_step(slpy[sz_i][iii], V[sz_i][iii + H_DNY], V[sz_i][iii], V[sz_i][iii - H_DNY]); | ||
1485 | minmod_step(slpz[sx_i][iii], V[sx_i][iii + H_DNZ], V[sx_i][iii], V[sx_i][iii - H_DNZ]); | ||
1486 | minmod_step(slpz[sy_i][iii], V[sy_i][iii + H_DNZ], V[sy_i][iii], V[sy_i][iii - H_DNZ]); | ||
1487 | |||
1488 | const real zx_lim = +(slpy[sz_i][iii] - slpz[sy_i][iii]) / 12.0; | ||
1489 | const real zy_lim = -(slpx[sz_i][iii] - slpz[sx_i][iii]) / 12.0; | ||
1490 | const real zz_lim = +(slpx[sy_i][iii] - slpy[sx_i][iii]) / 12.0; | ||
1491 | |||
1492 | const real Vzxi = V[zx_i][iii] - zx_lim * dx; | ||
1493 | const real Vzyi = V[zy_i][iii] - zy_lim * dx; | ||
1494 | const real Vzzi = V[zz_i][iii] - zz_lim * dx; | ||
1495 | |||
1496 | for (int face = 0; face != NFACE; ++face) { | ||
1497 | Uf[face][zx_i][iii] = Vzxi; | ||
1498 | Uf[face][zy_i][iii] = Vzyi; | ||
1499 | Uf[face][zz_i][iii] = Vzzi; | ||
1500 | } | ||
1501 | } | ||
1502 | } else { | ||
1503 | #pragma GCC ivdep | ||
1504 | for (integer iii = H_NX * H_NX; iii != H_N3 - H_NX * H_NX; ++iii) { | ||
1505 | const real Vzxi = V[zx_i][iii]; | ||
1506 | const real Vzyi = V[zy_i][iii]; | ||
1507 | const real Vzzi = V[zz_i][iii]; | ||
1508 | |||
1509 | for (int face = 0; face != NFACE; ++face) { | ||
1510 | Uf[face][zx_i][iii] = Vzxi; | ||
1511 | Uf[face][zy_i][iii] = Vzyi; | ||
1512 | Uf[face][zz_i][iii] = Vzzi; | ||
1513 | } | ||
1514 | } | ||
1515 | } | ||
1516 | for (integer field = 0; field != NF; ++field) { | ||
1517 | std::vector<real>& Vfield = V[field]; | ||
1518 | |||
1519 | std::vector<real>& UfFXPfield = Uf[FXP][field]; | ||
1520 | std::vector<real>& UfFXMfield = Uf[FXM][field]; | ||
1521 | std::vector<real> const& slpxfield = slpx[field]; | ||
1522 | |||
1523 | if (field >= zx_i && field <= zz_i) { | ||
1524 | continue; | ||
1525 | } | ||
1526 | if (!(field == sy_i || field == sz_i)) { | ||
1527 | #pragma GCC ivdep | ||
1528 | for (integer iii = 0; iii != H_N3 - H_NX * H_NX; ++iii) { | ||
1529 | const real& u0 = Vfield[iii]; | ||
1530 | UfFXPfield[iii] = UfFXMfield[iii + H_DNX] = (Vfield[iii + H_DNX] + u0) * HALF; | ||
1531 | } | ||
1532 | #pragma GCC ivdep | ||
1533 | for (integer iii = H_NX * H_NX; iii != H_N3 - H_NX * H_NX; ++iii) { | ||
1534 | const real& u0 = Vfield[iii]; | ||
1535 | const real& sx = slpxfield[iii]; | ||
1536 | UfFXPfield[iii] += (-(slpxfield[iii + H_DNX] - sx) / 3.0) * HALF; | ||
1537 | UfFXMfield[iii] += ( (slpxfield[iii - H_DNX] - sx) / 3.0) * HALF; | ||
1538 | limit_slope(UfFXMfield[iii], u0, UfFXPfield[iii]); | ||
1539 | } | ||
1540 | } else { | ||
1541 | #pragma GCC ivdep | ||
1542 | for (integer iii = H_NX * H_NX; iii != H_N3 - H_NX * H_NX; ++iii) { | ||
1543 | const real& u0 = Vfield[iii]; | ||
1544 | UfFXPfield[iii] = u0 + 0.5 * slpxfield[iii]; | ||
1545 | UfFXMfield[iii] = u0 - 0.5 * slpxfield[iii]; | ||
1546 | } | ||
1547 | } | ||
1548 | |||
1549 | std::vector<real>& UfFYPfield = Uf[FYP][field]; | ||
1550 | std::vector<real>& UfFYMfield = Uf[FYM][field]; | ||
1551 | std::vector<real> const& slpyfield = slpy[field]; | ||
1552 | |||
1553 | if (!(field == sx_i || field == sz_i)) { | ||
1554 | #pragma GCC ivdep | ||
1555 | for (integer iii = 0; iii != H_N3 - H_NX * H_NX; ++iii) { | ||
1556 | const real& u0 = Vfield[iii]; | ||
1557 | UfFYPfield[iii] = UfFYMfield[iii + H_DNY] = (Vfield[iii + H_DNY] + u0) * HALF; | ||
1558 | } | ||
1559 | #pragma GCC ivdep | ||
1560 | for (integer iii = H_NX * H_NX; iii != H_N3 - H_NX * H_NX; ++iii) { | ||
1561 | const real& u0 = Vfield[iii]; | ||
1562 | const real& sy = slpyfield[iii]; | ||
1563 | UfFYPfield[iii] += (-(slpyfield[iii + H_DNY] - sy) / 3.0) * HALF; | ||
1564 | UfFYMfield[iii] += ( (slpyfield[iii - H_DNY] - sy) / 3.0) * HALF; | ||
1565 | limit_slope(UfFYMfield[iii], u0, UfFYPfield[iii]); | ||
1566 | } | ||
1567 | } else { | ||
1568 | #pragma GCC ivdep | ||
1569 | for (integer iii = H_NX * H_NX; iii != H_N3 - H_NX * H_NX; ++iii) { | ||
1570 | const real& u0 = Vfield[iii]; | ||
1571 | UfFYPfield[iii] = u0 + 0.5 * slpyfield[iii]; | ||
1572 | UfFYMfield[iii] = u0 - 0.5 * slpyfield[iii]; | ||
1573 | } | ||
1574 | } | ||
1575 | |||
1576 | std::vector<real>& UfFZPfield = Uf[FZP][field]; | ||
1577 | std::vector<real>& UfFZMfield = Uf[FZM][field]; | ||
1578 | std::vector<real> const& slpzfield = slpz[field]; | ||
1579 | |||
1580 | if (!(field == sx_i || field == sy_i)) { | ||
1581 | #pragma GCC ivdep | ||
1582 | for (integer iii = 0; iii != H_N3 - H_NX * H_NX; ++iii) { | ||
1583 | const real& u0 = Vfield[iii]; | ||
1584 | UfFZPfield[iii] = UfFZMfield[iii + H_DNZ] = (Vfield[iii + H_DNZ] + u0) * HALF; | ||
1585 | } | ||
1586 | #pragma GCC ivdep | ||
1587 | for (integer iii = H_NX * H_NX; iii != H_N3 - H_NX * H_NX; ++iii) { | ||
1588 | const real& u0 = Vfield[iii]; | ||
1589 | const real& sz = slpzfield[iii]; | ||
1590 ![]() | 1.1% | UfFZPfield[iii] += (-(slpzfield[iii + H_DNZ] - sz) / 3.0) * HALF; | |
![]() ![]() | |||
1591 ![]() | 0.2% | UfFZMfield[iii] += ( (slpzfield[iii - H_DNZ] - sz) / 3.0) * HALF; | |
![]() ![]() | |||
1592 | limit_slope(UfFZMfield[iii], u0, UfFZPfield[iii]); | ||
1593 | } | ||
1594 | } else { | ||
1595 | #pragma GCC ivdep | ||
1596 | for (integer iii = H_NX * H_NX; iii != H_N3 - H_NX * H_NX; ++iii) { | ||
1597 | const real& u0 = Vfield[iii]; | ||
1598 | UfFZPfield[iii] = u0 + 0.5 * slpzfield[iii]; | ||
1599 | UfFZMfield[iii] = u0 - 0.5 * slpzfield[iii]; | ||
1600 | } | ||
1601 | } | ||
1602 | } | ||
1603 | |||
1604 | for (integer iii = 0; iii != H_N3; ++iii) { | ||
1605 | #pragma GCC ivdep | ||
1606 | for (integer face = 0; face != NFACE; ++face) { | ||
1607 | real w = 0.0; | ||
1608 | std::vector<std::vector<real> >& Ufface = Uf[face]; | ||
1609 | for (integer si = 0; si != NSPECIES; ++si) { | ||
1610 | w += Ufface[spc_i + si][iii]; | ||
1611 | } | ||
1612 | if (w > ZERO) { | ||
1613 | for (integer si = 0; si != NSPECIES; ++si) { | ||
1614 | Ufface[spc_i + si][iii] /= w; | ||
1615 | } | ||
1616 | } | ||
1617 | } | ||
1618 | } | ||
1619 | |||
1620 | if (node_server::is_gravity_on()) { | ||
1621 | //#pragma GCC ivdep | ||
1622 | std::vector<real>& UfFXMpot_i = Uf[FXM][pot_i]; | ||
1623 | std::vector<real>& UfFYMpot_i = Uf[FYM][pot_i]; | ||
1624 | std::vector<real>& UfFZMpot_i = Uf[FZM][pot_i]; | ||
1625 | |||
1626 | std::vector<real>& UfFXPpot_i = Uf[FXP][pot_i]; | ||
1627 | std::vector<real>& UfFYPpot_i = Uf[FYP][pot_i]; | ||
1628 | std::vector<real>& UfFZPpot_i = Uf[FZP][pot_i]; | ||
1629 | |||
1630 | for (integer iii = H_NX * H_NX; iii != H_N3 - H_NX * H_NX; ++iii) { | ||
1631 | const real phi_x = HALF * (UfFXMpot_i[iii] + UfFXPpot_i[iii - H_DNX]); | ||
1632 | const real phi_y = HALF * (UfFYMpot_i[iii] + UfFYPpot_i[iii - H_DNY]); | ||
1633 | const real phi_z = HALF * (UfFZMpot_i[iii] + UfFZPpot_i[iii - H_DNZ]); | ||
1634 | UfFXMpot_i[iii] = phi_x; | ||
1635 | UfFYMpot_i[iii] = phi_y; | ||
1636 | UfFZMpot_i[iii] = phi_z; | ||
1637 | UfFXPpot_i[iii - H_DNX] = phi_x; | ||
1638 | UfFYPpot_i[iii - H_DNY] = phi_y; | ||
1639 | UfFZPpot_i[iii - H_DNZ] = phi_z; | ||
1640 | } | ||
1641 | } | ||
1642 | for (integer field = 0; field != NF; ++field) { | ||
1643 | if (field != rho_i && field != tau_i) { | ||
1644 | #pragma GCC ivdep | ||
1645 | for (integer face = 0; face != NFACE; ++face) { | ||
1646 | std::vector<real>& Uffacefield = Uf[face][field]; | ||
1647 | std::vector<real> const& Uffacerho_i = Uf[face][rho_i]; | ||
1648 | for (integer iii = 0; iii != H_N3; ++iii) { | ||
1649 ![]() | 1.8% | Uffacefield[iii] *= Uffacerho_i[iii]; | |
![]() | |||
1650 | } | ||
1651 | } | ||
1652 | } | ||
1653 | } | ||
1654 | |||
1655 | for (integer i = H_BW - 1; i != H_NX - H_BW + 1; ++i) { | ||
1656 | for (integer j = H_BW - 1; j != H_NX - H_BW + 1; ++j) { | ||
1657 | #pragma GCC ivdep | ||
1658 | for (integer k = H_BW - 1; k != H_NX - H_BW + 1; ++k) { | ||
1659 | const integer iii = hindex(i, j, k); | ||
1660 | for (integer face = 0; face != NFACE; ++face) { | ||
1661 | std::vector<std::vector<real> >& Ufface = Uf[face]; | ||
1662 | real const Uffacerho_iii = Ufface[rho_i][iii]; | ||
1663 | |||
1664 | real x0 = ZERO; | ||
1665 | real y0 = ZERO; | ||
1666 | if (face == FXP) { | ||
1667 | x0 = +HALF * dx; | ||
1668 | } else if (face == FXM) { | ||
1669 | x0 = -HALF * dx; | ||
1670 | } else if (face == FYP) { | ||
1671 | y0 = +HALF * dx; | ||
1672 | } else if (face == FYM) { | ||
1673 | y0 = -HALF * dx; | ||
1674 | } | ||
1675 | |||
1676 | Ufface[sx_i][iii] -= omega * (X[YDIM][iii] + y0) * Uffacerho_iii; | ||
1677 | Ufface[sy_i][iii] += omega * (X[XDIM][iii] + x0) * Uffacerho_iii; | ||
1678 | Ufface[zz_i][iii] += sqr(dx) * omega * Uffacerho_iii / 6.0; | ||
1679 | Ufface[egas_i][iii] += HALF * sqr(Ufface[sx_i][iii]) / Uffacerho_iii; | ||
1680 | Ufface[egas_i][iii] += HALF * sqr(Ufface[sy_i][iii]) / Uffacerho_iii; | ||
1681 | Ufface[egas_i][iii] += HALF * sqr(Ufface[sz_i][iii]) / Uffacerho_iii; | ||
1682 | #ifdef WD_EOS | ||
1683 | Ufface[egas_i][iii] += ztwd_energy(Uffacerho_iii); | ||
1684 | #endif | ||
1685 | } | ||
1686 | } | ||
1687 | } | ||
1688 | } | ||
1689 | PROF_END; | ||
1690 | } | ||
1691 | |||
1692 | real grid::compute_fluxes() { | ||
1693 | PROF_BEGIN; | ||
1694 | const auto& Uf = TLS_Uf(); | ||
1695 | real max_lambda = ZERO; | ||
1696 | std::array<std::vector<real>, NF> ur, ul, f; | ||
1697 | std::vector<space_vector> x; | ||
1698 | |||
1699 | const integer line_sz = H_NX - 2 * H_BW + 1; | ||
1700 | for (integer field = 0; field != NF; ++field) { | ||
1701 | ur[field].resize(line_sz); | ||
1702 | ul[field].resize(line_sz); | ||
1703 | f[field].resize(line_sz); | ||
1704 | } | ||
1705 | x.resize(line_sz); | ||
1706 | |||
1707 | for (integer dim = 0; dim != NDIM; ++dim) { | ||
1708 | |||
1709 | const integer dx_i = dim; | ||
1710 | const integer dy_i = (dim == XDIM ? YDIM : XDIM); | ||
1711 | const integer dz_i = (dim == ZDIM ? YDIM : ZDIM); | ||
1712 | const integer face_p = 2 * dim + 1; | ||
1713 | const integer face_m = 2 * dim; | ||
1714 | |||
1715 | for (integer k = H_BW; k != H_NX - H_BW; ++k) { | ||
1716 | for (integer j = H_BW; j != H_NX - H_BW; ++j) { | ||
1717 | for (integer i = H_BW; i != H_NX - H_BW + 1; ++i) { | ||
1718 | const integer i0 = H_DN[dx_i] * i + H_DN[dy_i] * j + H_DN[dz_i] * k; | ||
1719 | const integer im = i0 - H_DN[dx_i]; | ||
1720 | for (integer field = 0; field != NF; ++field) { | ||
1721 | ur[field][i - H_BW] = Uf[face_m][field][i0]; | ||
1722 | ul[field][i - H_BW] = Uf[face_p][field][im]; | ||
1723 | } | ||
1724 | for (integer d = 0; d != NDIM; ++d) { | ||
1725 | x[i - H_BW][d] = (X[d][i0] + X[d][im]) * HALF; | ||
1726 | } | ||
1727 | } | ||
1728 | const real this_max_lambda = roe_fluxes(f, ul, ur, x, omega, dim, dx); | ||
1729 | max_lambda = std::max(max_lambda, this_max_lambda); | ||
1730 | for (integer field = 0; field != NF; ++field) { | ||
1731 | for (integer i = H_BW; i != H_NX - H_BW + 1; ++i) { | ||
1732 | const integer i0 = F_DN[dx_i] * (i - H_BW) + F_DN[dy_i] * (j - H_BW) + F_DN[dz_i] * (k - H_BW); | ||
1733 ![]() | 0.7% | F[dim][field][i0] = f[field][i - H_BW]; | |
![]() ![]() | |||
1734 | } | ||
1735 | } | ||
1736 | } | ||
1737 | } | ||
1738 | } | ||
1739 | |||
1740 | PROF_END; | ||
1741 | return max_lambda; | ||
1742 | } | ||
1743 | |||
1744 | void grid::set_max_level(integer l) { | ||
1745 | max_level = l; | ||
1746 | } | ||
1747 | |||
1748 | void grid::store() { | ||
1749 | for (integer field = 0; field != NF; ++field) { | ||
1750 | #pragma GCC ivdep | ||
1751 | for (integer i = 0; i != INX; ++i) { | ||
1752 | for (integer j = 0; j != INX; ++j) { | ||
1753 | for (integer k = 0; k != INX; ++k) { | ||
1754 | U0[field][h0index(i, j, k)] = U[field][hindex(i+H_BW,j+H_BW,k+H_BW)]; | ||
1755 | } | ||
1756 | } | ||
1757 | } | ||
1758 | } | ||
1759 | U_out0 = U_out; | ||
1760 | } | ||
1761 | |||
1762 | void grid::set_physical_boundaries(const geo::face& face, real t) { | ||
1763 | const auto dim = face.get_dimension(); | ||
1764 | const auto side = face.get_side(); | ||
1765 | const integer dni = dim == XDIM ? H_DNY : H_DNX; | ||
1766 | const integer dnj = dim == ZDIM ? H_DNY : H_DNZ; | ||
1767 | const integer dnk = H_DN[dim]; | ||
1768 | const integer klb = side == geo::MINUS ? 0 : H_NX - H_BW; | ||
1769 | const integer kub = side == geo::MINUS ? H_BW : H_NX; | ||
1770 | const integer ilb = 0; | ||
1771 | const integer iub = H_NX; | ||
1772 | const integer jlb = 0; | ||
1773 | const integer jub = H_NX; | ||
1774 | |||
1775 | if (opts.problem == SOD) { | ||
1776 | for (integer k = klb; k != kub; ++k) { | ||
1777 | for (integer j = jlb; j != jub; ++j) { | ||
1778 | for (integer i = ilb; i != iub; ++i) { | ||
1779 | const integer iii = i * dni + j * dnj + k * dnk; | ||
1780 | for (integer f = 0; f != NF; ++f) { | ||
1781 | U[f][iii] = 0.0; | ||
1782 | } | ||
1783 | sod_state_t s; | ||
1784 | real x = (X[XDIM][iii] + X[YDIM][iii] + X[ZDIM][iii]) / std::sqrt(3.0); | ||
1785 | exact_sod(&s, &sod_init, x, t); | ||
1786 | U[rho_i][iii] = s.rho; | ||
1787 | U[egas_i][iii] = s.p / (fgamma - 1.0); | ||
1788 | U[sx_i][iii] = s.rho * s.v / std::sqrt(3.0); | ||
1789 | U[sy_i][iii] = s.rho * s.v / std::sqrt(3.0); | ||
1790 | U[sz_i][iii] = s.rho * s.v / std::sqrt(3.0); | ||
1791 | U[tau_i][iii] = std::pow(U[egas_i][iii], 1.0 / fgamma); | ||
1792 | U[egas_i][iii] += s.rho * s.v * s.v / 2.0; | ||
1793 | U[spc_ac_i][iii] = s.rho; | ||
1794 | integer k0 = side == geo::MINUS ? H_BW : H_NX - H_BW - 1; | ||
1795 | U[zx_i][iii] = 0.0; | ||
1796 | U[zy_i][iii] = 0.0; | ||
1797 | U[zz_i][iii] = 0.0; | ||
1798 | } | ||
1799 | } | ||
1800 | } | ||
1801 | } else { | ||
1802 | for (integer field = 0; field != NF; ++field) { | ||
1803 | for (integer k = klb; k != kub; ++k) { | ||
1804 | for (integer j = jlb; j != jub; ++j) { | ||
1805 | for (integer i = ilb; i != iub; ++i) { | ||
1806 | integer k0; | ||
1807 | switch (boundary_types[face]) { | ||
1808 | case REFLECT: | ||
1809 | k0 = side == geo::MINUS ? (2 * H_BW - k - 1) : (2 * (H_NX - H_BW) - k - 1); | ||
1810 | break; | ||
1811 | case OUTFLOW: | ||
1812 | k0 = side == geo::MINUS ? H_BW : H_NX - H_BW - 1; | ||
1813 | break; | ||
1814 | default: | ||
1815 | k0 = -1; | ||
1816 | assert(false); | ||
1817 | abort(); | ||
1818 | } | ||
1819 | const real value = U[field][i * dni + j * dnj + k0 * dnk]; | ||
1820 | const integer iii = i * dni + j * dnj + k * dnk; | ||
1821 | real& ref = U[field][iii]; | ||
1822 | if (field == sx_i + dim) { | ||
1823 | real s0; | ||
1824 | if (field == sx_i) { | ||
1825 | s0 = -omega * X[YDIM][iii] * U[rho_i][iii]; | ||
1826 | } else if (field == sy_i) { | ||
1827 | s0 = +omega * X[XDIM][iii] * U[rho_i][iii]; | ||
1828 | } else { | ||
1829 | s0 = ZERO; | ||
1830 | } | ||
1831 | switch (boundary_types[face]) { | ||
1832 | case REFLECT: | ||
1833 | ref = -value; | ||
1834 | break; | ||
1835 | case OUTFLOW: | ||
1836 | const real before = value; | ||
1837 | if (side == geo::MINUS) { | ||
1838 | ref = s0 + std::min(value - s0, ZERO); | ||
1839 | } else { | ||
1840 | ref = s0 + std::max(value - s0, ZERO); | ||
1841 | } | ||
1842 | const real after = ref; | ||
1843 | assert(rho_i < field); | ||
1844 | assert(egas_i < field); | ||
1845 | real this_rho = U[rho_i][iii]; | ||
1846 | if (this_rho != ZERO) { | ||
1847 | U[egas_i][iii] += HALF * (after * after - before * before) / this_rho; | ||
1848 | } | ||
1849 | break; | ||
1850 | } | ||
1851 | } else { | ||
1852 | ref = +value; | ||
1853 | } | ||
1854 | } | ||
1855 | } | ||
1856 | } | ||
1857 | } | ||
1858 | } | ||
1859 | } | ||
1860 | |||
1861 | void grid::compute_sources(real t) { | ||
1862 | PROF_BEGIN; | ||
1863 | auto& src = dUdt; | ||
1864 | for (integer i = H_BW; i != H_NX - H_BW; ++i) { | ||
1865 | for (integer j = H_BW; j != H_NX - H_BW; ++j) { | ||
1866 | #pragma GCC ivdep | ||
1867 | for (integer k = H_BW; k != H_NX - H_BW; ++k) { | ||
1868 | const integer iii0 = h0index(i - H_BW, j - H_BW, k - H_BW); | ||
1869 | const integer iii = hindex(i, j, k); | ||
1870 | const integer iiif = findex(i - H_BW, j - H_BW, k - H_BW); | ||
1871 | const integer iiig = gindex(i - H_BW, j - H_BW, k - H_BW); | ||
1872 | for (integer field = 0; field != NF; ++field) { | ||
1873 | src[field][iii0] = ZERO; | ||
1874 | } | ||
1875 | const real rho = U[rho_i][iii]; | ||
1876 | src[zx_i][iii0] = (-(F[YDIM][sz_i][iiif + F_DNY] + F[YDIM][sz_i][iiif]) + (F[ZDIM][sy_i][iiif + F_DNZ] + F[ZDIM][sy_i][iiif])) * HALF; | ||
1877 | src[zy_i][iii0] = (+(F[XDIM][sz_i][iiif + F_DNX] + F[XDIM][sz_i][iiif]) - (F[ZDIM][sx_i][iiif + F_DNZ] + F[ZDIM][sx_i][iiif])) * HALF; | ||
1878 | src[zz_i][iii0] = (-(F[XDIM][sy_i][iiif + F_DNX] + F[XDIM][sy_i][iiif]) + (F[YDIM][sx_i][iiif + F_DNY] + F[YDIM][sx_i][iiif])) * HALF; | ||
1879 | if (node_server::is_gravity_on()) { | ||
1880 | src[sx_i][iii0] += rho * G[iiig][gx_i]; | ||
1881 | src[sy_i][iii0] += rho * G[iiig][gy_i]; | ||
1882 | src[sz_i][iii0] += rho * G[iiig][gz_i]; | ||
1883 | } | ||
1884 | src[sx_i][iii0] += omega * U[sy_i][iii]; | ||
1885 | src[sy_i][iii0] -= omega * U[sx_i][iii]; | ||
1886 | if (node_server::is_gravity_on()) { | ||
1887 | src[egas_i][iii0] -= omega * X[YDIM][iii] * rho * G[iiig][gx_i]; | ||
1888 | src[egas_i][iii0] += omega * X[XDIM][iii] * rho * G[iiig][gy_i]; | ||
1889 | } | ||
1890 | #ifdef USE_DRIVING | ||
1891 | const real period = (2.0 * M_PI / grid::omega); | ||
1892 | if (t < DRIVING_TIME * period) { | ||
1893 | const real ff = -DRIVING_RATE / period; | ||
1894 | const real rho = U[rho_i][iii]; | ||
1895 | const real sx = U[sx_i][iii]; | ||
1896 | const real sy = U[sy_i][iii]; | ||
1897 | const real zz = U[zz_i][iii]; | ||
1898 | const real x = X[XDIM][iii]; | ||
1899 | const real y = X[YDIM][iii]; | ||
1900 | const real R = std::sqrt(x * x + y * y); | ||
1901 | const real lz = (x * sy - y * sx); | ||
1902 | const real dsx = -y / R / R * lz * ff; | ||
1903 | const real dsy = +x / R / R * lz * ff; | ||
1904 | src[sx_i][iii0] += dsx; | ||
1905 | src[sy_i][iii0] += dsy; | ||
1906 | src[egas_i][iii0] += (sx * dsx + sy * dsy) / rho; | ||
1907 | src[zz_i][iii0] += ff * zz; | ||
1908 | |||
1909 | } | ||
1910 | #endif | ||
1911 | |||
1912 | } | ||
1913 | } | ||
1914 | } | ||
1915 | PROF_END; | ||
1916 | } | ||
1917 | |||
1918 | void grid::compute_dudt() { | ||
1919 | PROF_BEGIN; | ||
1920 | for (integer i = H_BW; i != H_NX - H_BW; ++i) { | ||
1921 | for (integer j = H_BW; j != H_NX - H_BW; ++j) { | ||
1922 | for (integer field = 0; field != NF; ++field) { | ||
1923 | #pragma GCC ivdep | ||
1924 | for (integer k = H_BW; k != H_NX - H_BW; ++k) { | ||
1925 | const integer iii0 = h0index(i - H_BW, j - H_BW, k - H_BW); | ||
1926 | const integer iiif = findex(i - H_BW, j - H_BW, k - H_BW); | ||
1927 | dUdt[field][iii0] -= (F[XDIM][field][iiif + F_DNX] - F[XDIM][field][iiif]) / dx; | ||
1928 | dUdt[field][iii0] -= (F[YDIM][field][iiif + F_DNY] - F[YDIM][field][iiif]) / dx; | ||
1929 | dUdt[field][iii0] -= (F[ZDIM][field][iiif + F_DNZ] - F[ZDIM][field][iiif]) / dx; | ||
1930 | } | ||
1931 | } | ||
1932 | if (node_server::is_gravity_on()) { | ||
1933 | |||
1934 | #pragma GCC ivdep | ||
1935 | for (integer k = H_BW; k != H_NX - H_BW; ++k) { | ||
1936 | const integer iii0 = h0index(i - H_BW, j - H_BW, k - H_BW); | ||
1937 | dUdt[egas_i][iii0] += dUdt[pot_i][iii0]; | ||
1938 | dUdt[pot_i][iii0] = ZERO; | ||
1939 | } | ||
1940 | } | ||
1941 | #pragma GCC ivdep | ||
1942 | for (integer k = H_BW; k != H_NX - H_BW; ++k) { | ||
1943 | const integer iii0 = h0index(i - H_BW, j - H_BW, k - H_BW); | ||
1944 | const integer iiig = gindex(i - H_BW, j - H_BW, k - H_BW); | ||
1945 | if (node_server::is_gravity_on()) { | ||
1946 | dUdt[egas_i][iii0] -= (dUdt[rho_i][iii0] * G[iiig][phi_i]) * HALF; | ||
1947 | } | ||
1948 | } | ||
1949 | } | ||
1950 | } | ||
1951 | PROF_END; | ||
1952 | // solve_gravity(DRHODT); | ||
1953 | } | ||
1954 | |||
1955 | void grid::egas_to_etot() { | ||
1956 | PROF_BEGIN; | ||
1957 | if (node_server::is_gravity_on()) { | ||
1958 | |||
1959 | for (integer i = H_BW; i != H_NX - H_BW; ++i) { | ||
1960 | for (integer j = H_BW; j != H_NX - H_BW; ++j) { | ||
1961 | #pragma GCC ivdep | ||
1962 | for (integer k = H_BW; k != H_NX - H_BW; ++k) { | ||
1963 | const integer iii = hindex(i, j, k); | ||
1964 | U[egas_i][iii] += U[pot_i][iii] * HALF; | ||
1965 | } | ||
1966 | } | ||
1967 | } | ||
1968 | } | ||
1969 | PROF_END; | ||
1970 | } | ||
1971 | |||
1972 | void grid::etot_to_egas() { | ||
1973 | PROF_BEGIN; | ||
1974 | if (node_server::is_gravity_on()) { | ||
1975 | |||
1976 | for (integer i = H_BW; i != H_NX - H_BW; ++i) { | ||
1977 | for (integer j = H_BW; j != H_NX - H_BW; ++j) { | ||
1978 | #pragma GCC ivdep | ||
1979 | for (integer k = H_BW; k != H_NX - H_BW; ++k) { | ||
1980 | const integer iii = hindex(i, j, k); | ||
1981 | U[egas_i][iii] -= U[pot_i][iii] * HALF; | ||
1982 | } | ||
1983 | } | ||
1984 | } | ||
1985 | } | ||
1986 | PROF_END; | ||
1987 | } | ||
1988 | |||
1989 | void grid::next_u(integer rk, real t, real dt) { | ||
1990 | PROF_BEGIN; | ||
1991 | for (integer i = H_BW; i != H_NX - H_BW; ++i) { | ||
1992 | for (integer j = H_BW; j != H_NX - H_BW; ++j) { | ||
1993 | #pragma GCC ivdep | ||
1994 | for (integer k = H_BW; k != H_NX - H_BW; ++k) { | ||
1995 | const integer iii0 = h0index(i - H_BW, j - H_BW, k - H_BW); | ||
1996 | const integer iii = hindex(i, j, k); | ||
1997 | dUdt[egas_i][iii0] += (dphi_dt[iii0] * U[rho_i][iii]) * HALF; | ||
1998 | dUdt[zx_i][iii0] -= omega * X[ZDIM][iii] * U[sx_i][iii]; | ||
1999 | dUdt[zy_i][iii0] -= omega * X[ZDIM][iii] * U[sy_i][iii]; | ||
2000 | dUdt[zz_i][iii0] += omega * (X[XDIM][iii] * U[sx_i][iii] + X[YDIM][iii] * U[sy_i][iii]); | ||
2001 | } | ||
2002 | } | ||
2003 | } | ||
2004 | |||
2005 | std::vector < real > du_out(NF, ZERO); | ||
2006 | |||
2007 | std::vector < real > ds(NDIM, ZERO); | ||
2008 | for (integer i = H_BW; i != H_NX - H_BW; ++i) { | ||
2009 | for (integer j = H_BW; j != H_NX - H_BW; ++j) { | ||
2010 | #pragma GCC ivdep | ||
2011 | for (integer k = H_BW; k != H_NX - H_BW; ++k) { | ||
2012 | const integer iii = hindex(i, j, k); | ||
2013 | const integer iii0 = h0index(i - H_BW, j - H_BW, k - H_BW); | ||
2014 | for (integer field = 0; field != NF; ++field) { | ||
2015 | const real u1 = U[field][iii] + dUdt[field][iii0] * dt; | ||
2016 | const real u0 = U0[field][h0index(i - H_BW, j - H_BW, k - H_BW)]; | ||
2017 | U[field][iii] = (ONE - rk_beta[rk]) * u0 + rk_beta[rk] * u1; | ||
2018 | } | ||
2019 | } | ||
2020 | } | ||
2021 | } | ||
2022 | |||
2023 | du_out[sx_i] += omega * U_out[sy_i] * dt; | ||
2024 | du_out[sy_i] -= omega * U_out[sx_i] * dt; | ||
2025 | |||
2026 | for (integer i = H_BW; i != H_NX - H_BW; ++i) { | ||
2027 | #pragma GCC ivdep | ||
2028 | for (integer j = H_BW; j != H_NX - H_BW; ++j) { | ||
2029 | const real dx2 = sqr(dx); | ||
2030 | const integer iii_p0 = findex(INX, i - H_BW, j - H_BW); | ||
2031 | const integer jjj_p0 = findex(j - H_BW, INX, i - H_BW); | ||
2032 | const integer kkk_p0 = findex(i - H_BW, j - H_BW, INX); | ||
2033 | const integer iii_m0 = findex(0, i - H_BW, j - H_BW); | ||
2034 | const integer jjj_m0 = findex(j - H_BW, 0, i - H_BW); | ||
2035 | const integer kkk_m0 = findex(i - H_BW, j - H_BW, 0); | ||
2036 | const integer iii_p = H_DNX * (H_NX - H_BW) + H_DNY * i + H_DNZ * j; | ||
2037 | const integer jjj_p = H_DNY * (H_NX - H_BW) + H_DNZ * i + H_DNX * j; | ||
2038 | const integer kkk_p = H_DNZ * (H_NX - H_BW) + H_DNX * i + H_DNY * j; | ||
2039 | const integer iii_m = H_DNX * (H_BW) + H_DNY * i + H_DNZ * j; | ||
2040 | const integer jjj_m = H_DNY * (H_BW) + H_DNZ * i + H_DNX * j; | ||
2041 | const integer kkk_m = H_DNZ * (H_BW) + H_DNX * i + H_DNY * j; | ||
2042 | std::vector < real > du(NF); | ||
2043 | for (integer field = 0; field != NF; ++field) { | ||
2044 | // if (field < zx_i || field > zz_i) { | ||
2045 | du[field] = ZERO; | ||
2046 | if (X[XDIM][iii_p] + pivot[XDIM] > scaling_factor) { | ||
2047 | du[field] += (F[XDIM][field][iii_p0]) * dx2; | ||
2048 | } | ||
2049 | if (X[YDIM][jjj_p] + pivot[YDIM] > scaling_factor) { | ||
2050 | du[field] += (F[YDIM][field][jjj_p0]) * dx2; | ||
2051 | } | ||
2052 | if (X[ZDIM][kkk_p] + pivot[ZDIM] > scaling_factor) { | ||
2053 | du[field] += (F[ZDIM][field][kkk_p0]) * dx2; | ||
2054 | } | ||
2055 | if (X[XDIM][iii_m] + pivot[XDIM] < -scaling_factor + dx) { | ||
2056 | du[field] += (-F[XDIM][field][iii_m0]) * dx2; | ||
2057 | } | ||
2058 | if (X[YDIM][jjj_m] + pivot[YDIM] < -scaling_factor + dx) { | ||
2059 | du[field] += (-F[YDIM][field][jjj_m0]) * dx2; | ||
2060 | } | ||
2061 | if (X[ZDIM][kkk_m] + pivot[ZDIM] < -scaling_factor + dx) { | ||
2062 | du[field] += (-F[ZDIM][field][kkk_m0]) * dx2; | ||
2063 | } | ||
2064 | // } | ||
2065 | } | ||
2066 | |||
2067 | if (X[XDIM][iii_p] + pivot[XDIM] > scaling_factor) { | ||
2068 | const real xp = X[XDIM][iii_p] - HALF * dx; | ||
2069 | du[zx_i] += (X[YDIM][iii_p] * F[XDIM][sz_i][iii_p0]) * dx2; | ||
2070 | du[zx_i] -= (X[ZDIM][iii_p] * F[XDIM][sy_i][iii_p0]) * dx2; | ||
2071 | du[zy_i] -= (xp * F[XDIM][sz_i][iii_p0]) * dx2; | ||
2072 | du[zy_i] += (X[ZDIM][iii_p] * F[XDIM][sx_i][iii_p0]) * dx2; | ||
2073 | du[zz_i] += (xp * F[XDIM][sy_i][iii_p0]) * dx2; | ||
2074 | du[zz_i] -= (X[YDIM][iii_p] * F[XDIM][sx_i][iii_p0]) * dx2; | ||
2075 | } | ||
2076 | if (X[YDIM][jjj_p] + pivot[YDIM] > scaling_factor) { | ||
2077 | const real yp = X[YDIM][jjj_p] - HALF * dx; | ||
2078 | du[zx_i] += (yp * F[YDIM][sz_i][jjj_p0]) * dx2; | ||
2079 | du[zx_i] -= (X[ZDIM][jjj_p] * F[YDIM][sy_i][jjj_p0]) * dx2; | ||
2080 | du[zy_i] -= (X[XDIM][jjj_p] * F[YDIM][sz_i][jjj_p0]) * dx2; | ||
2081 | du[zy_i] += (X[ZDIM][jjj_p] * F[YDIM][sx_i][jjj_p0]) * dx2; | ||
2082 | du[zz_i] += (X[XDIM][jjj_p] * F[YDIM][sy_i][jjj_p0]) * dx2; | ||
2083 | du[zz_i] -= (yp * F[YDIM][sx_i][jjj_p0]) * dx2; | ||
2084 | } | ||
2085 | if (X[ZDIM][kkk_p] + pivot[ZDIM] > scaling_factor) { | ||
2086 | const real zp = X[ZDIM][kkk_p] - HALF * dx; | ||
2087 | du[zx_i] -= (zp * F[ZDIM][sy_i][kkk_p0]) * dx2; | ||
2088 | du[zx_i] += (X[YDIM][kkk_p] * F[ZDIM][sz_i][kkk_p0]) * dx2; | ||
2089 | du[zy_i] += (zp * F[ZDIM][sx_i][kkk_p0]) * dx2; | ||
2090 | du[zy_i] -= (X[XDIM][kkk_p] * F[ZDIM][sz_i][kkk_p0]) * dx2; | ||
2091 | du[zz_i] += (X[XDIM][kkk_p] * F[ZDIM][sy_i][kkk_p0]) * dx2; | ||
2092 | du[zz_i] -= (X[YDIM][kkk_p] * F[ZDIM][sx_i][kkk_p0]) * dx2; | ||
2093 | } | ||
2094 | |||
2095 | if (X[XDIM][iii_m] + pivot[XDIM] < -scaling_factor + dx) { | ||
2096 | const real xm = X[XDIM][iii_m] - HALF * dx; | ||
2097 | du[zx_i] += (-X[YDIM][iii_m] * F[XDIM][sz_i][iii_m0]) * dx2; | ||
2098 | du[zx_i] -= (-X[ZDIM][iii_m] * F[XDIM][sy_i][iii_m0]) * dx2; | ||
2099 | du[zy_i] -= (-xm * F[XDIM][sz_i][iii_m0]) * dx2; | ||
2100 | du[zy_i] += (-X[ZDIM][iii_m] * F[XDIM][sx_i][iii_m0]) * dx2; | ||
2101 | du[zz_i] += (-xm * F[XDIM][sy_i][iii_m0]) * dx2; | ||
2102 | du[zz_i] -= (-X[YDIM][iii_m] * F[XDIM][sx_i][iii_m0]) * dx2; | ||
2103 | } | ||
2104 | if (X[YDIM][jjj_m] + pivot[YDIM] < -scaling_factor + dx) { | ||
2105 | const real ym = X[YDIM][jjj_m] - HALF * dx; | ||
2106 | du[zx_i] -= (-X[ZDIM][jjj_m] * F[YDIM][sy_i][jjj_m0]) * dx2; | ||
2107 | du[zx_i] += (-ym * F[YDIM][sz_i][jjj_m0]) * dx2; | ||
2108 | du[zy_i] -= (-X[XDIM][jjj_m] * F[YDIM][sz_i][jjj_m0]) * dx2; | ||
2109 | du[zy_i] += (-X[ZDIM][jjj_m] * F[YDIM][sx_i][jjj_m0]) * dx2; | ||
2110 | du[zz_i] += (-X[XDIM][jjj_m] * F[YDIM][sy_i][jjj_m0]) * dx2; | ||
2111 | du[zz_i] -= (-ym * F[YDIM][sx_i][jjj_m0]) * dx2; | ||
2112 | } | ||
2113 | if (X[ZDIM][kkk_m] + pivot[ZDIM] < -scaling_factor + dx) { | ||
2114 | const real zm = X[ZDIM][kkk_m] - HALF * dx; | ||
2115 | du[zx_i] -= (-zm * F[ZDIM][sy_i][kkk_m0]) * dx2; | ||
2116 | du[zx_i] += (-X[YDIM][kkk_m] * F[ZDIM][sz_i][kkk_m0]) * dx2; | ||
2117 | du[zy_i] += (-zm * F[ZDIM][sx_i][kkk_m0]) * dx2; | ||
2118 | du[zy_i] -= (-X[XDIM][kkk_m] * F[ZDIM][sz_i][kkk_m0]) * dx2; | ||
2119 | du[zz_i] += (-X[XDIM][kkk_m] * F[ZDIM][sy_i][kkk_m0]) * dx2; | ||
2120 | du[zz_i] -= (-X[YDIM][kkk_m] * F[ZDIM][sx_i][kkk_m0]) * dx2; | ||
2121 | } | ||
2122 | for (integer field = 0; field != NF; ++field) { | ||
2123 | du_out[field] += du[field] * dt; | ||
2124 | } | ||
2125 | } | ||
2126 | } | ||
2127 | #pragma GCC ivdep | ||
2128 | for (integer field = 0; field != NF; ++field) { | ||
2129 | const real out1 = U_out[field] + du_out[field]; | ||
2130 | const real out0 = U_out0[field]; | ||
2131 | U_out[field] = (ONE - rk_beta[rk]) * out0 + rk_beta[rk] * out1; | ||
2132 | } | ||
2133 | for (integer i = H_BW; i != H_NX - H_BW; ++i) { | ||
2134 | for (integer j = H_BW; j != H_NX - H_BW; ++j) { | ||
2135 | #pragma GCC ivdep | ||
2136 | for (integer k = H_BW; k != H_NX - H_BW; ++k) { | ||
2137 | const integer iii = hindex(i, j, k); | ||
2138 | if (opts.problem == SOD) { | ||
2139 | U[zx_i][iii] = U[zy_i][iii] = U[zz_i][iii] = 0.0; | ||
2140 | } | ||
2141 | U[rho_i][iii] = ZERO; | ||
2142 | for (integer si = 0; si != NSPECIES; ++si) { | ||
2143 | U[rho_i][iii] += U[spc_i + si][iii]; | ||
2144 | } | ||
2145 | if (U[tau_i][iii] < ZERO) { | ||
2146 | printf("Tau is negative- %e\n", double(U[tau_i][iii])); | ||
2147 | // abort(); | ||
2148 | } else if (U[rho_i][iii] <= ZERO) { | ||
2149 | printf("Rho is non-positive - %e %i %i %i\n", double(U[rho_i][iii]), int(i), int(j), int(k)); | ||
2150 | // abort(); | ||
2151 | } | ||
2152 | if (!opts.ang_con) { | ||
2153 | U[zx_i][iii] = U[zy_i][iii] = U[zx_i][iii] = 0.0; | ||
2154 | } | ||
2155 | } | ||
2156 | } | ||
2157 | } | ||
2158 | PROF_END; | ||
2159 | } | ||
2160 | |||
2161 | void grid::dual_energy_update() { | ||
2162 | PROF_BEGIN; | ||
2163 | // bool in_bnd; | ||
2164 | for (integer i = H_BW; i != H_NX - H_BW; ++i) { | ||
2165 | for (integer j = H_BW; j != H_NX - H_BW; ++j) { | ||
2166 | #pragma GCC ivdep | ||
2167 | for (integer k = H_BW; k != H_NX - H_BW; ++k) { | ||
2168 | const integer iii = hindex(i, j, k); | ||
2169 | real ek = ZERO; | ||
2170 | ek += HALF * pow(U[sx_i][iii], 2) / U[rho_i][iii]; | ||
2171 | ek += HALF * pow(U[sy_i][iii], 2) / U[rho_i][iii]; | ||
2172 | ek += HALF * pow(U[sz_i][iii], 2) / U[rho_i][iii]; | ||
2173 | real ei = U[egas_i][iii] - ek | ||
2174 | #ifdef WD_EOS | ||
2175 | - ztwd_energy(U[rho_i][iii]) | ||
2176 | #endif | ||
2177 | ; | ||
2178 | real et = U[egas_i][iii]; | ||
2179 | et = std::max(et, U[egas_i][iii + H_DNX]); | ||
2180 | et = std::max(et, U[egas_i][iii - H_DNX]); | ||
2181 | et = std::max(et, U[egas_i][iii + H_DNY]); | ||
2182 | et = std::max(et, U[egas_i][iii - H_DNY]); | ||
2183 | et = std::max(et, U[egas_i][iii + H_DNZ]); | ||
2184 | et = std::max(et, U[egas_i][iii - H_DNZ]); | ||
2185 | if (ei > de_switch1 * et) { | ||
2186 | U[tau_i][iii] = std::pow(ei, ONE / fgamma); | ||
2187 | } | ||
2188 | } | ||
2189 | } | ||
2190 | } | ||
2191 | PROF_END; | ||
2192 | } | ||
2193 | |||
2194 | std::vector<real> grid::conserved_outflows() const { | ||
2195 | auto Uret = U_out; | ||
2196 | if (node_server::is_gravity_on()) { | ||
2197 | Uret[egas_i] += Uret[pot_i]; | ||
2198 | } | ||
2199 | return Uret; | ||
2200 | } | ||
2201 |
Copyright (c) 2006-2012 Rogue Wave Software, Inc. All Rights Reserved.
Patents pending.