Skip to content

Commit

Permalink
Merge pull request #696 from beomki-yeo/adjust-text5
Browse files Browse the repository at this point in the history 
Revert the RIdders algorithm back to original form
  • Loading branch information
@beomki-yeo
beomki-yeo authored Mar 14, 2024
2 parents 50845bb + f8cd92c commit 7eb84c8
Show file tree
Hide file tree
Showing 2 changed files with 30 additions and 29 deletions.
21 changes: 7 additions & 14 deletions tests/integration_tests/cpu/propagator/jacobian_validation.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -56,14 +56,13 @@ namespace {
const vector3 B_z{0.f, 0.f, 1.996f * unit<scalar>::T};

// Initial delta for numerical differentiaion
const std::array<scalar, 5u> h_sizes_rect{2e0f, 2e0f, 2e-2f, 1e-3f, 1e-3f};
const std::array<scalar, 5u> h_sizes_wire{2e0f, 2e0f, 2e-2f, 1e-3f, 1e-3f};
const std::array<scalar, 5u> h_sizes_rect{3e0f, 3e0f, 2e-2f, 1e-3f, 1e-3f};
const std::array<scalar, 5u> h_sizes_wire{3e0f, 3e0f, 2e-2f, 1e-3f, 1e-3f};

// Ridders' algorithm setup
constexpr const unsigned int Nt = 50u;
const std::array<scalar, 5u> safe{3.0f, 3.0f, 3.0f, 3.0f, 3.0f};
const std::array<scalar, 5u> safe{5.0f, 5.0f, 5.0f, 5.0f, 5.0f};
const std::array<scalar, 5u> con{1.2f, 1.2f, 1.2f, 1.2f, 1.2f};
constexpr const scalar threshold_factor = 10.f;
constexpr const scalar big = std::numeric_limits<scalar>::max();

std::random_device rd;
Expand Down Expand Up @@ -161,24 +160,18 @@ struct ridders_derivative {
math::max(math::abs(Arr[j][q][p] - Arr[j][q - 1][p]),
math::abs(Arr[j][q][p] - Arr[j][q - 1][p - 1]));

const scalar V = getter::element(differentiated_jacobian, j, i);

if (errt[j] <= err[j]) {
if (complete[j] == false ||
math::abs(Arr[j][q][p]) >
threshold_factor * math::abs(V) ||
V * Arr[j][q][p] < 0.f) {

if (complete[j] == false) {
err[j] = errt[j];
getter::element(differentiated_jacobian, j, i) =
Arr[j][q][p];
/*
// Please leave this for debug
if (j == e_bound_theta && i == e_bound_phi) {
if (j == e_bound_theta && i == e_bound_loc1) {
std::cout << getter::element(
differentiated_jacobian, j, i)
<< " " << math::abs(Arr[j][q][p]) << " "
<< threshold_factor * math::abs(V) << " "
<< " " << math::abs(Arr[j][q][p])
<< std::endl;
}
*/
Expand All @@ -190,7 +183,7 @@ struct ridders_derivative {
for (unsigned int j = 0; j < 5u; j++) {
/*
// Please leave this for debug
if (j == e_bound_theta && i == e_bound_phi) {
if (j == e_bound_theta && i == e_bound_loc1) {
std::cout << getter::element(differentiated_jacobian, j, i)
<< " " << Arr[j][p][p] << " "
<< Arr[j][p - 1][p - 1] << " "
Expand Down
38 changes: 23 additions & 15 deletions tests/validation/root/rk_tolerance_comparison.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,9 +33,10 @@ double step_title_x = 0.175;
double step_title_y = 0.835;
double step_ygap = -0.05;
double x_label_offset = 0.015;
double x_margin = 1;
double y_label_offset = 0.015;
double rk_title_x_offset = 0.9;
double rk_title_x = -15.57;
double rk_title_offset_fraction = 0.0216;
double rk_title_y = 13.765;
double rk_ygap = -0.8;
double rk_header_text_size = 0.046;
Expand All @@ -56,7 +57,7 @@ double pad_y0 = 0.005f;
double pad_y1 = 1.f;
double ymin = -2;
double ymax = 4.;
std::array<double, 4u> ldim{0.189233, 0.509404, 0.907268, 0.954545};
std::array<double, 4u> ldim{0.189233, 0.510081, 0.907268, 0.954545};
} // namespace

std::vector<std::string> create_labels() {
Expand Down Expand Up @@ -238,16 +239,17 @@ void draw_graphs(const std::string header_title, const std::string geom_title,
mg->GetYaxis()->SetLabelSize(label_font_size_rk_tol);
mg->GetYaxis()->SetRangeUser(yaxis_min - yaxis_margin,
yaxis_max + yaxis_margin);
mg->GetYaxis()->SetNdivisions(406, kFALSE);
mg->GetYaxis()->SetLabelOffset(0.01);
mg->GetYaxis()->SetTitleFont(title_font);
mg->GetXaxis()->SetTitleFont(title_font);

mg->GetYaxis()->SetLabelSize(0);
mg->GetYaxis()->SetTickLength(0);
double x_min = x_vec.front();
double x_max = x_vec.back();

if (x_vec.size() > 10) {
mg->GetXaxis()->SetLimits(x_vec.front() - 1, x_vec.back() + 1);
x_min = x_min - x_margin;
x_max = x_max + x_margin;
mg->GetXaxis()->SetLimits(x_min, x_max);
mg->GetXaxis()->SetLabelSize(0);
mg->GetXaxis()->SetTickLength(0);

Expand All @@ -260,25 +262,31 @@ void draw_graphs(const std::string header_title, const std::string geom_title,
ga->SetLabelSize(label_font_size_rk_tol);
ga->SetLabelOffset(-0.0065);
ga->Draw();

mg->GetYaxis()->SetLabelSize(0);
mg->GetYaxis()->SetTickLength(0);
auto ga_y = new TGaxis(x_min, yaxis_min, x_min, yaxis_max, yaxis_min,
yaxis_max, 406, "N");
ga_y->SetLabelFont(label_font);
ga_y->SetLabelOffset(0.02);
ga_y->SetLabelSize(label_font_size_rk_tol);
ga_y->Draw();

} else {
mg->Draw("APL");
}

auto ga_y = new TGaxis(x_vec.front() - 1, yaxis_min, x_vec.front() - 1,
yaxis_max, yaxis_min, yaxis_max, 406, "N");
ga_y->SetLabelFont(label_font);
ga_y->SetLabelOffset(0.02);
ga_y->SetLabelSize(label_font_size_rk_tol);
ga_y->Draw();

TLegendEntry* header =
(TLegendEntry*)legend->GetListOfPrimitives()->First();
header->SetTextFont(22);
header->SetTextSize(.033);

double rk_title_deltaX =
(x_vec.back() - x_vec.front()) * rk_title_offset_fraction;

legend->Draw();
draw_text(rk_title_x, rk_title_y, rk_ygap, rk_header_text_size,
rk_geom_text_size, header_title, geom_title);
draw_text(rk_title_deltaX + x_vec.front(), rk_title_y, rk_ygap,
rk_header_text_size, rk_geom_text_size, header_title, geom_title);
}

void draw_mean_step_size(const std::string header_title,
Expand Down

0 comments on commit 7eb84c8

Please sign in to comment.