From 6c88778a8f8eab5ee0dac5b8d7154016dd660a89 Mon Sep 17 00:00:00 2001 From: Manfred Brath Date: Wed, 22 Apr 2026 16:18:20 +0200 Subject: [PATCH 1/3] Parallize run_cdisort --- src/disort.cc | 114 ++++++++++++++++++++++++++++++++++---------------- 1 file changed, 77 insertions(+), 37 deletions(-) diff --git a/src/disort.cc b/src/disort.cc index 0cc06a79e4..9488bdcb52 100644 --- a/src/disort.cc +++ b/src/disort.cc @@ -1006,8 +1006,6 @@ void run_cdisort(Workspace& ws, Numeric umu0 = 0.; //local azimuth angle of sun Numeric phi0 = 0.; - //Intensity of incident sun beam - Numeric fbeam = 0.; Index N_lev = p_grid.nelem(); @@ -1064,10 +1062,6 @@ void run_cdisort(Workspace& ws, ds.nphi = static_cast(nphi); Index Nlegendre = nstreams + 1; - /* Allocate memory */ - c_disort_state_alloc(&ds); - c_disort_out_alloc(&ds, &out); - // Looking direction of solar beam ds.bc.umu0 = umu0; ds.bc.phi0 = phi0; @@ -1075,17 +1069,6 @@ void run_cdisort(Workspace& ws, // Intensity of bottom-boundary isotropic illumination ds.bc.fluor = 0.; - // fill up azimuth angle and temperature array - for (Index i = 0; i < ds.nphi; i++) ds.phi[i] = aa_grid[i]; - - if (ds.flag.planck == TRUE) { - for (Index i = 0; i <= ds.nlyr; i++) - ds.temper[i] = temperature[ds.nlyr - i]; - } - - // Transform to mu, starting with negative values - for (Index i = 0; i < ds.numu; i++) ds.umu[i] = -cos(za_grid[i] * PI / 180); - //gas absorption Matrix ext_bulk_gas(nf, ds.nlyr + 1); get_gasoptprop(ws, @@ -1102,26 +1085,14 @@ void run_cdisort(Workspace& ws, get_angs(pfct_angs, scat_data, Npfct); nang = pfct_angs.nelem(); - //allocate + //allocate per-thread working matrices (firstprivate in OMP loop) Matrix ext_bulk_par(1, ds.nlyr + 1), abs_bulk_par(1, ds.nlyr + 1); Tensor3 pha_bulk_par(1, ds.nlyr + 1, nang); - Tensor3 pfct_bulk_par(1, ds.nlyr, nang); Tensor3 pmom(1, ds.nlyr, Nlegendre); - Vector f_grid_i(1); Matrix ssalb(1, ds.nlyr); Matrix ext_bulk_gas_i(1, ds.nlyr + 1); Matrix dtauc(1, ds.nlyr); - Matrix sca_coeff_gas_layer; - Matrix sca_bulk_par_layer; - Matrix sca_coeff_gas_level; - Matrix pmom_gas; - if (gas_scattering_do) { - sca_bulk_par_layer.resize(1, ds.nlyr); - sca_coeff_gas_layer.resize(1, ds.nlyr); - sca_coeff_gas_level.resize(1, ds.nlyr + 1); - pmom_gas.resize(ds.nlyr, Nlegendre); - } Matrix directbeam(nf, N_lev, 0); Matrix deltatau(nf, N_lev - 1, 0); Matrix snglsctalbedo(nf, N_lev - 1, 0); @@ -1134,8 +1105,49 @@ void run_cdisort(Workspace& ws, nlinspace(pfct_angs, 0, 180, nang); } + ArrayOfString fail_msg; + bool do_abort = false; + + WorkspaceOmpParallelCopyGuard wss{ws}; // loop over all frequencies +#pragma omp parallel for if (!arts_omp_in_parallel() && f_grid.nelem() > 1) \ + firstprivate(wss, \ + ds, \ + ext_bulk_gas_i, \ + ext_bulk_par, \ + abs_bulk_par, \ + pha_bulk_par, \ + umu0, \ + pmom, \ + ssalb, \ + dtauc, \ + out) for (Index f_index = 0; f_index < f_grid.nelem(); f_index++) { + + if (do_abort) { + continue; + } + + Vector f_grid_i(1); + + //Intensity of incident sun beam + Numeric fbeam = 0.; + + /* Allocate memory */ + c_disort_state_alloc(&ds); + c_disort_out_alloc(&ds, &out); + + // fill up azimuth angle and temperature array + for (Index i = 0; i < ds.nphi; i++) ds.phi[i] = aa_grid[i]; + + if (ds.flag.planck == TRUE) { + for (Index i = 0; i <= ds.nlyr; i++) + ds.temper[i] = temperature[ds.nlyr - i]; + } + + // Transform to mu, starting with negative values + for (Index i = 0; i < ds.numu; i++) ds.umu[i] = -cos(za_grid[i] * PI / 180); + f_grid_i = f_grid[f_index]; // Get particle bulk properties @@ -1179,6 +1191,7 @@ void run_cdisort(Workspace& ws, f_index); } + Tensor3 pfct_bulk_par(1, ds.nlyr, nang); get_pfct( pfct_bulk_par, pha_bulk_par, ext_bulk_par, abs_bulk_par, cboxlims); @@ -1195,10 +1208,18 @@ void run_cdisort(Workspace& ws, // gas scattering // layer averaged particle scattering coefficient + Matrix sca_coeff_gas_layer(1, ds.nlyr); + Matrix sca_bulk_par_layer(1, ds.nlyr); + Matrix sca_coeff_gas_level(1, ds.nlyr + 1); + Matrix pmom_gas(ds.nlyr, Nlegendre); get_scat_bulk_layer(sca_bulk_par_layer, ext_bulk_par, abs_bulk_par); // call gas_scattering_properties - get_gas_scattering_properties(ws, + sca_coeff_gas_layer = 0; + sca_coeff_gas_level = 0; + pmom_gas = 0; + + get_gas_scattering_properties(wss, sca_coeff_gas_layer, sca_coeff_gas_level, pmom_gas, @@ -1297,10 +1318,17 @@ void run_cdisort(Workspace& ws, ds.bc.umu0 = umu0; tries = Status::RETRY; - } else - throw e; + } else { +#pragma omp critical(run_cdisort_setabort) + do_abort = true; + + ostringstream os; + os << "failure at f_index " << f_index << ":\n" << e.what(); +#pragma omp critical(run_cdisort_push_fail_msg) + fail_msg.push_back(os.str()); + } } - } while (tries != Status::SUCCESS); + } while (tries != Status::SUCCESS and not do_abort); for (Index i = 0; i < ds.nphi; i++) { for (Index j = 0; j < ds.numu; j++) { @@ -1340,6 +1368,21 @@ void run_cdisort(Workspace& ws, exp(-dtauc(0, ds.nlyr - k + cboxlims[0]) / umu0); } } + + /* Free allocated memory */ + c_disort_out_free(&ds, &out); + c_disort_state_free(&ds); + } + + if (fail_msg.nelem()) { + ostringstream os; + + if (!do_abort) os << "\nError messages from failed run_cdisort:\n"; + for (ArrayOfString::const_iterator it = fail_msg.begin(); + it != fail_msg.end(); + it++) + os << *it << '\n'; + throw runtime_error(os.str()); } // Allocate aux data @@ -1372,9 +1415,6 @@ void run_cdisort(Workspace& ws, } } - /* Free allocated memory */ - c_disort_out_free(&ds, &out); - c_disort_state_free(&ds); #else ARTS_USER_ERROR("Did not compile with -DENABLE_ARTS_LGPL=0") From 01a4904a7eeb005d4f4d5fc923bb5e2eda7f51d2 Mon Sep 17 00:00:00 2001 From: Manfred Brath Date: Thu, 23 Apr 2026 11:47:11 +0200 Subject: [PATCH 2/3] Bugfix AngularGridsSetFluxCalc * Adjust za_grid_weights calculation to include PI/2 factor --- src/m_fluxes.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/m_fluxes.cc b/src/m_fluxes.cc index 0643b4ed2d..4abf2baa23 100644 --- a/src/m_fluxes.cc +++ b/src/m_fluxes.cc @@ -135,7 +135,7 @@ void AngularGridsSetFluxCalc(Vector& za_grid, // set the weights to the right component // by adjusting the domain, we also have to adjust the weights - za_grid_weights[i] = w[i] * sin(za_grid[i] * DEG2RAD); + za_grid_weights[i] = w[i] * sin(za_grid[i] * DEG2RAD)*PI/2; } } else if (za_grid_type == "linear_mu") { Vector x; From ce43ce46de7b16c7de12769bd797289bcce31fe1 Mon Sep 17 00:00:00 2001 From: Oliver Lemke Date: Fri, 24 Apr 2026 09:17:54 +0200 Subject: [PATCH 3/3] Minor fixes Unique name for critical section. Loop indentation. --- src/disort.cc | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/disort.cc b/src/disort.cc index 9488bdcb52..54b64025ad 100644 --- a/src/disort.cc +++ b/src/disort.cc @@ -1382,7 +1382,7 @@ void run_cdisort(Workspace& ws, it != fail_msg.end(); it++) os << *it << '\n'; - throw runtime_error(os.str()); + throw runtime_error(os.str()); } // Allocate aux data @@ -1818,7 +1818,7 @@ void run_cdisort_flux(Workspace& ws, ostringstream os; os << "failure at f_index " << f_index << ":\n" << e.what(); -#pragma omp critical(ybatchCalc_push_fail_msg) +#pragma omp critical(run_cdisort_flux_push_fail_msg) fail_msg.push_back(os.str()); } } @@ -1891,7 +1891,7 @@ void run_cdisort_flux(Workspace& ws, it != fail_msg.end(); it++) os << *it << '\n'; - throw runtime_error(os.str()); + throw runtime_error(os.str()); } // Allocate aux data