diff --git a/Source/InitialConditions/BosonStars/BosonStarSolution.impl.hpp b/Source/InitialConditions/BosonStars/BosonStarSolution.impl.hpp index b3e9b7c..0e3fd10 100644 --- a/Source/InitialConditions/BosonStars/BosonStarSolution.impl.hpp +++ b/Source/InitialConditions/BosonStars/BosonStarSolution.impl.hpp @@ -15,7 +15,7 @@ #include #include -// 1D boson-star solver +// 1D boson-star solver BosonStarSolution::BosonStarSolution() {} @@ -25,9 +25,9 @@ void BosonStarSolution::main() { if (BS_verbosity) { - pout() << "-----------------------------------------------" << endl; - pout() << "I am running iteration # " << iter << endl; - pout() << "-----------------------------------------------" << endl; + pout() << "-----------------------------------------------" << endl; + pout() << "I am running iteration # " << iter << endl; + pout() << "-----------------------------------------------" << endl; } // Set the initial conditions @@ -53,10 +53,11 @@ void BosonStarSolution::main() // Force the scalar field to zero after the point the amplitude diverges force_to_zero(matching_index); rk4_asymp(matching_index, false, omega_true); - rk4_asymp(matching_index, true, - omega_true); // (true) uses large radius adaptive stepsize to - // get asymptotics by integrating vacuum metric - // out to huge radius. + rk4_asymp( + matching_index, true, + omega_true); // (true) uses large radius adaptive stepsize to get + // asymptotics by integrating vacuum metric out to huge + // radius. PSI_INF = psi[gridsize - 1]; OM_INF = omega[gridsize - 1]; @@ -71,13 +72,12 @@ void BosonStarSolution::main() // Useful, when you may want to iterate through the for loop a few // times. if (fabs(PSI_INF - 1.) < 1e-05 && fabs(OM_INF - 1.) < 1e-05) - { + { if (BS_verbosity) { - pout() << "Found the desired solution... Central Density : " - << A[0] << ", PSI0 : " << psi[0] - << ", OM0 : " << omega[0] << ", w : " << sqrt(omega_true) - << endl; + pout() << "Found the desired solution... Central Density : " << A[0] << ", PSI0 : " << psi[0] + << ", OM0 : " << omega[0] << ", w : " << sqrt(omega_true) + << endl; } break; } @@ -88,26 +88,22 @@ void BosonStarSolution::main() if (BS_verbosity) { - pout() - << "For updated variables I am now using... Central Density : " - << A[0] << ", PSI0 : " << PSC << ", OM0 : " << OMC - << ", w : " << sqrt(omega_true) << endl; + pout() << "For updated variables I am now using... Central Density : " << A[0] << ", PSI0 : " << PSC + << ", OM0 : " << OMC << ", w : " << sqrt(omega_true) << endl; } - if (iter == niter - 1) + if (iter == niter-1) { - MayDay::Error("Ooopsies... I have reached the maximum number of " - "iterations and did not " - "find a BS solution ¯|_(ツ)_|¯. Have a look at " - "niter parameter and increase it, if needed."); + MayDay::Error("Ooopsies... I have reached the maximum number of iterations and did not " + "find a BS solution ¯|_(ツ)_|¯. Have a look at niter parameter and increase it, if needed."); } } if (BS_verbosity) { - pout() << "-----------------------------------------------" << endl; - pout() << "Computing the diagnostics of the solution" << endl; - pout() << "-----------------------------------------------" << endl; + pout() << "-----------------------------------------------" << endl; + pout() << "Computing the diagnostics of the solution" << endl; + pout() << "-----------------------------------------------" << endl; } rk4_match(matching_index, false, omega_true); @@ -120,13 +116,13 @@ void BosonStarSolution::main() if (BS_verbosity) { - pout() << "-----------------------------------------------" << endl; - pout() << "Central Density : " << A[0] << endl; - pout() << "ADM mass : " << adm_mass[gridsize - 1] << endl; - pout() << "Aspect mass : " << boson_mass[gridsize - 1] << endl; - pout() << "Radius : " << radius << endl; - pout() << "Compactness : " << compactness_value << endl; - pout() << "W : " << sqrt(omega_true) << endl; + pout() << "-----------------------------------------------" << endl; + pout() << "Central Density : " << A[0] << endl; + pout() << "ADM mass : " << adm_mass[gridsize - 1] << endl; + pout() << "Aspect mass : " << boson_mass[gridsize - 1] << endl; + pout() << "Radius : " << radius << endl; + pout() << "Compactness : " << compactness_value << endl; + pout() << "W : " << sqrt(omega_true) << endl; } } @@ -168,7 +164,7 @@ void BosonStarSolution::truncate_solution() } } -// Sets scalar field and gradient to zero after given point +// Sets scalar field and gradient to zero after given point void BosonStarSolution::force_to_zero(const int iter_crit) { for (int i = iter_crit + 1; i < gridsize; ++i) @@ -178,8 +174,8 @@ void BosonStarSolution::force_to_zero(const int iter_crit) } } -// Find the index (radius) where the amplitude starts to diverge, this will be -// used for matching in the asymptotic regime +// Find the index (radius) where the amplitude starts to diverge, this will be used +// for matching in the asymptotic regime int BosonStarSolution::find_matching_index() { int matching_index; @@ -307,11 +303,9 @@ int BosonStarSolution::found_zero_crossing() // asumptotics later in the main(). void BosonStarSolution::rk4(const double ww_) { - double k1 = 0, k2 = 0, k3 = 0, k4 = 0, q1 = 0, q2 = 0, q3 = 0, - q4 = 0; // for RK steps - double o1 = 0, o2 = 0, o3 = 0, o4 = 0, s1 = 0, s2 = 0, s3 = 0, s4 = 0, - r1 = 0, r2 = 0, r3 = 0, r4 = 0; // for RK steps - double x_ = 0., h = dx / 2.; // for step-size + double k1 = 0, k2 = 0, k3 = 0, k4 = 0, q1 = 0, q2 = 0, q3 = 0, q4 = 0; // for RK steps + double o1 = 0, o2 = 0, o3 = 0, o4 = 0, s1 = 0, s2 = 0, s3 = 0, s4 = 0, r1 = 0, r2 = 0, r3 = 0, r4 = 0; // for RK steps + double x_ = 0., h = dx / 2.; // for step-size const double DX = dx; double DX_ = DX; int index, jmax = 0; @@ -409,15 +403,15 @@ void BosonStarSolution::rk4(const double ww_) } } -// Intergartion starts at point (iter), enforcing scalar field to be in vacuum. -// The integral here is adaptive in the sense that it accelerates ar later -// radius in order to find correct asymptotic behaviour. It will give an error, -// if the radius reached is below 8e7. +// Intergartion starts at point (iter), enforcing scalar field to be in vacuum. The integral here is +// adaptive in the sense that it accelerates ar later radius in order to find correct +// asymptotic behaviour. It will give an error, if the radius reached is below 8e7. void BosonStarSolution::rk4_asymp(const int iter, const bool adaptive, const double ww_) { - double o1 = 0, o2 = 0, o3 = 0, o4 = 0, s1 = 0, s2 = 0, s3 = 0, s4 = 0, - r1 = 0, r2 = 0, r3 = 0, r4 = 0; // for RK steps + double o1 = 0, o2 = 0, o3 = 0, o4 = 0, s1 = 0, s2 = 0, s3 = 0, s4 = 0, r1 = 0, r2 = 0, r3 = 0, r4 = 0; // for RK steps + double k1 = 0, k2 = 0, k3 = 0, k4 = 0, q1 = 0, q2 = 0, q3 = 0, + q4 = 0; // for RK steps double x_ = iter * dx, h, delta = (double)gridsize; const double DX = dx; double DX_ = DX; @@ -496,18 +490,17 @@ void BosonStarSolution::rk4_asymp(const int iter, const bool adaptive, } if (adaptive and x_ < 8e7) - { + { pout() << "Radius is " << x_ << endl; - MayDay::Error("Asymptotic Radius Too Small"); + MayDay::Error("Asymptotic Radius Too Small"); } } -// Matches the solution to the coorect asymptotics at index iter +// Matches the solution to the coorect asymptotics at index iter void BosonStarSolution::rk4_match(const int iter, const bool adaptive, const double ww_) { - double o1 = 0, o2 = 0, o3 = 0, o = 0, s1 = 0, s2 = 0, s3 = 0, s4 = 0, - r1 = 0, r2 = 0, r3 = 0, r4 = 0; + double o1 = 0, o2 = 0, o3 = 0, o4 = 0, s1 = 0, s2 = 0, s3 = 0, s4 = 0, r1 = 0, r2 = 0, r3 = 0, r4 = 0; double r, dr; double c1, c2; double Amp, eta, mass, arealr, om0, epsilon; @@ -532,8 +525,6 @@ void BosonStarSolution::rk4_match(const int iter, const bool adaptive, { dr = radius_array[i] - radius_array[i - 1]; - h = DX_ / 2.; - // 1st RK step r = radius_array[i - 1]; om0 = pow(1 / omega[i - 1], 2); @@ -609,7 +600,7 @@ void BosonStarSolution::rk4_match(const int iter, const bool adaptive, } } -// RHS for BS amplitude +// RHS for BS amplitude double BosonStarSolution::A_RHS(const double x, const double A, const double DA, const double PSI, const double DPSI, const double OM, const double ww_) @@ -618,7 +609,7 @@ double BosonStarSolution::A_RHS(const double x, const double A, const double DA, return RHS; } -// RHS for the derivative of the BS amplitude +// RHS for the derivative of the BS amplitude double BosonStarSolution::DA_RHS(const double x, const double A, const double DA, const double PSI, const double DPSI, const double OM, @@ -630,7 +621,7 @@ double BosonStarSolution::DA_RHS(const double x, const double A, DA * (DOM / OM + DPSI / PSI + 2. / r); } -// RHS for the conformal factor +// RHS for the conformal factor double BosonStarSolution::PSI_RHS(const double x, const double A, const double DA, const double PSI, const double DPSI, const double OM, @@ -640,7 +631,7 @@ double BosonStarSolution::PSI_RHS(const double x, const double A, return RHS; } -// RHS for the derivative of the conformal factor +// RHS for the derivative of the conformal factor double BosonStarSolution::DPSI_RHS(const double x, const double A, const double DA, const double PSI, const double DPSI, const double OM, @@ -653,7 +644,7 @@ double BosonStarSolution::DPSI_RHS(const double x, const double A, ww_ * A * A * PSI * PSI / (OM * OM)); } -// RHS for the lapse +// RHS for the lapse double BosonStarSolution::OMEGA_RHS(const double x, const double A, const double DA, const double PSI, const double DPSI, const double OM, @@ -667,7 +658,7 @@ double BosonStarSolution::OMEGA_RHS(const double x, const double A, DPSI - 0.5 * x * DPSI * DPSI / PSI); } -// Boson star potential +// Boson star potential double BosonStarSolution::V(const double A) { if (!solitonic) @@ -680,7 +671,7 @@ double BosonStarSolution::V(const double A) } } -// Derivative of the potential +// Derivative of the potential double BosonStarSolution::DV(const double A) { if (!solitonic) @@ -693,7 +684,7 @@ double BosonStarSolution::DV(const double A) } } -// Find the aspect mass +// Find the aspect mass void BosonStarSolution::calculate_aspect_mass() { for (int i = 0; i < gridsize; ++i) @@ -702,7 +693,7 @@ void BosonStarSolution::calculate_aspect_mass() } } -// Find the ADM mass +// Find the ADM mass void BosonStarSolution::calculate_adm_mass() { for (int i = 0; i < gridsize; ++i) @@ -711,7 +702,7 @@ void BosonStarSolution::calculate_adm_mass() } } -// Find the 99% radius +// Find the 99% radius double BosonStarSolution::calculate_radius() { int i; @@ -723,7 +714,7 @@ double BosonStarSolution::calculate_radius() return radius_array[i + 1]; } -// 4th order error (cubic interpolation) for the amplitude +// 4th order error (cubic interpolation) for the amplitude double BosonStarSolution::get_A_interp(const double r) const { int iter = (int)floor( @@ -753,7 +744,7 @@ double BosonStarSolution::get_A_interp(const double r) const return interpolated_value; } -// 4th order error (cubic interpolation) for the derivative of the amplitude +// 4th order error (cubic interpolation) for the derivative of the amplitude double BosonStarSolution::get_dA_interp(const double r) const { int iter = (int)floor( @@ -762,14 +753,14 @@ double BosonStarSolution::get_dA_interp(const double r) const (r / dx) - floor(r / dx) - 0.5; // fraction from midpoint of two values, // a = +- 1/2 is the nearest gridpoints double interpolated_value = 0, f1, f2, f3, f4; - f1 = ((iter == 0) ? dA[1] : dA[iter - 1]); + f1 = ((iter == 0) ? dA[1] : dA[iter - 1]); f2 = dA[iter]; f3 = dA[iter + 1]; f4 = dA[iter + 2]; if (iter > gridsize - 3) { - MayDay::Error("FArrayBox domain exceeding star radius!"); + MayDay::Error("FArrayBox domain exceeding star radius!"); } // do the cubic spline, from mathematica script written by Robin @@ -781,7 +772,7 @@ double BosonStarSolution::get_dA_interp(const double r) const return interpolated_value; } -// 4th order error (cubic interpolation) for the lapse +// 4th order error (cubic interpolation) for the lapse double BosonStarSolution::get_lapse_interp(const double r) const { int iter = (int)floor( @@ -811,7 +802,7 @@ double BosonStarSolution::get_lapse_interp(const double r) const return interpolated_value; } -// 4th order error (cubic interpolation) for the conformal factor +// 4th order error (cubic interpolation) for the conformal factor double BosonStarSolution::get_psi_interp(const double r) const { int iter = (int)floor( @@ -820,7 +811,7 @@ double BosonStarSolution::get_psi_interp(const double r) const (r / dx) - floor(r / dx) - 0.5; // fraction from midpoint of two values, // a = +- 1/2 is the nearest gridpoints double interpolated_value = 0, f1, f2, f3, f4; - f1 = ((iter == 0) ? psi[1] : psi[iter - 1]); + f1 = ((iter == 0) ? psi[1] : psi[iter - 1]); f2 = psi[iter]; f3 = psi[iter + 1]; f4 = psi[iter + 2]; @@ -841,8 +832,7 @@ double BosonStarSolution::get_psi_interp(const double r) const return interpolated_value; } -// 4th order error (cubic interpolation) for the derivative of the conformal -// factor +// 4th order error (cubic interpolation) for the derivative of the conformal factor double BosonStarSolution::get_dpsi_interp(const double r) const { int iter = (int)floor( @@ -851,7 +841,7 @@ double BosonStarSolution::get_dpsi_interp(const double r) const (r / dx) - floor(r / dx) - 0.5; // fraction from midpoint of two values, // a = +- 1/2 is the nearest gridpoints double interpolated_value = 0, f1, f2, f3, f4; - f1 = ((iter == 0) ? dpsi[1] : dpsi[iter - 1]); + f1 = ((iter == 0) ? dpsi[1] : dpsi[iter - 1]); f2 = dpsi[iter]; f3 = dpsi[iter + 1]; f4 = dpsi[iter + 2]; @@ -870,7 +860,7 @@ double BosonStarSolution::get_dpsi_interp(const double r) const return interpolated_value; } -// 4th order error (cubic interpolation) for the derivative of the lapse +// 4th order error (cubic interpolation) for the derivative of the lapse double BosonStarSolution::get_dlapse_interp(const double r) const { int iter = (int)floor( @@ -898,25 +888,25 @@ double BosonStarSolution::get_dlapse_interp(const double r) const return interpolated_value; } -// Get the BS frequency (note that the right frequency is the sqaured root of -// omega_true) +// Get the BS frequency (note that the right frequency is the sqaured root of omega_true) double BosonStarSolution::get_BSfrequency() const { return sqrt(omega_true); } -// Initialise solver +// Initialise solver void BosonStarSolution::set_initialcondition_params( BosonStar_params_t m_params_BosonStar, Potential::params_t m_params_potential, const double max_r) { gridsize = m_params_BosonStar.gridpoints; - adaptive_buffer = 0.; // number of gridpoints to intergate more carefully - A.resize(gridsize); // scalar field modulus - dA.resize(gridsize); // scalar field modulus gradient - psi.resize(gridsize); // conformal factor - dpsi.resize(gridsize); // conformal factor gradient - omega.resize(gridsize); // lapse - radius_array.resize(gridsize); // radius array - boson_mass.resize(gridsize); // mass of the BS - adm_mass.resize(gridsize); // ADM BS mass + adaptive_buffer = + 0.; // number of gridpoints to intergate more carefully + A.resize(gridsize); // scalar field modulus + dA.resize(gridsize); // scalar field modulus gradient + psi.resize(gridsize); // conformal factor + dpsi.resize(gridsize); // conformal factor gradient + omega.resize(gridsize); // lapse + radius_array.resize(gridsize); // radius array + boson_mass.resize(gridsize); // mass of the BS + adm_mass.resize(gridsize); // ADM BS mass A0 = m_params_BosonStar.central_amplitude_CSF; mu = m_params_potential.scalar_mass * m_params_potential.scalar_mass; @@ -932,12 +922,12 @@ void BosonStarSolution::set_initialcondition_params( dx = L / double((gridsize - 1)); } -// IO routines +// IO routines void BosonStarSolution::output_csv() { std::ofstream A_file, dA_file, psi_file, dpsi_file, omega_file, r_file, mass_file; - + A_file.open("A.csv"); dA_file.open("dA.csv"); psi_file.open("psi.csv"); diff --git a/Source/utils/StarTracker.cpp b/Source/utils/StarTracker.cpp index 29cd0dc..e91faec 100644 --- a/Source/utils/StarTracker.cpp +++ b/Source/utils/StarTracker.cpp @@ -3,7 +3,7 @@ * Please refer to LICENSE in GRChombo's root directory. */ -#include "StarTracker.hpp" +#include "StarTracker.hpp" #include "DimensionDefinitions.hpp" #include "FittingMethod.hpp" #include "FittingModel.hpp"