diff --git a/tardis/montecarlo/src/cmontecarlo.c b/tardis/montecarlo/src/cmontecarlo.c index 750878e74a8..96fc89b377b 100644 --- a/tardis/montecarlo/src/cmontecarlo.c +++ b/tardis/montecarlo/src/cmontecarlo.c @@ -26,7 +26,7 @@ reverse_binary_search (const double *x, double x_insert, Have in mind that *x points to a reverse sorted array. That is large values will have small indices and small ones will have large indices. - */ + */ tardis_error_t ret_val = TARDIS_ERROR_OK; if (x_insert > x[imin] || x_insert < x[imax]) { @@ -34,34 +34,34 @@ reverse_binary_search (const double *x, double x_insert, } else { - int imid = (imin + imax)>>1; + int imid = (imin + imax) >> 1; while (imax - imin > 2) - { - if (x[imid] < x_insert) - { - imax = imid + 1; - } - else - { - imin = imid; - } - imid = (imin + imax)>>1; - } + { + if (x[imid] < x_insert) + { + imax = imid + 1; + } + else + { + imin = imid; + } + imid = (imin + imax) >> 1; + } if (imax - imin == 2 && x_insert < x[imin + 1]) - { - *result = imin + 1; - } + { + *result = imin + 1; + } else - { - *result = imin; - } + { + *result = imin; + } } return ret_val; } tardis_error_t line_search (const double *nu, double nu_insert, int64_t number_of_lines, - int64_t * result) + int64_t * result) { tardis_error_t ret_val = TARDIS_ERROR_OK; int64_t imin = 0; @@ -95,7 +95,7 @@ rpacket_doppler_factor (const rpacket_t *packet, const storage_model_t *storage) double bf_cross_section(const storage_model_t * storage, int64_t continuum_id, double comov_nu) { -// FIXME MR: this seems like it should not be used in production! + // FIXME MR: this seems like it should not be used in production! /* Temporary hardcoded values */ #define chi_bf_partial 0.25e-15 static const double cont_chi_bf[] = {chi_bf_partial, 0.0, 2.0 * chi_bf_partial, 0.3 * chi_bf_partial, 2.0 * chi_bf_partial}; @@ -123,16 +123,16 @@ void calculate_chi_bf(rpacket_t * packet, storage_model_t * storage) double bf_helper = 0; for(int64_t i = current_continuum_id; i < no_of_continuum_edges; i++) - { - // get the levelpopulation for the level ijk in the current shell: - double l_pop = storage->l_pop[shell_id * no_of_continuum_edges + i]; - // get the levelpopulation ratio \frac{n_{0,j+1,k}}{n_{i,j,k}} \frac{n_{i,j,k}}{n_{0,j+1,k}}^{*}: - double l_pop_r = storage->l_pop_r[shell_id * no_of_continuum_edges + i]; - bf_helper += l_pop * bf_cross_section(storage, i, comov_nu) * (1 - l_pop_r * boltzmann_factor); + { + // get the levelpopulation for the level ijk in the current shell: + double l_pop = storage->l_pop[shell_id * no_of_continuum_edges + i]; + // get the levelpopulation ratio \frac{n_{0,j+1,k}}{n_{i,j,k}} \frac{n_{i,j,k}}{n_{0,j+1,k}}^{*}: + double l_pop_r = storage->l_pop_r[shell_id * no_of_continuum_edges + i]; + bf_helper += l_pop * bf_cross_section(storage, i, comov_nu) * (1 - l_pop_r * boltzmann_factor); -// FIXME MR: Is this thread-safe? It doesn't look like it to me ... - storage->chi_bf_tmp_partial[i] = bf_helper; - } + // FIXME MR: Is this thread-safe? It doesn't look like it to me ... + storage->chi_bf_tmp_partial[i] = bf_helper; + } rpacket_set_chi_boundfree(packet, bf_helper * doppler_factor); } @@ -156,14 +156,14 @@ compute_distance2boundary (rpacket_t * packet, const storage_model_t * storage) { 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; - } + { + rpacket_set_next_shell_id (packet, 1); + return d_outer; + } else - { - d_inner = mu < 0.0 ? -r * mu - sqrt (check) : MISS_DISTANCE; - } + { + d_inner = mu < 0.0 ? -r * mu - sqrt (check) : MISS_DISTANCE; + } } if (d_inner < d_outer) { @@ -179,7 +179,7 @@ compute_distance2boundary (rpacket_t * packet, const storage_model_t * storage) tardis_error_t compute_distance2line (const rpacket_t * packet, const storage_model_t * storage, - double *result) + double *result) { tardis_error_t ret_val = TARDIS_ERROR_OK; if (rpacket_get_last_line (packet)) @@ -198,46 +198,46 @@ compute_distance2line (const rpacket_t * packet, const storage_model_t * storage double doppler_factor = 1.0 - mu * r * inverse_t_exp * INVERSE_C; double comov_nu = nu * doppler_factor; if (comov_nu < nu_line) - { - if (rpacket_get_next_line_id (packet) == storage->no_of_lines - 1) - { - fprintf (stderr, "last_line = %f\n", - storage-> - line_list_nu[rpacket_get_next_line_id (packet) - 1]); - fprintf (stderr, "Last line in line list reached!"); - } - else if (rpacket_get_next_line_id (packet) == 0) - { - fprintf (stderr, "First line in line list!"); - fprintf (stderr, "next_line = %f\n", - storage-> - line_list_nu[rpacket_get_next_line_id (packet) + 1]); - } - else - { - fprintf (stderr, "last_line = %f\n", - storage-> - line_list_nu[rpacket_get_next_line_id (packet) - 1]); - fprintf (stderr, "next_line = %f\n", - storage-> - line_list_nu[rpacket_get_next_line_id (packet) + 1]); - } - fprintf (stderr, "ERROR: Comoving nu less than nu_line!\n"); - fprintf (stderr, "comov_nu = %f\n", comov_nu); - fprintf (stderr, "nu_line = %f\n", nu_line); - fprintf (stderr, "(comov_nu - nu_line) / nu_line = %f\n", - (comov_nu - nu_line) / nu_line); - fprintf (stderr, "r = %f\n", r); - fprintf (stderr, "mu = %f\n", mu); - 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; - } + { + if (rpacket_get_next_line_id (packet) == storage->no_of_lines - 1) + { + fprintf (stderr, "last_line = %f\n", + storage-> + line_list_nu[rpacket_get_next_line_id (packet) - 1]); + fprintf (stderr, "Last line in line list reached!"); + } + else if (rpacket_get_next_line_id (packet) == 0) + { + fprintf (stderr, "First line in line list!"); + fprintf (stderr, "next_line = %f\n", + storage-> + line_list_nu[rpacket_get_next_line_id (packet) + 1]); + } + else + { + fprintf (stderr, "last_line = %f\n", + storage-> + line_list_nu[rpacket_get_next_line_id (packet) - 1]); + fprintf (stderr, "next_line = %f\n", + storage-> + line_list_nu[rpacket_get_next_line_id (packet) + 1]); + } + fprintf (stderr, "ERROR: Comoving nu less than nu_line!\n"); + fprintf (stderr, "comov_nu = %f\n", comov_nu); + fprintf (stderr, "nu_line = %f\n", nu_line); + fprintf (stderr, "(comov_nu - nu_line) / nu_line = %f\n", + (comov_nu - nu_line) / nu_line); + fprintf (stderr, "r = %f\n", r); + fprintf (stderr, "mu = %f\n", mu); + 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; - } + { + *result = ((comov_nu - nu_line) / nu) * C * t_exp; + } } return ret_val; } @@ -248,49 +248,49 @@ compute_distance2continuum(rpacket_t * packet, storage_model_t * storage) double chi_freefree, chi_electron, chi_continuum, d_continuum; if (storage->cont_status == CONTINUUM_ON) - { - calculate_chi_bf(packet, storage); - double chi_boundfree = rpacket_get_chi_boundfree(packet); - rpacket_set_chi_freefree(packet, 0.0); - chi_freefree = rpacket_get_chi_freefree(packet); // MR ?? this is always zero - chi_electron = storage->electron_densities[rpacket_get_current_shell_id(packet)] * storage->sigma_thomson * - rpacket_doppler_factor (packet, storage); - chi_continuum = chi_boundfree + chi_freefree + chi_electron; - d_continuum = rpacket_get_tau_event(packet) / chi_continuum; - } + { + calculate_chi_bf(packet, storage); + double chi_boundfree = rpacket_get_chi_boundfree(packet); + rpacket_set_chi_freefree(packet, 0.0); + chi_freefree = rpacket_get_chi_freefree(packet); // MR ?? this is always zero + chi_electron = storage->electron_densities[rpacket_get_current_shell_id(packet)] * storage->sigma_thomson * + rpacket_doppler_factor (packet, storage); + chi_continuum = chi_boundfree + chi_freefree + chi_electron; + d_continuum = rpacket_get_tau_event(packet) / chi_continuum; + } else - { -// FIXME MR: an assignment to chi_freefree seems to be missing here - chi_electron = storage->electron_densities[rpacket_get_current_shell_id(packet)] * storage->sigma_thomson; - chi_continuum = chi_electron; - d_continuum = storage->inverse_electron_densities[rpacket_get_current_shell_id (packet)] * - storage->inverse_sigma_thomson * rpacket_get_tau_event (packet); - } + { + // FIXME MR: an assignment to chi_freefree seems to be missing here + chi_electron = storage->electron_densities[rpacket_get_current_shell_id(packet)] * storage->sigma_thomson; + chi_continuum = chi_electron; + d_continuum = storage->inverse_electron_densities[rpacket_get_current_shell_id (packet)] * + storage->inverse_sigma_thomson * rpacket_get_tau_event (packet); + } if (rpacket_get_virtual_packet(packet) > 0) { - //Set all continuum distances to MISS_DISTANCE in case of an virtual_packet - rpacket_set_d_continuum(packet, MISS_DISTANCE); - rpacket_set_chi_boundfree(packet, 0.0); - rpacket_set_chi_electron(packet, chi_electron); + //Set all continuum distances to MISS_DISTANCE in case of an virtual_packet + rpacket_set_d_continuum(packet, MISS_DISTANCE); + rpacket_set_chi_boundfree(packet, 0.0); + rpacket_set_chi_electron(packet, chi_electron); rpacket_set_chi_freefree(packet, 0.0); rpacket_set_chi_continuum(packet, chi_continuum); - } - else - { - -// fprintf(stderr, "--------\n"); -// fprintf(stderr, "nu = %e \n", rpacket_get_nu(packet)); -// fprintf(stderr, "chi_electron = %e\n", chi_electron); -// fprintf(stderr, "chi_boundfree = %e\n", calculate_chi_bf(packet, storage)); -// fprintf(stderr, "chi_line = %e \n", rpacket_get_tau_event(packet) / rpacket_get_d_line(packet)); -// fprintf(stderr, "--------\n"); - - rpacket_set_chi_freefree(packet, chi_freefree); - rpacket_set_chi_electron(packet, chi_electron); - rpacket_set_chi_continuum(packet, chi_continuum); - rpacket_set_d_continuum(packet, d_continuum); - } + } + else + { + + // fprintf(stderr, "--------\n"); + // fprintf(stderr, "nu = %e \n", rpacket_get_nu(packet)); + // fprintf(stderr, "chi_electron = %e\n", chi_electron); + // fprintf(stderr, "chi_boundfree = %e\n", calculate_chi_bf(packet, storage)); + // fprintf(stderr, "chi_line = %e \n", rpacket_get_tau_event(packet) / rpacket_get_d_line(packet)); + // fprintf(stderr, "--------\n"); + + rpacket_set_chi_freefree(packet, chi_freefree); + rpacket_set_chi_electron(packet, chi_electron); + rpacket_set_chi_continuum(packet, chi_continuum); + rpacket_set_d_continuum(packet, d_continuum); + } } int64_t @@ -305,12 +305,12 @@ macro_atom (const rpacket_t * packet, const storage_model_t * storage, rk_state i = storage->macro_block_references[activate_level] - 1; double p = 0.0; do - { + { - probability_idx = ((++i) * storage->no_of_shells + - rpacket_get_current_shell_id (packet)); - p += storage->transition_probabilities[probability_idx]; - } + probability_idx = ((++i) * storage->no_of_shells + + rpacket_get_current_shell_id (packet)); + p += storage->transition_probabilities[probability_idx]; + } while (p <= event_random); emit = storage->transition_type[i]; activate_level = storage->destination_level_id[i]; @@ -326,38 +326,38 @@ move_packet (rpacket_t * packet, storage_model_t * storage, double distance) { double r = rpacket_get_r (packet); double new_r = - sqrt (r * r + distance * distance + - 2.0 * r * distance * rpacket_get_mu (packet)); + sqrt (r * r + distance * distance + + 2.0 * r * distance * rpacket_get_mu (packet)); rpacket_set_mu (packet, - (rpacket_get_mu (packet) * r + distance) / new_r); + (rpacket_get_mu (packet) * r + distance) / new_r); rpacket_set_r (packet, new_r); if (rpacket_get_virtual_packet (packet) <= 0) - { - double comov_energy = rpacket_get_energy (packet) * doppler_factor; - double comov_nu = rpacket_get_nu (packet) * doppler_factor; + { + double comov_energy = rpacket_get_energy (packet) * doppler_factor; + double comov_nu = rpacket_get_nu (packet) * doppler_factor; #ifdef WITHOPENMP #pragma omp atomic #endif - storage->js[rpacket_get_current_shell_id (packet)] += - comov_energy * distance; + storage->js[rpacket_get_current_shell_id (packet)] += + comov_energy * distance; #ifdef WITHOPENMP #pragma omp atomic #endif - storage->nubars[rpacket_get_current_shell_id (packet)] += - comov_energy * distance * comov_nu; - } + storage->nubars[rpacket_get_current_shell_id (packet)] += + comov_energy * distance * comov_nu; + } } return doppler_factor; } void increment_j_blue_estimator (const rpacket_t * packet, storage_model_t * storage, - double d_line, int64_t j_blue_idx) + double d_line, int64_t j_blue_idx) { double r = rpacket_get_r (packet); double r_interaction = sqrt (r * r + d_line * d_line + - 2.0 * r * d_line * rpacket_get_mu (packet)); + 2.0 * r * d_line * rpacket_get_mu (packet)); double mu_interaction = (rpacket_get_mu (packet) * r + d_line) / r_interaction; double doppler_factor = 1.0 - mu_interaction * r_interaction * storage->inverse_time_explosion * INVERSE_C; @@ -371,7 +371,7 @@ 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, rk_state *mt_state) + int64_t virtual_mode, rk_state *mt_state) { int64_t reabsorbed=-1; if (virtual_mode == 0) @@ -381,117 +381,117 @@ montecarlo_one_packet (storage_model_t * storage, rpacket_t * packet, else { if ((rpacket_get_nu (packet) > storage->spectrum_virt_start_nu) && (rpacket_get_nu(packet) < storage->spectrum_virt_end_nu)) - { - for (int64_t i = 0; i < rpacket_get_virtual_packet_flag (packet); i++) - { + { + for (int64_t i = 0; i < rpacket_get_virtual_packet_flag (packet); i++) + { double weight; rpacket_t virt_packet = *packet; double mu_min; - if (rpacket_get_r(&virt_packet) > storage->r_inner[0]) - { - mu_min = - -1.0 * sqrt (1.0 - - (storage->r_inner[0] / rpacket_get_r(&virt_packet)) * - (storage->r_inner[0] / rpacket_get_r(&virt_packet))); - } - else - { - 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); - switch (virtual_mode) - { - case -2: - weight = 1.0 / rpacket_get_virtual_packet_flag (packet); - break; - case -1: - weight = - 2.0 * rpacket_get_mu(&virt_packet) / - rpacket_get_virtual_packet_flag (packet); - break; - case 1: - weight = - (1.0 - - mu_min) / 2.0 / rpacket_get_virtual_packet_flag (packet); - break; - default: - fprintf (stderr, "Something has gone horribly wrong!\n"); + if (rpacket_get_r(&virt_packet) > storage->r_inner[0]) + { + mu_min = + -1.0 * sqrt (1.0 - + (storage->r_inner[0] / rpacket_get_r(&virt_packet)) * + (storage->r_inner[0] / rpacket_get_r(&virt_packet))); + } + else + { + 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); + switch (virtual_mode) + { + case -2: + weight = 1.0 / rpacket_get_virtual_packet_flag (packet); + break; + case -1: + weight = + 2.0 * rpacket_get_mu(&virt_packet) / + rpacket_get_virtual_packet_flag (packet); + break; + case 1: + weight = + (1.0 - + mu_min) / 2.0 / rpacket_get_virtual_packet_flag (packet); + break; + default: + fprintf (stderr, "Something has gone horribly wrong!\n"); // FIXME MR: we need to somehow signal an error here // I'm adding an exit() here to inform the compiler about the impossible path exit(1); - } - double doppler_factor_ratio = - rpacket_doppler_factor (packet, storage) / - rpacket_doppler_factor (&virt_packet, storage); - 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, mt_state); + } + double doppler_factor_ratio = + rpacket_doppler_factor (packet, storage) / + rpacket_doppler_factor (&virt_packet, storage); + 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, mt_state); #ifdef WITH_VPACKET_LOGGING #ifdef WITHOPENMP #pragma omp critical - { + { #endif // WITHOPENMP - if (storage->virt_packet_count >= storage->virt_array_size) - { - storage->virt_array_size *= 2; - storage->virt_packet_nus = safe_realloc(storage->virt_packet_nus, sizeof(double) * storage->virt_array_size); - storage->virt_packet_energies = safe_realloc(storage->virt_packet_energies, sizeof(double) * storage->virt_array_size); - storage->virt_packet_last_interaction_in_nu = safe_realloc(storage->virt_packet_last_interaction_in_nu, sizeof(double) * storage->virt_array_size); - storage->virt_packet_last_interaction_type = safe_realloc(storage->virt_packet_last_interaction_type, sizeof(int64_t) * storage->virt_array_size); - storage->virt_packet_last_line_interaction_in_id = safe_realloc(storage->virt_packet_last_line_interaction_in_id, sizeof(int64_t) * storage->virt_array_size); - storage->virt_packet_last_line_interaction_out_id = safe_realloc(storage->virt_packet_last_line_interaction_out_id, sizeof(int64_t) * storage->virt_array_size); - } - storage->virt_packet_nus[storage->virt_packet_count] = rpacket_get_nu(&virt_packet); - storage->virt_packet_energies[storage->virt_packet_count] = rpacket_get_energy(&virt_packet) * weight; - storage->virt_packet_last_interaction_in_nu[storage->virt_packet_count] = storage->last_interaction_in_nu[rpacket_get_id (packet)]; - storage->virt_packet_last_interaction_type[storage->virt_packet_count] = storage->last_interaction_type[rpacket_get_id (packet)]; - storage->virt_packet_last_line_interaction_in_id[storage->virt_packet_count] = storage->last_line_interaction_in_id[rpacket_get_id (packet)]; - storage->virt_packet_last_line_interaction_out_id[storage->virt_packet_count] = storage->last_line_interaction_out_id[rpacket_get_id (packet)]; - storage->virt_packet_count += 1; + if (storage->virt_packet_count >= storage->virt_array_size) + { + storage->virt_array_size *= 2; + storage->virt_packet_nus = safe_realloc(storage->virt_packet_nus, sizeof(double) * storage->virt_array_size); + storage->virt_packet_energies = safe_realloc(storage->virt_packet_energies, sizeof(double) * storage->virt_array_size); + storage->virt_packet_last_interaction_in_nu = safe_realloc(storage->virt_packet_last_interaction_in_nu, sizeof(double) * storage->virt_array_size); + storage->virt_packet_last_interaction_type = safe_realloc(storage->virt_packet_last_interaction_type, sizeof(int64_t) * storage->virt_array_size); + storage->virt_packet_last_line_interaction_in_id = safe_realloc(storage->virt_packet_last_line_interaction_in_id, sizeof(int64_t) * storage->virt_array_size); + storage->virt_packet_last_line_interaction_out_id = safe_realloc(storage->virt_packet_last_line_interaction_out_id, sizeof(int64_t) * storage->virt_array_size); + } + storage->virt_packet_nus[storage->virt_packet_count] = rpacket_get_nu(&virt_packet); + storage->virt_packet_energies[storage->virt_packet_count] = rpacket_get_energy(&virt_packet) * weight; + storage->virt_packet_last_interaction_in_nu[storage->virt_packet_count] = storage->last_interaction_in_nu[rpacket_get_id (packet)]; + storage->virt_packet_last_interaction_type[storage->virt_packet_count] = storage->last_interaction_type[rpacket_get_id (packet)]; + storage->virt_packet_last_line_interaction_in_id[storage->virt_packet_count] = storage->last_line_interaction_in_id[rpacket_get_id (packet)]; + storage->virt_packet_last_line_interaction_out_id[storage->virt_packet_count] = storage->last_line_interaction_out_id[rpacket_get_id (packet)]; + storage->virt_packet_count += 1; #ifdef WITHOPENMP - } + } #endif // WITHOPENMP #endif // WITH_VPACKET_LOGGING - if ((rpacket_get_nu(&virt_packet) < storage->spectrum_end_nu) && - (rpacket_get_nu(&virt_packet) > storage->spectrum_start_nu)) - { + if ((rpacket_get_nu(&virt_packet) < storage->spectrum_end_nu) && + (rpacket_get_nu(&virt_packet) > storage->spectrum_start_nu)) + { #ifdef WITHOPENMP #pragma omp critical - { + { #endif // WITHOPENMP - int64_t virt_id_nu = - floor ((rpacket_get_nu(&virt_packet) - - storage->spectrum_start_nu) / - storage->spectrum_delta_nu); - storage->spectrum_virt_nu[virt_id_nu] += - rpacket_get_energy(&virt_packet) * weight; + int64_t virt_id_nu = + floor ((rpacket_get_nu(&virt_packet) - + storage->spectrum_start_nu) / + storage->spectrum_delta_nu); + storage->spectrum_virt_nu[virt_id_nu] += + rpacket_get_energy(&virt_packet) * weight; #ifdef WITHOPENMP - } + } #endif // WITHOPENMP - } - } - } + } + } + } else - { - return 1; - } + { + return 1; + } } return reabsorbed; } void move_packet_across_shell_boundary (rpacket_t * packet, - storage_model_t * storage, double distance, rk_state *mt_state) + storage_model_t * storage, double distance, rk_state *mt_state) { move_packet (packet, storage, distance); if (rpacket_get_virtual_packet (packet) > 0) { double delta_tau_event = rpacket_get_chi_continuum(packet) * distance; rpacket_set_tau_event (packet, - rpacket_get_tau_event (packet) + - delta_tau_event); + rpacket_get_tau_event (packet) + + delta_tau_event); } else { @@ -500,21 +500,21 @@ move_packet_across_shell_boundary (rpacket_t * packet, if ((rpacket_get_current_shell_id (packet) < storage->no_of_shells - 1 && rpacket_get_next_shell_id (packet) == 1) || (rpacket_get_current_shell_id (packet) > 0 - && rpacket_get_next_shell_id (packet) == -1)) + && rpacket_get_next_shell_id (packet) == -1)) { rpacket_set_current_shell_id (packet, - rpacket_get_current_shell_id (packet) + - rpacket_get_next_shell_id (packet)); + rpacket_get_current_shell_id (packet) + + rpacket_get_next_shell_id (packet)); rpacket_set_recently_crossed_boundary (packet, - rpacket_get_next_shell_id - (packet)); + rpacket_get_next_shell_id + (packet)); } else if (rpacket_get_next_shell_id (packet) == 1) { 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); } @@ -529,15 +529,15 @@ move_packet_across_shell_boundary (rpacket_t * packet, 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, mt_state); - } + { + montecarlo_one_packet (storage, packet, -2, mt_state); + } } } void montecarlo_thomson_scatter (rpacket_t * packet, storage_model_t * storage, - double distance, rk_state *mt_state) + double distance, rk_state *mt_state) { double doppler_factor = move_packet (packet, storage, distance); double comov_nu = rpacket_get_nu (packet) * doppler_factor; @@ -571,29 +571,29 @@ montecarlo_bound_free_scatter (rpacket_t * packet, storage_model_t * storage, do int64_t ccontinuum = current_continuum_id; /* continuum_id of the continuum in which bf-absorption occurs */ while (storage->chi_bf_tmp_partial[ccontinuum] <= zrand_x_chibf) - { - ccontinuum++; - } -// Alternative way to choose a continuum for bf-absorption: -// error = -// binary_search(storage->chi_bf_tmp_partial, zrand_x_chibf, current_continuum_id,no_of_continuum_edges-1,&ccontinuum); -// if (error == TARDIS_ERROR_BOUNDS_ERROR) // x_insert < x[imin] -> set index equal to imin -// { -// ccontinuum = current_continuum_id; -// } + { + ccontinuum++; + } + // Alternative way to choose a continuum for bf-absorption: + // error = + // binary_search(storage->chi_bf_tmp_partial, zrand_x_chibf, current_continuum_id,no_of_continuum_edges-1,&ccontinuum); + // if (error == TARDIS_ERROR_BOUNDS_ERROR) // x_insert < x[imin] -> set index equal to imin + // { + // ccontinuum = current_continuum_id; + // } zrand = rk_double(mt_state); if (zrand < storage->continuum_list_nu[ccontinuum] / nu) - { - // go to ionization energy - rpacket_set_status (packet, TARDIS_PACKET_STATUS_REABSORBED); - } + { + // go to ionization energy + rpacket_set_status (packet, TARDIS_PACKET_STATUS_REABSORBED); + } else - { - //go to the thermal pool - //create_kpacket(packet); - rpacket_set_status (packet, TARDIS_PACKET_STATUS_REABSORBED); - } + { + //go to the thermal pool + //create_kpacket(packet); + rpacket_set_status (packet, TARDIS_PACKET_STATUS_REABSORBED); + } } void @@ -605,10 +605,10 @@ montecarlo_free_free_scatter(rpacket_t * packet, storage_model_t * storage, doub void montecarlo_line_scatter (rpacket_t * packet, storage_model_t * storage, - double distance, rk_state *mt_state) + 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); + * storage->no_of_shells + rpacket_get_current_shell_id (packet); if (rpacket_get_virtual_packet (packet) == 0) { increment_j_blue_estimator (packet, storage, distance, line2d_idx); @@ -626,7 +626,7 @@ montecarlo_line_scatter (rpacket_t * packet, storage_model_t * storage, if (rpacket_get_virtual_packet (packet) > 0) { rpacket_set_tau_event (packet, - rpacket_get_tau_event (packet) + tau_line); + rpacket_get_tau_event (packet) + tau_line); } else if (rpacket_get_tau_event (packet) < tau_combined) { @@ -636,57 +636,57 @@ montecarlo_line_scatter (rpacket_t * packet, storage_model_t * storage, double comov_energy = rpacket_get_energy (packet) * old_doppler_factor; rpacket_set_energy (packet, comov_energy * inverse_doppler_factor); storage->last_interaction_in_nu[rpacket_get_id (packet)] = - rpacket_get_nu (packet); + rpacket_get_nu (packet); storage->last_line_interaction_in_id[rpacket_get_id (packet)] = - rpacket_get_next_line_id (packet) - 1; + rpacket_get_next_line_id (packet) - 1; storage->last_line_interaction_shell_id[rpacket_get_id (packet)] = - rpacket_get_current_shell_id (packet); + rpacket_get_current_shell_id (packet); storage->last_interaction_type[rpacket_get_id (packet)] = 2; int64_t emission_line_id = 0; if (storage->line_interaction_id == 0) - { - emission_line_id = rpacket_get_next_line_id (packet) - 1; - } + { + emission_line_id = rpacket_get_next_line_id (packet) - 1; + } else if (storage->line_interaction_id >= 1) - { - emission_line_id = macro_atom (packet, storage, mt_state); - } + { + emission_line_id = macro_atom (packet, storage, mt_state); + } storage->last_line_interaction_out_id[rpacket_get_id (packet)] = - emission_line_id; + emission_line_id; rpacket_set_nu (packet, - storage->line_list_nu[emission_line_id] * - inverse_doppler_factor); + storage->line_list_nu[emission_line_id] * + 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, mt_state); rpacket_set_recently_crossed_boundary (packet, 0); if (rpacket_get_virtual_packet_flag (packet) > 0) - { - bool virtual_close_line = false; - if (!rpacket_get_last_line (packet) && - fabs (storage->line_list_nu[rpacket_get_next_line_id (packet)] - - rpacket_get_nu_line (packet)) < - (rpacket_get_nu_line (packet)* 1e-7)) - { - virtual_close_line = true; - } - // QUESTIONABLE!!! - bool old_close_line = rpacket_get_close_line (packet); - rpacket_set_close_line (packet, virtual_close_line); - montecarlo_one_packet (storage, packet, 1, mt_state); - rpacket_set_close_line (packet, old_close_line); - virtual_close_line = false; - } + { + bool virtual_close_line = false; + if (!rpacket_get_last_line (packet) && + fabs (storage->line_list_nu[rpacket_get_next_line_id (packet)] - + rpacket_get_nu_line (packet)) < + (rpacket_get_nu_line (packet)* 1e-7)) + { + virtual_close_line = true; + } + // QUESTIONABLE!!! + bool old_close_line = rpacket_get_close_line (packet); + rpacket_set_close_line (packet, virtual_close_line); + montecarlo_one_packet (storage, packet, 1, mt_state); + rpacket_set_close_line (packet, old_close_line); + virtual_close_line = false; + } } else { rpacket_set_tau_event (packet, - rpacket_get_tau_event (packet) - tau_line); + rpacket_get_tau_event (packet) - tau_line); } if (!rpacket_get_last_line (packet) && fabs (storage->line_list_nu[rpacket_get_next_line_id (packet)] - - rpacket_get_nu_line (packet)) < (rpacket_get_nu_line (packet)* - 1e-7)) + rpacket_get_nu_line (packet)) < (rpacket_get_nu_line (packet)* + 1e-7)) { rpacket_set_close_line (packet, true); } @@ -706,7 +706,7 @@ montecarlo_compute_distances (rpacket_t * packet, storage_model_t * storage) else { rpacket_set_d_boundary (packet, - compute_distance2boundary (packet, storage)); + compute_distance2boundary (packet, storage)); double d_line; compute_distance2line (packet, storage, &d_line); // FIXME MR: return status of compute_distance2line() is ignored @@ -717,7 +717,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, rk_state *mt_state) + double *distance, rk_state *mt_state) { montecarlo_compute_distances (packet, storage); double d_boundary = rpacket_get_d_boundary (packet); @@ -751,31 +751,31 @@ montecarlo_continuum_event_handler(rpacket_t * packet, storage_model_t * storage } else { - 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); - - if (zrand < normaliz_cont_th) - { - //Return the electron scatter event function - return &montecarlo_thomson_scatter; - } - else if (zrand < (normaliz_cont_th + normaliz_cont_bf)) - { - //Return the bound-free scatter event function - return &montecarlo_bound_free_scatter; - } - else - { - //Return the free-free scatter event function - return &montecarlo_free_free_scatter; - } + 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); + + if (zrand < normaliz_cont_th) + { + //Return the electron scatter event function + return &montecarlo_thomson_scatter; + } + else if (zrand < (normaliz_cont_th + normaliz_cont_bf)) + { + //Return the bound-free scatter event function + return &montecarlo_bound_free_scatter; + } + else + { + //Return the free-free scatter event function + return &montecarlo_free_free_scatter; + } } } int64_t montecarlo_one_packet_loop (storage_model_t * storage, rpacket_t * packet, - int64_t virtual_packet, rk_state *mt_state) + int64_t virtual_packet, rk_state *mt_state) { rpacket_set_tau_event (packet, 0.0); rpacket_set_nu_line (packet, 0.0); @@ -791,27 +791,27 @@ montecarlo_one_packet_loop (storage_model_t * storage, rpacket_t * packet, { // Check if we are at the end of line list. if (!rpacket_get_last_line (packet)) - { - rpacket_set_nu_line (packet, - storage-> - line_list_nu[rpacket_get_next_line_id - (packet)]); - } + { + rpacket_set_nu_line (packet, + storage-> + line_list_nu[rpacket_get_next_line_id + (packet)]); + } double distance; get_event_handler (packet, storage, &distance, mt_state) (packet, storage, - distance, mt_state); + distance, mt_state); if (virtual_packet > 0 && rpacket_get_tau_event (packet) > 10.0) - { - rpacket_set_tau_event (packet, 100.0); - rpacket_set_status (packet, TARDIS_PACKET_STATUS_EMITTED); - } + { + rpacket_set_tau_event (packet, 100.0); + rpacket_set_status (packet, TARDIS_PACKET_STATUS_EMITTED); + } } if (virtual_packet > 0) { rpacket_set_energy (packet, - rpacket_get_energy (packet) * exp (-1.0 * - rpacket_get_tau_event - (packet))); + rpacket_get_energy (packet) * exp (-1.0 * + rpacket_get_tau_event + (packet))); } return rpacket_get_status (packet) == TARDIS_PACKET_STATUS_REABSORBED ? 1 : 0; @@ -836,51 +836,51 @@ montecarlo_main_loop(storage_model_t * storage, int64_t virtual_packet_flag, int omp_set_num_threads(nthreads); #pragma omp parallel firstprivate(finished_packets) - { - rk_state mt_state; - rk_seed (seed + omp_get_thread_num(), &mt_state); + { + rk_state mt_state; + rk_seed (seed + omp_get_thread_num(), &mt_state); #pragma omp master - fprintf(stderr, "Running with OpenMP - %d threads\n", omp_get_num_threads()); + fprintf(stderr, "Running with OpenMP - %d threads\n", omp_get_num_threads()); #pragma omp for #else - rk_state mt_state; - rk_seed (seed, &mt_state); - fprintf(stderr, "Running without OpenMP\n"); + rk_state mt_state; + rk_seed (seed, &mt_state); + fprintf(stderr, "Running without OpenMP\n"); #endif - for (int64_t packet_index = 0; packet_index < storage->no_of_packets; packet_index++) - { - ++finished_packets; - if ( finished_packets%100 == 0 ) { + for (int64_t packet_index = 0; packet_index < storage->no_of_packets; packet_index++) + { + ++finished_packets; + if ( finished_packets%100 == 0 ) { #ifdef WITHOPENMP - // WARNING: This only works with a static sheduler and gives an approximation of progress. - // The alternative would be to have a shared variable but that could potentially decrease performance when using many threads. - if (omp_get_thread_num() == 0 ) - print_progress(finished_packets * omp_get_num_threads(), storage->no_of_packets); + // WARNING: This only works with a static sheduler and gives an approximation of progress. + // The alternative would be to have a shared variable but that could potentially decrease performance when using many threads. + if (omp_get_thread_num() == 0 ) + print_progress(finished_packets * omp_get_num_threads(), storage->no_of_packets); #else - print_progress(finished_packets, storage->no_of_packets); + print_progress(finished_packets, storage->no_of_packets); #endif - } - int reabsorbed = 0; - rpacket_t packet; - rpacket_set_id(&packet, packet_index); - rpacket_init(&packet, storage, packet_index, virtual_packet_flag); - if (virtual_packet_flag > 0) - { - reabsorbed = montecarlo_one_packet(storage, &packet, -1, &mt_state); - } - reabsorbed = montecarlo_one_packet(storage, &packet, 0, &mt_state); - storage->output_nus[packet_index] = rpacket_get_nu(&packet); - if (reabsorbed == 1) - { - storage->output_energies[packet_index] = -rpacket_get_energy(&packet); - } - else - { - storage->output_energies[packet_index] = rpacket_get_energy(&packet); - } - } + } + int reabsorbed = 0; + rpacket_t packet; + rpacket_set_id(&packet, packet_index); + rpacket_init(&packet, storage, packet_index, virtual_packet_flag); + if (virtual_packet_flag > 0) + { + reabsorbed = montecarlo_one_packet(storage, &packet, -1, &mt_state); + } + reabsorbed = montecarlo_one_packet(storage, &packet, 0, &mt_state); + storage->output_nus[packet_index] = rpacket_get_nu(&packet); + if (reabsorbed == 1) + { + storage->output_energies[packet_index] = -rpacket_get_energy(&packet); + } + else + { + storage->output_energies[packet_index] = rpacket_get_energy(&packet); + } + } #ifdef WITHOPENMP - } + } #endif print_progress(storage->no_of_packets, storage->no_of_packets); fprintf(stderr,"\n"); diff --git a/tardis/montecarlo/src/cmontecarlo.h b/tardis/montecarlo/src/cmontecarlo.h index 0ff73e141d3..b2970345f54 100644 --- a/tardis/montecarlo/src/cmontecarlo.h +++ b/tardis/montecarlo/src/cmontecarlo.h @@ -19,8 +19,8 @@ #endif typedef void (*montecarlo_event_handler_t) (rpacket_t * packet, - storage_model_t * storage, - double distance, rk_state *mt_state); + storage_model_t * storage, + double distance, rk_state *mt_state); void initialize_random_kit (unsigned long seed); @@ -34,7 +34,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, - const storage_model_t * storage); + const storage_model_t * storage); /** Calculate the distance the packet has to travel until it redshifts to the first spectral line. * @@ -44,11 +44,11 @@ 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); + const storage_model_t * storage, + double *result); /** Calculate the distance to the next continuum event, which can be a Thomson scattering, bound-free absorption or - free-free transition. + free-free transition. * * @param packet rpacket structure with packet information * @param storage storage model data @@ -60,23 +60,23 @@ void compute_distance2continuum (rpacket_t * packet, 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); + double distance); void increment_j_blue_estimator (const rpacket_t * packet, - storage_model_t * storage, - double d_line, int64_t j_blue_idx); + storage_model_t * storage, + double d_line, int64_t j_blue_idx); int64_t montecarlo_one_packet (storage_model_t * storage, rpacket_t * packet, - int64_t virtual_mode, rk_state *mt_state); + 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, rk_state *mt_state); + rpacket_t * packet, + int64_t virtual_packet, rk_state *mt_state); void montecarlo_main_loop(storage_model_t * storage, - int64_t virtual_packet_flag, - int nthreads, - unsigned long seed); + int64_t virtual_packet_flag, + int nthreads, + unsigned long seed); /* New handlers for continuum implementation */ diff --git a/tardis/montecarlo/src/cmontecarlo1.h b/tardis/montecarlo/src/cmontecarlo1.h index 563b8a2e699..be027c81c53 100644 --- a/tardis/montecarlo/src/cmontecarlo1.h +++ b/tardis/montecarlo/src/cmontecarlo1.h @@ -10,6 +10,6 @@ * @return index of the next line ot the red. If the key value is redder than the reddest line returns number_of_lines. */ tardis_error_t line_search (const double *nu, double nu_insert, - int64_t number_of_lines, int64_t * result); + int64_t number_of_lines, int64_t * result); #endif diff --git a/tardis/montecarlo/src/rpacket.c b/tardis/montecarlo/src/rpacket.c index 99abdec5f57..62f6fe4b2be 100644 --- a/tardis/montecarlo/src/rpacket.c +++ b/tardis/montecarlo/src/rpacket.c @@ -3,7 +3,7 @@ tardis_error_t rpacket_init (rpacket_t * packet, storage_model_t * storage, int packet_index, - int virtual_packet_flag) + int virtual_packet_flag) { int64_t current_line_id; tardis_error_t ret_val = TARDIS_ERROR_OK; @@ -15,16 +15,16 @@ rpacket_init (rpacket_t * packet, storage_model_t * storage, int packet_index, double current_r = storage->r_inner[0]; current_nu = current_nu / (1 - - (current_mu * current_r * storage->inverse_time_explosion * - INVERSE_C)); + (current_mu * current_r * storage->inverse_time_explosion * + INVERSE_C)); current_energy = current_energy / (1 - - (current_mu * current_r * - storage->inverse_time_explosion * INVERSE_C)); + (current_mu * current_r * + storage->inverse_time_explosion * INVERSE_C)); if ((ret_val = line_search (storage->line_list_nu, comov_current_nu, - storage->no_of_lines, - ¤t_line_id)) != TARDIS_ERROR_OK) + storage->no_of_lines, + ¤t_line_id)) != TARDIS_ERROR_OK) { return ret_val; } diff --git a/tardis/montecarlo/src/rpacket.h b/tardis/montecarlo/src/rpacket.h index c43c164e9b6..290ca211dd7 100644 --- a/tardis/montecarlo/src/rpacket.h +++ b/tardis/montecarlo/src/rpacket.h @@ -135,7 +135,7 @@ static inline unsigned int rpacket_get_current_shell_id (const rpacket_t * packe } static inline void rpacket_set_current_shell_id (rpacket_t * packet, - unsigned int current_shell_id) + unsigned int current_shell_id) { packet->current_shell_id = current_shell_id; } @@ -146,7 +146,7 @@ static inline unsigned int rpacket_get_next_line_id (const rpacket_t * packet) } static inline void rpacket_set_next_line_id (rpacket_t * packet, - unsigned int next_line_id) + unsigned int next_line_id) { packet->next_line_id = next_line_id; } @@ -177,7 +177,7 @@ static inline int rpacket_get_recently_crossed_boundary (const rpacket_t * packe } static inline void rpacket_set_recently_crossed_boundary (rpacket_t * packet, - int recently_crossed_boundary) + int recently_crossed_boundary) { packet->recently_crossed_boundary = recently_crossed_boundary; } @@ -188,7 +188,7 @@ static inline int rpacket_get_virtual_packet_flag (const rpacket_t * packet) } static inline void rpacket_set_virtual_packet_flag (rpacket_t * packet, - int virtual_packet_flag) + int virtual_packet_flag) { packet->virtual_packet_flag = virtual_packet_flag; } @@ -199,7 +199,7 @@ static inline int rpacket_get_virtual_packet (const rpacket_t * packet) } static inline void rpacket_set_virtual_packet (rpacket_t * packet, - int virtual_packet) + int virtual_packet) { packet->virtual_packet = virtual_packet; } @@ -270,7 +270,7 @@ static inline void rpacket_reset_tau_event (rpacket_t * packet, rk_state *mt_sta } tardis_error_t rpacket_init (rpacket_t * packet, storage_model_t * storage, - int packet_index, int virtual_packet_flag); + int packet_index, int virtual_packet_flag); /* New getter and setter methods for continuum implementation */