Skip to content

Commit

Permalink
Simplify distance2boundary calculation
Browse files Browse the repository at this point in the history 
The calculation that needs to be performed is a simple trigonomic
calculation with two possible results. For performance reason one can
introduce a previous check. For a detailed discussion see #497

This results in the following sheme:
 - For mu > 0 (outward propagation) the next boundary ist the outer
   boundary (case 1)
 - For mu < 0 we have to decide whether we hit the inner boundary or
   not and set the distance accordingly (case 2&3)

Additionally branch prediction optimization was performed in
distance2line. The compiler always assumes if statements to be true,
which means the default branch should always the one the if statement
checks for. This slightly reduces instruction misses.

Resolves: #497
  • Loading branch information
yeganer committed Mar 10, 2016
1 parent ca53169 commit 2dbbecf
Show file tree
Hide file tree
Showing 3 changed files with 50 additions and 70 deletions.
78 changes: 35 additions & 43 deletions tardis/montecarlo/src/cmontecarlo.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -137,67 +137,62 @@ void calculate_chi_bf(rpacket_t * packet, storage_model_t * storage)
rpacket_set_chi_boundfree(packet, bf_helper * doppler_factor);
}

double
void
compute_distance2boundary (rpacket_t * packet, const storage_model_t * storage)
{
double r = rpacket_get_r (packet);
double mu = rpacket_get_mu (packet);
double r_outer = storage->r_outer[rpacket_get_current_shell_id (packet)];
double r_inner = storage->r_inner[rpacket_get_current_shell_id (packet)];
double d_outer =
sqrt (r_outer * r_outer + ((mu * mu - 1.0) * r * r)) - (r * mu);
double d_inner;
double check, distance;
if (rpacket_get_recently_crossed_boundary (packet) == 1)
{
rpacket_set_next_shell_id (packet, 1);
return d_outer;
distance = sqrt (r_outer * r_outer + ((mu * mu - 1.0) * r * r)) - (r * mu);
}
else if (mu > 0.0)
{ // direction outward
rpacket_set_next_shell_id (packet, 1);
distance = sqrt (r_outer * r_outer + ((mu * mu - 1.0) * r * r)) - (r * mu);
}
else
{
double check = r_inner * r_inner + (r * r * (mu * mu - 1.0));
if (check < 0.0)
{
rpacket_set_next_shell_id (packet, 1);
return d_outer;
{ // going inward
if ( (check = r_inner * r_inner + (r * r * (mu * mu - 1.0)) )>= 0.0)
{ // hit inner boundary
rpacket_set_next_shell_id (packet, -1);
distance = - r * mu - sqrt (check);
}
else
{
d_inner = mu < 0.0 ? -r * mu - sqrt (check) : MISS_DISTANCE;
{ // miss inner boundary
rpacket_set_next_shell_id (packet, 1);
distance = sqrt (r_outer * r_outer + ((mu * mu - 1.0) * r * r)) - (r * mu);
}
}
if (d_inner < d_outer)
{
rpacket_set_next_shell_id (packet, -1);
return d_inner;
}
else
{
rpacket_set_next_shell_id (packet, 1);
return d_outer;
}
rpacket_set_d_boundary (packet, distance);
}

tardis_error_t
compute_distance2line (const rpacket_t * packet, const storage_model_t * storage,
double *result)
compute_distance2line (rpacket_t * packet, const storage_model_t * storage)
{
tardis_error_t ret_val = TARDIS_ERROR_OK;
if (rpacket_get_last_line (packet))
{
*result = MISS_DISTANCE;
}
else
if (!rpacket_get_last_line (packet))
{
double r = rpacket_get_r (packet);
double mu = rpacket_get_mu (packet);
double nu = rpacket_get_nu (packet);
double nu_line = rpacket_get_nu_line (packet);
double distance, nu_diff;
double t_exp = storage->time_explosion;
double inverse_t_exp = storage->inverse_time_explosion;
int64_t cur_zone_id = rpacket_get_current_shell_id (packet);
double doppler_factor = 1.0 - mu * r * inverse_t_exp * INVERSE_C;
double comov_nu = nu * doppler_factor;
if (comov_nu < nu_line)
if ( (nu_diff = comov_nu - nu_line) >= 0)
{
distance = (nu_diff / nu) * C * t_exp;
rpacket_set_d_line (packet, distance);
return TARDIS_ERROR_OK;
}
else
{
if (rpacket_get_next_line_id (packet) == storage->no_of_lines - 1)
{
Expand Down Expand Up @@ -232,14 +227,14 @@ compute_distance2line (const rpacket_t * packet, const storage_model_t * storage
fprintf (stderr, "nu = %f\n", nu);
fprintf (stderr, "doppler_factor = %f\n", doppler_factor);
fprintf (stderr, "cur_zone_id = %" PRIi64 "\n", cur_zone_id);
ret_val = TARDIS_ERROR_COMOV_NU_LESS_THAN_NU_LINE;
}
else
{
*result = ((comov_nu - nu_line) / nu) * C * t_exp;
return TARDIS_ERROR_COMOV_NU_LESS_THAN_NU_LINE;
}
}
return ret_val;
else
{
rpacket_set_d_line (packet, MISS_DISTANCE);
return TARDIS_ERROR_OK;
}
}

void
Expand Down Expand Up @@ -712,12 +707,9 @@ montecarlo_compute_distances (rpacket_t * packet, storage_model_t * storage)
}
else
{
rpacket_set_d_boundary (packet,
compute_distance2boundary (packet, storage));
double d_line;
compute_distance2line (packet, storage, &d_line);
compute_distance2boundary(packet, storage);
compute_distance2line (packet, storage);
// FIXME MR: return status of compute_distance2line() is ignored
rpacket_set_d_line (packet, d_line);
compute_distance2continuum (packet, storage);
}
}
Expand Down
7 changes: 3 additions & 4 deletions tardis/montecarlo/src/cmontecarlo.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,7 @@ double rpacket_doppler_factor(const rpacket_t *packet, const storage_model_t *st
*
* @return distance to shell boundary
*/
double compute_distance2boundary (rpacket_t * packet,
void compute_distance2boundary (rpacket_t * packet,
const storage_model_t * storage);

/** Calculate the distance the packet has to travel until it redshifts to the first spectral line.
Expand All @@ -43,9 +43,8 @@ double compute_distance2boundary (rpacket_t * packet,
*
* @return distance to the next spectral line
*/
tardis_error_t compute_distance2line (const rpacket_t * packet,
const storage_model_t * storage,
double *result);
tardis_error_t compute_distance2line (rpacket_t * packet,
const storage_model_t * storage);

/** Calculate the distance to the next continuum event, which can be a Thomson scattering, bound-free absorption or
free-free transition.
Expand Down
35 changes: 12 additions & 23 deletions tardis/montecarlo/src/test_cmontecarlo.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -220,23 +220,22 @@ test_compute_distance2boundary(void){
storage_model_t sm;
init_rpacket(&rp);
init_storage_model(&sm);
double D_BOUNDARY = compute_distance2boundary(&rp, &sm);
rpacket_set_d_boundary(&rp, D_BOUNDARY);
compute_distance2boundary(&rp, &sm);
double D_BOUNDARY = rpacket_get_d_boundary(&rp);
dealloc_storage_model(&sm);
return D_BOUNDARY;
return D_BOUNDARY;
}
double
test_compute_distance2line(void){
rpacket_t rp;
storage_model_t sm;
init_rpacket(&rp);
init_storage_model(&sm);
double D_LINE;
// FIXME MR: return status of compute_distance2line() is ignored
compute_distance2line(&rp, &sm, &D_LINE);
rpacket_set_d_line(&rp, D_LINE);
compute_distance2line(&rp, &sm);
double D_LINE = rpacket_get_d_line(&rp);
dealloc_storage_model(&sm);
return D_LINE;
return D_LINE;
}

double
Expand Down Expand Up @@ -279,17 +278,12 @@ test_increment_j_blue_estimator(void){
storage_model_t sm;
init_rpacket(&rp);
init_storage_model(&sm);
int64_t j_blue_idx = 0;
double D_BOUNDARY = compute_distance2boundary(&rp, &sm);
rpacket_set_d_boundary(&rp, D_BOUNDARY);
double D_LINE;
// FIXME MR: return status of compute_distance2line() is ignored
compute_distance2line(&rp, &sm, &D_LINE);
rpacket_set_d_line(&rp, D_LINE);
int64_t j_blue_idx = 0;
compute_distance2line(&rp, &sm);
move_packet(&rp, &sm, 1e13);
double d_line = rpacket_get_d_line(&rp);
increment_j_blue_estimator(&rp, &sm, d_line, j_blue_idx);
double res = sm.line_lists_j_blues[j_blue_idx];
double d_line = rpacket_get_d_line(&rp);
increment_j_blue_estimator(&rp, &sm, d_line, j_blue_idx);
double res = sm.line_lists_j_blues[j_blue_idx];
dealloc_storage_model(&sm);
return res;
}
Expand Down Expand Up @@ -379,12 +373,7 @@ test_calculate_chi_bf(void){
rk_state mt_state;
irandom(&mt_state);
int64_t j_blue_idx = 0;
double D_BOUNDARY = compute_distance2boundary(&rp, &sm);
rpacket_set_d_boundary(&rp, D_BOUNDARY);
double D_LINE;
// FIXME MR: return status of compute_distance2line() is ignored
compute_distance2line(&rp, &sm, &D_LINE);
rpacket_set_d_line(&rp, D_LINE);
compute_distance2line(&rp, &sm);
move_packet(&rp, &sm, 1e13);
double d_line = rpacket_get_d_line(&rp);
increment_j_blue_estimator(&rp, &sm, d_line, j_blue_idx);
Expand Down

0 comments on commit 2dbbecf

Please sign in to comment.