Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Alternative fix for issue #393 #401

Merged
merged 1 commit into from
Aug 26, 2015
Merged

Alternative fix for issue #393 #401

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
83 changes: 39 additions & 44 deletions tardis/montecarlo/src/cmontecarlo.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,14 +4,6 @@
#endif
#include "cmontecarlo.h"

rk_state mt_state;

static void
initialize_random_kit (unsigned long seed)
{
rk_seed (seed, &mt_state);
}

/** Look for a place to insert a value in an inversely sorted float array.
*
* @param x an inversely (largest to lowest) sorted float array
Expand Down Expand Up @@ -297,14 +289,14 @@ compute_distance2continuum(rpacket_t * packet, storage_model_t * storage)
}

int64_t
macro_atom (const rpacket_t * packet, const storage_model_t * storage)
macro_atom (const rpacket_t * packet, const storage_model_t * storage, rk_state *mt_state)
{
int emit = 0, i = 0, probability_idx = -1;
int activate_level =
storage->line2macro_level_upper[rpacket_get_next_line_id (packet) - 1];
while (emit != -1)
{
double event_random = rk_double (&mt_state);
double event_random = rk_double (mt_state);
i = storage->macro_block_references[activate_level] - 1;
double p = 0.0;
do
Expand Down Expand Up @@ -378,12 +370,12 @@ increment_j_blue_estimator (const rpacket_t * packet, storage_model_t * storage,

int64_t
montecarlo_one_packet (storage_model_t * storage, rpacket_t * packet,
int64_t virtual_mode)
int64_t virtual_mode, rk_state *mt_state)
{
int64_t reabsorbed=-1;
if (virtual_mode == 0)
{
reabsorbed = montecarlo_one_packet_loop (storage, packet, 0);
reabsorbed = montecarlo_one_packet_loop (storage, packet, 0, mt_state);
}
else
{
Expand All @@ -406,7 +398,7 @@ montecarlo_one_packet (storage_model_t * storage, rpacket_t * packet,
mu_min = 0.0;
}
double mu_bin = (1.0 - mu_min) / rpacket_get_virtual_packet_flag (packet);
rpacket_set_mu(&virt_packet,mu_min + (i + rk_double (&mt_state)) * mu_bin);
rpacket_set_mu(&virt_packet,mu_min + (i + rk_double (mt_state)) * mu_bin);
switch (virtual_mode)
{
case -2:
Expand Down Expand Up @@ -434,7 +426,7 @@ montecarlo_one_packet (storage_model_t * storage, rpacket_t * packet,
rpacket_set_energy(&virt_packet,
rpacket_get_energy (packet) * doppler_factor_ratio);
rpacket_set_nu(&virt_packet,rpacket_get_nu (packet) * doppler_factor_ratio);
reabsorbed = montecarlo_one_packet_loop (storage, &virt_packet, 1);
reabsorbed = montecarlo_one_packet_loop (storage, &virt_packet, 1, mt_state);
if ((rpacket_get_nu(&virt_packet) < storage->spectrum_end_nu) &&
(rpacket_get_nu(&virt_packet) > storage->spectrum_start_nu))
{
Expand Down Expand Up @@ -481,7 +473,7 @@ montecarlo_one_packet (storage_model_t * storage, rpacket_t * packet,

void
move_packet_across_shell_boundary (rpacket_t * packet,
storage_model_t * storage, double distance)
storage_model_t * storage, double distance, rk_state *mt_state)
{
move_packet (packet, storage, distance);
if (rpacket_get_virtual_packet (packet) > 0)
Expand All @@ -493,7 +485,7 @@ move_packet_across_shell_boundary (rpacket_t * packet,
}
else
{
rpacket_reset_tau_event (packet);
rpacket_reset_tau_event (packet, mt_state);
}
if ((rpacket_get_current_shell_id (packet) < storage->no_of_shells - 1
&& rpacket_get_next_shell_id (packet) == 1)
Expand All @@ -512,7 +504,7 @@ move_packet_across_shell_boundary (rpacket_t * packet,
rpacket_set_status (packet, TARDIS_PACKET_STATUS_EMITTED);
}
else if ((storage->reflective_inner_boundary == 0) ||
(rk_double (&mt_state) > storage->inner_boundary_albedo))
(rk_double (mt_state) > storage->inner_boundary_albedo))
{
rpacket_set_status (packet, TARDIS_PACKET_STATUS_REABSORBED);
}
Expand All @@ -521,40 +513,40 @@ move_packet_across_shell_boundary (rpacket_t * packet,
double doppler_factor = rpacket_doppler_factor (packet, storage);
double comov_nu = rpacket_get_nu (packet) * doppler_factor;
double comov_energy = rpacket_get_energy (packet) * doppler_factor;
rpacket_set_mu (packet, rk_double (&mt_state));
rpacket_set_mu (packet, rk_double (mt_state));
double inverse_doppler_factor = 1.0 / rpacket_doppler_factor (packet, storage);
rpacket_set_nu (packet, comov_nu * inverse_doppler_factor);
rpacket_set_energy (packet, comov_energy * inverse_doppler_factor);
rpacket_set_recently_crossed_boundary (packet, 1);
if (rpacket_get_virtual_packet_flag (packet) > 0)
{
montecarlo_one_packet (storage, packet, -2);
montecarlo_one_packet (storage, packet, -2, mt_state);
}
}
}

void
montecarlo_thomson_scatter (rpacket_t * packet, storage_model_t * storage,
double distance)
double distance, rk_state *mt_state)
{
double doppler_factor = move_packet (packet, storage, distance);
double comov_nu = rpacket_get_nu (packet) * doppler_factor;
double comov_energy = rpacket_get_energy (packet) * doppler_factor;
rpacket_set_mu (packet, 2.0 * rk_double (&mt_state) - 1.0);
rpacket_set_mu (packet, 2.0 * rk_double (mt_state) - 1.0);
double inverse_doppler_factor = 1.0 / rpacket_doppler_factor (packet, storage);
rpacket_set_nu (packet, comov_nu * inverse_doppler_factor);
rpacket_set_energy (packet, comov_energy * inverse_doppler_factor);
rpacket_reset_tau_event (packet);
rpacket_reset_tau_event (packet, mt_state);
rpacket_set_recently_crossed_boundary (packet, 0);
storage->last_interaction_type[rpacket_get_id (packet)] = 1;
if (rpacket_get_virtual_packet_flag (packet) > 0)
{
montecarlo_one_packet (storage, packet, 1);
montecarlo_one_packet (storage, packet, 1, mt_state);
}
}

void
montecarlo_bound_free_scatter (rpacket_t * packet, storage_model_t * storage, double distance)
montecarlo_bound_free_scatter (rpacket_t * packet, storage_model_t * storage, double distance, rk_state *mt_state)
{
/* current position in list of continuum edges -> indicates which bound-free processes are possible */
int64_t current_continuum_id = rpacket_get_current_continuum_id(packet);
Expand All @@ -563,7 +555,7 @@ montecarlo_bound_free_scatter (rpacket_t * packet, storage_model_t * storage, do
double nu = rpacket_get_nu(packet);
double chi_bf = rpacket_get_chi_boundfree(packet);
// get new zrand
double zrand = rk_double(&mt_state);
double zrand = rk_double(mt_state);
double zrand_x_chibf = zrand * chi_bf;

int64_t ccontinuum = current_continuum_id; /* continuum_id of the continuum in which bf-absorption occurs */
Expand All @@ -580,7 +572,7 @@ montecarlo_bound_free_scatter (rpacket_t * packet, storage_model_t * storage, do
// ccontinuum = current_continuum_id;
// }

zrand = rk_double(&mt_state);
zrand = rk_double(mt_state);
if (zrand < storage->continuum_list_nu[ccontinuum] / nu)
{
// go to ionization energy
Expand All @@ -595,15 +587,15 @@ montecarlo_bound_free_scatter (rpacket_t * packet, storage_model_t * storage, do
}

void
montecarlo_free_free_scatter(rpacket_t * packet, storage_model_t * storage, double distance)
montecarlo_free_free_scatter(rpacket_t * packet, storage_model_t * storage, double distance, rk_state *mt_state)
{
rpacket_set_status (packet, TARDIS_PACKET_STATUS_REABSORBED);
}


void
montecarlo_line_scatter (rpacket_t * packet, storage_model_t * storage,
double distance)
double distance, rk_state *mt_state)
{
int64_t line2d_idx = rpacket_get_next_line_id (packet)
* storage->no_of_shells + rpacket_get_current_shell_id (packet);
Expand All @@ -629,7 +621,7 @@ montecarlo_line_scatter (rpacket_t * packet, storage_model_t * storage,
else if (rpacket_get_tau_event (packet) < tau_combined)
{
double old_doppler_factor = move_packet (packet, storage, distance);
rpacket_set_mu (packet, 2.0 * rk_double (&mt_state) - 1.0);
rpacket_set_mu (packet, 2.0 * rk_double (mt_state) - 1.0);
double inverse_doppler_factor = 1.0 / rpacket_doppler_factor (packet, storage);
double comov_energy = rpacket_get_energy (packet) * old_doppler_factor;
rpacket_set_energy (packet, comov_energy * inverse_doppler_factor);
Expand All @@ -647,7 +639,7 @@ montecarlo_line_scatter (rpacket_t * packet, storage_model_t * storage,
}
else if (storage->line_interaction_id >= 1)
{
emission_line_id = macro_atom (packet, storage);
emission_line_id = macro_atom (packet, storage, mt_state);
}
storage->last_line_interaction_out_id[rpacket_get_id (packet)] =
emission_line_id;
Expand All @@ -656,7 +648,7 @@ montecarlo_line_scatter (rpacket_t * packet, storage_model_t * storage,
inverse_doppler_factor);
rpacket_set_nu_line (packet, storage->line_list_nu[emission_line_id]);
rpacket_set_next_line_id (packet, emission_line_id + 1);
rpacket_reset_tau_event (packet);
rpacket_reset_tau_event (packet, mt_state);
rpacket_set_recently_crossed_boundary (packet, 0);
if (rpacket_get_virtual_packet_flag (packet) > 0)
{
Expand All @@ -671,7 +663,7 @@ montecarlo_line_scatter (rpacket_t * packet, storage_model_t * storage,
// QUESTIONABLE!!!
bool old_close_line = rpacket_get_close_line (packet);
rpacket_set_close_line (packet, virtual_close_line);
montecarlo_one_packet (storage, packet, 1);
montecarlo_one_packet (storage, packet, 1, mt_state);
rpacket_set_close_line (packet, old_close_line);
virtual_close_line = false;
}
Expand Down Expand Up @@ -715,7 +707,7 @@ montecarlo_compute_distances (rpacket_t * packet, storage_model_t * storage)

static montecarlo_event_handler_t
get_event_handler (rpacket_t * packet, storage_model_t * storage,
double *distance)
double *distance, rk_state *mt_state)
{
montecarlo_compute_distances (packet, storage);
double d_boundary = rpacket_get_d_boundary (packet);
Expand All @@ -735,21 +727,21 @@ get_event_handler (rpacket_t * packet, storage_model_t * storage,
else
{
*distance = d_continuum;
handler = montecarlo_continuum_event_handler(packet, storage);
handler = montecarlo_continuum_event_handler(packet, storage, mt_state);
}
return handler;
}

montecarlo_event_handler_t
montecarlo_continuum_event_handler(rpacket_t * packet, storage_model_t * storage)
montecarlo_continuum_event_handler(rpacket_t * packet, storage_model_t * storage, rk_state *mt_state)
{
if (storage->cont_status == CONTINUUM_OFF)
{
return &montecarlo_thomson_scatter;
}
else
{
double zrand = (rk_double(&mt_state));
double zrand = (rk_double(mt_state));
double normaliz_cont_th = rpacket_get_chi_electron(packet)/rpacket_get_chi_continuum(packet);
double normaliz_cont_bf = rpacket_get_chi_boundfree(packet)/rpacket_get_chi_continuum(packet);

Expand All @@ -773,7 +765,7 @@ montecarlo_continuum_event_handler(rpacket_t * packet, storage_model_t * storage

int64_t
montecarlo_one_packet_loop (storage_model_t * storage, rpacket_t * packet,
int64_t virtual_packet)
int64_t virtual_packet, rk_state *mt_state)
{
rpacket_set_tau_event (packet, 0.0);
rpacket_set_nu_line (packet, 0.0);
Expand All @@ -782,7 +774,7 @@ montecarlo_one_packet_loop (storage_model_t * storage, rpacket_t * packet,
// Initializing tau_event if it's a real packet.
if (virtual_packet == 0)
{
rpacket_reset_tau_event (packet);
rpacket_reset_tau_event (packet,mt_state);
}
// For a virtual packet tau_event is the sum of all the tau's that the packet passes.
while (rpacket_get_status (packet) == TARDIS_PACKET_STATUS_IN_PROCESS)
Expand All @@ -796,8 +788,8 @@ montecarlo_one_packet_loop (storage_model_t * storage, rpacket_t * packet,
(packet)]);
}
double distance;
get_event_handler (packet, storage, &distance) (packet, storage,
distance);
get_event_handler (packet, storage, &distance, mt_state) (packet, storage,
distance, mt_state);
if (virtual_packet > 0 && rpacket_get_tau_event (packet) > 10.0)
{
rpacket_set_tau_event (packet, 100.0);
Expand Down Expand Up @@ -832,11 +824,14 @@ montecarlo_main_loop(storage_model_t * storage, int64_t virtual_packet_flag, int
omp_set_num_threads(nthreads);
#pragma omp parallel
{
initialize_random_kit(seed + omp_get_thread_num());
rk_state mt_state;
rk_seed (seed + omp_get_thread_num(), &mt_state);

#pragma omp for
#else
fprintf(stderr, "Running without OpenMP\n");
initialize_random_kit(seed);
rk_state mt_state;
rk_seed (seed, &mt_state);
#endif
for (int64_t packet_index = 0; packet_index < storage->no_of_packets; packet_index++)
{
Expand All @@ -846,9 +841,9 @@ montecarlo_main_loop(storage_model_t * storage, int64_t virtual_packet_flag, int
rpacket_init(&packet, storage, packet_index, virtual_packet_flag);
if (virtual_packet_flag > 0)
{
reabsorbed = montecarlo_one_packet(storage, &packet, -1);
reabsorbed = montecarlo_one_packet(storage, &packet, -1, &mt_state);
}
reabsorbed = montecarlo_one_packet(storage, &packet, 0);
reabsorbed = montecarlo_one_packet(storage, &packet, 0, &mt_state);
storage->output_nus[packet_index] = rpacket_get_nu(&packet);
if (reabsorbed == 1)
{
Expand Down
22 changes: 12 additions & 10 deletions tardis/montecarlo/src/cmontecarlo.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,9 @@

typedef void (*montecarlo_event_handler_t) (rpacket_t * packet,
storage_model_t * storage,
double distance);
double distance, rk_state *mt_state);

void initialize_random_kit (unsigned long seed);

double rpacket_doppler_factor(const rpacket_t *packet, const storage_model_t *storage);

Expand Down Expand Up @@ -49,7 +51,7 @@ tardis_error_t compute_distance2line (const rpacket_t * packet,
*/
void compute_distance2continuum (rpacket_t * packet, storage_model_t * storage);

int64_t macro_atom (const rpacket_t * packet, const storage_model_t * storage);
int64_t macro_atom (const rpacket_t * packet, const storage_model_t * storage, rk_state *mt_state);

double move_packet (rpacket_t * packet, storage_model_t * storage,
double distance);
Expand All @@ -59,11 +61,11 @@ void increment_j_blue_estimator (const rpacket_t * packet,
double d_line, int64_t j_blue_idx);

int64_t montecarlo_one_packet (storage_model_t * storage, rpacket_t * packet,
int64_t virtual_mode);
int64_t virtual_mode, rk_state *mt_state);

int64_t montecarlo_one_packet_loop (storage_model_t * storage,
rpacket_t * packet,
int64_t virtual_packet);
int64_t virtual_packet, rk_state *mt_state);

void montecarlo_main_loop(storage_model_t * storage,
int64_t virtual_packet_flag,
Expand All @@ -72,11 +74,11 @@ void montecarlo_main_loop(storage_model_t * storage,

/* New handlers for continuum implementation */

montecarlo_event_handler_t montecarlo_continuum_event_handler(rpacket_t * packet, storage_model_t * storage);
montecarlo_event_handler_t montecarlo_continuum_event_handler(rpacket_t * packet, storage_model_t * storage, rk_state *mt_state);

void montecarlo_free_free_scatter (rpacket_t * packet, storage_model_t * storage, double distance);
void montecarlo_free_free_scatter (rpacket_t * packet, storage_model_t * storage, double distance, rk_state *mt_state);

void montecarlo_bound_free_scatter (rpacket_t * packet, storage_model_t * storage, double distance);
void montecarlo_bound_free_scatter (rpacket_t * packet, storage_model_t * storage, double distance, rk_state *mt_state);

double
bf_cross_section(const storage_model_t * storage, int64_t continuum_id, double comov_nu);
Expand All @@ -85,14 +87,14 @@ void calculate_chi_bf(rpacket_t * packet, storage_model_t * storage);

void
move_packet_across_shell_boundary (rpacket_t * packet,
storage_model_t * storage, double distance);
storage_model_t * storage, double distance, rk_state *mt_state);

void
montecarlo_thomson_scatter (rpacket_t * packet, storage_model_t * storage,
double distance);
double distance, rk_state *mt_state);

void
montecarlo_line_scatter (rpacket_t * packet, storage_model_t * storage,
double distance);
double distance, rk_state *mt_state);

#endif // TARDIS_CMONTECARLO_H
2 changes: 0 additions & 2 deletions tardis/montecarlo/src/cmontecarlo1.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,4 @@
tardis_error_t line_search (const double *nu, double nu_insert,
int64_t number_of_lines, int64_t * result);

extern rk_state mt_state;

#endif
4 changes: 2 additions & 2 deletions tardis/montecarlo/src/rpacket.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -264,9 +264,9 @@ static inline void rpacket_set_id (rpacket_t * packet, int id)
packet->id = id;
}

static inline void rpacket_reset_tau_event (rpacket_t * packet)
static inline void rpacket_reset_tau_event (rpacket_t * packet, rk_state *mt_state)
{
rpacket_set_tau_event (packet, -log (rk_double (&mt_state)));
rpacket_set_tau_event (packet, -log (rk_double (mt_state)));
}

tardis_error_t rpacket_init (rpacket_t * packet, storage_model_t * storage,
Expand Down
Loading