Skip to content

Commit

Permalink
Merge branch 'imager_correction' of github.com:Magritte-code/Magritte…
Browse files Browse the repository at this point in the history 
… into imager_correction
  • Loading branch information
@ThomasCeulemans
ThomasCeulemans committed Oct 16, 2024
2 parents 6c486c6 + e3b007d commit 04e6878
Showing 1 changed file with 100 additions and 72 deletions.
172 changes: 100 additions & 72 deletions src/solver/solver.tpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1154,7 +1154,7 @@ accel inline void Solver ::get_eta_and_chi<None>(const Model& model, const Size
const Real freq, Real& eta, Real& chi) const {
// Initialize
eta = 0.0;
chi = model.parameters->min_opacity;
chi = 0.0; // model.parameters->min_opacity;

// Set line emissivity and opacity
for (Size l = 0; l < model.parameters->nlines(); l++) {
Expand Down Expand Up @@ -1182,7 +1182,7 @@ accel inline void Solver ::get_eta_and_chi_interpolated<None>(const Model& model
const Real freq, Real& eta, Real& chi) const {
// Initialize
eta = 0.0;
chi = model.parameters->min_opacity;
chi = 0.0; // model.parameters->min_opacity;

// Set line emissivity and opacity
for (Size l = 0; l < model.parameters->nlines(); l++) {
Expand Down Expand Up @@ -1216,7 +1216,7 @@ accel inline void Solver ::get_eta_and_chi<CloseLines>(const Model& model, const
const Real freq, Real& eta, Real& chi) const {
// Initialize
eta = 0.0;
chi = model.parameters->min_opacity;
chi = 0.0; // model.parameters->min_opacity;

const Real upper_bound_line_width =
model.parameters->max_distance_opacity_contribution
Expand Down Expand Up @@ -1263,7 +1263,7 @@ accel inline void Solver ::get_eta_and_chi_interpolated<CloseLines>(const Model&
const Real freq, Real& eta, Real& chi) const {
// Initialize
eta = 0.0;
chi = model.parameters->min_opacity;
chi = 0.0; // model.parameters->min_opacity;

const Real upper_bound_line_width =
model.parameters->max_distance_opacity_contribution
Expand Down Expand Up @@ -1318,7 +1318,7 @@ accel inline void Solver ::get_eta_and_chi<OneLine>(
const Real prof = gaussian(model.lines.inverse_width(p, l), diff);

eta = prof * model.lines.emissivity(p, l);
chi = prof * model.lines.opacity(p, l) + model.parameters->min_opacity;
chi = prof * model.lines.opacity(p, l); // + model.parameters->min_opacity;
}

/// Getter for the emissivity (eta) and the opacity (chi), with interpolation support
Expand Down Expand Up @@ -1347,7 +1347,7 @@ accel inline void Solver ::get_eta_and_chi_interpolated<OneLine>(const Model& mo
model.lines.opacity(p1, l), model.lines.opacity(p2, l), interp_factor);

eta = prof * interp_emissivity;
chi = prof * interp_opacity + model.parameters->min_opacity;
chi = prof * interp_opacity; // + model.parameters->min_opacity;
}

/// Apply trapezium rule to x_crt and x_nxt
Expand Down Expand Up @@ -1636,7 +1636,7 @@ accel inline Real Solver ::solve_shortchar_order_0_ray_forward(
Real eta, chi; // eta is dummy var
// chi needs to be recomputed, as we have no guarantee the for loop runs at least once
get_eta_and_chi<approx>(model, o, k, freq_line, eta, chi);
Real inverse_chi = 1.0 / chi;
Real inverse_chi = 1.0 / (chi + model.parameters->min_opacity);
Real phi = model.thermodynamics.profile(invr_mass, o, freq_line, freq);
// const Real lambda_factor =
// (dtau+expm1f(-dtau))/dtau;// If one wants to
Expand Down Expand Up @@ -1739,7 +1739,7 @@ accel inline Real Solver ::solve_shortchar_order_0_ray_backward(
Real eta, chi; // eta is dummy var
// chi has to be recomputed, as we have no guarantee the for loop runs at least once
get_eta_and_chi<approx>(model, o, k, freq_line, eta, chi);
Real inverse_chi = 1.0 / chi;
Real inverse_chi = 1.0 / (chi + model.parameters->min_opacity);
Real phi = model.thermodynamics.profile(invr_mass, o, freq_line, freq);
// const Real lambda_factor =
// (dtau+expm1f(-dtau))/dtau;// If one wants to
Expand Down Expand Up @@ -1796,7 +1796,7 @@ accel inline void Solver ::update_Lambda(Model& model, const Size rr, const Size
// * L_diag[centre] * inverse_opacity;
Real eta, chi; // eta is dummy var
get_eta_and_chi<approx>(model, nr[centre], k, freq_line, eta, chi);
Real inverse_chi = 1.0 / chi;
Real inverse_chi = 1.0 / (chi + model.parameters->min_opacity);

Real frq = freqs.nu(nr[centre], f) * shift[centre];
Real phi = thermodyn.profile(invr_mass, nr[centre], freq_line, frq);
Expand All @@ -1817,7 +1817,7 @@ accel inline void Solver ::update_Lambda(Model& model, const Size rr, const Size
// 1/HH_OVER_FOUR_PI L = constante * L_lower(m,n)
// * inverse_opacity;
get_eta_and_chi<approx>(model, nr[n], k, freq_line, eta, chi);
Real inverse_chi = 1.0 / chi;
Real inverse_chi = 1.0 / (chi + model.parameters->min_opacity);

frq = freqs.nu(nr[n], f) * shift[n];
phi = thermodyn.profile(invr_mass, nr[n], freq_line, frq);
Expand All @@ -1837,7 +1837,7 @@ accel inline void Solver ::update_Lambda(Model& model, const Size rr, const Size
// 1/HH_OVER_FOUR_PI L = constante * L_upper(m,n)
// * inverse_opacity;
get_eta_and_chi<approx>(model, nr[n], k, freq_line, eta, chi);
Real inverse_chi = 1.0 / chi;
Real inverse_chi = 1.0 / (chi + model.parameters->min_opacity);

frq = freqs.nu(nr[n], f) * shift[n];
phi = thermodyn.profile(invr_mass, nr[n], freq_line, frq);
Expand Down Expand Up @@ -1882,7 +1882,8 @@ inline void Solver ::compute_S_dtau_line_integrated<OneLine>(Model& model, Size
Real& Scurr, Real& Snext) {
// FIXME: line idx is wrong
dtau = compute_dtau_single_line(model, currpoint, nextpoint, currpoint_interp_idx,
nextpoint_interp_idx, lineidx, currfreq, nextfreq, curr_interp, next_interp, dZ);
nextpoint_interp_idx, lineidx, currfreq, nextfreq, curr_interp, next_interp, dZ)
+ model.parameters->min_dtau;
// const Real curr_opacity = interp_helper.interpolate_log(model.lines.opacity(currpoint,
// lineidx),
// model.lines.opacity(currpoint_interp_idx, lineidx), curr_interp);
Expand All @@ -1898,17 +1899,17 @@ inline void Solver ::compute_S_dtau_line_integrated<OneLine>(Model& model, Size
// Scurr = curr_emissivity / curr_opacity;
// Snext = next_emissivity / next_opacity;

Real Scurr_p =
model.lines.emissivity(currpoint, lineidx) / model.lines.opacity(currpoint, lineidx);
Real Snext_p =
model.lines.emissivity(nextpoint, lineidx) / model.lines.opacity(nextpoint, lineidx);
Real Scurr_interp_p = model.lines.emissivity(currpoint_interp_idx, lineidx)
/ model.lines.opacity(currpoint_interp_idx, lineidx);
Real Snext_interp_p = model.lines.emissivity(nextpoint_interp_idx, lineidx)
/ model.lines.opacity(nextpoint_interp_idx, lineidx);

Scurr = interp_helper.interpolate_log(Scurr_p, Scurr_interp_p, curr_interp);
Snext = interp_helper.interpolate_log(Snext_p, Snext_interp_p, next_interp);
// Real Scurr_p =
// model.lines.emissivity(currpoint, lineidx) / model.lines.opacity(currpoint, lineidx);
// Real Snext_p =
// model.lines.emissivity(nextpoint, lineidx) / model.lines.opacity(nextpoint, lineidx);
// Real Scurr_interp_p = model.lines.emissivity(currpoint_interp_idx, lineidx)
// / model.lines.opacity(currpoint_interp_idx, lineidx);
// Real Snext_interp_p = model.lines.emissivity(nextpoint_interp_idx, lineidx)
// / model.lines.opacity(nextpoint_interp_idx, lineidx);

// Scurr = interp_helper.interpolate_log(Scurr_p, Scurr_interp_p, curr_interp);
// Snext = interp_helper.interpolate_log(Snext_p, Snext_interp_p, next_interp);
}

/// Computer for the optical depth and source function when
Expand Down Expand Up @@ -1963,23 +1964,23 @@ inline void Solver ::compute_S_dtau_line_integrated<None>(Model& model, Size cur
// Real line_Scurr = curr_emissivity / curr_opacity;
// Real line_Snext = next_emissivity / next_opacity;

Real Scurr_p = model.lines.emissivity(currpoint, l) / model.lines.opacity(currpoint, l);
Real Snext_p = model.lines.emissivity(nextpoint, l) / model.lines.opacity(nextpoint, l);
Real Scurr_interp_p = model.lines.emissivity(currpoint_interp_idx, l)
/ model.lines.opacity(currpoint_interp_idx, l);
Real Snext_interp_p = model.lines.emissivity(nextpoint_interp_idx, l)
/ model.lines.opacity(nextpoint_interp_idx, l);
// Real Scurr_p = model.lines.emissivity(currpoint, l) / model.lines.opacity(currpoint, l);
// Real Snext_p = model.lines.emissivity(nextpoint, l) / model.lines.opacity(nextpoint, l);
// Real Scurr_interp_p = model.lines.emissivity(currpoint_interp_idx, l)
// / model.lines.opacity(currpoint_interp_idx, l);
// Real Snext_interp_p = model.lines.emissivity(nextpoint_interp_idx, l)
// / model.lines.opacity(nextpoint_interp_idx, l);

Real line_Scurr = interp_helper.interpolate_log(Scurr_p, Scurr_interp_p, curr_interp);
Real line_Snext = interp_helper.interpolate_log(Snext_p, Snext_interp_p, next_interp);
// Real line_Scurr = interp_helper.interpolate_log(Scurr_p, Scurr_interp_p, curr_interp);
// Real line_Snext = interp_helper.interpolate_log(Snext_p, Snext_interp_p, next_interp);

sum_dtau += line_dtau;
sum_dtau_times_Scurr += line_dtau * line_Scurr;
sum_dtau_times_Snext += line_dtau * line_Snext;
// sum_dtau += line_dtau;
// sum_dtau_times_Scurr += line_dtau * line_Scurr;
// sum_dtau_times_Snext += line_dtau * line_Snext;
}
dtau = sum_dtau;
Scurr = sum_dtau_times_Scurr / sum_dtau;
Snext = sum_dtau_times_Snext / sum_dtau;
dtau = sum_dtau + model.parameters->min_dtau;
// Scurr = sum_dtau_times_Scurr / sum_dtau;
// Snext = sum_dtau_times_Snext / sum_dtau;
}

/// Computer for the optical depth and source function when
Expand Down Expand Up @@ -2075,39 +2076,40 @@ inline void Solver ::compute_S_dtau_line_integrated<CloseLines>(Model& model, Si
// Real line_Scurr = curr_emissivity / curr_opacity;
// Real line_Snext = next_emissivity / next_opacity;

Real Scurr_p = model.lines.emissivity(currpoint, l) / model.lines.opacity(currpoint, l);
Real Snext_p = model.lines.emissivity(nextpoint, l) / model.lines.opacity(nextpoint, l);
Real Scurr_interp_p = model.lines.emissivity(currpoint_interp_idx, l)
/ model.lines.opacity(currpoint_interp_idx, l);
Real Snext_interp_p = model.lines.emissivity(nextpoint_interp_idx, l)
/ model.lines.opacity(nextpoint_interp_idx, l);
// Real Scurr_p = model.lines.emissivity(currpoint, l) / model.lines.opacity(currpoint, l);
// Real Snext_p = model.lines.emissivity(nextpoint, l) / model.lines.opacity(nextpoint, l);
// Real Scurr_interp_p = model.lines.emissivity(currpoint_interp_idx, l)
// / model.lines.opacity(currpoint_interp_idx, l);
// Real Snext_interp_p = model.lines.emissivity(nextpoint_interp_idx, l)
// / model.lines.opacity(nextpoint_interp_idx, l);

Real line_Scurr = interp_helper.interpolate_log(Scurr_p, Scurr_interp_p, curr_interp);
Real line_Snext = interp_helper.interpolate_log(Snext_p, Snext_interp_p, next_interp);
// Real line_Scurr = interp_helper.interpolate_log(Scurr_p, Scurr_interp_p, curr_interp);
// Real line_Snext = interp_helper.interpolate_log(Snext_p, Snext_interp_p, next_interp);

sum_dtau += line_dtau;
sum_dtau_times_Scurr += line_dtau * line_Scurr;
sum_dtau_times_Snext += line_dtau * line_Snext;
}
dtau = sum_dtau;
// needs extra bounding, as nothing may be added in the
// first place (above for loop may have looped over 0
// elements)
const Real bound_min_dtau = model.parameters->min_opacity * dZ;
// Correct way of bounding from below; should be able to
// deal with very minor computation errors around 0.
if (-bound_min_dtau < dtau) {
dtau = std::max(bound_min_dtau, dtau);
// sum_dtau_times_Scurr += line_dtau * line_Scurr;
// sum_dtau_times_Snext += line_dtau * line_Snext;
}
// Note: 0 source functions can be returned if no lines
// are nearby; but then the negligible lower bound gets
// returned
Scurr = sum_dtau_times_Scurr / dtau;
Snext = sum_dtau_times_Snext / dtau;

// note: due to interaction with dtau when computing all
// sources individually, we do need to recompute Scurr and
// Snext for all position increments
dtau = sum_dtau + model.parameters->min_dtau;

// // needs extra bounding, as nothing may be added in the
// // first place (above for loop may have looped over 0
// // elements)
// const Real bound_min_dtau = model.parameters->min_opacity * dZ;
// // Correct way of bounding from below; should be able to
// // deal with very minor computation errors around 0.
// if (-bound_min_dtau < dtau) {
// dtau = std::max(bound_min_dtau, dtau);
// }
// // Note: 0 source functions can be returned if no lines
// // are nearby; but then the negligible lower bound gets
// // returned
// Scurr = sum_dtau_times_Scurr / dtau;
// Snext = sum_dtau_times_Snext / dtau;

// // note: due to interaction with dtau when computing all
// // sources individually, we do need to recompute Scurr and
// // Snext for all position increments
}

/// Computes the source function and optical depth in a
Expand Down Expand Up @@ -2151,11 +2153,37 @@ accel inline void Solver ::compute_source_dtau(Model& model, Size currpoint, Siz

// fancy computation for large doppler shifts
if (using_large_shift) {
compute_curr_opacity = true;
compute_S_dtau_line_integrated<approx>(model, currpoint, nextpoint, currpoint_interp_idx,
nextpoint_interp_idx, line, curr_freq, next_freq, curr_interp, next_interp, dZ,
dtaunext, Scurr, Snext);
// OPACITY IS NOT COMPUTED IN THIS BRANCH!

// computing source function using exact same formula as in the other branch

// default computation using trapezoidal rule
if (compute_curr_opacity) // fancy computation does not
// compute the current
// opacity, so we might need
// to recompute it here
{
compute_curr_opacity = false;
Real eta_c = 0.0; // current emissivity
// also get previous opacity (emissivity does not
// matter)
// get_eta_and_chi<approx>(model, currpoint, line, curr_freq, eta_c, chicurr);
get_eta_and_chi_interpolated<approx>(model, currpoint, currpoint_interp_idx,
curr_interp, line, curr_freq, eta_c, chicurr);
Scurr = eta_c / (chicurr + model.parameters->min_opacity); // might as well compute the
// source function too
}

Real eta_n = 0.0;
// Get new radiative properties
get_eta_and_chi_interpolated<approx>(
model, nextpoint, nextpoint_interp_idx, next_interp, line, next_freq, eta_n, chinext);

Snext = eta_n / (chinext + model.parameters->min_opacity);

// OPACITY IS NOT COMPUTED IN THIS BRANCH! <before testing>
} else {
// default computation using trapezoidal rule
if (compute_curr_opacity) // fancy computation does not
Expand All @@ -2170,18 +2198,18 @@ accel inline void Solver ::compute_source_dtau(Model& model, Size currpoint, Siz
// get_eta_and_chi<approx>(model, currpoint, line, curr_freq, eta_c, chicurr);
get_eta_and_chi_interpolated<approx>(model, currpoint, currpoint_interp_idx,
curr_interp, line, curr_freq, eta_c, chicurr);
Scurr = eta_c / chicurr; // might as well compute the
// source function too
Scurr = eta_c / (chicurr + model.parameters->min_opacity); // might as well compute the
// source function too
}

Real eta_n = 0.0;
// Get new radiative properties
get_eta_and_chi_interpolated<approx>(
model, nextpoint, nextpoint_interp_idx, next_interp, line, next_freq, eta_n, chinext);

Snext = eta_n / chinext;
Snext = eta_n / (chinext + model.parameters->min_opacity);

dtaunext = half * (chicurr + chinext) * dZ;
dtaunext = half * (chicurr + chinext) * dZ + model.parameters->min_dtau;
}
}

Expand Down

0 comments on commit 04e6878

Please sign in to comment.