Skip to content

Commit f20c696

Browse files
willendwillend
authored
Merge pull request #1992 from mccode-dev/tobias_reso
Reso: Updates
2 parents 41cda38 + 1ccb3b7 commit f20c696

5 files changed

Lines changed: 492 additions & 320 deletions

File tree

‎mcstas-comps/monitors/Res_monitor.comp‎

Lines changed: 97 additions & 65 deletions
Original file line numberDiff line numberDiff line change
@@ -60,7 +60,7 @@
6060
*
6161
* DEFS: [struct] structure containing Monitor_nD Defines
6262
* Vars: [struct] structure containing Monitor_nD variables
63-
* buffer_index: [long] number of recorded events
63+
* num_events: [long] number of recorded events
6464
*
6565
* %E
6666
*******************************************************************************/
@@ -70,7 +70,8 @@ DEFINE COMPONENT Res_monitor
7070

7171
SETTING PARAMETERS (string res_sample_comp,
7272
string filename=0, string options=0, xwidth=.1, yheight=.1, zdepth=0, radius=0,
73-
xmin=0, xmax=0, ymin=0, ymax=0, zmin=0, zmax=0, bufsize=0, int restore_neutron=0, int live_calc=1)
73+
xmin=0, xmax=0, ymin=0, ymax=0, zmin=0, zmax=0, bufsize=0, int restore_neutron=0,
74+
int live_calc=1)
7475

7576
/* these are protected C variables */
7677

@@ -86,10 +87,11 @@ SHARE %{
8687
DECLARE%{
8788
MonitornD_Defines_type DEFS;
8889
MonitornD_Variables_type Vars;
89-
long buffer_index;
9090

91-
tl2_list_type *reso_events;
92-
tl2_list_type *reso_probabilities;
91+
unsigned long num_events;
92+
double* reso_events;
93+
double* reso_probabilities;
94+
9395
char res_pi_var[20];
9496
char res_ki_x_var[20];
9597
char res_ki_y_var[20];
@@ -105,17 +107,17 @@ DECLARE%{
105107

106108
INITIALIZE %{
107109

108-
// Use instance name for monitor output if no input was given
109-
if (!strcmp(filename,"\0")) sprintf(filename,"%s",NAME_CURRENT_COMP);
110-
111-
buffer_index = 0;
110+
/* Use instance name for monitor output if no input was given */
111+
if (!strcmp(filename, "\0"))
112+
sprintf(filename, "%s", NAME_CURRENT_COMP);
112113

114+
num_events = 0;
113115
reso_events = 0;
114116
reso_probabilities = 0;
115117

116118
if(live_calc) {
117-
reso_events = tl2_lst_create(0);
118-
reso_probabilities = tl2_lst_create(0);
119+
reso_events = malloc(mcget_ncount()*4*sizeof(double));
120+
reso_probabilities = malloc(mcget_ncount()*sizeof(double));
119121
}
120122

121123
int i;
@@ -130,16 +132,17 @@ INITIALIZE %{
130132
if (strstr(tmp, "list"))
131133
exit(fprintf(stderr, "Res_monitor: %s: Error: Only use geometry keywords (remove list from 'option').\n", NAME_CURRENT_COMP));
132134
if (!bufsize)
133-
sprintf(Vars.option, "%s borders list all, ud1 ud2 ud3 ud4 ud5 ud6 ud7 ud8 ud9 ud10", tmp);
135+
sprintf(Vars.option, "%s borders list all, ud1 ud2 ud3 ud4 ud5 ud6 ud7 ud8 ud9 ud10 n", tmp);
134136
else
135-
sprintf(Vars.option, "%s borders list=%f, ud1 ud2 ud3 ud4 ud5 ud6 ud7 ud8 ud9 ud10", tmp, bufsize);
137+
sprintf(Vars.option, "%s borders list=%f, ud1 ud2 ud3 ud4 ud5 ud6 ud7 ud8 ud9 ud10 n", tmp, bufsize);
136138

137-
if (radius) xwidth = 2*radius;
139+
if (radius)
140+
xwidth = 2*radius;
138141

139142
Monitor_nD_Init(&DEFS, &Vars, xwidth, yheight, zdepth, xmin, xmax, ymin, ymax, zmin, zmax, 0);
140143
Vars.Coord_Type[0] = DEFS.COORD_USERDOUBLE0; /* otherwise p is always the first variable */
141144

142-
if (Vars.Coord_Number != 10)
145+
if (Vars.Coord_Number != 11)
143146
exit(fprintf(stderr,"Res_monitor: %s: Error: Invalid number of variables to monitor (%li).\n", NAME_CURRENT_COMP, Vars.Coord_Number+1));
144147

145148
/* set the labels */
@@ -162,23 +165,27 @@ INITIALIZE %{
162165
/* Initialize uservar strings */
163166
int *index_ptr=COMP_GETPAR3(Res_sample, res_sample_comp, compindex);
164167
int index = *index_ptr;
165-
sprintf(res_pi_var,"res_pi_%i",index);
166-
sprintf(res_ki_x_var,"res_ki_x_%i",index);
167-
sprintf(res_ki_y_var,"res_ki_y_%i",index);
168-
sprintf(res_ki_z_var,"res_ki_z_%i",index);
169-
sprintf(res_kf_x_var,"res_kf_x_%i",index);
170-
sprintf(res_kf_y_var,"res_kf_y_%i",index);
171-
sprintf(res_kf_z_var,"res_kf_z_%i",index);
172-
sprintf(res_rx_var,"res_rx_%i",index);
173-
sprintf(res_ry_var,"res_ry_%i",index);
174-
sprintf(res_rz_var,"res_rz_%i",index);
168+
sprintf(res_pi_var, "res_pi_%i", index);
169+
sprintf(res_ki_x_var, "res_ki_x_%i", index);
170+
sprintf(res_ki_y_var, "res_ki_y_%i", index);
171+
sprintf(res_ki_z_var, "res_ki_z_%i", index);
172+
sprintf(res_kf_x_var, "res_kf_x_%i", index);
173+
sprintf(res_kf_y_var, "res_kf_y_%i", index);
174+
sprintf(res_kf_z_var, "res_kf_z_%i", index);
175+
sprintf(res_rx_var, "res_rx_%i", index);
176+
sprintf(res_ry_var, "res_ry_%i", index);
177+
sprintf(res_rz_var, "res_rz_%i", index);
175178
%}
176179

177180

178181
TRACE %{
179182
double t0 = 0;
180183
double t1 = 0;
181184
int intersect = 0;
185+
unsigned long event_idx = 0;
186+
double event_ki[3], event_kf[3];
187+
double event_pi, event_pf;
188+
double event_pos[3];
182189

183190
if (abs(Vars.Flag_Shape) == DEFS.SHAPE_SQUARE) /* square xy */
184191
{
@@ -224,49 +231,63 @@ TRACE %{
224231
}
225232

226233
/* Now fetch data from the Res_sample. */
227-
if(p && (!bufsize || buffer_index<bufsize))
234+
if(p > 0. && (!bufsize || num_events < bufsize))
228235
{
229236
/* old behaviour not supported by openacc:
230237
struct Res_sample_struct *s =
231238
(struct Res_sample_struct *)(COMP_GETPAR3(Res_sample, res_sample_comp, res_struct));*/
232239

233-
if(particle_getvar(_particle,res_pi_var,NULL))
240+
/* pi */
241+
event_pi = particle_getvar(_particle, res_pi_var, NULL);
242+
if(event_pi > 0.)
234243
{
244+
/* pf */
245+
event_pf = p / event_pi;
246+
235247
/* ki */
236-
Vars.UserDoubles[0] = particle_getvar(_particle,res_ki_x_var,NULL);
237-
Vars.UserDoubles[1] = particle_getvar(_particle,res_ki_y_var,NULL);
238-
Vars.UserDoubles[2] = particle_getvar(_particle,res_ki_z_var,NULL);
248+
event_ki[0] = particle_getvar(_particle, res_ki_x_var, NULL);
249+
event_ki[1] = particle_getvar(_particle, res_ki_y_var, NULL);
250+
event_ki[2] = particle_getvar(_particle, res_ki_z_var, NULL);
239251

240252
/* kf */
241-
Vars.UserDoubles[3] = particle_getvar(_particle,res_kf_x_var,NULL);
242-
Vars.UserDoubles[4] = particle_getvar(_particle,res_kf_y_var,NULL);
243-
Vars.UserDoubles[5] = particle_getvar(_particle,res_kf_z_var,NULL);
253+
event_kf[0] = particle_getvar(_particle, res_kf_x_var, NULL);
254+
event_kf[1] = particle_getvar(_particle, res_kf_y_var, NULL);
255+
event_kf[2] = particle_getvar(_particle, res_kf_z_var, NULL);
244256

245257
/* pos */
246-
Vars.UserDoubles[6] = particle_getvar(_particle,res_rx_var,NULL);
247-
Vars.UserDoubles[7] = particle_getvar(_particle,res_ry_var,NULL);
248-
Vars.UserDoubles[8] = particle_getvar(_particle,res_rz_var,NULL);
258+
event_pos[0] = particle_getvar(_particle, res_rx_var, NULL);
259+
event_pos[1] = particle_getvar(_particle, res_ry_var, NULL);
260+
event_pos[2] = particle_getvar(_particle, res_rz_var, NULL);
249261

250-
/* pi and pf */
251-
Vars.UserDoubles[9] = particle_getvar(_particle,res_pi_var,NULL);
252-
Vars.UserDoubles[10] = p/particle_getvar(_particle,res_pi_var,NULL);
262+
#pragma acc atomic capture
263+
{
264+
event_idx = num_events++;
265+
}
266+
267+
/* variables for Monitor_nD */
268+
Vars.UserDoubles[0] = event_ki[0];
269+
Vars.UserDoubles[1] = event_ki[1];
270+
Vars.UserDoubles[2] = event_ki[2];
271+
Vars.UserDoubles[3] = event_kf[0];
272+
Vars.UserDoubles[4] = event_kf[1];
273+
Vars.UserDoubles[5] = event_kf[2];
274+
Vars.UserDoubles[6] = event_pos[0];
275+
Vars.UserDoubles[7] = event_pos[1];
276+
Vars.UserDoubles[8] = event_pos[2];
277+
Vars.UserDoubles[9] = event_pi;
278+
Vars.UserDoubles[10] = event_pf;
253279

254280
Monitor_nD_Trace(&DEFS, &Vars, _particle);
255281

256282
/* live calculation */
257283
if(live_calc) {
258-
double *reso_event = (double*)malloc(4 * sizeof(double));
259-
reso_event[0] = Vars.UserDoubles[0] - Vars.UserDoubles[3];
260-
reso_event[1] = Vars.UserDoubles[1] - Vars.UserDoubles[4];
261-
reso_event[2] = Vars.UserDoubles[2] - Vars.UserDoubles[5];
262-
reso_event[3] = tl2_k_to_E(Vars.UserDoubles[0], Vars.UserDoubles[1], Vars.UserDoubles[2],
263-
Vars.UserDoubles[3], Vars.UserDoubles[4], Vars.UserDoubles[5]);
264-
265-
double *reso_prob = (double*)malloc(sizeof(double));
266-
*reso_prob = p;
267-
268-
tl2_lst_append(reso_events, reso_event);
269-
tl2_lst_append(reso_probabilities, reso_prob);
284+
reso_events[event_idx*4 + 0] = event_ki[0] - event_kf[0];
285+
reso_events[event_idx*4 + 1] = event_ki[1] - event_kf[1];
286+
reso_events[event_idx*4 + 2] = event_ki[2] - event_kf[2];
287+
reso_events[event_idx*4 + 3] = tl2_k_to_E(
288+
event_ki[0], event_ki[1], event_ki[2],
289+
event_kf[0], event_kf[1], event_kf[2]);
290+
reso_probabilities[event_idx] = event_pi * event_pf;
270291
}
271292
}
272293

@@ -281,22 +302,31 @@ TRACE %{
281302

282303

283304
SAVE %{
305+
int dbg_save_events = 0;
306+
284307
/* save results, but do not free pointers */
285308
Monitor_nD_Save(&DEFS, &Vars);
286309

287310
/* live calculation */
288-
if(live_calc) {
289-
double cov[4*4], reso[4*4];
290-
if(tl2_reso(reso_events->next, reso_probabilities->next, cov, reso)) {
291-
printf("Resolution matrix:"
292-
"\n\t%12.4e %12.4e %12.4e %12.4e"
293-
"\n\t%12.4e %12.4e %12.4e %12.4e"
294-
"\n\t%12.4e %12.4e %12.4e %12.4e"
295-
"\n\t%12.4e %12.4e %12.4e %12.4e\n",
296-
reso[0], reso[1], reso[2], reso[3],
297-
reso[4], reso[5], reso[6], reso[7],
298-
reso[8], reso[9], reso[10], reso[11],
299-
reso[12], reso[13], reso[14], reso[15]);
311+
if(live_calc && num_events > 0) {
312+
if(dbg_save_events) {
313+
/* save individual neutron events */
314+
#ifndef USE_MPI
315+
const char* event_filename = "reso_events.dat";
316+
#else
317+
char event_filename[256];
318+
sprintf(event_filename, "reso_events_%d.dat", mpi_node_rank);
319+
#endif
320+
321+
tl2_save_events(reso_events, reso_probabilities, event_filename, num_events);
322+
}
323+
324+
double cov[4*4], res[4*4];
325+
if(tl2_reso(reso_events, reso_probabilities, cov, res, num_events)) {
326+
printf("Resolution calculation used %d neutron events.\n", num_events);
327+
tl2_print_mat(cov, "Covariance matrix", 4, 4);
328+
tl2_print_mat(res, "Resolution matrix", 4, 4);
329+
printf("Please run \"mcresplot %s\" for a full analysis.\n", filename);
300330
}
301331
else {
302332
printf("Error: Resolution matrix could not be calculated.");
@@ -309,8 +339,10 @@ FINALLY %{
309339
/* free pointers */
310340
Monitor_nD_Finally(&DEFS, &Vars);
311341

312-
tl2_lst_free(reso_events);
313-
tl2_lst_free(reso_probabilities);
342+
if(live_calc) {
343+
free(reso_events);
344+
free(reso_probabilities);
345+
}
314346
%}
315347

316348

0 commit comments

Comments
 (0)