Skip to content

Commit

Permalink
clang formatting
Browse files Browse the repository at this point in the history
  • Loading branch information
@ShaunFell
ShaunFell committed Apr 17, 2024
1 parent 63f5165 commit 0d64bd1
Show file tree
Hide file tree
Showing 15 changed files with 423 additions and 349 deletions.
17 changes: 8 additions & 9 deletions Examples/ProcaKdS/EnhancedExcisionEvolution.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,21 +29,22 @@ class DumbExcisionEvolution

public:
DumbExcisionEvolution(const double a_dx,
const std::array<double, CH_SPACEDIM> a_center,
double inner_radius, double outer_radius, double a_horizon)
: m_dx(a_dx), m_center(a_center), m_inner_radius(inner_radius), m_outer_radius(outer_radius), m_outer_horizon(a_horizon)
const std::array<double, CH_SPACEDIM> a_center,
double inner_radius, double outer_radius,
double a_horizon)
: m_dx(a_dx), m_center(a_center), m_inner_radius(inner_radius),
m_outer_radius(outer_radius), m_outer_horizon(a_horizon)
{
}

void compute(const Cell<double> current_cell) const
{
// where are we?!
const Coordinates<double> coords(current_cell, m_dx, m_center);
const double radius { coords.get_radius() };
const double radius{coords.get_radius()};


bool is_excised { ( radius < m_inner_radius ) || ( radius > m_outer_radius ) };
bool is_Z_excised { radius < m_outer_horizon };
bool is_excised{(radius < m_inner_radius) || (radius > m_outer_radius)};
bool is_Z_excised{radius < m_outer_horizon};

if (is_excised)
{
Expand All @@ -63,9 +64,7 @@ class DumbExcisionEvolution
VarsTools::assign(vars, 0.0);
// assign values of rhs or vars in output box
current_cell.store_vars(vars);

}

}
};

Expand Down
16 changes: 9 additions & 7 deletions Examples/ProcaKdS/ProcaFieldLevel.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,8 @@ class ProcaFieldLevel : public BaseProcaFieldLevel<KerrdeSitter, ProcaField>

// Excise within horizon
DumbExcisionEvolution excisor(
m_dx, m_p.center, m_p.evolution_excision_inner, m_p.evolution_excision_outer, m_p.background_params.r_plus);
m_dx, m_p.center, m_p.evolution_excision_inner,
m_p.evolution_excision_outer, m_p.background_params.r_plus);

// Loop over box cells and excise cells within horizon
BoxLoops::loop(excisor, m_state_new, m_state_new, SKIP_GHOST_CELLS,
Expand All @@ -48,8 +49,9 @@ class ProcaFieldLevel : public BaseProcaFieldLevel<KerrdeSitter, ProcaField>
// Done!
}

//override default method for new excision routine
virtual void specificEvalRHS(GRLevelData &a_soln, GRLevelData &a_rhs, const double a_time) override
// override default method for new excision routine
virtual void specificEvalRHS(GRLevelData &a_soln, GRLevelData &a_rhs,
const double a_time) override
{
// Calculate right hand side
KerrdeSitter background_init{m_p.background_params, m_dx};
Expand All @@ -58,14 +60,14 @@ class ProcaFieldLevel : public BaseProcaFieldLevel<KerrdeSitter, ProcaField>
MatterEvolution<ProcaField, KerrdeSitter> matter_class(
proca_field, background_init, m_p.sigma, m_dx, m_p.center);
DumbExcisionEvolution excisor(
m_dx, m_p.center, m_p.evolution_excision_inner, m_p.evolution_excision_outer, m_p.background_params.r_plus);
m_dx, m_p.center, m_p.evolution_excision_inner,
m_p.evolution_excision_outer, m_p.background_params.r_plus);

BoxLoops::loop(matter_class, a_soln, a_rhs, SKIP_GHOST_CELLS);
BoxLoops::loop(excisor, a_soln, a_rhs, SKIP_GHOST_CELLS, disable_simd());
BoxLoops::loop(excisor, a_soln, a_rhs, SKIP_GHOST_CELLS,
disable_simd());
}



/* //override method to add black hole evolution
virtual void additionalPostTimeStep() override
{
Expand Down
37 changes: 24 additions & 13 deletions Examples/ProcaKdS/SimulationParameters.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,9 +7,9 @@
#include "ProcaField.hpp"
#include "ProcaSimulationParameters.hpp"

//For interpolating the BH horizon
#include "SimpleDataReader.hpp"
// For interpolating the BH horizon
#include "DataManipulation.hpp"
#include "SimpleDataReader.hpp"

class SimulationParameters : public ProcaSimulationParameters
{
Expand Down Expand Up @@ -51,7 +51,7 @@ class SimulationParameters : public ProcaSimulationParameters
// turn on black hole evolution
pp.load("bh_evolution", evolve_bh, false);

//excisoin
// excisoin
pp.load("evolution_excision_inner", evolution_excision_inner, 0.0);
pp.load("evolution_excision_outer", evolution_excision_outer, 999.0);
}
Expand All @@ -71,25 +71,36 @@ class SimulationParameters : public ProcaSimulationParameters

void compute_bh_horizon()
{
SimpleDataReader<double> reader {"KerrdeSitter_rPlus.dat"};
SimpleDataReader<double> reader{"KerrdeSitter_rPlus.dat"};
DataContainer<double> data = reader.get_data();

pout() << "Interpolating BH horizon" << std::endl;

std::vector<std::vector<double>> coords {data.get_coords()};
std::vector<double> horizon_values {data.get_data()};
std::vector<std::vector<double>> coords{data.get_coords()};
std::vector<double> horizon_values{data.get_data()};

std::vector<double> query_point { background_params.cosmo_constant, background_params.spin };
std::vector<double> query_point{background_params.cosmo_constant,
background_params.spin};

std::pair< std::vector<double>, std::vector<std::vector<double>> > nearest_3_neighbors { DataManipulation::find_nearest_neighbors(coords, query_point, 3) };
std::pair<std::vector<double>, std::vector<std::vector<double>>>
nearest_3_neighbors{DataManipulation::find_nearest_neighbors(
coords, query_point, 3)};

std::vector<double> point1 { nearest_3_neighbors.second[0][0], nearest_3_neighbors.second[0][1] , horizon_values[nearest_3_neighbors.first[0]] };
std::vector<double> point2 { nearest_3_neighbors.second[1][0], nearest_3_neighbors.second[1][1] , horizon_values[nearest_3_neighbors.first[1]] };
std::vector<double> point3 { nearest_3_neighbors.second[2][0], nearest_3_neighbors.second[2][1] , horizon_values[nearest_3_neighbors.first[2]] };
std::vector<double> point1{
nearest_3_neighbors.second[0][0], nearest_3_neighbors.second[0][1],
horizon_values[nearest_3_neighbors.first[0]]};
std::vector<double> point2{
nearest_3_neighbors.second[1][0], nearest_3_neighbors.second[1][1],
horizon_values[nearest_3_neighbors.first[1]]};
std::vector<double> point3{
nearest_3_neighbors.second[2][0], nearest_3_neighbors.second[2][1],
horizon_values[nearest_3_neighbors.first[2]]};

double interpolated_horizon { DataManipulation::lin_interp_2d(point1, point2, point3, query_point) };
double interpolated_horizon{DataManipulation::lin_interp_2d(
point1, point2, point3, query_point)};

pout() << "Interpolated horizon value: " << interpolated_horizon << std::endl;
pout() << "Interpolated horizon value: " << interpolated_horizon
<< std::endl;
background_params.r_plus = interpolated_horizon;
}

Expand Down
132 changes: 94 additions & 38 deletions Source/Background/KerrdeSitter.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -96,28 +96,38 @@ class KerrdeSitter
data_t Lambda_r{1 - CC / 3 * r2};
data_t Gamma{2 * M * r * delta_theta + rho2 * Theta * Lambda_r};
data_t Upsilon{delta_theta * (r2 + spin2) - delta_r};
data_t spin2r2 { spin2 + r2};
data_t spin2r2{spin2 + r2};

// Compute the 3+1 metric variables

// lapse
data_t lapse_2{ delta_theta * Lambda_r * Lambda_r * rho2 / Gamma };
// lapse
data_t lapse_2{delta_theta * Lambda_r * Lambda_r * rho2 / Gamma};
metric_vars.lapse = sqrt(lapse_2);

// shift
Tensor<1, data_t> spherical_shift;
data_t denominator { 2 * spin2 * M * r * sin_theta2 * Lambda_r + Theta * ( 2 * M * r + (spin2 + r2) * Lambda_r ) * rho2 };
spherical_shift[0] = 2 * M * r * ( spin2 + r2 ) * delta_theta * Lambda_r / denominator;
data_t denominator{2 * spin2 * M * r * sin_theta2 * Lambda_r +
Theta * (2 * M * r + (spin2 + r2) * Lambda_r) *
rho2};
spherical_shift[0] =
2 * M * r * (spin2 + r2) * delta_theta * Lambda_r / denominator;
spherical_shift[1] = 0.;
spherical_shift[2] = - 2 * spin * M * r * delta_theta * Lambda_r / denominator;
spherical_shift[2] =
-2 * spin * M * r * delta_theta * Lambda_r / denominator;

// spatial metric
Tensor<2, data_t> spherical_gamma;
spherical_gamma[0][0] = rho2 * ( 2 * M * r + spin2r2 * Lambda_r) / (Lambda_r * Lambda_r * spin2r2 * spin2r2);
spherical_gamma[0][0] = rho2 * (2 * M * r + spin2r2 * Lambda_r) /
(Lambda_r * Lambda_r * spin2r2 * spin2r2);
spherical_gamma[1][1] = rho2 / delta_theta;
spherical_gamma[2][2] = sin_theta2 * ( spin2r2 * Theta + 2 * spin2 * M * r * sin_theta2 / rho2 ) / (Theta * Theta) ;
spherical_gamma[2][2] =
sin_theta2 *
(spin2r2 * Theta + 2 * spin2 * M * r * sin_theta2 / rho2) /
(Theta * Theta);
spherical_gamma[0][1] = 0.;
spherical_gamma[0][2] = - 2 * spin * M * r * sin_theta2 / ( spin2 * Theta * Lambda_r + r2 * Theta * Lambda_r );
spherical_gamma[0][2] =
-2 * spin * M * r * sin_theta2 /
(spin2 * Theta * Lambda_r + r2 * Theta * Lambda_r);
spherical_gamma[1][0] = 0.;
spherical_gamma[1][2] = 0.;
spherical_gamma[2][0] = spherical_gamma[0][2];
Expand All @@ -126,59 +136,109 @@ class KerrdeSitter
// Now the 1st derivatives

// first derivative of the metric variables
data_t delta_r_dr{-2. / 3. * (3. * M + r * (- 3. + spin2 * CC + 2. * r2 * CC))};
data_t delta_theta_dtheta{- 2. / 3. * spin2 * CC * cos_theta * sin_theta};
data_t delta_r_dr{-2. / 3. *
(3. * M + r * (-3. + spin2 * CC + 2. * r2 * CC))};
data_t delta_theta_dtheta{-2. / 3. * spin2 * CC * cos_theta *
sin_theta};
data_t rho2_dr{2. * r};
data_t rho2_dtheta{-2. * spin2 * cos_theta * sin_theta};
data_t Lambda_r_dr{-2. / 3. * r * CC};
data_t Gamma_dr { 2. * M * delta_theta + rho2_dr * Theta * Lambda_r + rho2 * Theta * Lambda_r_dr };
data_t Gamma_dtheta { 2. * M * r * delta_theta_dtheta + rho2_dtheta * Theta * Lambda_r };
data_t Gamma_dr{2. * M * delta_theta + rho2_dr * Theta * Lambda_r +
rho2 * Theta * Lambda_r_dr};
data_t Gamma_dtheta{2. * M * r * delta_theta_dtheta +
rho2_dtheta * Theta * Lambda_r};

// First derivatives of lapse
Tensor<1, data_t> dlapse_dr;
// derivative of lapse w.r.t r
dlapse_dr[0] = delta_theta * Lambda_r / (2. * Gamma * Gamma * metric_vars.lapse) * ( - Lambda_r * rho2 * Gamma_dr + Gamma * ( 2. * rho2 * Lambda_r_dr + Lambda_r * rho2_dr));
dlapse_dr[0] = delta_theta * Lambda_r /
(2. * Gamma * Gamma * metric_vars.lapse) *
(-Lambda_r * rho2 * Gamma_dr +
Gamma * (2. * rho2 * Lambda_r_dr + Lambda_r * rho2_dr));
// derivative of lapse w.r.t theta
dlapse_dr[1] = Lambda_r * Lambda_r / (2. * Gamma * Gamma * metric_vars.lapse) * ( - delta_theta * rho2 * Gamma_dtheta + Gamma * ( rho2 * delta_theta_dtheta + delta_theta * rho2_dtheta));
dlapse_dr[1] =
Lambda_r * Lambda_r / (2. * Gamma * Gamma * metric_vars.lapse) *
(-delta_theta * rho2 * Gamma_dtheta +
Gamma * (rho2 * delta_theta_dtheta + delta_theta * rho2_dtheta));
// derivative of lapse w.r.t. phi
dlapse_dr[2] = 0.;


// First derivatives of the shift
Tensor<1, Tensor<1, data_t>> dshift_dr;
// Zero initialize
FOR2(i, j) { dshift_dr[i][j] = 0; };
// Set the only 2 non-zero components
// derivative w.r.t. r
dshift_dr[0][0] = 2 * M * delta_theta * ( 2. * M * r2 * spin2r2 * Theta * rho2 * Lambda_r_dr +
2. * M * r2 * Theta * Lambda_r * ( 2. * r * rho2 - spin2r2 * rho2_dr) +
Lambda_r * Lambda_r * (4. * spin2 * M * r2 * r * sin_theta2 + spin2r2 * spin2r2 * Theta * (rho2 - r * rho2_dr)));
dshift_dr[0][0] =
2 * M * delta_theta *
(2. * M * r2 * spin2r2 * Theta * rho2 * Lambda_r_dr +
2. * M * r2 * Theta * Lambda_r *
(2. * r * rho2 - spin2r2 * rho2_dr) +
Lambda_r * Lambda_r *
(4. * spin2 * M * r2 * r * sin_theta2 +
spin2r2 * spin2r2 * Theta * (rho2 - r * rho2_dr)));
dshift_dr[0][0] *= 1. / denominator / denominator;
dshift_dr[0][1] = 2. * M * r * spin2r2 * Lambda_r * ( delta_theta_dtheta * ( Theta * rho2 * (2. * M * r + spin2r2 * Lambda_r) + 2 * spin2 * M * r * Lambda_r * sin_theta2) + delta_theta * ( - 2. * spin2 * M * r * Lambda_r * 2. * sin_theta * cos_theta - Theta * rho2_dtheta * (2 * M * r + Lambda_r * spin2r2 ) ) );
dshift_dr[0][1] =
2. * M * r * spin2r2 * Lambda_r *
(delta_theta_dtheta *
(Theta * rho2 * (2. * M * r + spin2r2 * Lambda_r) +
2 * spin2 * M * r * Lambda_r * sin_theta2) +
delta_theta *
(-2. * spin2 * M * r * Lambda_r * 2. * sin_theta * cos_theta -
Theta * rho2_dtheta * (2 * M * r + Lambda_r * spin2r2)));
dshift_dr[0][1] *= 1. / denominator / denominator;
dshift_dr[2][0] = 2. * spin * M * delta_theta * Theta * ( rho2 * ( Lambda_r * Lambda_r * (r2 - spin2) - 2. * M * r2 * Lambda_r_dr ) + r * Lambda_r * rho2_dr * (2. * M * r + Lambda_r * spin2r2 ) );
dshift_dr[2][0] =
2. * spin * M * delta_theta * Theta *
(rho2 * (Lambda_r * Lambda_r * (r2 - spin2) -
2. * M * r2 * Lambda_r_dr) +
r * Lambda_r * rho2_dr * (2. * M * r + Lambda_r * spin2r2));
dshift_dr[2][0] *= 1. / denominator / denominator;
dshift_dr[2][1] = 2. * spin * M * r * Lambda_r * ( - delta_theta_dtheta * ( Theta * rho2 * ( 2. * M * r + Lambda_r * spin2r2) + 2. * spin2 * M * r * Lambda_r * sin_theta2 ) + delta_theta * ( 4. * spin2 * M * r * Lambda_r * cos_theta * sin_theta + Theta * rho2_dtheta * ( 2 * M * r + Lambda_r * spin2r2 ) ) );
dshift_dr[2][1] =
2. * spin * M * r * Lambda_r *
(-delta_theta_dtheta *
(Theta * rho2 * (2. * M * r + Lambda_r * spin2r2) +
2. * spin2 * M * r * Lambda_r * sin_theta2) +
delta_theta *
(4. * spin2 * M * r * Lambda_r * cos_theta * sin_theta +
Theta * rho2_dtheta * (2 * M * r + Lambda_r * spin2r2)));
dshift_dr[2][1] *= 1. / denominator / denominator;


// First derivatives of the spatial metric
Tensor<2, Tensor<1, data_t>> dgamma_dr;
// Zero initialize
FOR3(i, j, k) { dgamma_dr[i][j][k] = 0; }
// derivative w.r.t. r
dgamma_dr[0][0][0] = rho2 * ( 2 * Lambda_r * ( M * (spin2 - 3 * r2) - r * Lambda_r * spin2r2 ) - spin2r2 * Lambda_r_dr * ( 4 * M * r + spin2r2 * Lambda_r ) ) + spin2r2 * Lambda_r * rho2_dr * (2 * M * r + spin2r2 * Lambda_r);
dgamma_dr[0][0][0] *= 1. / (Lambda_r * Lambda_r * Lambda_r * spin2r2 * spin2r2 * spin2r2);
dgamma_dr[0][2][0] = 2 * spin * M * sin_theta2 * ( Lambda_r * (r2 - spin2) + r * Lambda_r_dr * spin2r2) / ( Lambda_r * Lambda_r * Theta * spin2r2 * spin2r2 );
dgamma_dr[0][0][0] =
rho2 * (2 * Lambda_r *
(M * (spin2 - 3 * r2) - r * Lambda_r * spin2r2) -
spin2r2 * Lambda_r_dr * (4 * M * r + spin2r2 * Lambda_r)) +
spin2r2 * Lambda_r * rho2_dr * (2 * M * r + spin2r2 * Lambda_r);
dgamma_dr[0][0][0] *=
1. / (Lambda_r * Lambda_r * Lambda_r * spin2r2 * spin2r2 * spin2r2);
dgamma_dr[0][2][0] =
2 * spin * M * sin_theta2 *
(Lambda_r * (r2 - spin2) + r * Lambda_r_dr * spin2r2) /
(Lambda_r * Lambda_r * Theta * spin2r2 * spin2r2);
dgamma_dr[1][1][0] = rho2_dr / delta_theta;
dgamma_dr[2][0][0] = dgamma_dr[0][2][0];
dgamma_dr[2][2][0] = (sin_theta2 / Theta / Theta ) * ( 2 * r * Theta + (2 * spin2 * M * sin_theta2 * ( rho2 - r * rho2_dr) ) / rho2 / rho2 );
//derivative w.r.t. theta
dgamma_dr[0][0][1] = rho2_dtheta * (2 * M * r + Lambda_r * spin2r2 ) / (Lambda_r * Lambda_r * spin2r2 * spin2r2 ) ;
dgamma_dr[0][2][1] = - 4 * spin * M * r * cos_theta * sin_theta / ( spin2 * Theta * Lambda_r + r2 * Theta * Lambda_r ) ;
dgamma_dr[1][1][1] = ( -rho2 * delta_theta_dtheta + delta_theta * rho2_dtheta ) / delta_theta / delta_theta;
dgamma_dr[2][2][0] =
(sin_theta2 / Theta / Theta) *
(2 * r * Theta +
(2 * spin2 * M * sin_theta2 * (rho2 - r * rho2_dr)) / rho2 / rho2);
// derivative w.r.t. theta
dgamma_dr[0][0][1] = rho2_dtheta * (2 * M * r + Lambda_r * spin2r2) /
(Lambda_r * Lambda_r * spin2r2 * spin2r2);
dgamma_dr[0][2][1] = -4 * spin * M * r * cos_theta * sin_theta /
(spin2 * Theta * Lambda_r + r2 * Theta * Lambda_r);
dgamma_dr[1][1][1] =
(-rho2 * delta_theta_dtheta + delta_theta * rho2_dtheta) /
delta_theta / delta_theta;
dgamma_dr[2][0][1] = dgamma_dr[0][2][1];
dgamma_dr[2][2][1] = ( 8 * spin2 * M * r * rho2 * cos_theta * sin_theta * sin_theta2 + spin2r2 * Theta * rho2 * rho2 * 2 * sin_theta * cos_theta - 2 * spin2 * M * r * sin_theta2 * sin_theta2 * rho2_dtheta ) / (Theta * Theta * rho2 * rho2 );
dgamma_dr[2][2][1] =
(8 * spin2 * M * r * rho2 * cos_theta * sin_theta * sin_theta2 +
spin2r2 * Theta * rho2 * rho2 * 2 * sin_theta * cos_theta -
2 * spin2 * M * r * sin_theta2 * sin_theta2 * rho2_dtheta) /
(Theta * Theta * rho2 * rho2);

// Now transform to the rectangular coordinates

Expand Down Expand Up @@ -239,26 +299,23 @@ class KerrdeSitter
{
metric_vars.d1_gamma[i][j][k] = 0.;

FOR2(a,b)
FOR2(a, b)
{
metric_vars.d1_gamma[i][j][k] +=
spherical_gamma[a][b] *
(jacobian_BL_to_Cart_deriv[a][i][k] *
jacobian_BL_to_Cart[b][j] +
jacobian_BL_to_Cart[a][i] *
jacobian_BL_to_Cart_deriv[b][j][k]);

}

FOR3(a, b,n)
FOR3(a, b, n)
{
metric_vars.d1_gamma[i][j][k] +=
dgamma_dr[a][b][n] * jacobian_BL_to_Cart[n][k] *
jacobian_BL_to_Cart[a][i] *
jacobian_BL_to_Cart[b][j];

}

}
}
}
Expand Down Expand Up @@ -286,7 +343,6 @@ class KerrdeSitter
{
metric_vars.K += gamma_UU[i][j] * metric_vars.K_tensor[i][j];
}

}

// Ultimate values referenced from arXiv:1011.0479v1
Expand Down
Loading

0 comments on commit 0d64bd1

Please sign in to comment.