20283049329

Committed 16 Dec 2025 09:19PM UTC coverage: 88.351% (-0.8%) from 89.184%

Build # 20283049329

Build Type

Pull #91

github

Committed by

web-flow

Commit Message

Merge d2f36f095 into a73a760a4

Pull Request Pull Request #91: Implement multi-threading in NativeRunner

Run Details

366 of 468 new or added lines in 11 files covered. (78.21%)

34 existing lines in 6 files now uncovered.

6007 of 6799 relevant lines covered (88.35%)

7475791.74 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

44.58

/coretrace/simulation_runner/native_runner/pt_optimizations.cpp


#include "pt_optimizations.hpp"

#include <algorithm>

#include "constants.hpp"
#include "simulation_data_export.hpp"

namespace SolTrace::NativeRunner
{

        void SetupPTOptimizations(
                // system info
                TSystem *System, const bool AsPowerTower,
                // outputs
                st_hash_tree &sun_hash,
                st_hash_tree &rec_hash,
                double (&reccm_helio)[3])
        {
                // Calculate the center of mass of the receiver stage (StageList[1]) in
                // heliostat stage coordinates.
                double reccm[] = {0., 0., 0.};
                int nelrec = 0;
                if (AsPowerTower)
                {
                        for (uint_fast64_t j = 0; j < System->StageList[1]->ElementList.size(); j++)
                        {
                                TElement *el = System->StageList[1]->ElementList.at(j).get();

                                // if (!el->Enabled)
                                //         continue;

                                nelrec++;

                                for (int jj = 0; jj < 3; jj++)
                                        reccm[jj] += el->Origin[jj];
                        }
                        for (int jj = 0; jj < 3; jj++)
                                reccm[jj] /= (double)nelrec; // average

                        // Transform to reference
                        double dum1[] = {0., 0., 1.};
                        double dum2[3];
                        double reccm_global[3];
                        TransformToReference(reccm, dum1, System->StageList[1]->Origin,
                                                                 System->StageList[1]->RLocToRef, reccm_global, dum2);

                        // Transform to local (heliostat). reccm_helio is the x,y,z position
                        // of the receiver centroid in heliostat stage coordinates.
                        TransformToLocal(reccm_global, dum1, System->StageList[0]->Origin,
                                                         System->StageList[0]->RRefToLoc, reccm_helio, dum2);
                }
                // Create an array that stores the element address and the projected size
                // in polar coordinates
                std::vector<eprojdat> el_proj_dat;
                el_proj_dat.reserve(System->StageList[0]->ElementList.size());

                // calculate the smallest zone size. This should be on the order of the
                // largest element in the stage. load stage 0 elements into the mesh
                double d_elm_max = -9.e9;
                double d_elm;

                for (uint_fast64_t i = 0; i < System->StageList[0]->ElementList.size(); i++)
                {
                        TElement *el = System->StageList[0]->ElementList.at(i).get();

                        // el->element_number = i + 1; // use index for element number

                        d_elm = el->aperture->diameter_circumscribed_circle();

                        // double d_elm;

                        // switch (el->ShapeIndex)
                        // {
                        // // circular aperture
                        // case 'c':
                        // case 'C':
                        //         // hexagonal aperture
                        // case 'h':
                        // case 'H':
                        //         // triangular aperture
                        // case 't':
                        // case 'T':
                        //         d_elm = el->ParameterA;
                        //         break;
                        //         // rectangular aperture
                        // case 'r':
                        // case 'R':
                        //         d_elm = sqrt(el->ParameterA * el->ParameterA + el->ParameterB * el->ParameterB);
                        //         break;
                        //         // annular aperture
                        // case 'a':
                        // case 'A':
                        //         d_elm = el->ParameterB;
                        //         break;
                        // case 'l':
                        // case 'L':
                        //         // off axis aperture section of line focus trough  or cylinder
                        //         d_elm = sqrt(el->ParameterB * el->ParameterB * 4. + el->ParameterC * el->ParameterC);
                        //         break;
                        //         // Irregular triangle
                        // case 'i':
                        // case 'I':
                        //         // irregular quadrilateral
                        // case 'q':
                        // case 'Q':
                        // {
                        //         double xmax = fmax(el->ParameterA, fmax(el->ParameterC, el->ParameterE));
                        //         double xmin = fmin(el->ParameterA, fmin(el->ParameterC, el->ParameterE));
                        //         double ymax = fmax(el->ParameterB, fmax(el->ParameterD, el->ParameterF));
                        //         double ymin = fmin(el->ParameterB, fmin(el->ParameterD, el->ParameterF));

                        //         if (el->ShapeIndex == 'q' || el->ShapeIndex == 'Q')
                        //         {
                        //                 xmax = fmax(xmax, el->ParameterG);
                        //                 xmin = fmin(xmin, el->ParameterG);
                        //                 ymax = fmax(ymax, el->ParameterH);
                        //                 ymin = fmin(ymin, el->ParameterH);
                        //         }

                        //         double dx = xmax - xmin;
                        //         double dy = ymax - ymin;

                        //         d_elm = sqrt(dx * dx + dy * dy);

                        //         break;
                        // }
                        // default:
                        //         break;
                        // }

                        d_elm_max = fmax(d_elm_max, d_elm);

                        if (AsPowerTower)
                        {
                                // Calculate the distance from the receiver to the element and the max projected size
                                double dX[3];
                                for (int jj = 0; jj < 3; jj++)
                                        dX[jj] = el->Origin[jj] - reccm_helio[jj]; // vector from receiver to heliostat (not unitized)
                                double r_elm = 0.;
                                for (int jj = 0; jj < 3; jj++)
                                        r_elm += dX[jj] * dX[jj];
                                r_elm = sqrt(r_elm);                           // vector length
                                double d_elm_proj = d_elm / r_elm; // Projected size of the element from the view of the receiver (radians)

                                // calculate az,zen coordinate
                                double az, zen;
                                az = atan2(dX[0] / r_elm, dX[1] / r_elm); // Az coordinate of the heliostat from the receiver's perspective
                                zen = asin(dX[2] / r_elm);                                  // Zen coordinate """"

                                el_proj_dat.push_back(eprojdat(el, d_elm_proj, az, zen));
                        }
                }

                if (AsPowerTower)
                {
                        // Sort the polar projections by size, largest to smallest
                        std::sort(el_proj_dat.begin(), el_proj_dat.end(), eprojdat_compare_refactored);
                }

                // set up the layout data object that provides configuration details for
                // the hash tree
                KDLayoutData sun_ld;
                sun_ld.xlim[0] = System->Sun.MinXSun;
                sun_ld.xlim[1] = System->Sun.MaxXSun;
                sun_ld.ylim[0] = System->Sun.MinYSun;
                sun_ld.ylim[1] = System->Sun.MaxYSun;
                sun_ld.min_unit_dx = d_elm_max;
                sun_ld.min_unit_dy = d_elm_max;

                sun_hash.create_mesh(sun_ld);

                // load stage 0 elements into the mesh
                for (uint_fast64_t i = 0; i < System->StageList[0]->ElementList.size(); i++)
                {
                        TElement *el = System->StageList[0]->ElementList.at(i).get();
                        sun_hash.add_object((void *)el, el->PosSunCoords[0], el->PosSunCoords[1]);
                }

                // calculate and associate neighbors with each zone
                sun_hash.add_neighborhood_data();

                if (AsPowerTower)
                {
                        // Set things up for the polar coordinate tree
                        KDLayoutData rec_ld;
                        rec_ld.xlim[0] = -PI;
                        rec_ld.xlim[1] = PI;
                        rec_ld.ylim[0] = -PI / 2.;
                        rec_ld.ylim[1] = PI / 2.;
                        // use smallest element to set the minimum size
                        rec_ld.min_unit_dx = rec_ld.min_unit_dy = el_proj_dat.back().d_proj; // radians at equator

                        rec_hash.create_mesh(rec_ld);

                        // load stage 0 elements into the receiver mesh in the order of largest projection to smallest
                        for (int i = 0; i < el_proj_dat.size(); i++)
                        {
                                eprojdat *D = &el_proj_dat.at(i);

                                // Calculate the angular span of the element
                                double angspan[2];
                                double adjmult = 1.5;
                                angspan[0] = D->d_proj / cos(fabs(D->zen)) * adjmult; // azimuthal span
                                angspan[0] = fmin(angspan[0], 2. * PI);                                  // limit to circumference
                                angspan[1] = D->d_proj / PI * adjmult;                                  // zenithal span
                                rec_hash.add_object((void *)D->el_addr, D->az, D->zen, angspan);
                        }

                        // associate neighbors with each zone
                        rec_hash.add_neighborhood_data();
                }
        }

        uint_fast64_t GetPTElements(
                // system info
                const bool AsPowerTower,

                // Stage info
                // const TStage *Stage
                const tstage_ptr Stage,
                const int i,

                // Ray info
                const bool in_multi_hit_loop,
                const double (&PosRayStage)[3],
                const double (&reccm_helio)[3],
                st_hash_tree *rec_hash,

                const std::vector<void *> &sunint_elements,

                // Outputs
                std::vector<void *> &reflint_elements,
                bool &has_elements)
        {
                uint_fast64_t nintelements = 0;

                if (i == 0)
                {
                        if (in_multi_hit_loop)
                        {
                                if (AsPowerTower)
                                {
                                        //>=Second time through - checking for first stage multiple element interactions

                                        // get ray position in receiver polar coordinates
                                        double raypvec[3];
                                        for (int jj = 0; jj < 3; jj++)
                                                raypvec[jj] = PosRayStage[jj] - reccm_helio[jj];
                                        double raypvecmag = sqrt(raypvec[0] * raypvec[0] + raypvec[1] * raypvec[1] + raypvec[2] * raypvec[2]);
                                        double raypol[2];
                                        raypol[0] = atan2(raypvec[0], raypvec[1]);
                                        raypol[1] = asin(raypvec[2] / raypvecmag);
                                        // get elements in the vicinity of the ray's polar coordinates
                                        reflint_elements.clear();
                                        rec_hash->get_all_data_at_loc(reflint_elements, raypol[0], raypol[1]);
                                        nintelements = reflint_elements.size();
                                        has_elements = nintelements > 0;
                                }
                                else
                                {
                                        nintelements = Stage->ElementList.size();
                                }
                        }
                        else
                        {
                                // First time through - checking for sun ray intersections
                                if (has_elements)
                                        nintelements = sunint_elements.size();
                                else
                                        nintelements = 0;
                        }
                }
                else
                        nintelements = Stage->ElementList.size();

                return nintelements;
        }

} // namespace SolTrace::NativeRunner

1
2	#include "pt_optimizations.hpp"
3
4	#include <algorithm>
5
6	#include "constants.hpp"
7	#include "simulation_data_export.hpp"
8
9	namespace SolTrace::NativeRunner
10	{
11
12	void SetupPTOptimizations(	4✔
13	// system info
14	TSystem *System, const bool AsPowerTower,
15	// outputs
16	st_hash_tree &sun_hash,
17	st_hash_tree &rec_hash,
18	double (&reccm_helio)[3])
19	{
20	// Calculate the center of mass of the receiver stage (StageList[1]) in
21	// heliostat stage coordinates.
22	double reccm[] = {0., 0., 0.};	4✔
23	int nelrec = 0;	4✔
24	if (AsPowerTower)	4✔
25	{
NEW 26	for (uint_fast64_t j = 0; j < System->StageList[1]->ElementList.size(); j++)	×
27	{
NEW 28	TElement *el = System->StageList[1]->ElementList.at(j).get();	×
29
30	// if (!el->Enabled)
31	// continue;
32
NEW 33	nelrec++;	×
34
NEW 35	for (int jj = 0; jj < 3; jj++)	×
NEW 36	reccm[jj] += el->Origin[jj];	×
37	}
UNCOV 38	for (int jj = 0; jj < 3; jj++)	×
NEW 39	reccm[jj] /= (double)nelrec; // average	×
40
41	// Transform to reference
NEW 42	double dum1[] = {0., 0., 1.};	×
43	double dum2[3];
44	double reccm_global[3];
NEW 45	TransformToReference(reccm, dum1, System->StageList[1]->Origin,	×
NEW 46	System->StageList[1]->RLocToRef, reccm_global, dum2);	×
47
48	// Transform to local (heliostat). reccm_helio is the x,y,z position
49	// of the receiver centroid in heliostat stage coordinates.
NEW 50	TransformToLocal(reccm_global, dum1, System->StageList[0]->Origin,	×
NEW 51	System->StageList[0]->RRefToLoc, reccm_helio, dum2);	×
52	}
53	// Create an array that stores the element address and the projected size
54	// in polar coordinates
55	std::vector<eprojdat> el_proj_dat;	4✔
56	el_proj_dat.reserve(System->StageList[0]->ElementList.size());	4✔
57
58	// calculate the smallest zone size. This should be on the order of the
59	// largest element in the stage. load stage 0 elements into the mesh
60	double d_elm_max = -9.e9;	4✔
61	double d_elm;
62
63	for (uint_fast64_t i = 0; i < System->StageList[0]->ElementList.size(); i++)	12,590✔
64	{
65	TElement *el = System->StageList[0]->ElementList.at(i).get();	12,586✔
66
67	// el->element_number = i + 1; // use index for element number
68
69	d_elm = el->aperture->diameter_circumscribed_circle();	12,586✔
70
71	// double d_elm;
72
73	// switch (el->ShapeIndex)
74	// {
75	// // circular aperture
76	// case 'c':
77	// case 'C':
78	// // hexagonal aperture
79	// case 'h':
80	// case 'H':
81	// // triangular aperture
82	// case 't':
83	// case 'T':
84	// d_elm = el->ParameterA;
85	// break;
86	// // rectangular aperture
87	// case 'r':
88	// case 'R':
89	// d_elm = sqrt(el->ParameterA * el->ParameterA + el->ParameterB * el->ParameterB);
90	// break;
91	// // annular aperture
92	// case 'a':
93	// case 'A':
94	// d_elm = el->ParameterB;
95	// break;
96	// case 'l':
97	// case 'L':
98	// // off axis aperture section of line focus trough or cylinder
99	// d_elm = sqrt(el->ParameterB * el->ParameterB * 4. + el->ParameterC * el->ParameterC);
100	// break;
101	// // Irregular triangle
102	// case 'i':
103	// case 'I':
104	// // irregular quadrilateral
105	// case 'q':
106	// case 'Q':
107	// {
108	// double xmax = fmax(el->ParameterA, fmax(el->ParameterC, el->ParameterE));
109	// double xmin = fmin(el->ParameterA, fmin(el->ParameterC, el->ParameterE));
110	// double ymax = fmax(el->ParameterB, fmax(el->ParameterD, el->ParameterF));
111	// double ymin = fmin(el->ParameterB, fmin(el->ParameterD, el->ParameterF));
112
113	// if (el->ShapeIndex == 'q' \|\| el->ShapeIndex == 'Q')
114	// {
115	// xmax = fmax(xmax, el->ParameterG);
116	// xmin = fmin(xmin, el->ParameterG);
117	// ymax = fmax(ymax, el->ParameterH);
118	// ymin = fmin(ymin, el->ParameterH);
119	// }
120
121	// double dx = xmax - xmin;
122	// double dy = ymax - ymin;
123
124	// d_elm = sqrt(dx * dx + dy * dy);
125
126	// break;
127	// }
128	// default:
129	// break;
130	// }
131
132	d_elm_max = fmax(d_elm_max, d_elm);	12,586✔
133
134	if (AsPowerTower)	12,586✔
135	{
136	// Calculate the distance from the receiver to the element and the max projected size
137	double dX[3];
NEW 138	for (int jj = 0; jj < 3; jj++)	×
NEW 139	dX[jj] = el->Origin[jj] - reccm_helio[jj]; // vector from receiver to heliostat (not unitized)	×
NEW 140	double r_elm = 0.;	×
NEW 141	for (int jj = 0; jj < 3; jj++)	×
NEW 142	r_elm += dX[jj] * dX[jj];	×
NEW 143	r_elm = sqrt(r_elm); // vector length	×
NEW 144	double d_elm_proj = d_elm / r_elm; // Projected size of the element from the view of the receiver (radians)	×
145
146	// calculate az,zen coordinate
147	double az, zen;
NEW 148	az = atan2(dX[0] / r_elm, dX[1] / r_elm); // Az coordinate of the heliostat from the receiver's perspective	×
NEW 149	zen = asin(dX[2] / r_elm); // Zen coordinate """"	×
150
NEW 151	el_proj_dat.push_back(eprojdat(el, d_elm_proj, az, zen));	×
152	}
153	}
154
155	if (AsPowerTower)	4✔
156	{
157	// Sort the polar projections by size, largest to smallest
NEW 158	std::sort(el_proj_dat.begin(), el_proj_dat.end(), eprojdat_compare_refactored);	×
159	}
160
161	// set up the layout data object that provides configuration details for
162	// the hash tree
163	KDLayoutData sun_ld;
164	sun_ld.xlim[0] = System->Sun.MinXSun;	4✔
165	sun_ld.xlim[1] = System->Sun.MaxXSun;	4✔
166	sun_ld.ylim[0] = System->Sun.MinYSun;	4✔
167	sun_ld.ylim[1] = System->Sun.MaxYSun;	4✔
168	sun_ld.min_unit_dx = d_elm_max;	4✔
169	sun_ld.min_unit_dy = d_elm_max;	4✔
170
171	sun_hash.create_mesh(sun_ld);	4✔
172
173	// load stage 0 elements into the mesh
174	for (uint_fast64_t i = 0; i < System->StageList[0]->ElementList.size(); i++)	12,590✔
175	{
176	TElement *el = System->StageList[0]->ElementList.at(i).get();	12,586✔
177	sun_hash.add_object((void *)el, el->PosSunCoords[0], el->PosSunCoords[1]);	12,586✔
178	}
179
180	// calculate and associate neighbors with each zone
181	sun_hash.add_neighborhood_data();	4✔
182
183	if (AsPowerTower)	4✔
184	{
185	// Set things up for the polar coordinate tree
186	KDLayoutData rec_ld;
NEW 187	rec_ld.xlim[0] = -PI;	×
NEW 188	rec_ld.xlim[1] = PI;	×
NEW 189	rec_ld.ylim[0] = -PI / 2.;	×
NEW 190	rec_ld.ylim[1] = PI / 2.;	×
191	// use smallest element to set the minimum size
NEW 192	rec_ld.min_unit_dx = rec_ld.min_unit_dy = el_proj_dat.back().d_proj; // radians at equator	×
193
NEW 194	rec_hash.create_mesh(rec_ld);	×
195
196	// load stage 0 elements into the receiver mesh in the order of largest projection to smallest
NEW 197	for (int i = 0; i < el_proj_dat.size(); i++)	×
198	{
NEW 199	eprojdat *D = &el_proj_dat.at(i);	×
200
201	// Calculate the angular span of the element
202	double angspan[2];
NEW 203	double adjmult = 1.5;	×
NEW 204	angspan[0] = D->d_proj / cos(fabs(D->zen)) * adjmult; // azimuthal span	×
NEW 205	angspan[0] = fmin(angspan[0], 2. * PI); // limit to circumference	×
NEW 206	angspan[1] = D->d_proj / PI * adjmult; // zenithal span	×
NEW 207	rec_hash.add_object((void *)D->el_addr, D->az, D->zen, angspan);	×
208	}
209
210	// associate neighbors with each zone
NEW 211	rec_hash.add_neighborhood_data();	×
212	}
213	}	4✔
214
215	uint_fast64_t GetPTElements(	532,205✔
216	// system info
217	const bool AsPowerTower,
218
219	// Stage info
220	// const TStage *Stage
221	const tstage_ptr Stage,
222	const int i,
223
224	// Ray info
225	const bool in_multi_hit_loop,
226	const double (&PosRayStage)[3],
227	const double (&reccm_helio)[3],
228	st_hash_tree *rec_hash,
229
230	const std::vector<void *> &sunint_elements,
231
232	// Outputs
233	std::vector<void *> &reflint_elements,
234	bool &has_elements)
235	{
236	uint_fast64_t nintelements = 0;	532,205✔
237
238	if (i == 0)	532,205✔
239	{
240	if (in_multi_hit_loop)	466,326✔
241	{
242	if (AsPowerTower)	64,255✔
243	{
244	//>=Second time through - checking for first stage multiple element interactions
245
246	// get ray position in receiver polar coordinates
247	double raypvec[3];
NEW 248	for (int jj = 0; jj < 3; jj++)	×
NEW 249	raypvec[jj] = PosRayStage[jj] - reccm_helio[jj];	×
NEW 250	double raypvecmag = sqrt(raypvec[0] * raypvec[0] + raypvec[1] * raypvec[1] + raypvec[2] * raypvec[2]);	×
251	double raypol[2];
NEW 252	raypol[0] = atan2(raypvec[0], raypvec[1]);	×
NEW 253	raypol[1] = asin(raypvec[2] / raypvecmag);	×
254	// get elements in the vicinity of the ray's polar coordinates
NEW 255	reflint_elements.clear();	×
NEW 256	rec_hash->get_all_data_at_loc(reflint_elements, raypol[0], raypol[1]);	×
NEW 257	nintelements = reflint_elements.size();	×
NEW 258	has_elements = nintelements > 0;	×
259	}
260	else
261	{
262	nintelements = Stage->ElementList.size();	64,255✔
263	}
264	}
265	else
266	{
267	// First time through - checking for sun ray intersections
268	if (has_elements)	402,071✔
269	nintelements = sunint_elements.size();	334,240✔
270	else
271	nintelements = 0;	67,831✔
272	}
273	}
274	else
275	nintelements = Stage->ElementList.size();	65,879✔
276
277	return nintelements;	532,205✔
278	}
279
280	} // namespace SolTrace::NativeRunner

NREL / SolTrace / 20283049329

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous