@@ -16,6 +16,8 @@ LevelCache::LevelCache(const PolarGrid& grid, const DensityProfileCoefficients&
1616 , art_(" art" , cache_domain_geometry ? grid.numberOfNodes() : 0)
1717 , detDF_(" detDF" , cache_domain_geometry ? grid.numberOfNodes() : 0)
1818{
19+ // Pre-compute and store alpha/beta coefficients at all grid nodes to avoid
20+ // repeated expensive evaluations during runtime computations
1921 if (cache_density_profile_coefficients_) {
2022#pragma omp parallel for
2123 for (int i_r = 0 ; i_r < grid.nr (); i_r++) {
@@ -31,15 +33,20 @@ LevelCache::LevelCache(const PolarGrid& grid, const DensityProfileCoefficients&
3133 }
3234 }
3335
36+ // Pre-compute and store Jacobian matrix elements (arr, att, art, detDF) at all grid nodes
37+ // to avoid repeated coordinate transformation calculations during domain operations
3438 if (cache_domain_geometry_) {
39+ // We split the loops into two regions to better respect the
40+ // access patterns of the smoother and improve cache locality
41+
42+ // Circular Indexing Section
3543#pragma omp parallel for
3644 for (int i_r = 0 ; i_r < grid.numberSmootherCircles (); i_r++) {
3745 const double r = grid.radius (i_r);
3846 for (int i_theta = 0 ; i_theta < grid.ntheta (); i_theta++) {
39- const double theta = grid.theta (i_theta);
40- const int index = grid.index (i_r, i_theta);
41-
42- double coeff_alpha = density_profile_coefficients.alpha (r, theta);
47+ const double theta = grid.theta (i_theta);
48+ const int index = grid.index (i_r, i_theta);
49+ const double coeff_alpha = density_profile_coefficients.alpha (r, theta);
4350
4451 double arr, att, art, detDF;
4552 compute_jacobian_elements (domain_geometry_, r, theta, coeff_alpha, arr, att, art, detDF);
@@ -49,21 +56,14 @@ LevelCache::LevelCache(const PolarGrid& grid, const DensityProfileCoefficients&
4956 art_ (index) = art;
5057 }
5158 }
52-
59+ // Radial Indexing Section
5360#pragma omp parallel for
5461 for (int i_theta = 0 ; i_theta < grid.ntheta (); i_theta++) {
5562 const double theta = grid.theta (i_theta);
5663 for (int i_r = grid.numberSmootherCircles (); i_r < grid.nr (); i_r++) {
57- const double r = grid.radius (i_r);
58- const int index = grid.index (i_r, i_theta);
59-
60- double coeff_alpha;
61- if (cache_density_profile_coefficients_ && !cache_domain_geometry_) {
62- coeff_alpha = coeff_alpha_ (index);
63- }
64- else {
65- coeff_alpha = density_profile_coefficients.alpha (r, theta);
66- }
64+ const double r = grid.radius (i_r);
65+ const int index = grid.index (i_r, i_theta);
66+ const double coeff_alpha = density_profile_coefficients.alpha (r, theta);
6767
6868 double arr, att, art, detDF;
6969 compute_jacobian_elements (domain_geometry_, r, theta, coeff_alpha, arr, att, art, detDF);
@@ -76,62 +76,6 @@ LevelCache::LevelCache(const PolarGrid& grid, const DensityProfileCoefficients&
7676 }
7777}
7878
79- LevelCache::LevelCache (const Level& previous_level, const PolarGrid& current_grid)
80- : domain_geometry_(previous_level.levelCache().domainGeometry())
81- , density_profile_coefficients_(previous_level.levelCache().densityProfileCoefficients())
82- , cache_density_profile_coefficients_(previous_level.levelCache().cacheDensityProfileCoefficients())
83- , coeff_alpha_(" coeff_alpha" ,
84- previous_level.levelCache().coeff_alpha().size() > 0 ? current_grid.numberOfNodes() : 0)
85- , coeff_beta_(" coeff_beta" , previous_level.levelCache().coeff_beta().size() > 0 ? current_grid.numberOfNodes() : 0)
86- , cache_domain_geometry_(previous_level.levelCache().cacheDomainGeometry())
87- , arr_(" arr" , previous_level.levelCache().arr().size() > 0 ? current_grid.numberOfNodes() : 0)
88- , att_(" att" , previous_level.levelCache().att().size() > 0 ? current_grid.numberOfNodes() : 0)
89- , art_(" art" , previous_level.levelCache().art().size() > 0 ? current_grid.numberOfNodes() : 0)
90- , detDF_(" detDF" , previous_level.levelCache().detDF().size() > 0 ? current_grid.numberOfNodes() : 0)
91- {
92- const auto & previous_level_cache = previous_level.levelCache ();
93-
94- if (previous_level_cache.cacheDensityProfileCoefficients ()) {
95- #pragma omp parallel for
96- for (int i_r = 0 ; i_r < current_grid.nr (); i_r++) {
97- for (int i_theta = 0 ; i_theta < current_grid.ntheta (); i_theta++) {
98- const int current_index = current_grid.index (i_r, i_theta);
99- const int previous_index = previous_level.grid ().index (2 * i_r, 2 * i_theta);
100-
101- if (!previous_level_cache.cacheDomainGeometry ()) {
102- coeff_alpha_[current_index] = previous_level_cache.coeff_alpha ()[previous_index];
103- }
104- coeff_beta_[current_index] = previous_level_cache.coeff_beta ()[previous_index];
105- }
106- }
107- }
108-
109- if (previous_level_cache.cacheDomainGeometry ()) {
110- #pragma omp parallel for
111- for (int i_r = 0 ; i_r < current_grid.numberSmootherCircles (); i_r++) {
112- for (int i_theta = 0 ; i_theta < current_grid.ntheta (); i_theta++) {
113- const int current_index = current_grid.index (i_r, i_theta);
114- const int previous_index = previous_level.grid ().index (2 * i_r, 2 * i_theta);
115- arr_[current_index] = previous_level_cache.arr ()[previous_index];
116- att_[current_index] = previous_level_cache.att ()[previous_index];
117- art_[current_index] = previous_level_cache.art ()[previous_index];
118- detDF_[current_index] = previous_level_cache.detDF ()[previous_index];
119- }
120- }
121- #pragma omp parallel for
122- for (int i_theta = 0 ; i_theta < current_grid.ntheta (); i_theta++) {
123- for (int i_r = current_grid.numberSmootherCircles (); i_r < current_grid.nr (); i_r++) {
124- const int current_index = current_grid.index (i_r, i_theta);
125- const int previous_index = previous_level.grid ().index (2 * i_r, 2 * i_theta);
126- arr_[current_index] = previous_level_cache.arr ()[previous_index];
127- att_[current_index] = previous_level_cache.att ()[previous_index];
128- art_[current_index] = previous_level_cache.art ()[previous_index];
129- detDF_[current_index] = previous_level_cache.detDF ()[previous_index];
130- }
131- }
132- }
133- }
134-
13579const DensityProfileCoefficients& LevelCache::densityProfileCoefficients () const
13680{
13781 return density_profile_coefficients_;
0 commit comments