From 74dc408029cdd492cde5ec8aa19df71a606ed01d Mon Sep 17 00:00:00 2001 From: Peter Willendrup Date: Fri, 20 Dec 2024 09:46:51 +0100 Subject: [PATCH 1/3] Minor updates to Union comps/lib from Mads, work for NMX at ESS --- mcstas-comps/share/union-lib.c | 1 + .../union/Union_abs_logger_event.comp | 131 +++- mcstas-comps/union/Union_abs_logger_nD.comp | 564 ++++++++++++++++++ mcstas-comps/union/Union_master.comp | 15 + 4 files changed, 701 insertions(+), 10 deletions(-) create mode 100644 mcstas-comps/union/Union_abs_logger_nD.comp diff --git a/mcstas-comps/share/union-lib.c b/mcstas-comps/share/union-lib.c index 5ca190853..7a8b7edae 100755 --- a/mcstas-comps/share/union-lib.c +++ b/mcstas-comps/share/union-lib.c @@ -219,6 +219,7 @@ union abs_logger_data_union{ struct a_1D_time_to_lambda_abs_storage_struct *p_1D_time_to_lambda_abs_storage; struct a_event_abs_storage_struct *p_event_abs_storage; struct a_1D_event_abs_storage_struct *p_1D_event_abs_storage; + struct a_nD_abs_storage_struct *p_nD_abs_storage; // Additional logger storage structs to be addedd }; diff --git a/mcstas-comps/union/Union_abs_logger_event.comp b/mcstas-comps/union/Union_abs_logger_event.comp index d3e5fb6ac..6452ae515 100644 --- a/mcstas-comps/union/Union_abs_logger_event.comp +++ b/mcstas-comps/union/Union_abs_logger_event.comp @@ -74,6 +74,7 @@ DEFINE COMPONENT Union_abs_logger_event SETTING PARAMETERS(string target_geometry="NULL", string filename="NULL", + xwidth=0, xbins=0, yheight=0, ybins=0, zdepth=0, zbins=0, order_total=-1, order_volume=-1, logger_conditional_extend_index=-1, string init="init") @@ -120,19 +121,40 @@ struct a_event_abs_storage_struct { int order; int order_in_this_volume; int order_process_in_this_volume; + + double stored_xwidth; + double stored_yheight; + double stored_zdepth; + + int stored_xbins; + int stored_ybins; + int stored_zbins; Coords position; Rotation rotation; Rotation t_rotation; }; +double round_to_nearest_bin(double value, double w, int N) { + // Calculate the width of each bin + double bin_width = w / N; + + // Calculate the index of the nearest bin + int bin_index = round((value + w / 2) / bin_width); + + // Calculate the center of the nearest bin + double bin_center = (bin_index * bin_width) - (w / 2) + (bin_width / 2); + + return bin_center; +} + // record_to_temp // Would be nice if x y z, k_new and k_old were all coords void record_to_temp_abs_event(Coords *position, double *k, double p, double time, int scattered_in_this_volume, int total_number_of_scattering_events, struct abs_logger_struct *abs_logger, struct abs_logger_with_data_struct *abs_logger_with_data_array) { struct a_event_abs_storage_struct *storage; storage = abs_logger->data_union.p_event_abs_storage; - + int add_point = 1; if (storage->order != -1) { @@ -141,7 +163,7 @@ void record_to_temp_abs_event(Coords *position, double *k, double p, double time else add_point = 0; } - + if (storage->order_in_this_volume != -1) { if (storage->order_in_this_volume == scattered_in_this_volume) add_point = 1; @@ -149,11 +171,51 @@ void record_to_temp_abs_event(Coords *position, double *k, double p, double time add_point = 0; } + if (storage->stored_xwidth != 0) { + if (position->x < 0.5*storage->stored_xwidth && position->x > -0.5*storage->stored_xwidth) + add_point = 1; + else + add_point = 0; + } + + if (storage->stored_yheight != 0) { + if (position->y < 0.5*storage->stored_yheight && position->y > -0.5*storage->stored_yheight) + add_point = 1; + else + add_point = 0; + } + + if (storage->stored_zdepth != 0) { + if (position->z < 0.5*storage->stored_zdepth && position->z > -0.5*storage->stored_xwidth) + add_point = 1; + else + add_point = 0; + } + if (add_point == 1) { int i; double given_x_pos, given_y_pos, given_z_pos; + double used_x_pos, used_y_pos, used_z_pos; coords_get(*position, &given_x_pos, &given_y_pos, &given_z_pos); + + if (storage->stored_xbins != 0) { + used_x_pos = round_to_nearest_bin(given_x_pos, storage->stored_xwidth, storage->stored_xbins); + } else { + used_x_pos = given_x_pos; + } + + if (storage->stored_ybins != 0) { + used_y_pos = round_to_nearest_bin(given_y_pos, storage->stored_yheight, storage->stored_ybins); + } else { + used_y_pos = given_y_pos; + } + + if (storage->stored_zbins != 0) { + used_z_pos = round_to_nearest_bin(given_z_pos, storage->stored_zdepth, storage->stored_zbins); + } else { + used_z_pos = given_z_pos; + } double given_x_vel, given_y_vel, given_z_vel; given_x_vel = k[0]*K2V; @@ -162,9 +224,9 @@ void record_to_temp_abs_event(Coords *position, double *k, double p, double time if (storage->temp_abs_event_data.num_elements < storage->temp_abs_event_data.allocated_elements) { - storage->temp_abs_event_data.elements[storage->temp_abs_event_data.num_elements].x_pos = given_x_pos; - storage->temp_abs_event_data.elements[storage->temp_abs_event_data.num_elements].y_pos = given_y_pos; - storage->temp_abs_event_data.elements[storage->temp_abs_event_data.num_elements].z_pos = given_z_pos; + storage->temp_abs_event_data.elements[storage->temp_abs_event_data.num_elements].x_pos = used_x_pos; + storage->temp_abs_event_data.elements[storage->temp_abs_event_data.num_elements].y_pos = used_y_pos; + storage->temp_abs_event_data.elements[storage->temp_abs_event_data.num_elements].z_pos = used_z_pos; storage->temp_abs_event_data.elements[storage->temp_abs_event_data.num_elements].x_vel = given_x_vel; storage->temp_abs_event_data.elements[storage->temp_abs_event_data.num_elements].y_vel = given_y_vel; storage->temp_abs_event_data.elements[storage->temp_abs_event_data.num_elements].z_vel = given_z_vel; @@ -270,12 +332,52 @@ void record_to_perm_abs_event(Coords *position, double *k, double p, double time add_point = 0; } + if (storage->stored_xwidth != 0) { + if (position->x < 0.5*storage->stored_xwidth && position->x > -0.5*storage->stored_xwidth) + add_point = 1; + else + add_point = 0; + } + + if (storage->stored_yheight != 0) { + if (position->y < 0.5*storage->stored_yheight && position->y > -0.5*storage->stored_yheight) + add_point = 1; + else + add_point = 0; + } + + if (storage->stored_zdepth != 0) { + if (position->z < 0.5*storage->stored_zdepth && position->z > -0.5*storage->stored_xwidth) + add_point = 1; + else + add_point = 0; + } + if (add_point == 1) { //printf("storage was set \n"); double given_x_pos, given_y_pos, given_z_pos; + double used_x_pos, used_y_pos, used_z_pos; coords_get(*position, &given_x_pos, &given_y_pos, &given_z_pos); - + + if (storage->stored_xbins != 0) { + used_x_pos = round_to_nearest_bin(given_x_pos, storage->stored_xwidth, storage->stored_xbins); + } else { + used_x_pos = given_x_pos; + } + + if (storage->stored_ybins != 0) { + used_y_pos = round_to_nearest_bin(given_y_pos, storage->stored_yheight, storage->stored_ybins); + } else { + used_y_pos = given_y_pos; + } + + if (storage->stored_zbins != 0) { + used_z_pos = round_to_nearest_bin(given_z_pos, storage->stored_zdepth, storage->stored_zbins); + } else { + used_z_pos = given_z_pos; + } + double given_x_vel, given_y_vel, given_z_vel; given_x_vel = k[0]*K2V; given_y_vel = k[1]*K2V; @@ -283,9 +385,9 @@ void record_to_perm_abs_event(Coords *position, double *k, double p, double time _class_particle _localparticle; - _localparticle.x = given_x_pos; - _localparticle.y = given_y_pos; - _localparticle.z = given_z_pos; + _localparticle.x = used_x_pos; + _localparticle.y = used_y_pos; + _localparticle.z = used_z_pos; _localparticle.vx = given_x_vel; _localparticle.vy = given_y_vel; @@ -478,8 +580,17 @@ INITIALIZE this_abs_storage.order = order_total; this_abs_storage.order_in_this_volume = order_volume; + this_abs_storage.stored_xwidth = xwidth; + this_abs_storage.stored_xbins = (int) xbins; + + this_abs_storage.stored_yheight = yheight; + this_abs_storage.stored_ybins = (int) ybins; + + this_abs_storage.stored_zdepth = zdepth; + this_abs_storage.stored_zbins = (int) zbins; + this_abs_storage.temp_abs_event_data.num_elements=0; - + this_abs_storage.temp_abs_event_data.allocated_elements = 10; this_abs_storage.temp_abs_event_data.elements = malloc(this_abs_storage.temp_abs_event_data.allocated_elements*sizeof(struct temp_abs_event_data_element_struct)); diff --git a/mcstas-comps/union/Union_abs_logger_nD.comp b/mcstas-comps/union/Union_abs_logger_nD.comp new file mode 100644 index 000000000..bd013fadc --- /dev/null +++ b/mcstas-comps/union/Union_abs_logger_nD.comp @@ -0,0 +1,564 @@ +/******************************************************************************* +* +* McStas, neutron ray-tracing package +* Copyright(C) 2007 Risoe National Laboratory. +* +* %I +* Written by: Mads Bertelsen +* Date: 19.06.20 +* Version: $Revision: 0.1 $ +* Origin: ESS DMSC +* +* A logger of absorption as events with Monitor_nD options +* +* %D +* Part of the Union components, a set of components that work together and thus +* separates geometry and physics within McStas. +* The use of this component requires other components to be used. +* +* 1) One specifies a number of processes using process components +* 2) These are gathered into material definitions using Union_make_material +* 3) Geometries are placed using Union_box/cylinder/sphere, assigned a material +* 4) Logger and conditional components can be placed which will record what happens +* 5) A Union_master component placed after all of the above +* +* Only in step 5 will any simulation happen, and per default all geometries +* defined before this master, but after the previous will be simulated here. +* +* There is a dedicated manual available for the Union_components +* +* This component is an absorption logger, and thus placed in point 4) above. +* +* An absorption logger will log something for each absorption event happening +* in the geometry or geometries on which it is attached. These are specified +* in the target_geometry string. By leaving it blank, all geometries are +* logged, even the ones not defined at this point in the instrument file. +* Multiple geometries are specified as a comma separated list. +* +* This absorption logger stores absorption as events, with position, velocity, +* time and weight. The Monitor_nD libraries are used to write the event files. +* This version is a close copy of Monitor_nD, having the same interface, though +* the user must be aware that no propagation happens for rays to hit the +* detector pixels, instead it uses the position of absorbed. Use the previous +* keyword to tell Monitor_nD that this is going on. +* +* This absorption logger needs to be placed in space, the position and velocity +* is recorded in the coordinate system of the logger component. +* +* It is possible to attach one or more conditional components to this absorption +* logger. Such a conditional component would impose a condition on the state of +* the neutron after the Union_master component that executes the simulation, +* and the absorption logger will only record the event if this condition is true. +* +* To use the logger_conditional_extend function, set it to some integer value n +* and make and extend section to the master component that runs the geometry. +* In this extend function, logger_conditional_extend[n] is 1 if the conditional +* stack evaluated to true, 0 if not. This way one can check what rays is logged +* using regular McStas monitors. Only works if a conditional is applied to this +* logger. +* +* %P +* INPUT PARAMETERS: +* target_geometry: [string] Comma separated list of geometry names that will be logged, leave empty for all volumes (even not defined yet) +* filename: [string] Filename of produced data file +* order_total: [1] Only log rays that have scattered n times, -1 for all orders +* order_volume: [1] Only log rays that have scattered n times in the same geometry, -1 for all orders +* logger_conditional_extend_index: [1] If a conditional is used with this logger, the result of each conditional calculation can be made available in extend as a array called "logger_conditional_extend", and one would then access logger_conditional_extend[n] if logger_conditional_extend_index is set to n +* init: [string] Name of Union_init component (typically "init", default) +* +* OUTPUT PARAMETERS: +* +* GLOBAL PARAMETERS: +* +* %L +* +* %E +******************************************************************************/ + +DEFINE COMPONENT Union_abs_logger_nD +DEFINITION PARAMETERS () +SETTING PARAMETERS(string target_geometry="NULL", + order_total=-1, order_volume=-1, logger_conditional_extend_index=-1, string init="init", + string user1="", string user2="", string user3="", + xwidth=0, yheight=0, zdepth=0, + xmin=0, xmax=0, ymin=0, ymax=0, zmin=0, zmax=0, + int bins=0, min=-1e40, max=1e40, int restore_neutron=0, radius=0, + string options="NULL", string filename="NULL",string geometry="NULL", int nowritefile=0, + string username1="NULL", string username2="NULL", string username3="NULL") +OUTPUT PARAMETERS () + +/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */ + +SHARE +%{ +#ifndef Union +#error "The Union_init component must be included before this Union_abs_logger_nD component" +#endif + +%include "monitor_nd-lib" // Internally avoids double import + +struct temp_abs_nD_data_element_struct { + double x_pos; + double y_pos; + double z_pos; + + double x_vel; + double y_vel; + double z_vel; + + double x_pol; + double y_pol; + double z_pol; + + double time; + double weight; +}; + +struct temp_abs_nD_data_struct { + int num_elements; + int allocated_elements; + struct temp_abs_nD_data_element_struct *elements; +}; + +struct a_nD_abs_storage_struct { + MonitornD_Defines_type *p_DEFS; + MonitornD_Variables_type *p_Vars; + + struct temp_abs_nD_data_struct temp_abs_nD_data; + //some type + int order; + int order_in_this_volume; + int order_process_in_this_volume; + + Coords position; + Rotation rotation; + Rotation t_rotation; +}; + +// record_to_temp +// Would be nice if x y z, k_new and k_old were all coords +void record_to_temp_abs_nD(Coords *position, double *k, double p, double time, int scattered_in_this_volume, int total_number_of_scattering_events, struct abs_logger_struct *abs_logger, struct abs_logger_with_data_struct *abs_logger_with_data_array) { + + struct a_nD_abs_storage_struct *storage; + storage = abs_logger->data_union.p_nD_abs_storage; + + int add_point = 1; + + if (storage->order != -1) { + if (storage->order == total_number_of_scattering_events) + add_point = 1; + else + add_point = 0; + } + + if (storage->order_in_this_volume != -1) { + if (storage->order_in_this_volume == scattered_in_this_volume) + add_point = 1; + else + add_point = 0; + } + + if (add_point == 1) { + + int i; + double given_x_pos, given_y_pos, given_z_pos; + coords_get(*position, &given_x_pos, &given_y_pos, &given_z_pos); + + double given_x_vel, given_y_vel, given_z_vel; + given_x_vel = k[0]*K2V; + given_y_vel = k[1]*K2V; + given_z_vel = k[2]*K2V; + + if (storage->temp_abs_nD_data.num_elements < storage->temp_abs_nD_data.allocated_elements) { + + storage->temp_abs_nD_data.elements[storage->temp_abs_nD_data.num_elements].x_pos = given_x_pos; + storage->temp_abs_nD_data.elements[storage->temp_abs_nD_data.num_elements].y_pos = given_y_pos; + storage->temp_abs_nD_data.elements[storage->temp_abs_nD_data.num_elements].z_pos = given_z_pos; + storage->temp_abs_nD_data.elements[storage->temp_abs_nD_data.num_elements].x_vel = given_x_vel; + storage->temp_abs_nD_data.elements[storage->temp_abs_nD_data.num_elements].y_vel = given_y_vel; + storage->temp_abs_nD_data.elements[storage->temp_abs_nD_data.num_elements].z_vel = given_z_vel; + storage->temp_abs_nD_data.elements[storage->temp_abs_nD_data.num_elements].x_pol = 0.0; + storage->temp_abs_nD_data.elements[storage->temp_abs_nD_data.num_elements].y_pol = 0.0; + storage->temp_abs_nD_data.elements[storage->temp_abs_nD_data.num_elements].z_pol = 1.0; + storage->temp_abs_nD_data.elements[storage->temp_abs_nD_data.num_elements].time = time; + storage->temp_abs_nD_data.elements[storage->temp_abs_nD_data.num_elements++].weight = p; + } else { + // No more space, need to allocate a larger buffer for this logger. Wish I had generics. + + // copy current data to temp + struct temp_abs_nD_data_struct temporary_storage; + temporary_storage.num_elements = storage->temp_abs_nD_data.num_elements; + temporary_storage.elements = malloc(temporary_storage.num_elements*sizeof(struct temp_abs_nD_data_element_struct)); + + int index; + for (index=0;indextemp_abs_nD_data.num_elements;index++) { + temporary_storage.elements[index].x_pos = storage->temp_abs_nD_data.elements[index].x_pos; + temporary_storage.elements[index].y_pos = storage->temp_abs_nD_data.elements[index].y_pos; + temporary_storage.elements[index].z_pos = storage->temp_abs_nD_data.elements[index].z_pos; + temporary_storage.elements[index].x_vel = storage->temp_abs_nD_data.elements[index].x_vel; + temporary_storage.elements[index].y_vel = storage->temp_abs_nD_data.elements[index].y_vel; + temporary_storage.elements[index].z_vel = storage->temp_abs_nD_data.elements[index].z_vel; + temporary_storage.elements[index].x_pol = storage->temp_abs_nD_data.elements[index].x_pol; + temporary_storage.elements[index].y_pol = storage->temp_abs_nD_data.elements[index].y_pol; + temporary_storage.elements[index].z_pol = storage->temp_abs_nD_data.elements[index].z_pol; + temporary_storage.elements[index].time = storage->temp_abs_nD_data.elements[index].time; + temporary_storage.elements[index].weight = storage->temp_abs_nD_data.elements[index].weight; + } + + // free current data + free(storage->temp_abs_nD_data.elements); + + // allocate larger array (10 larger) + storage->temp_abs_nD_data.allocated_elements = 10 + storage->temp_abs_nD_data.num_elements; + storage->temp_abs_nD_data.elements = malloc(storage->temp_abs_nD_data.allocated_elements*sizeof(struct temp_abs_nD_data_element_struct)); + + // copy back from temp + for (index=0;indextemp_abs_nD_data.num_elements;index++) { + storage->temp_abs_nD_data.elements[index].x_pos = temporary_storage.elements[index].x_pos; + storage->temp_abs_nD_data.elements[index].y_pos = temporary_storage.elements[index].y_pos; + storage->temp_abs_nD_data.elements[index].z_pos = temporary_storage.elements[index].z_pos; + storage->temp_abs_nD_data.elements[index].x_vel = temporary_storage.elements[index].x_vel; + storage->temp_abs_nD_data.elements[index].y_vel = temporary_storage.elements[index].y_vel; + storage->temp_abs_nD_data.elements[index].z_vel = temporary_storage.elements[index].z_vel; + storage->temp_abs_nD_data.elements[index].x_pol = temporary_storage.elements[index].x_pol; + storage->temp_abs_nD_data.elements[index].y_pol = temporary_storage.elements[index].y_pol; + storage->temp_abs_nD_data.elements[index].z_pol = temporary_storage.elements[index].z_pol; + storage->temp_abs_nD_data.elements[index].time = temporary_storage.elements[index].time; + storage->temp_abs_nD_data.elements[index].weight = temporary_storage.elements[index].weight; + } + + // free temporary data + free(temporary_storage.elements); + + // add new data point + storage->temp_abs_nD_data.elements[storage->temp_abs_nD_data.num_elements].x_pos = given_x_pos; + storage->temp_abs_nD_data.elements[storage->temp_abs_nD_data.num_elements].y_pos = given_y_pos; + storage->temp_abs_nD_data.elements[storage->temp_abs_nD_data.num_elements].z_pos = given_z_pos; + storage->temp_abs_nD_data.elements[storage->temp_abs_nD_data.num_elements].x_vel = given_x_vel; + storage->temp_abs_nD_data.elements[storage->temp_abs_nD_data.num_elements].y_vel = given_y_vel; + storage->temp_abs_nD_data.elements[storage->temp_abs_nD_data.num_elements].z_vel = given_z_vel; + storage->temp_abs_nD_data.elements[storage->temp_abs_nD_data.num_elements].x_pol = 0.0; + storage->temp_abs_nD_data.elements[storage->temp_abs_nD_data.num_elements].y_pol = 0.0; + storage->temp_abs_nD_data.elements[storage->temp_abs_nD_data.num_elements].z_pol = 1.0; + storage->temp_abs_nD_data.elements[storage->temp_abs_nD_data.num_elements].time = time; + storage->temp_abs_nD_data.elements[storage->temp_abs_nD_data.num_elements++].weight = p; + } + + // If this is the first time this ray is being recorded in this logger, add it to the list of loggers that write to temp and may get it moved to perm + if (storage->temp_abs_nD_data.num_elements == 1) + add_to_abs_logger_with_data(abs_logger_with_data_array,abs_logger); + + } +} + +// clear_temp +void clear_temp_abs_nD(union abs_logger_data_union *data_union) { + data_union->p_nD_abs_storage->temp_abs_nD_data.num_elements = 0; +} + +// record_to_perm +void record_to_perm_abs_nD(Coords *position, double *k, double p, double time, int scattered_in_this_volume, int total_number_of_scattering_events, struct abs_logger_struct *abs_logger, struct abs_logger_with_data_struct *abs_logger_with_data_array) { + + //printf("In record to permanent \n"); + struct a_nD_abs_storage_struct *storage; + storage = abs_logger->data_union.p_nD_abs_storage; + + int add_point = 1; + + if (storage->order != -1) { + if (storage->order == total_number_of_scattering_events) + add_point = 1; + else + add_point = 0; + } + + if (storage->order_in_this_volume != -1) { + if (storage->order_in_this_volume == scattered_in_this_volume) + add_point = 1; + else + add_point = 0; + } + + if (add_point == 1) { + //printf("storage was set \n"); + + double given_x_pos, given_y_pos, given_z_pos; + coords_get(*position, &given_x_pos, &given_y_pos, &given_z_pos); + + double given_x_vel, given_y_vel, given_z_vel; + given_x_vel = k[0]*K2V; + given_y_vel = k[1]*K2V; + given_z_vel = k[2]*K2V; + + _class_particle _localparticle; + + _localparticle.x = given_x_pos; + _localparticle.y = given_y_pos; + _localparticle.z = given_z_pos; + + _localparticle.vx = given_x_vel; + _localparticle.vy = given_y_vel; + _localparticle.vz = given_z_vel; + + _localparticle.sx = 0.0; + _localparticle.sy = 0.0; + _localparticle.sz = 1.0; + + _localparticle.p = p; + _localparticle.t = time; + + int pp; + pp = Monitor_nD_Trace(storage->p_DEFS, storage->p_Vars, &_localparticle); + + } + +} + +// write_temp_to_perm +void write_temp_to_perm_abs_nD(union abs_logger_data_union *data_union) { + + struct a_nD_abs_storage_struct *storage; + storage = data_union->p_nD_abs_storage; + + int index; + // Add all data points to the historgram, they are saved as index / weight combinations + for (index=0;indextemp_abs_nD_data.num_elements;index++) { + + _class_particle _localparticle; + + _localparticle.x = storage->temp_abs_nD_data.elements[index].x_pos; + _localparticle.y = storage->temp_abs_nD_data.elements[index].y_pos; + _localparticle.z = storage->temp_abs_nD_data.elements[index].z_pos; + + _localparticle.vx = storage->temp_abs_nD_data.elements[index].x_vel; + _localparticle.vy = storage->temp_abs_nD_data.elements[index].y_vel; + _localparticle.vz = storage->temp_abs_nD_data.elements[index].z_vel; + + _localparticle.sx = 0.0; + _localparticle.sy = 0.0; + _localparticle.sz = 1.0; + + _localparticle.p = storage->temp_abs_nD_data.elements[index].weight; + _localparticle.t = storage->temp_abs_nD_data.elements[index].time; + + int pp; + pp = Monitor_nD_Trace(storage->p_DEFS, storage->p_Vars, &_localparticle); + + } + clear_temp_abs_nD(data_union); +} + +void write_temp_to_perm_final_p_abs_nD(union abs_logger_data_union *data_union, double final_weight) { + + struct a_nD_abs_storage_struct *storage; + storage = data_union->p_nD_abs_storage; + + int index; + // Add all data points to the historgram, they are saved as index / weight combinations + for (index=0;indextemp_abs_nD_data.num_elements;index++) { + + _class_particle _localparticle; + + _localparticle.x = storage->temp_abs_nD_data.elements[index].x_pos; + _localparticle.y = storage->temp_abs_nD_data.elements[index].y_pos; + _localparticle.z = storage->temp_abs_nD_data.elements[index].z_pos; + + _localparticle.vx = storage->temp_abs_nD_data.elements[index].x_vel; + _localparticle.vy = storage->temp_abs_nD_data.elements[index].y_vel; + _localparticle.vz = storage->temp_abs_nD_data.elements[index].z_vel; + + _localparticle.sx = 0.0; + _localparticle.sy = 0.0; + _localparticle.sz = 1.0; + + _localparticle.p = final_weight; + _localparticle.t = storage->temp_abs_nD_data.elements[index].time; + + int pp; + pp = Monitor_nD_Trace(storage->p_DEFS, storage->p_Vars, &_localparticle); + } + clear_temp_abs_nD(data_union); +} + +// Only need to define linking function for loggers once. +#ifndef UNION_ABS_LOGGER +#define UNION_ABS_LOGGER Dummy +// Linking function for loggers, finds the indicies of the specified geometries on the global_geometry_list +void manual_linking_function_abs_logger_volumes(char *input_string, struct pointer_to_global_geometry_list *global_geometry_list, struct pointer_to_1d_int_list *accepted_volumes, char *component_name) { + // Need to check a input_string of text for an occurance of name. If it is in the inputstring, yes return 1, otherwise 0. + char *token; + int loop_index; + char local_string[512]; + + strcpy(local_string,input_string); + // get the first token + token = strtok(local_string,","); + + // walk through other tokens + while( token != NULL ) + { + //printf( " %s\n", token ); + for (loop_index=0;loop_indexnum_elements;loop_index++) { + if (strcmp(token,global_geometry_list->elements[loop_index].name) == 0) { + add_element_to_int_list(accepted_volumes,loop_index); + break; + } + + if (loop_index == global_geometry_list->num_elements - 1) { + // All possible geometry names have been looked through, and the break was not executed. + // Alert the user to this problem by showing the geometry name that was not found and the currently available geometires + printf("\n"); + printf("ERROR: The target_geometry string \"%s\" in Union logger component \"%s\" had an entry that did not match a specified geometry. \n",input_string,component_name); + printf(" The unrecoignized geometry name was: \"%s\" \n",token); + printf(" The geometries available at this point (need to be defined before the logger): \n"); + for (loop_index=0;loop_indexnum_elements;loop_index++) + printf(" %s\n",global_geometry_list->elements[loop_index].name); + exit(1); + } + } + + // Updates the token + token = strtok(NULL,","); + } +} + +#endif + +%} + +DECLARE +%{ +int loop_index; +int found_process; +int specified_processes; +char local_string[256]; + +// Reused for logger +struct pointer_to_1d_int_list accepted_processes; + +struct global_abs_logger_element_struct abs_logger_list_element; + +struct pointer_to_1d_int_list accepted_volumes; + +struct abs_logger_struct this_abs_logger; +struct a_nD_abs_storage_struct this_abs_storage; + +struct abs_loggers_struct *abs_loggers_on_target_volume; +struct Volume_struct *target_volume; + +MonitornD_Defines_type DEFS; +MonitornD_Variables_type Vars; +MCDETECTOR detector; +off_struct offdata; +%} + +INITIALIZE COPY Monitor_nD EXTEND +%{ + + printf("Completed Monitor_nD initialization, continuing Union_abs_logger_nD initialization."); + + // input sanitation, provide error if commands used in options that are not supported in Union context + // capture, 3He_pressure, slit, absorb + if (strstr(Vars.option, "capture") || strstr(Vars.option, "3He_pressure") || strstr(Vars.option, "slit") || strstr(Vars.option, "absorb")) + exit(printf("Abs_logger_nD: %s Union version does not support following options. Aborting (capture, 3He_pressure, slit, absorb).\n", NAME_CURRENT_COMP)); + + // May need to manipulate options string for: min pixel ID (if auto_pixel_ID) and previous + + this_abs_storage.p_DEFS = &DEFS; // Grab DEFS prepared by Monitor_nD initialization + this_abs_storage.p_Vars = &Vars; // Grab Vars prepared by Monitor_nD initialization + + this_abs_storage.order = order_total; + this_abs_storage.order_in_this_volume = order_volume; + + this_abs_storage.temp_abs_nD_data.num_elements=0; + + this_abs_storage.temp_abs_nD_data.allocated_elements = 10; + this_abs_storage.temp_abs_nD_data.elements = malloc(this_abs_storage.temp_abs_nD_data.allocated_elements*sizeof(struct temp_abs_nD_data_element_struct)); + + if (_getcomp_index(init) < 0) { + fprintf(stderr,"Union_abs_logger_nD:%s: Error identifying Union_init component, %s is not a known component name.\n", + NAME_CURRENT_COMP, init); + exit(-1); + } + + struct global_positions_to_transform_list_struct *global_positions_to_transform_list = COMP_GETPAR3(Union_init, init, global_positions_to_transform_list); + struct global_rotations_to_transform_list_struct *global_rotations_to_transform_list = COMP_GETPAR3(Union_init, init, global_rotations_to_transform_list); + // Test position and rotation stored in a data storage, and pointers assigned to transform lists + this_abs_logger.position = POS_A_CURRENT_COMP; + add_position_pointer_to_list(global_positions_to_transform_list, &this_abs_logger.position); + + rot_copy(this_abs_logger.rotation,ROT_A_CURRENT_COMP); + add_rotation_pointer_to_list(global_rotations_to_transform_list, &this_abs_logger.rotation); + + rot_transpose(ROT_A_CURRENT_COMP,this_abs_logger.t_rotation); + add_rotation_pointer_to_list(global_rotations_to_transform_list, &this_abs_logger.t_rotation); + + + // Book keeping + this_abs_logger.abs_logger_extend_index = logger_conditional_extend_index; + this_abs_logger.function_pointers.active_record_function = &record_to_perm_abs_nD; // Assume no conditional + this_abs_logger.function_pointers.inactive_record_function = &record_to_temp_abs_nD; // If an assume is present, these two pointers are switched + // Temp to perm functions, and standard identifier + //this_abs_logger.function_pointers.select_t_to_p = 1; // 1: temp_to_perm, 2: temp_to_perm_final_p // Not relevant for abs + this_abs_logger.function_pointers.temp_to_perm = &write_temp_to_perm_abs_nD; + this_abs_logger.function_pointers.temp_to_perm_final_p = &write_temp_to_perm_final_p_abs_nD; + this_abs_logger.function_pointers.clear_temp = &clear_temp_abs_nD; + + // Initializing for conditional + this_abs_logger.conditional_list.num_elements = 0; + + this_abs_logger.data_union.p_nD_abs_storage = &this_abs_storage; + sprintf(this_abs_logger.name,"%s",NAME_CURRENT_COMP); + sprintf(abs_logger_list_element.name,"%s",NAME_CURRENT_COMP); + abs_logger_list_element.component_index = INDEX_CURRENT_COMP; + abs_logger_list_element.abs_logger = &this_abs_logger; + + struct pointer_to_global_geometry_list *global_geometry_list = COMP_GETPAR3(Union_init, init, global_geometry_list); + struct pointer_to_global_abs_logger_list *global_specific_volumes_abs_logger_list = COMP_GETPAR3(Union_init, init, global_specific_volumes_abs_logger_list); + // In order to run the logger at the right times, pointers to this logger is stored in each volume it logs, + // and additionally for each avaiable process. If a process is not logged, the pointer is simply not stored. + int process_index; + // Need to find the volumes for which the processes should have a reference to this logger + if (target_geometry && strlen(target_geometry) && strcmp(target_geometry,"NULL") && strcmp(target_geometry, "0")) { + // Certain volumes were selected, find the indicies in the global_geometry_list + manual_linking_function_abs_logger_volumes(target_geometry, global_geometry_list, &accepted_volumes, NAME_CURRENT_COMP); + // Add this logger to the global_specific_volumes_logger_list (so that conditionals can affect it) + add_element_to_abs_logger_list(global_specific_volumes_abs_logger_list, abs_logger_list_element); + + for (loop_index=0;loop_indexelements[accepted_volumes.elements[loop_index]].Volume; + // Add an element to its abs_logger list + + add_initialized_abs_logger_in_volume(&target_volume->abs_loggers); + target_volume->abs_loggers.p_abs_logger[target_volume->abs_loggers.num_elements-1] = &this_abs_logger; + + } + } else { + // Send to global_all_volumes_logger_list + // Here there is no system for selecting processes as well + struct pointer_to_global_abs_logger_list *global_all_volume_abs_logger_list = COMP_GETPAR3(Union_init, init, global_all_volume_abs_logger_list); + add_element_to_abs_logger_list(global_all_volume_abs_logger_list, abs_logger_list_element); + } + + + %} + +TRACE +%{ +%} + +SAVE COPY Monitor_nD + +FINALLY COPY Monitor_nD EXTEND +%{ +// Remember to clean up allocated lists +if (this_abs_storage.temp_abs_nD_data.allocated_elements>0) free(this_abs_storage.temp_abs_nD_data.elements); + +if (accepted_processes.num_elements > 0) free(accepted_processes.elements); +if (accepted_volumes.num_elements > 0) free(accepted_volumes.elements); +%} + +MCDISPLAY COPY Monitor_nD + +END + diff --git a/mcstas-comps/union/Union_master.comp b/mcstas-comps/union/Union_master.comp index 5a2c9fac9..97c94a753 100755 --- a/mcstas-comps/union/Union_master.comp +++ b/mcstas-comps/union/Union_master.comp @@ -58,6 +58,7 @@ SETTING PARAMETERS(verbal = 1, history_limit=300000, enable_conditionals=1, inherit_number_of_scattering_events=0, + weight_ratio_limit=0.0000001, string init="init") NOACC @@ -987,6 +988,12 @@ TRACE printf("\n\n\n\n\n----------- NEW RAY -------------------------------------------------\n"); printf("Union_master component name: %s \n \n",NAME_CURRENT_COMP); #endif + + double start_weight; + start_weight = p; + + double weight_limit; + weight_limit = p*weight_ratio_limit; // Initialize logic done = 0; @@ -1107,6 +1114,14 @@ TRACE print_intersection_table(&intersection_time_table); printf("current_volume = %d \n",current_volume); #endif + + if (weight_ratio_limit && p < weight_limit) { + #ifdef Union_trace_verbal_setting + printf("Weight reduced more than ratio_limit p=%lf, p0=%lf, (p_limit=%lf, limit=%d)\n", p, start_weight, weight_limit, limit); + #endif + //printf("Weight reduced more than ratio_limit p=%30.28lf, p0=%15.13lf, (p_limit=%15.13lf, scatter=%d)\n", p, start_weight, weight_limit, 100000-limit); + ABSORB; + } // Calculating intersections with the necessary volumes. The relevant set of volumes depend on the current volume and the mask status array. // First the volumes on the current volumes intersect list is checked, then its mask interset list. Before checking the volume itself, it is From e5976ee40f85ac482a462630c4ed384cc68040aa Mon Sep 17 00:00:00 2001 From: Peter Willendrup Date: Fri, 20 Dec 2024 09:49:17 +0100 Subject: [PATCH 2/3] Common return statement for slaves / master alike --- common/lib/share/mccode-r.c | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/common/lib/share/mccode-r.c b/common/lib/share/mccode-r.c index 9ff689ed4..e531ace83 100644 --- a/common/lib/share/mccode-r.c +++ b/common/lib/share/mccode-r.c @@ -2037,8 +2037,9 @@ int mcdetector_out_data_nexus(NXhandle f, MCDETECTOR detector) mcdetector_out_array_nexus(f, "ncount", detector.p0, detector); } else if (strcasestr(detector.format, "pixels")) { mcdetector_out_array_nexus( f, "pixels", detector.p1, detector); - } else + } else { mcdetector_out_array_nexus( f, "events", detector.p1, detector); + } NXclosegroup(f); NXopengroup(f, data_name, "NXdata"); NXgetgroupID(nxhandle, &pLink); @@ -2089,7 +2090,6 @@ MCDETECTOR mcdetector_out_list_slaves(MCDETECTOR detector) || mc_MPI_Send(detector.p1, mnp[0]*mnp[1]*mnp[2], MPI_DOUBLE, mpi_node_root) != MPI_SUCCESS) fprintf(stderr, "Warning: proc %i to master: MPI_Send p1 list error: mnp=%i (mcdetector_out_list_slaves)\n", mpi_node_rank, abs(mnp[0]*mnp[1]*mnp[2])); /* slaves are done: sent mnp and p1 */ - return (detector); } /* end slaves */ /* MPI master: receive data from slaves sequentially: 2 MPI_Recv calls */ @@ -2122,6 +2122,7 @@ MCDETECTOR mcdetector_out_list_slaves(MCDETECTOR detector) printf("\n** Done ** \n"); ); } + // Common return statement for slaves / master alike return(detector); } #endif From f30e837ba8dc2c623f44a7fe7f30ace29e402e36 Mon Sep 17 00:00:00 2001 From: Peter Willendrup Date: Fri, 20 Dec 2024 09:50:23 +0100 Subject: [PATCH 3/3] Increase MPI_REDUCE_BLOCKSIZE 100 fold. Improves NeXus output stability in massively parallel (>500 core) scenarios. --- common/lib/share/mccode-r.h.in | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/common/lib/share/mccode-r.h.in b/common/lib/share/mccode-r.h.in index 879ce0b66..d5348e208 100644 --- a/common/lib/share/mccode-r.h.in +++ b/common/lib/share/mccode-r.h.in @@ -376,7 +376,7 @@ void destroy_darr3d(DArray3d a); } #ifndef MPI_REDUCE_BLOCKSIZE -#define MPI_REDUCE_BLOCKSIZE 1000 +#define MPI_REDUCE_BLOCKSIZE 100000 #endif int mc_MPI_Sum(double* buf, long count);