diff --git a/Examples/Fluid_Kerr_merged/ConservedQuantities.hpp b/Examples/Fluid_Kerr_merged/ConservedQuantities.hpp index 47dfb2574..50ddbf4cc 100644 --- a/Examples/Fluid_Kerr_merged/ConservedQuantities.hpp +++ b/Examples/Fluid_Kerr_merged/ConservedQuantities.hpp @@ -31,11 +31,12 @@ void PtoC(const data_t P_of_rho, vars_t &vars) data_t WW = 1. / (1. - v2); data_t hh = 1. + vars.eps + P_of_rho / vars.rho; - vars.D = vars.rho * sqrt(WW); - vars.tau = vars.rho * hh * WW - P_of_rho - vars.D; + data_t rho_conformal = vars.rho / pow(chi_regularised, 1.5); + vars.D = rho_conformal * sqrt(WW); + vars.tau = rho_conformal * hh * WW - P_of_rho - vars.D; // S_j (note lower index) = - n^a T_ai - FOR(i) { vars.Sj[i] = vars.rho * hh * WW * vi_D[i]; } + FOR(i) { vars.Sj[i] = rho_conformal * hh * WW * vi_D[i]; } } } // namespace ConservedQuantities diff --git a/Examples/Fluid_Kerr_merged/EoS.hpp b/Examples/Fluid_Kerr_merged/EoS.hpp index 8aa3a2303..bd7e394dc 100644 --- a/Examples/Fluid_Kerr_merged/EoS.hpp +++ b/Examples/Fluid_Kerr_merged/EoS.hpp @@ -19,9 +19,9 @@ class EoS void compute_eos(data_t &P_of_rho, const vars_t &vars) const { // The pressure value as a function of rho - const double K = 1./3.;//100.; + const double K = 1. / 3.; // 100.; const double n = 1.; - const double Gamma = 1.;// + 1. / n; + const double Gamma = 1.; // + 1. / n; P_of_rho = K * pow(vars.rho, Gamma); } }; diff --git a/Examples/Fluid_Kerr_merged/Fluxes.hpp b/Examples/Fluid_Kerr_merged/Fluxes.hpp index 9114bb7f8..fe7af5e10 100644 --- a/Examples/Fluid_Kerr_merged/Fluxes.hpp +++ b/Examples/Fluid_Kerr_merged/Fluxes.hpp @@ -41,10 +41,13 @@ vars_t compute_flux(const data_t P_of_rho, const vars_t &vars, data_t WW = 1. / (1. - v2); data_t hh = 1. + vars.eps + P_of_rho / vars.rho; + data_t rho_conformal = vars.rho / pow(vars.chi, 1.5); + FOR(j) { - out.Sj[j] = vars.lapse * vars.rho * hh * WW * vars.vi[idir] * vi_D[j] - - vars.shift[idir] * vars.Sj[j]; + out.Sj[j] = + vars.lapse * rho_conformal * hh * WW * vars.vi[idir] * vi_D[j] - + vars.shift[idir] * vars.Sj[j]; FOR(k) out.Sj[j] += vars.lapse * P_of_rho * h_UU[idir][k] * vars.h[j][k]; } diff --git a/Examples/Fluid_Kerr_merged/InitialFluidData.impl.hpp b/Examples/Fluid_Kerr_merged/InitialFluidData.impl.hpp index 8f1c575e4..e0290796b 100644 --- a/Examples/Fluid_Kerr_merged/InitialFluidData.impl.hpp +++ b/Examples/Fluid_Kerr_merged/InitialFluidData.impl.hpp @@ -32,7 +32,8 @@ void InitialFluidData::compute(Cell current_cell) const // calculate the field value matter_vars.rho = - m_params.rho0;// * (exp(-pow(rr / 2. / m_params.awidth, 2.0))) + m_params.delta; + m_params.rho0; // * (exp(-pow(rr / 2. / m_params.awidth, 2.0))) + + // m_params.delta; data_t v2 = 0.; FOR(i, j) v2 += metric_vars.h[i][j] * matter_vars.vi[i] * matter_vars.vi[j] / @@ -44,13 +45,15 @@ void InitialFluidData::compute(Cell current_cell) const data_t WW = 1. / (1. - v2); data_t hh = 1. + matter_vars.eps + P_of_rho / matter_vars.rho; - matter_vars.D = matter_vars.rho * sqrt(WW); - matter_vars.tau = matter_vars.rho * hh * WW - P_of_rho - matter_vars.D; + data_t rho_conformal = rho / pow(chi_regularised, 1.5); + + matter_vars.D = rho_conformal * sqrt(WW); + matter_vars.tau = rho_conformal * hh * WW - P_of_rho - matter_vars.D; FOR(i) { matter_vars.Sj[i] = 0.; FOR(j) - matter_vars.Sj[i] += matter_vars.rho * hh * WW * metric_vars.h[i][j] * + matter_vars.Sj[i] += rho_conformal * hh * WW * metric_vars.h[i][j] * matter_vars.vi[j] / chi_regularised; } @@ -58,7 +61,7 @@ void InitialFluidData::compute(Cell current_cell) const current_cell.store_vars(matter_vars); using namespace TensorAlgebra; - metric_vars.K += 24. * M_PI * matter_vars.rho; + metric_vars.K += 24. * M_PI * rho_conformal; current_cell.store_vars(metric_vars); } diff --git a/Examples/Fluid_Kerr_merged/PerfectFluid.impl.hpp b/Examples/Fluid_Kerr_merged/PerfectFluid.impl.hpp index 8ef31ec0b..519bf3f57 100644 --- a/Examples/Fluid_Kerr_merged/PerfectFluid.impl.hpp +++ b/Examples/Fluid_Kerr_merged/PerfectFluid.impl.hpp @@ -80,14 +80,21 @@ void PerfectFluid::add_matter_rhs( 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; - 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); + rhs.D = source.D /*+ vars.D * (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)*/ + ; 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)*/ + ; + + /* FOR(i) { vars_t flux = Fluxes::compute_flux(P_of_rho, vars, i); rhs.D += GR_SPACEDIM / 2. * d1.chi[i] / chi_regularised * flux.D; @@ -95,7 +102,7 @@ void PerfectFluid::add_matter_rhs( 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) { diff --git a/Examples/Fluid_Kerr_merged/PerfectFluidLevel.cpp b/Examples/Fluid_Kerr_merged/PerfectFluidLevel.cpp index 4d2f177a4..2ce8b0e83 100644 --- a/Examples/Fluid_Kerr_merged/PerfectFluidLevel.cpp +++ b/Examples/Fluid_Kerr_merged/PerfectFluidLevel.cpp @@ -39,7 +39,8 @@ void PerfectFluidLevel::specificAdvance() // Enforce trace free A_ij and positive chi and alpha BoxLoops::loop(make_compute_pack(TraceARemoval(), PositiveChiAndAlpha(), PositiveDensity()), - m_state_new, m_state_new, INCLUDE_GHOST_CELLS, disable_simd()); + m_state_new, m_state_new, INCLUDE_GHOST_CELLS, + disable_simd()); // Check for nan's if (m_p.nan_check) @@ -73,10 +74,10 @@ void PerfectFluidLevel::prePlotLevel() fillAllGhosts(); EoS eos; PerfectFluidEoS perfect_fluid(m_dx, m_p.lambda, eos); - BoxLoops::loop(MatterConstraints(perfect_fluid, m_dx, - m_p.G_Newton, c_Ham, - Interval(c_Mom1, c_Mom3)), - m_state_new, m_state_diagnostics, EXCLUDE_GHOST_CELLS, disable_simd()); + BoxLoops::loop( + MatterConstraints(perfect_fluid, m_dx, m_p.G_Newton, + c_Ham, Interval(c_Mom1, c_Mom3)), + m_state_new, m_state_diagnostics, EXCLUDE_GHOST_CELLS, disable_simd()); } #endif @@ -88,7 +89,7 @@ void PerfectFluidLevel::specificEvalRHS(GRLevelData &a_soln, GRLevelData &a_rhs, BoxLoops::loop(make_compute_pack(TraceARemoval(), PositiveDensity(), PositiveChiAndAlpha(), PrimitiveRecovery()), - a_soln, a_soln, INCLUDE_GHOST_CELLS/*, disable_simd()*/); + a_soln, a_soln, INCLUDE_GHOST_CELLS /*, disable_simd()*/); // Calculate MatterCCZ4 right hand side with matter_t = ScalarField EoS eos; @@ -100,7 +101,8 @@ void PerfectFluidLevel::specificEvalRHS(GRLevelData &a_soln, GRLevelData &a_rhs, // FourthOrderDerivatives> my_ccz4_matter(perfect_fluid, m_p.ccz4_params, m_dx, m_p.sigma, m_p.formulation, m_p.G_Newton); - BoxLoops::loop(my_ccz4_matter, a_soln, a_rhs, EXCLUDE_GHOST_CELLS, disable_simd()); + BoxLoops::loop(my_ccz4_matter, a_soln, a_rhs, EXCLUDE_GHOST_CELLS, + disable_simd()); } else if (m_p.max_spatial_derivative_order == 6) { @@ -109,7 +111,8 @@ void PerfectFluidLevel::specificEvalRHS(GRLevelData &a_soln, GRLevelData &a_rhs, // FourthOrderDerivatives> my_ccz4_matter(perfect_fluid, m_p.ccz4_params, m_dx, m_p.sigma, m_p.formulation, m_p.G_Newton); - BoxLoops::loop(my_ccz4_matter, a_soln, a_rhs, EXCLUDE_GHOST_CELLS, disable_simd()); + BoxLoops::loop(my_ccz4_matter, a_soln, a_rhs, EXCLUDE_GHOST_CELLS, + disable_simd()); } } @@ -121,7 +124,7 @@ void PerfectFluidLevel::specificUpdateODE(GRLevelData &a_soln, BoxLoops::loop(make_compute_pack(TraceARemoval(), PositiveDensity(), PositiveChiAndAlpha(), PrimitiveRecovery()), - a_soln, a_soln, INCLUDE_GHOST_CELLS/*, disable_simd()*/); + a_soln, a_soln, INCLUDE_GHOST_CELLS /*, disable_simd()*/); } void PerfectFluidLevel::preTagCells() @@ -134,5 +137,6 @@ void PerfectFluidLevel::computeTaggingCriterion( FArrayBox &tagging_criterion, const FArrayBox ¤t_state, const FArrayBox ¤t_state_diagnostics) { - BoxLoops::loop(ChiTaggingCriterion(m_dx), current_state, tagging_criterion, disable_simd()); + BoxLoops::loop(ChiTaggingCriterion(m_dx), current_state, tagging_criterion, + disable_simd()); } diff --git a/Examples/Fluid_Kerr_merged/PerfectFluidLevel.hpp b/Examples/Fluid_Kerr_merged/PerfectFluidLevel.hpp index 40cffd497..0adcb438f 100644 --- a/Examples/Fluid_Kerr_merged/PerfectFluidLevel.hpp +++ b/Examples/Fluid_Kerr_merged/PerfectFluidLevel.hpp @@ -47,7 +47,8 @@ class PerfectFluidLevel : public GRAMRLevel //! Things to do in UpdateODE step, after soln + rhs update virtual void specificUpdateODE(GRLevelData &a_soln, - const GRLevelData &a_rhs, Real a_dt) override; + const GRLevelData &a_rhs, + Real a_dt) override; /// Things to do before tagging cells (i.e. filling ghosts) virtual void preTagCells() override; diff --git a/Examples/Fluid_Kerr_merged/PositiveDensity.hpp b/Examples/Fluid_Kerr_merged/PositiveDensity.hpp index 46fc960d4..80c469878 100644 --- a/Examples/Fluid_Kerr_merged/PositiveDensity.hpp +++ b/Examples/Fluid_Kerr_merged/PositiveDensity.hpp @@ -17,24 +17,26 @@ class PositiveDensity template struct Vars { data_t chi; - data_t D; - data_t tau; - Tensor<1, data_t> Sj; - Tensor<2, data_t> h; + data_t D; + data_t tau; + Tensor<1, data_t> Sj; + Tensor<2, data_t> h; template void enum_mapping(mapping_function_t mapping_function) { - using namespace VarsTools; //define_enum_mapping is part of VarsTools - define_enum_mapping(mapping_function, c_chi, chi); - define_enum_mapping(mapping_function, c_D, D); - define_enum_mapping(mapping_function, c_tau, tau); - define_enum_mapping( - mapping_function, GRInterval(), Sj); - define_symmetric_enum_mapping( - mapping_function, GRInterval(), h); + using namespace VarsTools; // define_enum_mapping is part of + // VarsTools + define_enum_mapping(mapping_function, c_chi, chi); + define_enum_mapping(mapping_function, c_D, D); + define_enum_mapping(mapping_function, c_tau, tau); + define_enum_mapping(mapping_function, GRInterval(), + Sj); + define_symmetric_enum_mapping(mapping_function, + GRInterval(), h); } }; + private: const double m_min_D; const double m_min_v; @@ -49,19 +51,21 @@ class PositiveDensity template void compute(Cell current_cell) const { auto vars = current_cell.template load_vars(); - const data_t chi_regularised = simd_max(1e-6, vars.chi); + const data_t chi_regularised = simd_max(1e-6, vars.chi); + + // Take into account that conservative variables are conformally + // rescaled + + // D >= min_D + vars.D = simd_max(vars.D, m_min_D / pow(chi_regularised, 1.5)); - // Take into account that conservative variables are conformally rescaled - - //D >= min_D - vars.D = simd_max(vars.D, m_min_D / pow(chi_regularised, 1.5)); + data_t S2_over_chi = 0.; + FOR(i, j) S2_over_chi += vars.h[i][j] * vars.Sj[i] * vars.Sj[j]; - data_t S2_over_chi = 0.; - FOR(i,j) S2_over_chi += vars.h[i][j] * vars.Sj[i] * vars.Sj[j]; - - // S^2 <= (D + tau)^2 - data_t factor = simd_min(1., vars.chi * (vars.D + vars.tau) / sqrt(S2_over_chi)); - FOR(i) vars.Sj[i] *= factor; + // S^2 <= (D + tau)^2 + data_t factor = + simd_min(1., vars.chi * (vars.D + vars.tau) / sqrt(S2_over_chi)); + FOR(i) vars.Sj[i] *= factor; current_cell.store_vars(vars.D, c_D); current_cell.store_vars(vars.Sj, GRInterval()); diff --git a/Examples/Fluid_Kerr_merged/PrimitiveRecovery.hpp b/Examples/Fluid_Kerr_merged/PrimitiveRecovery.hpp index b521e5073..7f80f6f94 100644 --- a/Examples/Fluid_Kerr_merged/PrimitiveRecovery.hpp +++ b/Examples/Fluid_Kerr_merged/PrimitiveRecovery.hpp @@ -11,12 +11,12 @@ #include "Cell.hpp" #include "DefaultEoS.hpp" #include "EoS.hpp" -#include "simd.hpp" #include "Tensor.hpp" #include "TensorAlgebra.hpp" #include "UserVariables.hpp" #include "UsingNamespace.H" #include "VarsTools.hpp" +#include "simd.hpp" // template class PrimitiveRecovery : public EoS @@ -51,7 +51,7 @@ class PrimitiveRecovery : public EoS data_t xa = 1.5 * (1. + q); Wa = sqrt(pow(xa, 2.) / (pow(xa, 2.) - r)); - vars.rho = vars.D / Wa; + vars.rho = pow(vars.chi, 1.5) * vars.D / Wa; vars.eps = -1. + xa / Wa * (1. - Wa * Wa) + Wa * (1. + q); EoS::compute_eos(P_of_rho, vars); xn = Wa * (1. + vars.eps + P_of_rho / vars.rho); @@ -59,16 +59,16 @@ class PrimitiveRecovery : public EoS int i = 0; - data_t empty; - while (!simd_all_false(simd_compare_lt(tolerance, diff), empty)) + data_t empty; + while (!simd_all_false(simd_compare_lt(tolerance, diff), empty)) { i++; Wa = sqrt(pow(xa, 2.) / (pow(xa, 2.) - r)); - vars.rho = vars.D / Wa; + vars.rho = pow(vars.chi, 1.5) * vars.D / Wa; vars.eps = -1. + xa / Wa * (1. - Wa * Wa) + Wa * (1. + q); - EoS::compute_eos(P_of_rho, vars); - xn = Wa * (1. + vars.eps + P_of_rho / vars.rho); + EoS::compute_eos(P_of_rho, vars); + xn = Wa * (1. + vars.eps + P_of_rho / vars.rho); diff = abs(xn - xa); xa = xn; if (i >= 100) diff --git a/Examples/Fluid_Kerr_merged/SimulationParameters.hpp b/Examples/Fluid_Kerr_merged/SimulationParameters.hpp index 9c591037e..0383a5415 100644 --- a/Examples/Fluid_Kerr_merged/SimulationParameters.hpp +++ b/Examples/Fluid_Kerr_merged/SimulationParameters.hpp @@ -37,21 +37,6 @@ class SimulationParameters : public SimulationParametersBase pp.load("fluid_delta", initial_params.delta, 0.2); pp.load("lambda", lambda, 1.); // eigenvalue for numerical flux - // Reading data - // pp.load("spacing", initial_params.spacing); - //pp.load("lines", lines); - //double rho_1D[lines]; - - //double tmp_data; - //ifstream read_file("rho_vs_r.csv"); - - //for (int i = 0; i < lines; ++i) - //{ - // read_file >> tmp_data; - // rho_1D[i] = tmp_data; - //} - //initial_params.rho_1D = rho_1D; - // Initial Kerr data pp.load("kerr_mass", kerr_params.mass); pp.load("kerr_spin", kerr_params.spin); @@ -80,7 +65,6 @@ class SimulationParameters : public SimulationParametersBase // Initial data for matter and potential and BH double G_Newton; double lambda; - //int lines; InitialFluidData::params_t initial_params; KerrBH::params_t kerr_params; }; diff --git a/Examples/Fluid_Kerr_merged/Sources.hpp b/Examples/Fluid_Kerr_merged/Sources.hpp index 55c3e8a43..26250e586 100644 --- a/Examples/Fluid_Kerr_merged/Sources.hpp +++ b/Examples/Fluid_Kerr_merged/Sources.hpp @@ -28,6 +28,8 @@ vars_t compute_source(const data_t P_of_rho, const vars_t &vars, data_t WW = 1. / (1. - v2); data_t hh = 1. + vars.eps + P_of_rho / vars.rho; + data_t rho_conformal = vars.rho / pow(chi_regularised, 1.5); + out.D = 0.; FOR(j) { @@ -40,8 +42,8 @@ vars_t compute_source(const data_t P_of_rho, const vars_t &vars, out.Sj[j] += vars.lapse / 2. * (d1.h[i][k][j] - vars.h[i][k] * d1.chi[j] / chi_regularised) * - (vars.rho * hh * WW * vars.vi[i] * vars.vi[k] / - chi_regularised + + (rho_conformal * hh * WW * vars.vi[i] * + vars.vi[k] / chi_regularised + P_of_rho * h_UU[i][k]); } } @@ -50,7 +52,7 @@ vars_t compute_source(const data_t P_of_rho, const vars_t &vars, FOR(i, j) { out.tau += vars.lapse * (vars.A[i][j] + vars.h[i][j] / 3. * vars.K) * - (vars.rho * hh * WW * vars.vi[i] * vars.vi[j] / + (rho_conformal * hh * WW * vars.vi[i] * vars.vi[j] / chi_regularised + P_of_rho * h_UU[i][j]) - vars.chi * h_UU[i][j] * vars.Sj[i] * d1.lapse[j];