/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.