diff --git a/src/bindings/pybindings.cpp b/src/bindings/pybindings.cpp index 95c26cab..5befbbdc 100644 --- a/src/bindings/pybindings.cpp +++ b/src/bindings/pybindings.cpp @@ -192,12 +192,9 @@ PYBIND11_MODULE(core, module) { "Feautrier solver.") /// Solver is bugged, so removed from the api, as the shortchar solver can /// replace it - // .def ( - // "compute_radiation_field_feautrier_order_2_uv", - // &Model::compute_radiation_field_feautrier_order_2_uv, - // "Compute the radiation field for the modle using the 2nd-order - // Feautrier solver." - // ) + .def("compute_radiation_field_feautrier_order_2_uv", + &Model::compute_radiation_field_feautrier_order_2_uv, + "Compute the radiation field for the modle using the 2nd-order Feautrier solver.") .def("compute_radiation_field_feautrier_order_2_anis", &Model::compute_radiation_field_feautrier_order_2_anis, "Compute the radiation field for the modle using the 2nd-order " diff --git a/src/model/model.cpp b/src/model/model.cpp index 5d0d7fac..5c44bce3 100644 --- a/src/model/model.cpp +++ b/src/model/model.cpp @@ -355,28 +355,25 @@ int Model ::compute_radiation_field_feautrier_order_2() { /// BUGGED: v computation is incorrect // /// Computer for the radiation field // ///////////////////////////////////// -// int Model :: compute_radiation_field_feautrier_order_2_uv () -// { -// cout << "Computing radiation field..." << endl; -// -// Solver solver; -// solver.setup (*this); -// -// if (parameters->one_line_approximation) -// { -// solver.solve_feautrier_order_2_uv (*this); -// return (0); -// } -// -// if (parameters->sum_opacity_emissivity_over_all_lines) -// { -// solver.solve_feautrier_order_2_uv (*this); -// return (0); -// } -// -// solver.solve_feautrier_order_2_uv (*this); -// return (0); -// } +int Model ::compute_radiation_field_feautrier_order_2_uv() { + cout << "Computing radiation field..." << endl; + + Solver solver; + solver.setup(*this); + + if (parameters->one_line_approximation) { + solver.solve_feautrier_order_2_uv(*this); + return (0); + } + + if (parameters->sum_opacity_emissivity_over_all_lines) { + solver.solve_feautrier_order_2_uv(*this); + return (0); + } + + solver.solve_feautrier_order_2_uv(*this); + return (0); +} /// Computer for the radiation field ///////////////////////////////////// diff --git a/src/model/model.hpp b/src/model/model.hpp index 087798fd..4327b569 100644 --- a/src/model/model.hpp +++ b/src/model/model.hpp @@ -106,7 +106,7 @@ struct Model { int compute_image_for_point(const Size ray_nr, const Size p); - // int compute_radiation_field_feautrier_order_2_uv (); + int compute_radiation_field_feautrier_order_2_uv(); int compute_radiation_field_feautrier_order_2_anis(); int compute_radiation_field_feautrier_order_2_sparse(); diff --git a/src/solver/solver.hpp b/src/solver/solver.hpp index 8c9ad656..b3476da3 100644 --- a/src/solver/solver.hpp +++ b/src/solver/solver.hpp @@ -187,13 +187,10 @@ struct Solver { accel inline void solve_shortchar_order_0(Model& model, const Size o, const Size r); /// BUGGED: v computation is incorrect - // template - // accel inline void solve_feautrier_order_2_uv (Model& - // model); - // - // template - // accel inline void solve_feautrier_order_2_uv (Model& - // model, const Size o, const Size f); + template accel inline void solve_feautrier_order_2_uv(Model& model); + + template + accel inline void solve_feautrier_order_2_uv(Model& model, const Size o, const Size f); // Getters for emissivities, opacities, and boundary // conditions diff --git a/src/solver/solver.tpp b/src/solver/solver.tpp index 261e614c..85cbeb2c 100644 --- a/src/solver/solver.tpp +++ b/src/solver/solver.tpp @@ -234,81 +234,56 @@ template inline void Solver ::solve_shortchar_order_0 } /// BUGGED: v computation is incorrect -// template -// inline void Solver :: solve_feautrier_order_2_uv (Model& -// model) -// { -// // Allocate memory if not pre-allocated -// if (!model.parameters->store_intensities) -// { -// model.radiation.u.resize -// (model.parameters->hnrays(), -// model.parameters->npoints(), -// model.parameters->nfreqs()); -// model.radiation.v.resize -// (model.parameters->hnrays(), -// model.parameters->npoints(), -// model.parameters->nfreqs()); -// } -// -// -// // For each ray, solve transfer equation -// for (Size rr = 0; rr < model.parameters->hnrays(); -// rr++) -// { -// const Size ar = model.geometry.rays.antipod[rr]; -// -// cout << "--- rr = " << rr << endl; -// -// accelerated_for (o, model.parameters->npoints(), -// { -// const Real dshift_max = get_dshift_max -// (model, o); -// -// nr_ ()[centre] = o; -// shift_()[centre] = 1.0; -// -// first_() = trace_ray -// (model.geometry, o, rr, dshift_max, -1, -// centre-1, centre-1) + 1; last_ () = trace_ray -// (model.geometry, o, ar, -// dshift_max, +1, centre+1, centre ) - 1; -// n_tot_() = (last_()+1) - first_(); -// -// if (n_tot_() > 1) -// { -// for (Size f = 0; f < -// model.parameters->nfreqs(); f++) -// { -// solve_feautrier_order_2_uv -// (model, o, f); -// -// model.radiation.u(rr,o,f) = -// Su_()[centre]; -// model.radiation.v(rr,o,f) = -// Sv_()[centre]; -// } -// } -// else -// { -// for (Size f = 0; f < -// model.parameters->nfreqs(); f++) -// { -// model.radiation.u(rr,o,f) = -// boundary_intensity(model, o, -// model.radiation.frequencies.nu(o, -// f)); model.radiation.v(rr,o,f) = -// 0.0; -// } -// } -// }) -// -// pc::accelerator::synchronize(); -// } -// -// model.radiation.u.copy_ptr_to_vec(); -// model.radiation.v.copy_ptr_to_vec(); -// } +template inline void Solver ::solve_feautrier_order_2_uv(Model& model) { + // Allocate memory if not pre-allocated + if (!model.parameters->store_intensities) { + model.radiation.u.resize( + model.parameters->hnrays(), model.parameters->npoints(), model.parameters->nfreqs()); + model.radiation.v.resize( + model.parameters->hnrays(), model.parameters->npoints(), model.parameters->nfreqs()); + } + + // For each ray, solve transfer equation + for (Size rr = 0; rr < model.parameters->hnrays(); rr++) { + const Size ar = model.geometry.rays.antipod[rr]; + + cout << "--- rr = " << rr << endl; + + accelerated_for(o, model.parameters->npoints(), { + const Real dshift_max = get_dshift_max(model, o); + + nr_()[centre] = o; + shift_()[centre] = 1.0; + + first_() = + trace_ray(model.geometry, o, rr, dshift_max, -1, centre - 1, centre - 1) + + 1; + last_() = + trace_ray(model.geometry, o, ar, dshift_max, +1, centre + 1, centre) - 1; + n_tot_() = (last_() + 1) - first_(); + + if (n_tot_() > 1) { + for (Size f = 0; f < model.parameters->nfreqs(); f++) { + solve_feautrier_order_2_uv(model, o, f); + + model.radiation.u(rr, o, f) = Su_()[centre]; + model.radiation.v(rr, o, f) = Sv_()[centre]; + } + } else { + for (Size f = 0; f < model.parameters->nfreqs(); f++) { + model.radiation.u(rr, o, f) = + boundary_intensity(model, o, model.radiation.frequencies.nu(o, f)); + model.radiation.v(rr, o, f) = 0.0; + } + } + }) + + pc::accelerator::synchronize(); + } + + model.radiation.u.copy_ptr_to_vec(); + model.radiation.v.copy_ptr_to_vec(); +} template inline void Solver ::solve_feautrier_order_2_sparse(Model& model) { @@ -2356,144 +2331,161 @@ accel inline void Solver ::image_optical_depth(Model& model, const Size o, const optical_depth_() = tau; } -/// BUGGED: v computation is incorrect -// /// Solver for Feautrier equation along ray pairs using -// the (ordinary) -// /// 2nd-order solver, without adaptive optical depth -// increments -// /////////////////////////////////////////////////////////////////////// -// template -// accel inline void Solver :: solve_feautrier_order_2_uv -// (Model& model, const Size o, const Size f) -// { -// const Real freq = model.radiation.frequencies.nu(o, -// f); const Size l = -// model.radiation.frequencies.corresponding_line[f]; -// -// Real eta_c, chi_c, dtau_c, term_c; -// Real eta_n, chi_n, dtau_n, term_n; -// -// const Size first = first_(); -// const Size last = last_ (); -// const Size n_tot = n_tot_(); -// -// Vector& dZ = dZ_ (); -// Vector& nr = nr_ (); -// Vector& shift = shift_(); -// -// Vector& inverse_chi = inverse_chi_(); -// -// Vector& Su = Su_(); -// Vector& Sv = Sv_(); -// -// Vector& A = A_ (); -// Vector& C = C_ (); -// Vector& inverse_A = inverse_A_(); -// Vector& inverse_C = inverse_C_(); -// -// Vector& FF = FF_(); -// Vector& FI = FI_(); -// -// bool compute_curr_opacity = true; // for the first -// point, we need to compute both the curr and next -// opacity (and source) -// -// compute_source_dtau(model, nr[first], -// nr[first+1], l, freq*shift[first], -// freq*shift[first+1], shift[first], shift[first+1], -// dZ[first], compute_curr_opacity, dtau_n, chi_c, -// chi_n, term_c, term_n); -// -// // Set boundary conditions -// const Real inverse_dtau_f = one / dtau_n; -// -// C[first] = two * inverse_dtau_f * -// inverse_dtau_f; -// inverse_C[first] = one / C[first]; // Required for -// Lambda_diag -// -// const Real Bf_min_Cf = one + two * inverse_dtau_f; -// const Real Bf = Bf_min_Cf + C[first]; -// const Real I_bdy_f = boundary_intensity (model, -// nr[first], freq*shift[first]); -// -// Su[first] = term_c + two * I_bdy_f * inverse_dtau_f; -// Sv[first] = two * inverse_dtau_f * (I_bdy_f - -// term_c); -// -// Su[first] /= Bf; -// Sv[first] /= Bf; -// -// /// Write economically: F[first] = (B[first] - -// C[first]) / C[first]; FF[first] = half * Bf_min_Cf * -// dtau_n * dtau_n; FI[first] = one / (one + FF[first]); -// -// -// /// Set body of Feautrier matrix -// for (Size n = first+1; n < last; n++) -// { -// term_c = term_n; -// dtau_c = dtau_n; -// eta_c = eta_n; -// chi_c = chi_n; -// -// compute_source_dtau(model, nr[n], -// nr[n+1], l, freq*shift[n], freq*shift[n+1], -// shift[n], shift[n+1], dZ[n], -// compute_curr_opacity, dtau_n, chi_c, chi_n, -// term_c, term_n); -// -// const Real dtau_avg = half * (dtau_c + dtau_n); -// inverse_A[n] = dtau_avg * dtau_c; -// inverse_C[n] = dtau_avg * dtau_n; -// -// A[n] = one / inverse_A[n]; -// C[n] = one / inverse_C[n]; -// -// /// Use the previously stored value of the source -// function Su[n] = term_c; -// -// FF[n] = (A[n] * FF[n-1] * FI[n-1] + one) * -// inverse_C[n]; FI[n] = one / (one + FF[n]); Su[n] -// = (A[n] * Su[n-1] + Su[n]) * FI[n] * -// inverse_C[n]; Sv[n] = (A[n] * Sv[n-1] ) * -// FI[n] * inverse_C[n]; -// } -// -// -// /// Set boundary conditions -// const Real inverse_dtau_l = one / dtau_n; -// -// A[last] = two * inverse_dtau_l * inverse_dtau_l; -// -// const Real Bl_min_Al = one + two * inverse_dtau_l; -// const Real Bl = Bl_min_Al + A[last]; -// -// const Real denominator = one / (Bl * FF[last-1] + -// Bl_min_Al); -// -// const Real I_bdy_l = boundary_intensity (model, -// nr[last], freq*shift[last]); -// -// Su[last] = term_n + two * I_bdy_l * inverse_dtau_l; -// Sv[last] = two * inverse_dtau_l * (I_bdy_l - term_n); -// -// Su[last] = (A[last] * Su[last-1] + Su[last]) * (one + -// FF[last-1]) * denominator; Sv[last] = (A[last] * -// Sv[last-1] ) * (one + FF[last-1]) * -// denominator; -// -// if (centre < last) -// { -// for (long n = last-1; n >= centre; n--) // use -// long in reverse loops! -// { -// Su[n] += Su[n+1] * FI[n]; -// Sv[n] += Sv[n+1] * FI[n]; -// } -// } -// -// } +/// Solver for Feautrier equation along ray pairs using the (ordinary) +/// 2nd-order solver, without adaptive optical depth increments +/// Computes both the mean intensity u and the flux v +////////////////////////////////////////////////////////////////////////// +template +accel inline void Solver ::solve_feautrier_order_2_uv(Model& model, const Size o, const Size f) { + // Note: in comments, we have another option for computing v (Sv at every location), by + // mirroring the solution steps for u with slightly different data. This might be slower, and a + // bit more work to make second-order accurate. + // Currently, only sv at the center position gets computed using du/dtau = -v + const Real freq = model.radiation.frequencies.nu(o, f); + const Size l = model.radiation.frequencies.corresponding_line[f]; + + Real eta_c, chi_c, dtau_c, term_c; + Real eta_n, chi_n, dtau_n, term_n; + + const Size first = first_(); + const Size last = last_(); + const Size n_tot = n_tot_(); + + Vector& dZ = dZ_(); + Vector& nr = nr_(); + Vector& shift = shift_(); + + Vector& inverse_chi = inverse_chi_(); + + Vector& Su = Su_(); + Vector& Sv = Sv_(); + + Vector& A = A_(); + Vector& C = C_(); + Vector& inverse_A = inverse_A_(); + Vector& inverse_C = inverse_C_(); + + Vector& FF = FF_(); + Vector& FI = FI_(); + + bool compute_curr_opacity = + true; // for the first point, we need to compute both the curr and next opacity(and source) + + compute_source_dtau(model, nr[first], nr[first + 1], l, freq * shift[first], + freq * shift[first + 1], shift[first], shift[first + 1], dZ[first], compute_curr_opacity, + dtau_n, chi_c, chi_n, term_c, term_n); + + // Set boundary conditions + const Real inverse_dtau_f = one / dtau_n; + + C[first] = two * inverse_dtau_f * inverse_dtau_f; + inverse_C[first] = one / C[first]; // Required for Lambda_diag + + const Real Bf_min_Cf = one + two * inverse_dtau_f; + const Real Bf = Bf_min_Cf + C[first]; + const Real I_bdy_f = boundary_intensity(model, nr[first], freq * shift[first]); + // TODO: if re-implementing other option, make this dSdtau second order accurate + // Real dSdtau = (term_n - term_c) / dtau_n; + // only first order accurate dSdtau, as otherwise the solver needs to be rewritten + // Current algorithm only allows us to access the data at current and next point + + Su[first] = term_c + two * I_bdy_f * inverse_dtau_f; + // Sv[first] = -dSdtau + two * inverse_dtau_f * (I_bdy_f - term_c); + + Su[first] /= Bf; + // Sv[first] /= Bf; + + /// Write economically: F[first] = (B[first] - C[first]) / C[first]; + FF[first] = half * Bf_min_Cf * dtau_n * dtau_n; + FI[first] = one / (one + FF[first]); + + /// Set body of Feautrier matrix + for (Size n = first + 1; n < last; n++) { + term_c = term_n; + dtau_c = dtau_n; + eta_c = eta_n; + chi_c = chi_n; + + compute_source_dtau(model, nr[n], nr[n + 1], l, freq * shift[n], + freq * shift[n + 1], shift[n], shift[n + 1], dZ[n], compute_curr_opacity, dtau_n, chi_c, + chi_n, term_c, term_n); + + const Real dtau_avg = half * (dtau_c + dtau_n); + // TODO: if re-implementing other option, make this dSdtau second order accurate + // dSdtau = (term_n - term_c) / dtau_n; + inverse_A[n] = dtau_avg * dtau_c; + inverse_C[n] = dtau_avg * dtau_n; + + A[n] = one / inverse_A[n]; + C[n] = one / inverse_C[n]; + + /// Use the previously stored value of the source function + Su[n] = term_c; + + FF[n] = (A[n] * FF[n - 1] * FI[n - 1] + one) * inverse_C[n]; + FI[n] = one / (one + FF[n]); + Su[n] = (A[n] * Su[n - 1] + Su[n]) * FI[n] * inverse_C[n]; + // Sv[n] = (A[n] * Sv[n - 1] - dSdtau) * FI[n] * inverse_C[n]; + } + + /// Set boundary conditions + const Real inverse_dtau_l = one / dtau_n; + + A[last] = two * inverse_dtau_l * inverse_dtau_l; + + const Real Bl_min_Al = one + two * inverse_dtau_l; + const Real Bl = Bl_min_Al + A[last]; + + const Real denominator = one / (Bl * FF[last - 1] + Bl_min_Al); + + const Real I_bdy_l = boundary_intensity(model, nr[last], freq * shift[last]); + + Su[last] = term_n + two * I_bdy_l * inverse_dtau_l; + // Sv[last] = -dSdtau - two * inverse_dtau_l * (I_bdy_l - term_n); + // Different sign for Sv last boundary condition extra term (2/dtau(I-S))! (should be + // assymetric) + + Su[last] = (A[last] * Su[last - 1] + Su[last]) * (one + FF[last - 1]) * denominator; + // Sv[last] = (A[last] * Sv[last - 1] + Sv[last]) * (one + FF[last - 1]) * denominator; + + if (centre < last) { + for (long n = last - 1; n >= centre; n--) // use long in reverse loops ! + { + Su[n] += Su[n + 1] * FI[n]; + // Sv[n] += Sv[n + 1] * FI[n]; + } + } else { + // Compute v using boundary at the end + Sv[last] = Su[last] - I_bdy_l; + return; + } + // Do one extra step for computing the derivative of the mean intensity + if (centre > first) { + Su[centre - 1] += Su[centre] * FI[centre - 1]; + // Recompute the optical depth increments, as we forgot to save them; FIXME: add + // threadprivate dtau_() to solver.hpp + compute_curr_opacity = true; + compute_source_dtau(model, nr[centre - 1], nr[centre], l, freq * shift[centre - 1], + freq * shift[centre], shift[centre - 1], shift[centre], dZ[centre - 1], + compute_curr_opacity, dtau_n, chi_c, chi_n, term_c, term_n); + const Real dtaumin = dtau_n; + compute_source_dtau(model, nr[centre], nr[centre + 1], l, freq * shift[centre], + freq * shift[centre + 1], shift[centre], shift[centre + 1], dZ[centre], + compute_curr_opacity, dtau_n, chi_c, chi_n, term_c, term_n); + const Real dtauplus = dtau_n; + // TODO: optimize this calculation + const Real coeffmin = -dtauplus / (dtaumin * dtaumin + dtaumin * dtauplus); + const Real coeffplus = dtaumin / (dtaumin * dtauplus + dtauplus * dtauplus); + const Real coeffzero = -coeffmin - coeffplus; + + Sv[centre] = + -(coeffmin * Su[centre - 1] + coeffplus * Su[centre + 1] + coeffzero * Su[centre]); + } else { + // Compute v using boundary at the start + Sv[first] = I_bdy_f - Su[first]; + return; + } +} accel inline void Solver ::set_eta_and_chi(Model& model, const Size rr) const { model.eta.resize(model.parameters->npoints(), model.parameters->nfreqs()); diff --git a/tests/benchmarks/analytic/all_constant_single_ray.py b/tests/benchmarks/analytic/all_constant_single_ray.py index f7bfc2bc..cce4200d 100644 --- a/tests/benchmarks/analytic/all_constant_single_ray.py +++ b/tests/benchmarks/analytic/all_constant_single_ray.py @@ -99,6 +99,13 @@ def run_model (nosave=False): timer4.stop() u_2f = np.array(model.radiation.u) + timer5 = tools.Timer('feautrier 2 uv') + timer5.start() + model.compute_radiation_field_feautrier_order_2_uv () + timer5.stop() + u_2f_uv = np.array(model.radiation.u) + v_2f_uv = np.array(model.radiation.v) + x = np.array(model.geometry.points.position)[:,0] nu = np.array(model.radiation.frequencies.nu) @@ -125,6 +132,8 @@ def u_ (x): error_u_0s = np.abs(tools.relative_error (u_(x), u_0s[0,:,0])) error_u_2f = np.abs(tools.relative_error (u_(x), u_2f[0,:,0])) + error_I_0_2f_uv = np.abs(tools.relative_error (I_0(x), u_2f_uv[0,:,0]+v_2f_uv[0,:,0])) + error_I_1_2f_uv = np.abs(tools.relative_error (I_1(x), u_2f_uv[0,:,0]-v_2f_uv[0,:,0])) result = f'--- Benchmark name ----------------------------\n' result += f'{modelName }\n' @@ -136,6 +145,8 @@ def u_ (x): result += f'--- Accuracy ----------------------------------\n' result += f'max error in shortchar 0 = {np.max(error_u_0s)}\n' result += f'max error in feautrier 2 = {np.max(error_u_2f)}\n' + result += f'max error in I_0 2f uv = {np.max(error_I_0_2f_uv)}\n' + result += f'max error in I_1 2f uv = {np.max(error_I_1_2f_uv)}\n' result += f'--- Timers ------------------------------------\n' result += f'{timer1.print() }\n' result += f'{timer2.print() }\n' @@ -168,14 +179,17 @@ def u_ (x): #error bounds are chosen somewhat arbitrarily, based on previously obtained results; this should prevent serious regressions. FEAUTRIER_AS_EXPECTED=(np.max(error_u_2f)<1.7e-4) FIRSTORDER_AS_EXPECTED=(np.max(error_u_0s)<5.1e-8) + FEAUTRIER_UV_AS_EXPECTED=(np.max(error_I_0_2f_uv)<2.0e-4) and (np.max(error_I_1_2f_uv)<2.0e-4) if not FIRSTORDER_AS_EXPECTED: print("First order solver max error too large: ", np.max(error_u_0s)) if not FEAUTRIER_AS_EXPECTED: print("Feautrier solver max error too large: ", np.max(error_u_2f)) + if not FEAUTRIER_UV_AS_EXPECTED: + print("Feautrier solver with uv max error too large: ", np.max(error_I_0_2f_uv), np.max(error_I_1_2f_uv)) - return (FEAUTRIER_AS_EXPECTED&FIRSTORDER_AS_EXPECTED) + return (FEAUTRIER_AS_EXPECTED&FIRSTORDER_AS_EXPECTED&FEAUTRIER_UV_AS_EXPECTED) def run_test (nosave=False): diff --git a/tests/benchmarks/analytic/all_constant_single_ray_testuv.py b/tests/benchmarks/analytic/all_constant_single_ray_testuv.py new file mode 100644 index 00000000..d61c2dd6 --- /dev/null +++ b/tests/benchmarks/analytic/all_constant_single_ray_testuv.py @@ -0,0 +1,212 @@ +import os +import sys + +curdir = os.path.dirname(os.path.realpath(__file__)) +datdir = f'{curdir}/../../data/' +moddir = f'{curdir}/../../models/' +resdir = f'{curdir}/../../results/' + +import numpy as np +import matplotlib.pyplot as plt +import magritte.tools as tools +import magritte.setup as setup +import magritte.core as magritte + + +dimension = 1 +npoints = 50 +nrays = 2 +nspecs = 3 +nlspecs = 1 +nquads = 1 + +nH2 = 1.0E+12 # [m^-3] +nTT = 1.0E+03 # [m^-3] +temp = 4.5E+00 # [K] +turb = 0.0E+00 # [m/s] +dx = 1.0E+12 # [m] +dv = 0.0E+00 / magritte.CC # [fraction of speed of light] + + +def create_model (): + """ + Create a model file for the all_constant benchmark, single ray. + """ + + modelName = f'all_constant_single_ray_testuv' + modelFile = f'{moddir}{modelName}.hdf5' + lamdaFile = f'{datdir}test.txt' + + + model = magritte.Model () + model.parameters.set_spherical_symmetry(False) + model.parameters.set_model_name (modelFile) + model.parameters.set_dimension (dimension) + model.parameters.set_npoints (npoints) + model.parameters.set_nrays (nrays) + model.parameters.set_nspecs (nspecs) + model.parameters.set_nlspecs (nlspecs) + model.parameters.set_nquads (nquads) + + model.geometry.points.position.set([[i*dx, 0, 0] for i in range(npoints)]) + model.geometry.points.velocity.set([[i*dv, 0, 0] for i in range(npoints)]) + + model.chemistry.species.abundance = [[nTT, nH2, 0.0] for _ in range(npoints)] + model.chemistry.species.symbol = ['test', 'H2', 'e-'] + + # model.thermodynamics.temperature.gas .set( temp * np.ones(npoints)) + #Changing the temperature such that the source function is not constant + model.thermodynamics.temperature.gas .set( temp * (1.0+10*np.arange(npoints)/npoints)) + model.thermodynamics.turbulence.vturb2.set((turb/magritte.CC)**2 * np.ones(npoints)) + + model = setup.set_Delaunay_neighbor_lists (model) + model = setup.set_Delaunay_boundary (model) + model = setup.set_boundary_condition_CMB (model) + model = setup.set_uniform_rays (model) + model = setup.set_linedata_from_LAMDA_file(model, lamdaFile) + model = setup.set_quadrature (model) + + model.write() + + return #magritte.Model (modelFile) + + +def run_model (nosave=False): + + modelName = f'all_constant_single_ray_testuv' + modelFile = f'{moddir}{modelName}.hdf5' + timestamp = tools.timestamp() + + timer1 = tools.Timer('reading model') + timer1.start() + model = magritte.Model (modelFile) + timer1.stop() + + timer2 = tools.Timer('setting model') + timer2.start() + model.compute_spectral_discretisation () + model.compute_inverse_line_widths () + model.compute_LTE_level_populations () + timer2.stop() + + timer3 = tools.Timer('shortchar 0 ') + timer3.start() + model.compute_radiation_field_shortchar_order_0 () + timer3.stop() + u_0s = np.array(model.radiation.u) + I_0s = np.array(model.radiation.I) + + timer4 = tools.Timer('feautrier 2 ') + timer4.start() + model.compute_radiation_field_feautrier_order_2 () + timer4.stop() + u_2f = np.array(model.radiation.u) + + timer5 = tools.Timer('feautrier 2 uv') + timer5.start() + model.compute_radiation_field_feautrier_order_2_uv () + timer5.stop() + u_2f_uv = np.array(model.radiation.u) + v_2f_uv = np.array(model.radiation.v) + + x = np.array(model.geometry.points.position)[:,0] + nu = np.array(model.radiation.frequencies.nu) + + ld = model.lines.lineProducingSpecies[0].linedata + + k = 0 + + frq = ld.frequency[k] + pop = tools.LTEpop (ld, temp) * nTT + phi = tools.profile (ld, k, temp, (turb/magritte.CC)**2, frq) + eta = tools.lineEmissivity (ld, pop)[k] * phi + chi = tools.lineOpacity (ld, pop)[k] * phi + src = tools.lineSource (ld, pop)[k] + bdy = tools.I_CMB (frq) + + # def I_0 (x): + # return src + (bdy-src)*np.exp(-chi*x) + + # def I_1 (x): + # return src + (bdy-src)*np.exp(-chi*(x[-1]-x)) + + # def u_ (x): + # return 0.5 * (I_0(x) + I_1(x)) + + # error_u_0s = np.abs(tools.relative_error (u_(x), u_0s[0,:,0])) + # error_u_2f = np.abs(tools.relative_error (u_(x), u_2f[0,:,0])) + error_I_0_2f_uv = np.abs(tools.relative_error (I_0s[0,:,0], u_2f_uv[0,:,0]+v_2f_uv[0,:,0])) + error_I_1_2f_uv = np.abs(tools.relative_error (I_0s[1,:,0], u_2f_uv[0,:,0]-v_2f_uv[0,:,0])) + + result = f'--- Benchmark name ----------------------------\n' + result += f'{modelName }\n' + result += f'--- Parameters --------------------------------\n' + result += f'dimension = {model.parameters.dimension() }\n' + result += f'npoints = {model.parameters.npoints () }\n' + result += f'nrays = {model.parameters.nrays () }\n' + result += f'nquads = {model.parameters.nquads () }\n' + result += f'--- Accuracy ----------------------------------\n' + result += f'max error in I_0 2f uv = {np.max(error_I_0_2f_uv)}\n' + result += f'max error in I_1 2f uv = {np.max(error_I_1_2f_uv)}\n' + result += f'--- Timers ------------------------------------\n' + result += f'{timer1.print() }\n' + result += f'{timer2.print() }\n' + result += f'{timer3.print() }\n' + result += f'{timer4.print() }\n' + result += f'-----------------------------------------------\n' + + print(result) + + if not nosave: + with open(f'{resdir}{modelName}-{timestamp}.log' ,'w') as log: + log.write(result) + + fig = plt.figure(dpi=150) + plt.title(modelName) + plt.scatter(x, u_0s[0,:,0], s=0.5, label='0s', zorder=1) + plt.scatter(x, u_2f[0,:,0], s=0.5, label='2f', zorder=1) + # plt.plot(x, u_(x), c='lightgray', zorder=0) + plt.legend() + # plt.xscale('log') + plt.xlabel('r [m]') + plt.ylabel('Mean intensity [W/m$^{2}$]') + # plt.savefig(f'{resdir}{modelName}-{timestamp}.png', dpi=150) + plt.figure() + plt.plot(x, I_0s[0,:,0], label='I shortchar 0') + plt.plot(x, u_2f_uv[0,:,0]+v_2f_uv[0,:,0], label='I feautrier uv') + plt.legend() + plt.figure() + plt.plot(x, I_0s[1,:,0], label='I shortchar 1') + plt.plot(x, u_2f_uv[0,:,0]-v_2f_uv[0,:,0], label='I feautrier uv') + plt.legend() + plt.figure() + plt.plot(x, I_0s[0,:,0]-(u_2f_uv[0,:,0]+v_2f_uv[0,:,0]), label='diff') + plt.figure() + plt.plot(x, (I_0s[0,:,0]-I_0s[1,:,0])/2.0, label='shortchar diff') + plt.plot(x, v_2f_uv[0,:,0], label='feautrier v') + plt.legend() + plt.show() + + #error bounds are chosen somewhat arbitrarily, based on previously obtained results; this should prevent serious regressions. + FEAUTRIER_UV_AS_EXPECTED=(np.mean(error_I_0_2f_uv)<0.021) and (np.mean(error_I_1_2f_uv)<0.01) + + if not FEAUTRIER_UV_AS_EXPECTED: + print("Feautrier solver with uv mean error too large: ", np.mean(error_I_0_2f_uv), np.mean(error_I_1_2f_uv)) + + + return (FEAUTRIER_UV_AS_EXPECTED) + + +def run_test (nosave=False): + + create_model () + run_model (nosave) + + return + + +if __name__ == '__main__': + + nosave = (len(sys.argv) > 1) and (sys.argv[1] == 'nosave') + + run_test (nosave) diff --git a/tests/benchmarks/analytic/test_benchmarks.py b/tests/benchmarks/analytic/test_benchmarks.py index 37c036f2..ff31f56e 100644 --- a/tests/benchmarks/analytic/test_benchmarks.py +++ b/tests/benchmarks/analytic/test_benchmarks.py @@ -1,6 +1,8 @@ from reject_nonexistent_species import create_model as reject_nonexistent_species_setup from all_constant_single_ray import create_model as all_constant_setup from all_constant_single_ray import run_model as all_constant_run +from all_constant_single_ray_testuv import create_model as all_constant_testuv_setup +from all_constant_single_ray_testuv import run_model as all_constant_testuv_run from all_zero_single_ray import create_model as all_zero_setup from all_zero_single_ray import run_model as all_zero_run from density_distribution_1D import create_model as density_dist_setup @@ -47,6 +49,13 @@ def test_all_constant_setup(self): def test_all_constant_run(self): assert all_constant_run(nosave=True) + @pytest.mark.incremental + class TestUVFeautrier: + def test_uv_feautrier_setup(self): + all_constant_testuv_setup() + + def test_uv_feautrier_run(self): + assert all_constant_testuv_run(nosave=True) @pytest.mark.incremental class TestAllZero: @@ -125,5 +134,3 @@ def test_velocity_gradient_3D_new_imager_setup(self): def test_velocity_gradient_3D_new_imager_run(self): assert velocity_gradient_new_imager_run(nosave=True, use_widgets=False) - -#TODO ADD ALL ANALYTIC BENCHMARKS, get some manner of performance somewhere else