Skip to content

Commit

Permalink
Change to conformal definition of conservative variables
Browse files Browse the repository at this point in the history
  • Loading branch information
@dinatraykova
dinatraykova committed Aug 1, 2024
1 parent a4295cb commit 4e002a6
Show file tree
Hide file tree
Showing 11 changed files with 89 additions and 80 deletions.
7 changes: 4 additions & 3 deletions Examples/Fluid_Kerr_merged/ConservedQuantities.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,11 +31,12 @@ void PtoC(const data_t P_of_rho, vars_t<data_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
Expand Down
4 changes: 2 additions & 2 deletions Examples/Fluid_Kerr_merged/EoS.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,9 +19,9 @@ class EoS
void compute_eos(data_t &P_of_rho, const vars_t<data_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);
}
};
Expand Down
7 changes: 5 additions & 2 deletions Examples/Fluid_Kerr_merged/Fluxes.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,10 +41,13 @@ vars_t<data_t> compute_flux(const data_t P_of_rho, const vars_t<data_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];
}
Expand Down
13 changes: 8 additions & 5 deletions Examples/Fluid_Kerr_merged/InitialFluidData.impl.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,8 @@ void InitialFluidData::compute(Cell<data_t> 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] /
Expand All @@ -44,21 +45,23 @@ void InitialFluidData::compute(Cell<data_t> 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;
}

// store the vars
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);
}

Expand Down
23 changes: 15 additions & 8 deletions Examples/Fluid_Kerr_merged/PerfectFluid.impl.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -80,22 +80,29 @@ void PerfectFluid<eos_t>::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<data_t> flux = Fluxes::compute_flux(P_of_rho, 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)
{
Expand Down
24 changes: 14 additions & 10 deletions Examples/Fluid_Kerr_merged/PerfectFluidLevel.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -73,10 +74,10 @@ void PerfectFluidLevel::prePlotLevel()
fillAllGhosts();
EoS eos;
PerfectFluidEoS perfect_fluid(m_dx, m_p.lambda, eos);
BoxLoops::loop(MatterConstraints<PerfectFluidEoS>(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<PerfectFluidEoS>(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

Expand All @@ -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;
Expand All @@ -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)
{
Expand All @@ -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());
}
}

Expand All @@ -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()
Expand All @@ -134,5 +137,6 @@ void PerfectFluidLevel::computeTaggingCriterion(
FArrayBox &tagging_criterion, const FArrayBox &current_state,
const FArrayBox &current_state_diagnostics)
{
BoxLoops::loop(ChiTaggingCriterion(m_dx), current_state, tagging_criterion, disable_simd());
BoxLoops::loop(ChiTaggingCriterion(m_dx), current_state, tagging_criterion,
disable_simd());
}
3 changes: 2 additions & 1 deletion Examples/Fluid_Kerr_merged/PerfectFluidLevel.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down
50 changes: 27 additions & 23 deletions Examples/Fluid_Kerr_merged/PositiveDensity.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,24 +17,26 @@ class PositiveDensity
template <class data_t> 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 <typename mapping_function_t>
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<c_Sj1, c_Sj3>(), Sj);
define_symmetric_enum_mapping(
mapping_function, GRInterval<c_h11, c_h33>(), 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<c_Sj1, c_Sj3>(),
Sj);
define_symmetric_enum_mapping(mapping_function,
GRInterval<c_h11, c_h33>(), h);
}
};

private:
const double m_min_D;
const double m_min_v;
Expand All @@ -49,19 +51,21 @@ class PositiveDensity
template <class data_t> void compute(Cell<data_t> current_cell) const
{
auto vars = current_cell.template load_vars<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<c_Sj1, c_Sj3>());
Expand Down
14 changes: 7 additions & 7 deletions Examples/Fluid_Kerr_merged/PrimitiveRecovery.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 eos_t = DefaultEoS>
class PrimitiveRecovery : public EoS
Expand Down Expand Up @@ -51,24 +51,24 @@ 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);
diff = abs(xn - xa);

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)
Expand Down
16 changes: 0 additions & 16 deletions Examples/Fluid_Kerr_merged/SimulationParameters.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down Expand Up @@ -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;
};
Expand Down
8 changes: 5 additions & 3 deletions Examples/Fluid_Kerr_merged/Sources.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,8 @@ vars_t<data_t> compute_source(const data_t P_of_rho, const vars_t<data_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)
{
Expand All @@ -40,8 +42,8 @@ vars_t<data_t> compute_source(const data_t P_of_rho, const vars_t<data_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]);
}
}
Expand All @@ -50,7 +52,7 @@ vars_t<data_t> compute_source(const data_t P_of_rho, const vars_t<data_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];
Expand Down

0 comments on commit 4e002a6

Please sign in to comment.