/home/users/khuck/src/octotiger/src/grid.cpp

Line% of fetchesSource
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.