diff --git a/mcstas-comps/share/union-lib.c b/mcstas-comps/share/union-lib.c index e7f62da10..2ff1696fc 100644 --- a/mcstas-comps/share/union-lib.c +++ b/mcstas-comps/share/union-lib.c @@ -25,7 +25,6 @@ void scattering_process_struct_init(struct scattering_process_struct *sps) sps->scattering_function = NULL; sps->non_isotropic_rot_index = -1; sps->needs_cross_section_focus = -1; - sps->needs_numerical_integration = -1; sps->sampling_points = -1; } @@ -8257,6 +8256,283 @@ void overwrite_if_empty(char *input_string, char *overwrite) { } } +//============================================================================== +//=========== Helper functions to Union_master ================================= +//============================================================================== + int + volume_is_only_absorber (struct Volume_struct* Volume) { + // This function returns true if a volume does not have any physical processes + // and if the volume is not a vacuum. + if (!Volume->p_physics->number_of_processes && !Volume->p_physics->is_vacuum) + return 1; + return 0; + } + + int + process_needs_inhomogenous_sampling (struct physics_struct* physics, struct scattering_process_struct* process) { + if (physics->sampling_points != -1) { + if (process->needs_cross_section_focus || process->sampling_points != -1) + return 1; + } + return 0; + } + + void + adjust_abs_weight_factor (struct Volume_struct* Volume, double* my_sum_plus_abs, + double* length_to_boundary, double v_length, double time_to_boundary, + double* abs_weight_factor, int* abs_weight_factor_set) { + *my_sum_plus_abs = Volume->p_physics->my_a * (2200 / v_length); + *length_to_boundary = time_to_boundary * v_length; + + *abs_weight_factor = exp (-Volume->p_physics->my_a * 2200 * time_to_boundary); + *abs_weight_factor_set = 1; + + #ifdef Union_trace_verbal_setting + printf ("name of material: %s \n", Volume->name); + printf ("length to boundery = %f\n", *length_to_boundary); + printf ("absorption cross section = %f\n", Volume->p_physics->my_a); + printf ("chance to get through this length of absorber: %f %%\n", 100 * exp (-Volume->p_physics->my_a * *length_to_boundary)); + #endif + } + + + void + transform_wavevector_into_local_coord_system (struct Volume_struct* Volume, Coords* wavevector_rotated, double* k_rotated, int p_index, + Coords* wavevector, Coords* ray_position_geometry) { + + int non_isotropic_rot_index = Volume->p_physics->p_scattering_array[p_index].non_isotropic_rot_index; + *wavevector_rotated = rot_apply (Volume->geometry.process_rot_matrix_array[non_isotropic_rot_index], *wavevector); + coords_get (*wavevector_rotated, &k_rotated[0], &k_rotated[1], &k_rotated[2]); + + if (Volume->p_physics->p_scattering_array[p_index].needs_cross_section_focus == 1) { + // Prepare focus data using ray_position_geometry of forced scattering point which will be prepared if any process needs cross_section time + // focusing + Coords ray_position_geometry_rotated = rot_apply (Volume->geometry.process_rot_matrix_array[non_isotropic_rot_index], *ray_position_geometry); + + int focus_data_index = Volume->geometry.focus_array_indices.elements[p_index]; + struct focus_data_struct* this_focus_data = &Volume->geometry.focus_data_array.elements[focus_data_index]; + this_focus_data->RayAim = coords_sub (this_focus_data->Aim, ray_position_geometry_rotated); // Aim vector for this ray + + #ifdef Union_trace_verbal_setting + printf ("Checking process number : %d, it was not isotropic, so RayAim updated \n", p_index); + print_position (*ray_position_geometry, "ray_position_geometry"); + print_position (ray_position_geometry_rotated, "ray_position_geometry_rotated"); + print_position (this_focus_data->RayAim, "this_focus_data->RayAim"); + #endif + } + } + + + + void + move_and_aim_neutron (struct physics_struct* physics, int i, struct scattering_process_struct* process, _class_particle* _particle, + struct Volume_struct* Volume, int p_index, Coords* ray_velocity, Coords* ray_position, struct focus_data_struct* this_focus_data) { + // Transport neutron to place inside geometry + *ray_velocity = coords_set (_particle->vx, _particle->vy, _particle->vz); + // Find location of scattering point in master coordinate system without changing main position / velocity variables + Coords direction = coords_scalar_mult (*ray_velocity, 1.0 / length_of_position_vector (*ray_velocity)); + Coords sampling_displacement = coords_scalar_mult (direction, physics->cumul_dists[i]); + Coords sampling_point = coords_add (*ray_position, sampling_displacement); + Coords sampling_point_geometry = coords_sub (sampling_point, Volume->geometry.center); + // Also focus the ray at this point, if the component needs focusing + if (process->needs_cross_section_focus) { + if (Volume->p_physics->p_scattering_array[p_index].non_isotropic_rot_index != -1) { + int non_isotropic_rot_index = Volume->p_physics->p_scattering_array[p_index].non_isotropic_rot_index; + sampling_point_geometry = rot_apply (Volume->geometry.process_rot_matrix_array[non_isotropic_rot_index], sampling_point_geometry); + } + this_focus_data->RayAim = coords_sub (this_focus_data->Aim, sampling_point_geometry); + } + // Calculate mu and probability + coords_get (sampling_point_geometry, &_particle->x, &_particle->y, &_particle->z); + } + + int + mu_and_intersect_dist_safeguard (double mu_sum, double length_to_boundary, double safety_distance2, int* scattering_event) { + if (mu_sum < 1E-18) { + *scattering_event = 0; + return 0; + } + if (length_to_boundary < safety_distance2) { + *scattering_event = 0; + return 0; + } + return 1; + } + + + void + focus_in_cross_section_set_forced_point_and_dir (double* forced_length_to_scattering, double safety_distance, double safety_distance2, double length_to_boundary, + _class_particle* _particle, Coords* ray_velocity, Coords* ray_position_geometry, Coords* ray_position, + struct Volume_struct* Volume, struct focus_data_struct* this_focus_data) { + // Sample length_to_scattering in linear manner + *forced_length_to_scattering = safety_distance + rand01 () * (length_to_boundary - safety_distance2); + + *ray_velocity = coords_set (_particle->vx, _particle->vy, _particle->vz); // Test for root cause + // Find location of scattering point in master coordinate system without changing main position / velocity variables + Coords direction = coords_scalar_mult (*ray_velocity, 1.0 / length_of_position_vector (*ray_velocity)); + Coords scattering_displacement = coords_scalar_mult (direction, *forced_length_to_scattering); + Coords forced_ray_scattering_point = coords_add (*ray_position, scattering_displacement); + *ray_position_geometry = coords_sub (forced_ray_scattering_point, Volume->geometry.center); // ray_position relative to geometry center + + // Calculate the aim for non isotropic processes + this_focus_data = &Volume->geometry.focus_data_array.elements[0]; + this_focus_data->RayAim = coords_sub (this_focus_data->Aim, *ray_position_geometry); // Aim vector for this ray + + #ifdef Union_trace_verbal_setting + printf ("Prepared for focus in cross section calculation in volume: %s \n", Volume->name); + printf ("forced_length_to_scattering =%lf \n", forced_length_to_scattering); + print_position (*ray_position, "ray_position"); + print_position (direction, "direction"); + print_position (scattering_displacement, "scattering_displacement"); + print_position (forced_ray_scattering_point, "forced_ray_scattering_point"); + print_position (*ray_position_geometry, "ray_position_geometry"); + printf ("for isotropic processes this RayAim is used \n"); + print_position (this_focus_data->RayAim, "this_focus_data->RayAim"); + #endif + } + + void + focus_in_cross_section_set_scat_length (double* length_to_scattering, double forced_length_to_scattering, double* weight, double length_to_boundary, + double my_sum_plus_abs) { + // Respect forced length to scattering chosen by process + *length_to_scattering = forced_length_to_scattering; + // Drawing between 0 and L from constant s = 1/L and should have been q = A*exp(-kz). + // Normalizing A*exp(-kz) over 0 to L: A = k/(1-exp(-k*L)) + // Weight correction is ratio between s and q, L*A*exp(-kz) = L*k*exp(-kz)/(1-exp(-Lk)) + *weight *= length_to_boundary * my_sum_plus_abs * exp (-*length_to_scattering * my_sum_plus_abs) / (1.0 - exp (-length_to_boundary * my_sum_plus_abs)); + #ifdef Union_trace_verbal_setting + printf ("Used forced length to scattering, %lf \n", *length_to_scattering); + #endif + } + + + void + inhomogenous_set_cumul_dist_array (struct physics_struct* physics, int i) { + physics->cumul_dists[i] = (i > 0) ? physics->cumul_dists[i - 1] + physics->dist : physics->dist / 2; + } + + void + inhomogenous_sample_transmission_probability (struct physics_struct* physics, struct Volume_struct* Volume, double* real_transmission_probability, + double v_length) { + // Calculate the probabilities and then add them cumulatively + memset (physics->total_mus, 0, sizeof (double) * physics->sampling_points); + + for (int i = 0; i < Volume->p_physics->number_of_processes; i++) { + struct scattering_process_struct* process_i = &Volume->p_physics->p_scattering_array[i]; + if (process_i->sampling_points == -1){ + for (int j = 0; j < physics->sampling_points; j++) { + physics->total_mus[j] += physics->mus[i][0]; + } + } else { + for (int j = 0; j < physics->sampling_points; j++) { + physics->total_mus[j] += physics->mus[i][j]; + } + } + } + double mu_abs_at_speed = Volume->p_physics->my_a * (2200 / v_length); + for (int j = 0; j < physics->sampling_points; j++) { + physics->total_mus[j] += mu_abs_at_speed; + } + double trans_prob; + for (int i = 0; i < physics->sampling_points; i++) { + trans_prob = exp (-physics->total_mus[i] * physics->dist); + if (i == 0) + physics->cumul_transmission_prob[i] = trans_prob; + else + physics->cumul_transmission_prob[i] = physics->cumul_transmission_prob[i - 1] * trans_prob; + } + + *real_transmission_probability = physics->cumul_transmission_prob[physics->sampling_points - 1]; + } + + double + inhomogenous_sample_scattering_point (struct Volume_struct* Volume, struct physics_struct* physics, _class_particle* _particle, double* abs_weight_factor, double v_length, + double safety_distance, int* selected_sampling) { + + // Numerical integration happens, and therefore we must choose between the different samples + // We do this by drawing a random number between 0 and max cumul prob, + // and then seeing which cumul prob is the first to include it. + *abs_weight_factor = 1; + double mu_abs_at_speed = Volume->p_physics->my_a * (2200 / v_length); + double pseudo_rand = rand01 () * (1 - physics->cumul_transmission_prob[physics->sampling_points - 1]); + for (int i = 0; i < physics->sampling_points; i++) { + if (pseudo_rand < 1 - physics->cumul_transmission_prob[i]){ + *selected_sampling = i; + break; + } + } + *abs_weight_factor + *= (physics->total_mus[*selected_sampling] - mu_abs_at_speed ) / physics->total_mus[*selected_sampling]; + + double sampled_dist = safety_distance + - log (1.0 - rand01 () * (1.0 - exp (-physics->total_mus[*selected_sampling] * physics->dist))) + / physics->total_mus[*selected_sampling]; + return physics->cumul_dists[*selected_sampling] - physics->dist / 2 + sampled_dist; + } + + + void + inhomogenous_choose_process (struct physics_struct* physics, struct Volume_struct* Volume, double* culmative_probability, double mc_prop, + double my_sum, int selected_sampling, int* selected_process) { + for (int i = 0; i < Volume->p_physics->number_of_processes; i++) { + *culmative_probability += physics->mus[i][selected_sampling] / my_sum; + if (*culmative_probability > mc_prop) { + *selected_process = i; + break; + } + } + } + + + int + p_interact_is_set (struct Volume_struct* Volume) { + if (Volume->geometry.geometry_p_interact != 0) { + return 1; + } + return 0; + } + + void + p_interact_check_scattering_event (struct Volume_struct* Volume, _class_particle* _particle, int* scattering_event, double* weight, double real_transmission_prob) { + double mc_transmission_prob = 1 - Volume->geometry.geometry_p_interact; + *scattering_event = rand01 () > mc_transmission_prob; + if (*scattering_event) { + // Scattering event happens, this is the correction for the weight + *weight *= (1.0 - real_transmission_prob) / (1.0 - mc_transmission_prob); + } else { + // Scattering event does not happen, this is the appropriate correction + *weight *= real_transmission_prob / mc_transmission_prob; + } + } + + void + p_interact_select_process (struct Volume_struct* Volume, _class_particle* _particle, double* my_trace_fraction_control, double* my_trace, double* total_process_interact, + double* culmative_probability, double* mc_prop, double* weight, double my_sum, int* selected_process) { + // Interact_fraction is used to influence the choice of process in this material + *mc_prop = rand01 (); + *culmative_probability = 0; + *total_process_interact = 1.0; + + // If any of the processes have probability 0, they are excluded from the selection + for (int i = 0; i < Volume->p_physics->number_of_processes; i++) { + if (my_trace[i] < 1E-18) { + // When this happens, the total force probability is corrected and the probability for this particular instance is set to 0 + *total_process_interact -= Volume->p_physics->p_scattering_array[i].process_p_interact; + my_trace_fraction_control[i] = 0; + // In cases where my_trace is not zero, the forced fraction is still used. + } else + my_trace_fraction_control[i] = Volume->p_physics->p_scattering_array[i].process_p_interact; + } + // Randomly select a process using the weights stored in my_trace_fraction_control divided by total_process_interact + for (int i = 0; i < Volume->p_physics->number_of_processes; i++) { + *culmative_probability += my_trace_fraction_control[i] / *total_process_interact; + if (*culmative_probability > *mc_prop) { + *selected_process = i; + *weight *= (my_trace[i] / my_sum) * (*total_process_interact / my_trace_fraction_control[i]); + break; + } + } + } struct union_state_struct *union_acquire_state(int master_present) { (void)master_present; union_state.refcount++; diff --git a/mcstas-comps/share/union-lib.h b/mcstas-comps/share/union-lib.h index 45e70d4c4..82d126b0c 100755 --- a/mcstas-comps/share/union-lib.h +++ b/mcstas-comps/share/union-lib.h @@ -637,7 +637,6 @@ struct scattering_process_struct double process_p_interact; // double between 0 and 1 that describes the fraction of events forced to undergo this process. -1 for disable int non_isotropic_rot_index; // -1 if process is isotrpic, otherwise is the index of the process rotation matrix in the volume int needs_cross_section_focus; // 1 if physics_my needs to call focus functions, otherwise -1 - int needs_numerical_integration; // 1 if the process is inhomogenous and therefore needs numerical integration, otherwise -1. Rotation rotation_matrix; // rotation matrix of process, reported by component in local frame, transformed and moved to volume struct in main double *inhomogenous_cumul_prob; // The cumulative probability of a process in case of inhomogenous processes double *inhomogenous_distances; // The distance of each step in which the cumulative probabilities will be calculated. @@ -1467,6 +1466,64 @@ void fill_surface_stack(char *input_string, struct pointer_to_global_surface_lis void overwrite_if_empty(char *input_string, char *overwrite); +// ------------- Union_master scattering helpers ------------------------------------------- +int volume_is_only_absorber (struct Volume_struct *Volume); + +int process_needs_inhomogenous_sampling (struct physics_struct *current_p_physics, + struct scattering_process_struct *process); + +void adjust_abs_weight_factor (struct Volume_struct *Volume, double *my_sum_plus_abs, + double *length_to_boundary, double v_length, double time_to_boundary, + double *abs_weight_factor, int *abs_weight_factor_set); + +void transform_wavevector_into_local_coord_system (struct Volume_struct *Volume, Coords *wavevector_rotated, + double *k_rotated, int p_index, Coords *wavevector, + Coords *ray_position_geometry); + +void move_and_aim_neutron (struct physics_struct *current_p_physics, int i, + struct scattering_process_struct *process, _class_particle *_particle, + struct Volume_struct *Volume, int p_index, Coords *ray_velocity, + Coords *ray_position, struct focus_data_struct *this_focus_data); + +int mu_and_intersect_dist_safeguard (double mu_sum, double length_to_boundary, double safety_distance2, + int *scattering_event); + +void focus_in_cross_section_set_forced_point_and_dir (double *forced_length_to_scattering, + double safety_distance, double safety_distance2, + double length_to_boundary, _class_particle *_particle, + Coords *ray_velocity, Coords *ray_position_geometry, + Coords *ray_position, struct Volume_struct *Volume, + struct focus_data_struct *this_focus_data); + +void focus_in_cross_section_set_scat_length (double *length_to_scattering, double forced_length_to_scattering, + double *weight, double length_to_boundary, double my_sum_plus_abs); + +void inhomogenous_set_cumul_dist_array (struct physics_struct *current_p_physics, int i); + +void inhomogenous_sample_transmission_probability (struct physics_struct *current_p_physics, + struct Volume_struct *Volume, + double *real_transmission_probability, double v_length); + +double inhomogenous_sample_scattering_point (struct Volume_struct *Volume, + struct physics_struct *current_p_physics, + _class_particle* _particle, + double *abs_weight_factor, double v_length, + double safety_distance, int *selected_sampling); + +void inhomogenous_choose_process (struct physics_struct *current_p_physics, struct Volume_struct *Volume, + double *culmative_probability, double mc_prop, double my_sum, + int selected_sampling, int *selected_process); + +int p_interact_is_set (struct Volume_struct *Volume); + +void p_interact_check_scattering_event (struct Volume_struct *Volume, _class_particle* _particle, int *scattering_event, + double *weight, double real_transmission_prob); + +void p_interact_select_process (struct Volume_struct *Volume, _class_particle* _particle, double *my_trace_fraction_control, + double *my_trace, double *total_process_interact, + double *culmative_probability, double *mc_prop, double *weight, + double my_sum, int *selected_process); + // ------------- Shared Union state ---------------------------------------------------------- // Union components communicate through lists that every component appends to // in INITIALIZE and each Union_master consumes. The lists live in one instance @@ -1521,4 +1578,3 @@ void union_free_state(void); void union_release(const char *comp_name); #endif /* UNION_LIB_H */ - diff --git a/mcstas-comps/union/Inhomogenous_incoherent_process.comp b/mcstas-comps/union/Inhomogenous_incoherent_process.comp index a796d94f3..bdab50263 100755 --- a/mcstas-comps/union/Inhomogenous_incoherent_process.comp +++ b/mcstas-comps/union/Inhomogenous_incoherent_process.comp @@ -369,7 +369,6 @@ INITIALIZE This_process.data_transfer.Inhomogenous_incoherent_struct = &Inhomogenous_storage; This_process.probability_for_scattering_function = &Inhomogenous_incoherent_physics_my; This_process.scattering_function = &Inhomogenous_incoherent_physics_scattering; - This_process.needs_numerical_integration = 1; This_process.sampling_points = number_of_sample_points; // This will be the same for all process's, and can thus be moved to an include. diff --git a/mcstas-comps/union/Union_master.comp b/mcstas-comps/union/Union_master.comp index 5ec11608c..78e86c8d1 100755 --- a/mcstas-comps/union/Union_master.comp +++ b/mcstas-comps/union/Union_master.comp @@ -84,6 +84,8 @@ SHARE #ifndef MASTER_DETECTOR #define MASTER_DETECTOR dummy #endif + + %} DECLARE @@ -851,10 +853,9 @@ INITIALIZE for (int mat_idx = 0; mat_idx < global_material_list_master->num_elements; mat_idx++) { struct physics_struct* physics = global_material_list_master->elements[mat_idx].physics; + physics->sampling_points = -1; // Initialize n points to avoid undefined behaviour for (int proc_idx = 0; proc_idx < physics->number_of_processes; proc_idx++) { struct scattering_process_struct spec_process = physics->p_scattering_array[proc_idx]; - if (spec_process.needs_numerical_integration != 1) - continue; if (physics->sampling_points < spec_process.sampling_points) physics->sampling_points = spec_process.sampling_points; } @@ -1541,28 +1542,16 @@ TRACE scattering_event = 0; // Assume a scattering event will not occur // Check if a scattering event should occur - if (current_volume != 0) { // Volume 0 is always vacuum, and if this is the current volume, an event will not occur - if (Volumes[current_volume]->p_physics->number_of_processes == 0) { // If there are no processes, the volume could be vacuum or an absorber - if (Volumes[current_volume]->p_physics->is_vacuum == 0) { - // This volume does not have physical processes but does have an absorption cross section, so the ray weight is reduced accordingly - - my_sum_plus_abs = Volumes[current_volume]->p_physics->my_a * (2200 / v_length); - length_to_boundary = time_to_boundery * v_length; - - abs_weight_factor = exp (-Volumes[current_volume]->p_physics->my_a * 2200 * time_to_boundery); - abs_weight_factor_set = 1; - - #ifdef Union_trace_verbal_setting - printf ("name of material: %s \n", Volumes[current_volume]->name); - printf ("length to boundery = %f\n", length_to_boundary); - printf ("absorption cross section = %f\n", Volumes[current_volume]->p_physics->my_a); - printf ("chance to get through this length of absorber: %f %%\n", 100 * exp (-Volumes[current_volume]->p_physics->my_a * length_to_boundary)); - #endif - } + if (current_volume != 0) { // Volume 0 is always vacuum, and if this is the current volume, an event will not occur + if (volume_is_only_absorber (Volumes[current_volume])) { // If there are no processes, check if the volume is a pure absorber and not a vacuum + adjust_abs_weight_factor (Volumes[current_volume], &my_sum_plus_abs, &length_to_boundary, v_length, time_to_boundery, &abs_weight_factor, + &abs_weight_factor_set); } else { // Since there is a non-zero number of processes in this material, all the scattering cross section for these are calculated struct physics_struct* current_p_physics = Volumes[current_volume]->p_physics; int selected_sampling = -1; + double forced_length_to_scattering; + Coords ray_position_geometry; my_sum = 0; k[0] = V2K * vx; k[1] = V2K * vy; @@ -1570,56 +1559,10 @@ TRACE p_my_trace = my_trace; wavevector = coords_set (k[0], k[1], k[2]); length_to_boundary = time_to_boundery * v_length; - - double forced_length_to_scattering; - Coords ray_position_geometry; - // If any process in this material needs focusing, sample scattering position and update focus_data accordingly if (current_p_physics->any_process_needs_cross_section_focus == 1) { - // Sample length_to_scattering in linear manner - forced_length_to_scattering = safety_distance + rand01 () * (length_to_boundary - safety_distance2); - - ray_velocity = coords_set (vx, vy, vz); // Test for root cause - // Find location of scattering point in master coordinate system without changing main position / velocity variables - Coords direction = coords_scalar_mult (ray_velocity, 1.0 / length_of_position_vector (ray_velocity)); - Coords scattering_displacement = coords_scalar_mult (direction, forced_length_to_scattering); - Coords forced_ray_scattering_point = coords_add (ray_position, scattering_displacement); - ray_position_geometry - = coords_sub (forced_ray_scattering_point, Volumes[current_volume]->geometry.center); // ray_position relative to geometry center - - // Calculate the aim for non isotropic processes - this_focus_data = &Volumes[current_volume]->geometry.focus_data_array.elements[0]; - this_focus_data->RayAim = coords_sub (this_focus_data->Aim, ray_position_geometry); // Aim vector for this ray - - #ifdef Union_trace_verbal_setting - printf ("Prepared for focus in cross section calculation in volume: %s \n", Volumes[current_volume]->name); - printf ("forced_length_to_scattering =%lf \n", forced_length_to_scattering); - print_position (ray_position, "ray_position"); - print_position (direction, "direction"); - print_position (scattering_displacement, "scattering_displacement"); - print_position (forced_ray_scattering_point, "forced_ray_scattering_point"); - print_position (ray_position_geometry, "ray_position_geometry"); - printf ("for isotropic processes this RayAim is used \n"); - print_position (this_focus_data->RayAim, "this_focus_data->RayAim"); - #endif - - /* - // update focus data for this ray (could limit this to only update the necessary focus_data element, but there are typically very few) - int f_index; - for (f_index=0; f_index < Volumes[current_volume]->geometry.focus_data_array.num_elements; f_index++) { - this_focus_data = &Volumes[current_volume]->geometry.focus_data_array.elements[f_index]; - // Coords ray_position_geometry_rotated = rot_apply(this_focus_data.absolute_rotation, ray_position_geometry); - this_focus_data->RayAim = coords_sub(this_focus_data->Aim, ray_position_geometry); // Aim vector for this ray - } - - printf("calculated forced_length_to_scattering = %lf, new RayAim \n", forced_length_to_scattering); - print_position(direction, "direction"); - print_position(scattering_displacement, "scattering_displacement"); - print_position(forced_ray_scattering_point, "forced_ray_scattering_point"); - print_position(ray_position_geometry, "ray_position_geometry"); - print_position(this_focus_data->RayAim, "this_focus_data->RayAim"); - */ - + focus_in_cross_section_set_forced_point_and_dir (&forced_length_to_scattering, safety_distance, safety_distance2, length_to_boundary, _particle, + &ray_velocity, &ray_position_geometry, &ray_position, Volumes[current_volume], this_focus_data); } else { forced_length_to_scattering = -1.0; // Signals that no forcing needed, could also if on the selected process struct } @@ -1633,24 +1576,8 @@ TRACE if (Volumes[current_volume]->p_physics->p_scattering_array[p_index].non_isotropic_rot_index != -1) { // If the process is not isotropic, the wavevector is transformed into the local coordinate system of the process - int non_isotropic_rot_index = Volumes[current_volume]->p_physics->p_scattering_array[p_index].non_isotropic_rot_index; - wavevector_rotated = rot_apply (Volumes[current_volume]->geometry.process_rot_matrix_array[non_isotropic_rot_index], wavevector); - coords_get (wavevector_rotated, &k_rotated[0], &k_rotated[1], &k_rotated[2]); - - if (Volumes[current_volume]->p_physics->p_scattering_array[p_index].needs_cross_section_focus == 1) { - // Prepare focus data using ray_position_geometry of forced scattering point which will be prepared if any process needs cross_section time - // focusing - Coords ray_position_geometry_rotated - = rot_apply (Volumes[current_volume]->geometry.process_rot_matrix_array[non_isotropic_rot_index], ray_position_geometry); - this_focus_data->RayAim = coords_sub (this_focus_data->Aim, ray_position_geometry_rotated); // Aim vector for this ray - #ifdef Union_trace_verbal_setting - printf ("Checking process number : %d, it was not isotropic, so RayAim updated \n", p_index); - print_position (ray_position_geometry, "ray_position_geometry"); - print_position (ray_position_geometry_rotated, "ray_position_geometry_rotated"); - print_position (this_focus_data->RayAim, "this_focus_data->RayAim"); - #endif - } - + transform_wavevector_into_local_coord_system (Volumes[current_volume], &wavevector_rotated, k_rotated, p_index, &wavevector, + &ray_position_geometry); } else { k_rotated[0] = k[0]; k_rotated[1] = k[1]; @@ -1665,41 +1592,25 @@ TRACE double mu; current_p_physics->dist = length_to_boundary / current_p_physics->sampling_points - safety_distance; - if (process->sampling_points != -1 || (process->needs_cross_section_focus == 1 && current_p_physics->sampling_points != 0)) { - // Populate length and probability arrays: + if (process_needs_inhomogenous_sampling (current_p_physics, process)) { + // If mulitple sampling is required for this process, populate the neutron distance and probability arrays: Coords original_position = coords_set (x, y, z); + double mu; + *p_my_trace = 0; + this_focus_data = &Volumes[current_volume]->geometry.focus_data_array.elements[0]; for (int i = 0; i < current_p_physics->sampling_points; i++) { - current_p_physics->cumul_dists[i] = (i > 0) ? current_p_physics->cumul_dists[i - 1] + current_p_physics->dist : current_p_physics->dist / 2; - // Transport neutron to place inside geometry - ray_velocity = coords_set (vx, vy, vz); - // Find location of scattering point in master coordinate system without changing main position / velocity variables - Coords direction = coords_scalar_mult (ray_velocity, 1.0 / length_of_position_vector (ray_velocity)); - Coords sampling_displacement = coords_scalar_mult (direction, current_p_physics->cumul_dists[i]); - Coords sampling_point = coords_add (ray_position, sampling_displacement); - Coords sampling_point_geometry = coords_sub (sampling_point, Volumes[current_volume]->geometry.center); - // Also focus the ray at this point, if the component needs focusing - if (process->needs_cross_section_focus) { - this_focus_data = &Volumes[current_volume]->geometry.focus_data_array.elements[0]; - if (Volumes[current_volume]->p_physics->p_scattering_array[p_index].non_isotropic_rot_index != -1) { - int non_isotropic_rot_index = Volumes[current_volume]->p_physics->p_scattering_array[p_index].non_isotropic_rot_index; - sampling_point_geometry - = rot_apply (Volumes[current_volume]->geometry.process_rot_matrix_array[non_isotropic_rot_index], sampling_point_geometry); - } - this_focus_data->RayAim = coords_sub (this_focus_data->Aim, sampling_point_geometry); - } + inhomogenous_set_cumul_dist_array (current_p_physics, i); + move_and_aim_neutron (current_p_physics, i, process, _particle, Volumes[current_volume], p_index, &ray_velocity, &ray_position, this_focus_data); // Calculate mu and probability - coords_get (sampling_point_geometry, &x, &y, &z); UNION_PHYSICS_MY (physics_output, process, &mu, k_rotated, this_focus_data, _particle); current_p_physics->mus[p_index][i] = mu; + *p_my_trace += mu / current_p_physics->sampling_points; } - coords_get (original_position, &x, &y, &z); - *p_my_trace = 0; - for (int i = 0; i < current_p_physics->sampling_points; i++) - *p_my_trace += current_p_physics->mus[p_index][i] / current_p_physics->sampling_points; - } else { + } + else {// process does not need inhomogenous sampling UNION_PHYSICS_MY (physics_output, process, p_my_trace, k_rotated, this_focus_data, _particle); - if (current_p_physics->sampling_points != 0) { + if (current_p_physics->sampling_points != -1) { current_p_physics->mus[p_index][0] = *p_my_trace; } } @@ -1724,59 +1635,19 @@ TRACE // New flow:length_to_boundary // Calculate if scattering happens based on my_sub_plus_abs - if (my_sum < 1E-18) { - scattering_event = 0; - } else if (length_to_boundary < safety_distance2) { - scattering_event = 0; - } else { - if (current_p_physics->sampling_points == 0) { + if (mu_and_intersect_dist_safeguard (my_sum, length_to_boundary, safety_distance2, &scattering_event)) { + // First calculate the transmission probability + if (current_p_physics->sampling_points == -1) { real_transmission_probability = exp (-length_to_boundary * my_sum_plus_abs); - } else if (current_p_physics->sampling_points != 0) { - // Calculate the probabilities and then add them cumulatively - memset (current_p_physics->total_mus, 0, sizeof (double) * current_p_physics->sampling_points); - - for (int i = 0; i < Volumes[current_volume]->p_physics->number_of_processes; i++) { - struct scattering_process_struct* process_i = &Volumes[current_volume]->p_physics->p_scattering_array[i]; - if (process_i->needs_numerical_integration != 1) - for (int j = 0; j < current_p_physics->sampling_points; j++) { - current_p_physics->total_mus[j] += current_p_physics->mus[i][0] * current_p_physics->dist; - } - else - for (int j = 0; j < current_p_physics->sampling_points; j++) { - current_p_physics->total_mus[j] += current_p_physics->mus[i][j] * current_p_physics->dist; - } - } - // for (int i =0;isampling_points;i++){ - // printf("\nTotalmu=%g\tinteger=%d\n", current_p_physics->total_mus[i], i); - // } - double mu_at_speed = Volumes[current_volume]->p_physics->my_a * (2200 / v_length); - for (int j = 0; j < current_p_physics->sampling_points; j++) { - current_p_physics->total_mus[j] += mu_at_speed * current_p_physics->dist; - } - double trans_prob; - for (int i = 0; i < current_p_physics->sampling_points; i++) { - trans_prob = exp (-current_p_physics->total_mus[i]); - if (i == 0) - current_p_physics->cumul_transmission_prob[i] = trans_prob; - else - current_p_physics->cumul_transmission_prob[i] = current_p_physics->cumul_transmission_prob[i - 1] * trans_prob; - } - - real_transmission_probability = current_p_physics->cumul_transmission_prob[current_p_physics->sampling_points - 1]; + } else { + inhomogenous_sample_transmission_probability (current_p_physics, Volumes[current_volume], &real_transmission_probability, v_length); } - // printf("Trans prop = %g\n", real_transmission_probability); - if (Volumes[current_volume]->geometry.geometry_p_interact != 0) { - mc_transmission_probability = (1.0 - Volumes[current_volume]->geometry.geometry_p_interact); - if ((scattering_event = (rand01 () > mc_transmission_probability))) { - // Scattering event happens, this is the correction for the weight - p *= (1.0 - real_transmission_probability) / (1.0 - mc_transmission_probability); - } else { - // Scattering event does not happen, this is the appropriate correction - p *= real_transmission_probability / mc_transmission_probability; - } + + // Then check if we scatter. + if (p_interact_is_set (Volumes[current_volume])) { + p_interact_check_scattering_event (Volumes[current_volume], _particle, &scattering_event, &p, real_transmission_probability); } else { // probability to scatter is the natural value - // printf("Real transmission prob %g\n", real_transmission_probability); scattering_event = rand01 () > real_transmission_probability; } } @@ -1792,67 +1663,20 @@ TRACE if (my_sum / my_sum_plus_abs > 1.0) printf ("WARNING: Absorption weight factor above 1! Should not happen! \n"); // Select distance to scattering position - if (current_p_physics->sampling_points != 0) { - // Numerical integration happens, and therefore we must choose between the different samples - // We do this by drawing a random number between 0 and max cumul prob, - // and then seeing which cumul prob is the first to include it. - abs_weight_factor = 1; - double mu_at_speed = Volumes[current_volume]->p_physics->my_a * (2200 / v_length); - double pseudo_rand = rand01 () * (1 - current_p_physics->cumul_transmission_prob[current_p_physics->sampling_points - 1]); - for (int i = 0; i < current_p_physics->sampling_points; i++) { - // printf("\nCumul trans prob = %g\t pseudo rand = %g\n", current_p_physics->cumul_transmission_prob[i], pseudo_rand); - if (pseudo_rand >= 1 - current_p_physics->cumul_transmission_prob[i]) - continue; - selected_sampling = i; - break; - } - abs_weight_factor - *= (current_p_physics->total_mus[selected_sampling] - mu_at_speed * current_p_physics->dist) / current_p_physics->total_mus[selected_sampling]; - - // printf("\nSelected_sampling = %d\n", selected_sampling); - - // printf("dist i = %g\tdist=%g\n", dist_i, dist); - double sampled_dist = safety_distance - - log (1.0 - rand01 () * (1.0 - exp (-current_p_physics->total_mus[selected_sampling]))) - / current_p_physics->total_mus[selected_sampling] * current_p_physics->dist; - length_to_scattering = current_p_physics->cumul_dists[selected_sampling] - current_p_physics->dist / 2 + sampled_dist; + if (current_p_physics->sampling_points != -1) { + length_to_scattering = inhomogenous_sample_scattering_point (Volumes[current_volume], current_p_physics, _particle, &abs_weight_factor, v_length, + safety_distance, &selected_sampling); } // Select process if (Volumes[current_volume]->p_physics->number_of_processes == 1) { // trivial case // Select the only available process, which will always have index 0 selected_process = 0; } else { - if (Volumes[current_volume]->p_physics->interact_control == 1) { - // Interact_fraction is used to influence the choice of process in this material - mc_prop = rand01 (); - culmative_probability = 0; - total_process_interact = 1.0; - - // If any of the processes have probability 0, they are excluded from the selection - for (iterator = 0; iterator < Volumes[current_volume]->p_physics->number_of_processes; iterator++) { - if (my_trace[iterator] < 1E-18) { - // When this happens, the total force probability is corrected and the probability for this particular instance is set to 0 - total_process_interact -= Volumes[current_volume]->p_physics->p_scattering_array[iterator].process_p_interact; - my_trace_fraction_control[iterator] = 0; - // In cases where my_trace is not zero, the forced fraction is still used. - } else - my_trace_fraction_control[iterator] = Volumes[current_volume]->p_physics->p_scattering_array[iterator].process_p_interact; - } - // Randomly select a process using the weights stored in my_trace_fraction_control divided by total_process_interact - for (iterator = 0; iterator < Volumes[current_volume]->p_physics->number_of_processes; iterator++) { - culmative_probability += my_trace_fraction_control[iterator] / total_process_interact; - if (culmative_probability > mc_prop) { - selected_process = iterator; - p *= (my_trace[iterator] / my_sum) * (total_process_interact / my_trace_fraction_control[iterator]); - break; - } - } - - } else { + if (Volumes[current_volume]->p_physics->interact_control != 1) { // Select a process based on their relative attenuations factors mc_prop = rand01 (); culmative_probability = 0; - if (current_p_physics->sampling_points == 0) { + if (current_p_physics->sampling_points == -1) { for (iterator = 0; iterator < Volumes[current_volume]->p_physics->number_of_processes; iterator++) { culmative_probability += my_trace[iterator] / my_sum; if (culmative_probability > mc_prop) { @@ -1861,31 +1685,20 @@ TRACE } } } else { - for (iterator = 0; iterator < Volumes[current_volume]->p_physics->number_of_processes; iterator++) { - culmative_probability += current_p_physics->mus[iterator][selected_sampling] / my_sum; - if (culmative_probability > mc_prop) { - selected_process = iterator; - break; - } - } + inhomogenous_choose_process (current_p_physics, Volumes[current_volume], &culmative_probability, mc_prop, my_sum, selected_sampling, + &selected_process); } + } else { + p_interact_select_process (Volumes[current_volume], _particle, my_trace_fraction_control, my_trace, &total_process_interact, &culmative_probability, + &mc_prop, &p, my_sum, &selected_process); } } process = &Volumes[current_volume]->p_physics->p_scattering_array[selected_process]; - if (current_p_physics->sampling_points == 0) { + if (current_p_physics->sampling_points == -1) { // No numerical integration is necessary. if (process->needs_cross_section_focus == 1) { - // Respect forced length to scattering chosen by process - length_to_scattering = forced_length_to_scattering; - // Drawing between 0 and L from constant s = 1/L and should have been q = A*exp(-kz). - // Normalizing A*exp(-kz) over 0 to L: A = k/(1-exp(-k*L)) - // Weight correction is ratio between s and q, L*A*exp(-kz) = L*k*exp(-kz)/(1-exp(-Lk)) - p *= length_to_boundary * my_sum_plus_abs * exp (-length_to_scattering * my_sum_plus_abs) / (1.0 - exp (-length_to_boundary * my_sum_plus_abs)); - #ifdef Union_trace_verbal_setting - printf ("Used forced length to scattering, %lf \n", length_to_scattering); - #endif - + focus_in_cross_section_set_scat_length (&length_to_scattering, forced_length_to_scattering, &p, length_to_boundary, my_sum_plus_abs); } else { // Decided the ray scatters, choose where on truncated exponential from safety_distance to length_to_boundary - safety_distance length_to_scattering