Skip to content

Commit

Permalink
Pre-pull request update
Browse files Browse the repository at this point in the history
  • Loading branch information
melley95 committed Oct 2, 2024
1 parent 6b7a777 commit de568b3
Show file tree
Hide file tree
Showing 5 changed files with 69 additions and 87 deletions.
18 changes: 7 additions & 11 deletions Examples/FRWScalarField/FRWScalarLevel.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -16,8 +16,8 @@
#include "FixedGridsTaggingCriterion.hpp"

// Problem specific includes
#include "InitialScalarData.hpp"
#include "FRW.hpp"
#include "InitialScalarData.hpp"
#include "MatterEvolution.hpp"
#include "ScalarField.hpp"
#include "ScalarPotential.hpp"
Expand All @@ -32,18 +32,17 @@ void FRWScalarLevel::initialData()
// First set everything to zero ... we don't want undefined values in
// constraints etc, then initial conditions for fields
SetValue set_zero(0.0);
FRW frw_bg(m_p.bg_params, m_dx, m_time);
FRW frw_bg(m_p.bg_params, m_dx, m_time);
InitialScalarData initial_sf(m_p.initial_params, m_dx);
auto compute_pack = make_compute_pack(set_zero, frw_bg);
BoxLoops::loop(compute_pack, m_state_diagnostics, m_state_diagnostics,
SKIP_GHOST_CELLS);
BoxLoops::loop(initial_sf, m_state_new, m_state_new, FILL_GHOST_CELLS);

}

// Things to do in RHS update, at each RK4 step
void FRWScalarLevel::specificEvalRHS(GRLevelData &a_soln, GRLevelData &a_rhs,
const double a_time)
const double a_time)
{
// Calculate MatterCCZ4 right hand side with matter_t = ScalarField
// We don't want undefined values floating around in the constraints so
Expand All @@ -54,7 +53,6 @@ void FRWScalarLevel::specificEvalRHS(GRLevelData &a_soln, GRLevelData &a_rhs,
MatterEvolution<ScalarFieldWithPotential, FRW> my_matter(
scalar_field, frw_bg, m_p.sigma, m_dx, m_p.center);
BoxLoops::loop(my_matter, a_soln, a_rhs, SKIP_GHOST_CELLS);

}

void FRWScalarLevel::specificPostTimeStep()
Expand All @@ -65,25 +63,23 @@ void FRWScalarLevel::specificPostTimeStep()
disable_simd());

// At any level, but after the min_level timestep
int min_level = 0;
int min_level = 0;
bool calculate_diagnostics = at_level_timestep_multiple(min_level);
if (calculate_diagnostics)
{
fillAllGhosts();
ScalarPotential potential(m_p.initial_params);
ScalarFieldWithPotential scalar_field(potential);
FRW frw_bg(m_p.bg_params, m_dx, m_time);
BoxLoops::loop(frw_bg, m_state_new,
m_state_diagnostics, SKIP_GHOST_CELLS);

BoxLoops::loop(frw_bg, m_state_new, m_state_diagnostics,
SKIP_GHOST_CELLS);
}

// write out the integral after each timestep on minimum level

}

void FRWScalarLevel::computeTaggingCriterion(FArrayBox &tagging_criterion,
const FArrayBox &current_state)
const FArrayBox &current_state)
{
BoxLoops::loop(FixedGridsTaggingCriterion(m_dx, m_level, m_p.L, m_p.center),
current_state, tagging_criterion);
Expand Down
2 changes: 1 addition & 1 deletion Examples/FRWScalarField/Main_FRWScalar.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,7 @@ int runGRChombo(int argc, char *argv[])
// and an associated LevelFactory)
GRAMR gr_amr;
DefaultLevelFactory<FRWScalarLevel> scalar_field_level_fact(gr_amr,
sim_params);
sim_params);
setupAMRObject(gr_amr, scalar_field_level_fact);

// call this after amr object setup so grids known
Expand Down
5 changes: 2 additions & 3 deletions Examples/FRWScalarField/SimulationParameters.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -11,8 +11,8 @@
#include "GRParmParse.hpp"

// Problem specific includes:
#include "InitialScalarData.hpp"
#include "FRW.hpp"
#include "InitialScalarData.hpp"

class SimulationParameters : public FixedBGSimulationParametersBase
{
Expand All @@ -35,8 +35,8 @@ class SimulationParameters : public FixedBGSimulationParametersBase
// Initial and Kerr data
pp.load("frw_rho0", bg_params.rho0, 0.0);
pp.load("frw_omega", bg_params.omega, 0.0);
pp.load("frw_alpha", bg_params.alpha, 0.0);
pp.load("center", bg_params.center, center);

}

void check_params()
Expand All @@ -45,7 +45,6 @@ class SimulationParameters : public FixedBGSimulationParametersBase
initial_params.mass < 0.2 / coarsest_dx / dt_multiplier,
"oscillations of scalar field do not appear to be "
"resolved on coarsest level");

}

// Problem specific parameters
Expand Down
3 changes: 2 additions & 1 deletion Examples/FRWScalarField/params.txt
Original file line number Diff line number Diff line change
Expand Up @@ -42,7 +42,8 @@ scalar_amplitude = 0.1


frw_rho0 = 1.0 # Initial energy density
frw_omega = 0.0 # Equation of state parameter
frw_omega = 0.0 # Background equation of state
frw_alpha = 0.0 # = 0 (cosmic time), 1 (conformal time)



Expand Down
128 changes: 57 additions & 71 deletions Source/Background/FRW.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -17,123 +17,109 @@

//! Class which computes an FRW background


class FRW
{
public:
//! Struct for the params of the BH
struct params_t
{
double rho0; // initial energy density
double rho0; // initial energy density
double omega; // equation of state parameter = 0 (MD), 1/3 (RD)
std::array<double, CH_SPACEDIM> center;

double alpha; // can set to conformal time
std::array<double, CH_SPACEDIM> center;
};

template <class data_t> using Vars = ADMFixedBGVars::Vars<data_t>;


const params_t m_params;

const double m_time;
const double m_dx;


FRW(params_t a_params, double a_dx, double a_time) : m_params(a_params), m_dx(a_dx), m_time(a_time){}

FRW(params_t a_params, double a_dx, double a_time)
: m_params(a_params), m_dx(a_dx), m_time(a_time)
{
}

template <class data_t> void compute(Cell<data_t> current_cell) const
{



const Coordinates<data_t> coords(current_cell, m_dx, m_params.center);
Vars<data_t> metric_vars;


compute_metric_background(metric_vars, coords);




//Reconstruct FRW quantities to check as diagnostics
data_t scalefac = sqrt(metric_vars.gamma[0][0]);
data_t hubbleparam = -metric_vars.K_tensor[0][0];
// Reconstruct FRW quantities to check as diagnostics
data_t scalefac = sqrt(metric_vars.gamma[0][0]);
data_t hubbleparam = -metric_vars.K_tensor[0][0];

current_cell.store_vars(scalefac, c_scalefac);
current_cell.store_vars(hubbleparam, c_hubbleparam);



current_cell.store_vars(scalefac, c_scalefac);
current_cell.store_vars(hubbleparam, c_hubbleparam);
}


//Analytical background
template <class data_t, template <typename> class vars_t>
// Analytical background
template <class data_t, template <typename> class vars_t>
void compute_metric_background(vars_t<data_t> &vars,
const Coordinates<data_t> &coords) const
{

data_t scalefac;
data_t hubparam;

// Zero components - to do this in level?
vars.K = 0.0;

FOR(i) vars.shift[i] = 0.0;
FOR(i,j) vars.gamma[i][j] = 0.0;
FOR(i,j) vars.K_tensor[i][j] = 0.0;


FOR(i,j, k) vars.d1_gamma[i][j][k] = 0.0;
FOR(i) vars.d1_lapse[i] = 0.0;
FOR(i,j) vars.d1_shift[i][j] = 0.0;
data_t scalefac;
data_t hubparam;

// ds^2 = -dt^2 + a^2 dx_i dx^i assumin flat - TODO: curvature, conformal time
// ds^2 = -dt^2 + a^2 dx_i dx^i

vars.lapse = 1.0; // Conformal time - a^2

if (m_params.rho0 == 0.0){
if (m_params.rho0 == 0.0)
{

scalefac = 1.0;

FOR(i) vars.gamma[i][i] = scalefac; // Minkowski

vars.lapse = 1.0;
}


else{

data_t K0 = -sqrt(24.0*M_PI*m_params.rho0); //Ham constraint

data_t t0 = -1.0/((1.0+m_params.omega)*K0); //Defining t0 from time evolution of scale factor



scalefac = pow((m_time/t0)+1, 2/(3*(1+m_params.omega))); //Scale factor evolution

hubparam = 2.0/(3.0*(1+m_params.omega)*(m_time+t0)); //Hubble parameter evolution




vars.K = -3.0*hubparam; //K = -3H (not used but included for completeness)

FOR(i) vars.gamma[i][i] = pow(scalefac, 2.0); //gamma_ij = a^2

FOR(i) vars.K_tensor[i][i] = -hubparam; //Kij = -H gamma_ij

}


}
else
{

data_t K0 = -sqrt(24.0 * M_PI * m_params.rho0); // Ham constraint

data_t t0 =
-1.0 / ((1.0 + m_params.omega) *
K0); // Defining t0 from time evolution of scale factor

if (m_params.omega == -1)
{ // Cosmo constant

hubparam = -K0 / 3.0;



scalefac = exp(hubparam * m_time);
}

else
{

scalefac = pow(
(m_time / t0) + 1,
2 / (3 * (1 + m_params.omega))); // Scale factor evolution

hubparam = 2.0 / (3.0 * (1 + m_params.omega) *
(m_time + t0)); // Hubble parameter evolution
}

vars.lapse = pow(
scalefac,
m_params
.alpha); // Cosmic time alpha = 0, conformal time alpha = 1

vars.K =
-3.0 *
hubparam; // K = -3H (not used but included for completeness)

FOR(i) vars.gamma[i][i] = pow(scalefac, 2.0); // gamma_ij = a^2

FOR(i) vars.K_tensor[i][i] = -hubparam; // Kij = -H gamma_ij
}
}
};

#endif /* FRW_HPP_ */

0 comments on commit de568b3

Please sign in to comment.