diff --git a/Examples/PerfectFluid/Fluxes.hpp b/Examples/PerfectFluid/Fluxes.hpp index 385574335..3c7dcfffe 100644 --- a/Examples/PerfectFluid/Fluxes.hpp +++ b/Examples/PerfectFluid/Fluxes.hpp @@ -13,50 +13,57 @@ namespace Fluxes { -// Primitive to conservative variables template class vars_t> -vars_t compute_flux(const vars_t &vars, const int idir, - const double lambda) +vars_t compute_flux(const vars_t &vars, const int idir) { const auto h_UU = TensorAlgebra::compute_inverse_sym(vars.h); vars_t out; - out.D = (vars.lapse * vars.vi[idir] - vars.shift[idir]) * vars.D - - lambda * vars.D; + out.D = (vars.lapse * vars.vi[idir] - vars.shift[idir]) * vars.D; Tensor<1, data_t> Sj_U; - FOR(i) - { - Sj_U[i] = 0; - FOR(j) Sj_U[i] += vars.chi * h_UU[i][j] * vars.Sj[j]; + FOR(i) { + Sj_U[i] = 0; + FOR(j) Sj_U[i] += vars.chi * h_UU[i][j] * vars.Sj[j]; } - + data_t chi_regularised = simd_max(1e-6, vars.chi); Tensor<1, data_t> vi_D; FOR(i) { vi_D[i] = 0.; - FOR(j) { vi_D[i] += vars.h[i][j] * vars.vi[j] / chi_regularised; } + FOR(j) + { + vi_D[i] += vars.h[i][j] * vars.vi[j] / chi_regularised; + } } - + data_t v2 = 0.; FOR(i) v2 += vars.vi[i] * vi_D[i]; - - data_t P = - vars.rho * (1. + vars.eps) / 3.; // for now we assume a conformal fluid + + data_t P = vars.rho * (1. + vars.eps) / 3.; //for now we assume a conformal fluid data_t WW = 1. / (1. - v2); data_t hh = 1. + vars.eps + P / vars.rho; - - FOR(j) + + FOR(j) { - out.Sj[j] = vars.lapse * vars.rho * hh * WW * vars.vi[idir] * vi_D[j] - - vars.shift[idir] * vars.Sj[j] - lambda * vars.Sj[j]; - FOR(k) out.Sj[j] += vars.lapse * P * h_UU[idir][k] * vars.h[j][k]; + out.Sj[j] = vars.lapse * vars.rho * hh * WW * vars.vi[idir] * vi_D[j] - + vars.shift[idir] * vars.Sj[j]; + FOR(k) out.Sj[j] += vars.lapse * P * h_UU[idir][k] * vars.h[j][k]; } - + out.tau = vars.lapse * (Sj_U[idir] - vars.D * vars.vi[idir]) - - vars.shift[idir] * vars.tau - lambda * vars.tau; - + vars.shift[idir] * vars.tau; + + return out; +} +template class vars_t> +vars_t compute_num_flux(const vars_t &vars, const int idir, const double lambda) { + vars_t out = compute_flux(vars, idir); + out.D -= lambda * vars.D; + FOR(j) out.Sj[j] -= lambda * vars.Sj[j]; + out.tau -= lambda * vars.tau; return out; + } -} // namespace Fluxes +} #endif /* FLUXES_HPP_ */ diff --git a/Examples/PerfectFluid/InitialFluidData.hpp b/Examples/PerfectFluid/InitialFluidData.hpp index 7955396b2..41e27d35a 100644 --- a/Examples/PerfectFluid/InitialFluidData.hpp +++ b/Examples/PerfectFluid/InitialFluidData.hpp @@ -11,19 +11,19 @@ #include "Coordinates.hpp" #include "FluidCCZ4RHS.hpp" #include "PerfectFluid.hpp" -#include "PrimitiveRecovery.hpp" #include "Tensor.hpp" #include "UserVariables.hpp" //This files needs NUM_VARS - total no. components #include "VarsTools.hpp" #include "simd.hpp" +#include "PrimitiveRecovery.hpp" //! Class which sets the initial fluid matter config // aka Kevin data class InitialFluidData { - // Now the non grid ADM vars - template - using MetricVars = ADMConformalVars::VarsWithGauge; + // Now the non grid ADM vars + template + using MetricVars = ADMConformalVars::VarsWithGauge; public: //! A structure for the input params for fluid properties and initial @@ -31,18 +31,19 @@ class InitialFluidData struct params_t { std::array - center; //!< Centre of perturbation in initial SF bubble + center; //!< Centre of perturbation in initial SF bubble double rho0; double uflow; double amplitude; double awidth; double sigma; - std::array ycenter; + std::array + ycenter; }; //! The constructor - InitialFluidData(params_t a_params, double a_dx) - : m_dx(a_dx), m_params(a_params) + InitialFluidData(params_t a_params, double a_dx) + : m_dx(a_dx), m_params(a_params) { } @@ -52,59 +53,47 @@ class InitialFluidData // where am i? Coordinates coords(current_cell, m_dx, m_params.center); - MetricVars metric_vars; - - data_t x = coords.x; - double y = coords.y; - double z = coords.z; - - data_t vx, vy, vz, eps; - data_t nn; - - vx = m_params.uflow * - (tanh((y - m_params.ycenter[0]) / m_params.awidth) - - tanh((y - m_params.ycenter[1]) / m_params.awidth) - 1.); - vy = m_params.amplitude * sin(2. * M_PI * x) * - (exp(-pow((y - m_params.ycenter[0]) / m_params.sigma, 2)) + - exp(-pow((y - m_params.ycenter[1]) / m_params.sigma, 2))); - vz = 0.; - - // DT: Thought eps was meant to be 0? - nn = 1. + 0.5 * (tanh((y - m_params.ycenter[0]) / m_params.awidth) - - tanh((y - m_params.ycenter[1]) / m_params.awidth)); - eps = 0.; - // gamma_fac = 1./sqrt(1. - vx*vx - vy*vy - vz*vz); - // ux = gamma_fac * vx; - // uy = gamma_fac * vy; - // uz = gamma_fac * vz; + MetricVars metric_vars; + + data_t x = coords.x; + double y = coords.y; + double z = coords.z; + + Tensor<1, data_t> vi, Sj; + data_t chi_regularised = simd_max(metric_vars.chi, 1e-6); + + vi[0] = m_params.uflow * (tanh((y - m_params.ycenter[0]) / m_params.awidth) + - tanh((y - m_params.ycenter[1]) / m_params.awidth) - 1.); + vi[1] = m_params.amplitude * sin(2. * M_PI * x) + * (exp(-pow((y - m_params.ycenter[0]) / m_params.sigma,2)) + + exp(-pow((y - m_params.ycenter[1]) / m_params.sigma,2))); + vi[2] = 0.; + + data_t nn = 1. + 0.5 * (tanh((y - m_params.ycenter[0]) / m_params.awidth) + - tanh((y - m_params.ycenter[1]) / m_params.awidth)); + data_t eps = 0.; data_t rho = m_params.rho0; - data_t v2 = vx * vx + vy * vy + vz * vz; - data_t P = rho * (1. + eps) / 3.; - data_t WW = 1. / (1. - v2); + data_t v2 = 0.; + FOR(i,j) v2 += metric_vars.h[i][j] * vi[i] * vi[j] / chi_regularised; + data_t P = rho * (1. + eps) / 3.; + data_t WW = 1. / (1. - v2); data_t hh = 1. + eps + P / rho; - - data_t D = rho * sqrt(WW); - data_t tau = rho * hh * WW - P - D; - data_t Sj1 = rho * hh * WW * vx; - data_t Sj2 = rho * hh * WW * vy; - data_t Sj3 = rho * hh * WW * vz; - + + data_t D = rho * sqrt(WW); + data_t tau = rho * hh * WW - P - D; + FOR(i) Sj[i] = rho * hh * WW * vi[i]; + // store the vars - current_cell.store_vars(rho, c_rho); - current_cell.store_vars(vx, c_vi1); - current_cell.store_vars(vy, c_vi2); - current_cell.store_vars(vz, c_vi3); + current_cell.store_vars(rho, c_rho); + current_cell.store_vars(vi, GRInterval()); current_cell.store_vars(D, c_D); - current_cell.store_vars(Sj1, c_Sj1); - current_cell.store_vars(Sj2, c_Sj2); - current_cell.store_vars(Sj3, c_Sj3); - current_cell.store_vars(tau, c_tau); + current_cell.store_vars(Sj, GRInterval()); + current_cell.store_vars(tau, c_tau); } - - protected: - double m_dx; - const params_t m_params; +protected: + double m_dx; + const params_t m_params; }; #endif /* INITIALFLUIDDATA_HPP_ */ diff --git a/Examples/PerfectFluid/PerfectFluid.impl.hpp b/Examples/PerfectFluid/PerfectFluid.impl.hpp index ee00e3ac4..beaf910a2 100644 --- a/Examples/PerfectFluid/PerfectFluid.impl.hpp +++ b/Examples/PerfectFluid/PerfectFluid.impl.hpp @@ -78,9 +78,25 @@ void PerfectFluid::add_matter_rhs( vars_t source = Sources::compute_source(vars, d1); // evolution equations for the fluid conservative variables - rhs.D = source.D; - FOR(i) rhs.Sj[i] = source.Sj[i]; - rhs.tau = source.tau; + data_t divshift = TensorAlgebra::compute_trace(d1.shift); + data_t chi_regularised = simd_max(vars.chi, 1e-6); + data_t advec_chi = 0.; + FOR(i) advec_chi += vars.shift[i] * d1.chi[i] / chi_regularised; + // source - vars * \partial_t\sqrt{\gamma}/\sqrt{\gamma} + rhs.D = source.D - vars.D * (vars.lapse * vars.K + - divshift + GR_SPACEDIM / 2. * advec_chi); + FOR(i) rhs.Sj[i] = source.Sj[i] - vars.Sj[i] * (vars.lapse * + vars.K - divshift + GR_SPACEDIM / 2. * advec_chi); + rhs.tau = source.tau - vars.tau * (vars.lapse * vars.K + - divshift + GR_SPACEDIM / 2. * advec_chi); + + // - F^i\partial_i\sqrt{\gamma}/\sqrt{\gamma} + FOR(i) { + vars_t flux = Fluxes::compute_flux(vars, i); + rhs.D -= GR_SPACEDIM / 2. * d1.chi[i] / chi_regularised * flux.D; + FOR(j) rhs.Sj[j] -= GR_SPACEDIM / 2. * d1.chi[i] / chi_regularised * flux.Sj[j]; + rhs.tau -= GR_SPACEDIM / 2. * d1.chi[i] / chi_regularised * flux.tau; + } FOR(idir) { @@ -90,7 +106,7 @@ void PerfectFluid::add_matter_rhs( FOR(j) { vars_right_p.vi[j] = rp.vi[j][idir]; } PrimitiveRecovery::PtoC(vars_right_p); vars_t flux_right_p = - Fluxes::compute_flux(vars_right_p, idir, m_lambda); + Fluxes::compute_num_flux(vars_right_p, idir, m_lambda); vars_t vars_right_m; vars_right_m.rho = rm.rho[idir]; @@ -98,7 +114,7 @@ void PerfectFluid::add_matter_rhs( FOR(j) vars_right_m.vi[j] = rm.vi[j][idir]; PrimitiveRecovery::PtoC(vars_right_m); vars_t flux_right_m = - Fluxes::compute_flux(vars_right_m, idir, m_lambda); + Fluxes::compute_num_flux(vars_right_m, idir, m_lambda); rhs.D += -1. / (2. * m_dx) * (flux_right_p.D + flux_right_m.D); FOR(j) @@ -117,7 +133,7 @@ void PerfectFluid::add_matter_rhs( FOR(j) { vars_left_p.vi[j] = lp.vi[j][idir]; } PrimitiveRecovery::PtoC(vars_left_p); vars_t flux_left_p = - Fluxes::compute_flux(vars_left_p, idir, m_lambda); + Fluxes::compute_num_flux(vars_left_p, idir, m_lambda); vars_t vars_left_m; vars_left_m.rho = lm.rho[idir]; @@ -125,7 +141,7 @@ void PerfectFluid::add_matter_rhs( FOR(j) { vars_left_m.vi[j] = lm.vi[j][idir]; } PrimitiveRecovery::PtoC(vars_left_m); vars_t flux_left_m = - Fluxes::compute_flux(vars_left_m, idir, m_lambda); + Fluxes::compute_num_flux(vars_left_m, idir, m_lambda); rhs.D += 1. / (2. * m_dx) * (flux_left_p.D + flux_left_m.D); FOR(j)