From 194f5676f16f105caa0dd65dbaea1fe1ce0031dd Mon Sep 17 00:00:00 2001 From: geogunow Date: Wed, 23 Mar 2016 22:42:25 -0400 Subject: [PATCH 1/3] Updated solver to include semi-actual ray tracing by z-stack --- src/init.c | 5 +- src/solver.c | 301 +++++++++++++++++++++++++++++++-------------------- 2 files changed, 187 insertions(+), 119 deletions(-) diff --git a/src/init.c b/src/init.c index 61b7f5b..22ceede 100644 --- a/src/init.c +++ b/src/init.c @@ -20,7 +20,10 @@ void calculate_derived_inputs( Input * I ) I->ntracks_2D = 2 * ( I->ntracks_2D / 2 ); - I->z_stacked = (int) ( I->height / (I->axial_z_sep * I->decomp_assemblies_ax)); + // calculate the number of tracks + // note that tracks extend roughly twice the height of the domain + I->z_stacked = (int) ( 2 * I->height / + (I->axial_z_sep * I->decomp_assemblies_ax)); I->ntracks = I->ntracks_2D * I->n_polar_angles * I->z_stacked; I->domain_height = I->height / I->decomp_assemblies_ax; diff --git a/src/solver.c b/src/solver.c index 0027db1..ea6c6cb 100644 --- a/src/solver.c +++ b/src/solver.c @@ -371,119 +371,198 @@ void transport_sweep( Params * params, Input * I ) pos_z_dir = false; float p_angle = params->polar_angles[j]; float mu = cos(p_angle); + float sin_theta = sin(p_angle) - // start with all z stacked rays - int begin_stacked = 0; - int end_stacked = I->z_stacked; + // estimate the upper and lower points of the first track in the + // z-stack + float first_track_lower_z = -I->assembly_width / sin_theta; + float first_track_upper_z = 0; + // loop over all 2D segments in the 2D track for( int n = 0; n < params->tracks_2D[i].n_segments; n++) { - // calculate distance traveled in cell if segment completed - float s_full = params->tracks_2D[i].segments[n].length - / sin(p_angle); - - // allocate varaible for distance traveled in an FSR - float ds = 0; - - // loop over remaining z-stacked rays - for( int k = begin_stacked; k < end_stacked; k++) - { - // initialize s to full length - float s = s_full; - - // select current track - Track * track = ¶ms->tracks[i][j][k]; - - // set flag for completeion of segment - bool seg_complete = false; - - // calculate interval - int curr_interval; - if( pos_z_dir) - curr_interval = get_pos_interval(track->z_height, - fine_delta_z); - else - curr_interval = get_neg_interval(track->z_height, - fine_delta_z); - - while( !seg_complete ) - { - // flag to reset z position - bool reset = false; - - - /* calculate new height based on s - * (distance traveled in FSR) */ - float z = track->z_height + s * cos(p_angle); - - // check if still in same FSR (fine axial interval) - int new_interval; - if( pos_z_dir ) - new_interval = get_pos_interval(z, - fine_delta_z); - else - new_interval = get_neg_interval(z, - fine_delta_z); - - if( new_interval == curr_interval ) + // determine the start and end points in the first track + float first_start_z; + float first_end_z; + if( pos_z_dir ) + { + first_start_z = first_track_lower_z; + first_end_z = first_track_upper_z; + } + else + { + first_start_z = first_track_upper_z; + first_end_z = first_track_lower_z; + } + + // loop over all 3D SRs associated with the 2D segment + int num_ax_srs = I->cai*I->fai; + for( int z_iter = 0; z_iter < num_ax_srs; z_iter++) + { + // pick a random SR (cache miss expected) + #ifdef OPENMP + long QSR_id = rand_r(&seed) % + I->n_source_regions_per_node; + #else + long QSR_id = rand() % + I->n_source_regions_per_node; + #endif + + + // If traveling in the negative direction, loop thorugh + // SRs from the top + int z_ind; + if ( pos_z_dir ) + z_ind = z_iter; + else + z_ind = num_ax_srs - z_iter - 1; + + // calculate the boundaries of the source region + float z_min = fine_delta_z * z_ind; + float z_max = z_min + fine_delta_z; + + // compute track indexes that cross the source region + int start_track = ceil((z_min - first_track_upper_z) / + z_spacing); + int start_full = ceil((z_min - first_track_lower_z) / + z_spacing); + int end_full = ceil((z_max - first_track_upper_z) / + z_spacing); + int end_track = ceil((z_max - first_track_lower_z) / + z_spacing); + + // Treat lower tracks that do not cross the entire 2D + // length + int min_lower = min(start_full, end_full); + for( int k = start_track; k < min_lower; k++) { + + // select current track + Track * track = ¶ms->tracks[i][j][k]; + + // calculate the distance traveled in the SR + float end_z = first_track_upper_z + k*fine_delta_z; + float ds = (end_z - z_min) / abs(mu); + + /* update sources and fluxes from attenuation + * over SR */ + if( I->axial_exp == 2 ) { - seg_complete = true; - ds = s; + attenuate_fluxes( track, true, + ¶ms->sources[QSR_id], + I, params, ds, mu, + params->tracks_2D[i].az_weight, &A ); + + segments_processed++; } + else if( I->axial_exp == 0 ) + { + attenuate_FSR_fluxes( track, true, + ¶ms->sources[QSR_id], + I, params, ds, mu, + params->tracks_2D[i].az_weight, &A ); - // otherwise, we need to recalculate distances + segments_processed++; + } else { - // correct z - if( pos_z_dir ) - { - curr_interval++; - z = fine_delta_z * (float) curr_interval; - } - else{ - curr_interval--; - z = fine_delta_z * (float) curr_interval; - } - - // calculate distance travelled in FSR (ds) - ds = (z - track->z_height) / cos(p_angle); - - // update track length remaining - s -= ds; - - /* check remaining track length to protect - * against potential roundoff errors */ - if( s <= 0 ) - seg_complete = true; - - // check if out of bounds or track complete - if( z <= 0 || z >= node_delta_z ) - { - // mark segment as completed - seg_complete = true; - - // remember to no longer treat this track - if ( pos_z_dir ) - end_stacked--; - else - begin_stacked++; - - // reset z height - reset = true; - } + printf("Error: invalid axial expansion order"); + printf("\n Please input 0 or 2\n"); + exit(1); } - - // pick a random FSR (cache miss expected) - #ifdef OPENMP - long QSR_id = rand_r(&seed) % - I->n_source_regions_per_node; - #else - long QSR_id = rand() % - I->n_source_regions_per_node; - #endif - - /* update sources and fluxes from attenuation - * over FSR */ + } + + // Find if there are tracks traversing the entire 2D length + if( end_full > start_full ) { + float ds = params->tracks_2D[i].segments[n].length + / sin(p_angle); + for (int k = start_full; k < end_full; k++) { + + // select current track + Track * track = ¶ms->tracks[i][j][k]; + + /* update sources and fluxes from attenuation + * over SR */ + if( I->axial_exp == 2 ) + { + attenuate_fluxes( track, true, + ¶ms->sources[QSR_id], + I, params, ds, mu, + params->tracks_2D[i].az_weight, &A ); + + segments_processed++; + } + else if( I->axial_exp == 0 ) + { + attenuate_FSR_fluxes( track, true, + ¶ms->sources[QSR_id], + I, params, ds, mu, + params->tracks_2D[i].az_weight, &A ); + + segments_processed++; + } + else + { + printf("Error: invalid axial expansion order"); + printf("\n Please input 0 or 2\n"); + exit(1); + } + } + } + + // determine if there are tracks that cross the entire + // z-height of the source region (mutually exclusive + // with any traveling the entire 2D distance) + else if( start_full > end_full ) + { + float ds = (z_max - z_min) / abs(mu); + + for (int k = end_full; k < start_full; k++) { + + // select current track + Track * track = ¶ms->tracks[i][j][k]; + + /* update sources and fluxes from attenuation + * over SR */ + if( I->axial_exp == 2 ) + { + attenuate_fluxes( track, true, + ¶ms->sources[QSR_id], + I, params, ds, mu, + params->tracks_2D[i].az_weight, &A ); + + segments_processed++; + } + else if( I->axial_exp == 0 ) + { + attenuate_FSR_fluxes( track, true, + ¶ms->sources[QSR_id], + I, params, ds, mu, + params->tracks_2D[i].az_weight, &A ); + + segments_processed++; + } + else + { + printf("Error: invalid axial expansion order"); + printf("\n Please input 0 or 2\n"); + exit(1); + } + } + } + // Treat upper tracks that do not cross the entire 2D + // length + int min_upper = max(start_full, end_full); + for( int k = min_upper; k < end_track; k++) { + + // select current track + Track * track = ¶ms->tracks[i][j][k]; + + // calculate the distance traveled in the SR + float start_z = first_track_lower_z + k*fine_delta_z; + float ds = (z_max - start_z) / abs(mu); + + /* update sources and fluxes from attenuation + * over SR */ if( I->axial_exp == 2 ) { attenuate_fluxes( track, true, @@ -493,7 +572,6 @@ void transport_sweep( Params * params, Input * I ) segments_processed++; } - else if( I->axial_exp == 0 ) { attenuate_FSR_fluxes( track, true, @@ -509,20 +587,7 @@ void transport_sweep( Params * params, Input * I ) printf("\n Please input 0 or 2\n"); exit(1); } - - // update with new z height or reset if finished - if( n == params->tracks_2D[i].n_segments - 1 - || reset) - { - if( pos_z_dir) - track->z_height = I->axial_z_sep * k; - else - track->z_height = I->axial_z_sep * (k+1); - } - else - track->z_height = z; - - } + } } } } From 82fb72f6f45ecc65157e0c8dad7aa07ba932a7c1 Mon Sep 17 00:00:00 2001 From: geogunow Date: Wed, 6 Apr 2016 20:15:07 -0400 Subject: [PATCH 2/3] Fixed bugs in tracking update --- src/Makefile | 2 +- src/SimpleMOC_header.h | 8 ++++++++ src/solver.c | 10 +++++----- 3 files changed, 14 insertions(+), 6 deletions(-) diff --git a/src/Makefile b/src/Makefile index 8500cf1..9b59376 100644 --- a/src/Makefile +++ b/src/Makefile @@ -2,7 +2,7 @@ # User Options #=============================================================================== -COMPILER = intel +COMPILER = gcc MPI = no OPENMP = yes OPTIMIZE = yes diff --git a/src/SimpleMOC_header.h b/src/SimpleMOC_header.h index 55e2fff..e2e9df5 100644 --- a/src/SimpleMOC_header.h +++ b/src/SimpleMOC_header.h @@ -251,4 +251,12 @@ void fast_transfer_boundary_fluxes( Params params, Input I, CommGrid grid); void transfer_boundary_fluxes( Params params, Input I, CommGrid grid); #endif +// min/max functions +#define min(a,b) ({ __typeof__ (a) _a = (a); \ + __typeof__ (b) _b = (b); \ + _a < _b ? _a : _b; }) +#define max(a,b) ({ __typeof__ (a) _a = (a); \ + __typeof__ (b) _b = (b); \ + _a > _b ? _a : _b; }) + #endif diff --git a/src/solver.c b/src/solver.c index ea6c6cb..b005026 100644 --- a/src/solver.c +++ b/src/solver.c @@ -371,7 +371,7 @@ void transport_sweep( Params * params, Input * I ) pos_z_dir = false; float p_angle = params->polar_angles[j]; float mu = cos(p_angle); - float sin_theta = sin(p_angle) + float sin_theta = sin(p_angle); // estimate the upper and lower points of the first track in the // z-stack @@ -423,13 +423,13 @@ void transport_sweep( Params * params, Input * I ) // compute track indexes that cross the source region int start_track = ceil((z_min - first_track_upper_z) / - z_spacing); + fine_delta_z); int start_full = ceil((z_min - first_track_lower_z) / - z_spacing); + fine_delta_z); int end_full = ceil((z_max - first_track_upper_z) / - z_spacing); + fine_delta_z); int end_track = ceil((z_max - first_track_lower_z) / - z_spacing); + fine_delta_z); // Treat lower tracks that do not cross the entire 2D // length From 07f036a182c0bf0275df42cdb546e73d0aacdbf2 Mon Sep 17 00:00:00 2001 From: geogunow Date: Fri, 8 Apr 2016 13:09:39 -0400 Subject: [PATCH 3/3] Updated absolute value for floating point --- src/solver.c | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/solver.c b/src/solver.c index b005026..e05f99b 100644 --- a/src/solver.c +++ b/src/solver.c @@ -441,7 +441,7 @@ void transport_sweep( Params * params, Input * I ) // calculate the distance traveled in the SR float end_z = first_track_upper_z + k*fine_delta_z; - float ds = (end_z - z_min) / abs(mu); + float ds = (end_z - z_min) / fabs(mu); /* update sources and fluxes from attenuation * over SR */ @@ -514,7 +514,7 @@ void transport_sweep( Params * params, Input * I ) // with any traveling the entire 2D distance) else if( start_full > end_full ) { - float ds = (z_max - z_min) / abs(mu); + float ds = (z_max - z_min) / fabs(mu); for (int k = end_full; k < start_full; k++) { @@ -559,7 +559,7 @@ void transport_sweep( Params * params, Input * I ) // calculate the distance traveled in the SR float start_z = first_track_lower_z + k*fine_delta_z; - float ds = (z_max - start_z) / abs(mu); + float ds = (z_max - start_z) / fabs(mu); /* update sources and fluxes from attenuation * over SR */