Replace printout macro with function (to prepare for C++23 std::format)

.clang-tidy

@@ -16,6 +16,7 @@ Checks: >
     -bugprone-narrowing-conversions,
     -boost-use-ranges,
     -llvmlibc-*,
+    -cert-dcl50-cpp,
     -cert-err33-c,
     -cert-err34-c,
     -cert-err58-cpp,

atomic.h

@@ -230,9 +230,6 @@ inline auto get_elementindex(const int Z) -> int {
     return static_cast<int>(elem - globals::elements.begin());
   }
 
-  // printout("[debug] get_elementindex: element Z=%d was not found in atomic data ... skip readin of cross sections
-  // for this element\n",Z); printout("[fatal] get_elementindex: element Z=%d was not found in atomic data ...
-  // abort\n"); abort();;
   return -100;
 }
 

decay.cc

@@ -973,9 +973,6 @@ void init_nuclides(const std::vector<int> &custom_zlist, const std::vector<int>
         nuclides[alphanucindex].endecay_q[DECAYTYPE_BETAMINUS] = Q_total_betadec * MEV;
         nuclides[alphanucindex].branchprobs[DECAYTYPE_ALPHA] = branch_alpha;
         nuclides[alphanucindex].endecay_q[DECAYTYPE_ALPHA] = Q_total_alphadec * MEV;
-
-        // printout("alphadecay file: Adding (Z=%d)%s%d endecay_alpha %g endecay_gamma %g tau_s %g\n",
-        //          z, get_elname(z), a, e_alpha_mev, e_gamma_mev, tau_sec);
       }
     }
   }

exspec.cc

@@ -82,7 +82,6 @@ void do_angle_bin(const int a, Packet *pkts, bool load_allrank_packets, Spectra
     int nesc_gamma = 0;
     int nesc_rpkt = 0;
     for (int ii = 0; ii < globals::npkts; ii++) {
-      // printout("packet %d escape_type %d type %d", ii, pkts[ii].escape_type, pkts[ii].type);
       if (pkts_start[ii].type == TYPE_ESCAPE) {
         nesc_tot++;
         if (pkts_start[ii].escape_type == TYPE_RPKT) {

gammapkt.cc

@@ -296,13 +296,10 @@ auto choose_f(const double xx, const double zrand) -> double
   int count = 0;
   double err = 1e20;
 
-  // printout("new\n");
-
   double ftry = (f_max + f_min) / 2;
   while ((err > 1.e-4) && (count < 1000)) {
     ftry = (f_max + f_min) / 2;
     const double sigma_try = sigma_compton_partial(xx, ftry);
-    // printout("ftry %g %g %g %g %g\n",ftry, f_min, f_max, try, norm);
     if (sigma_try > norm) {
       f_max = ftry;
       err = (sigma_try - norm) / norm;
@@ -780,8 +777,6 @@ void transport_gamma(Packet &pkt, const double t2) {
 
   assert_always(tdist >= 0);
 
-  // printout("sdist, tdist, edist %g %g %g\n",sdist, tdist, edist);
-
   if ((sdist < tdist) && (sdist < edist)) {
     move_pkt_withtime(pkt, sdist / 2.);
 
@@ -1056,9 +1051,6 @@ __host__ __device__ void pellet_gamma_decay(Packet &pkt) {
   }
 
   pkt.pol_dir = vec_norm(pkt.pol_dir);
-  // printout("initialise pol state of packet %g, %g, %g, %g,
-  // %g\n",pkt.stokes_qu[0],pkt.stokes_qu[1],pkt.pol_dir[0],pkt.pol_dir[1],pkt.pol_dir[2]);
-  // printout("pkt direction %g, %g, %g\n",pkt.dir[0],pkt.dir[1],pkt.dir[2]);
 }
 
 __host__ __device__ void do_gamma(Packet &pkt, const int nts, const double t2) {

grid.cc

@@ -111,9 +111,6 @@ void read_possible_yefile() {
     const int mgi = mgiplusone - 1;
     if (mgi >= 0 && mgi < get_npts_model()) {
       set_initelectronfrac(mgi, initelecfrac);
-      // printout("Ye.txt: setting mgi %d init_ye %g\n", mgi, initelecfrac);
-    } else {
-      // printout("Ye.txt: ignoring mgi %d init_ye %g\n", mgi, initelecfrac);
     }
   }
   fclose(filein);
@@ -1021,8 +1018,6 @@ void read_3d_model() {
     std::istringstream ssline(line);
 
     assert_always(ssline >> cellnumberin >> cellpos_in[0] >> cellpos_in[1] >> cellpos_in[2] >> rho_model);
-    // printout("cell %d, posz %g, posy %g, posx %g, rho %g, rho_init %g\n",dum1,dum3,dum4,dum5,rho_model,rho_model*
-    // pow( (t_model/globals::tmin), 3.));
 
     if (mgi == 0) {
       first_cellindex = cellnumberin;
@@ -2379,7 +2374,6 @@ void grid_init(const int my_rank) {
 
       if (totmassradionuclide_actual > 0.) {
         const double ratio = totmassradionuclide[nucindex] / totmassradionuclide_actual;
-        // printout("nuclide %d ratio %g\n", nucindex, ratio);
         for (int nonemptymgi = 0; nonemptymgi < get_nonempty_npts_model(); nonemptymgi++) {
           const int mgi = grid::get_mgi_of_nonemptymgi(nonemptymgi);
           const double prev_abund = get_modelinitnucmassfrac(mgi, nucindex);
@@ -2470,7 +2464,6 @@ auto get_totmassradionuclide(const int z, const int a) -> double {
   const auto negdirections = std::array<enum cell_boundary, 3>{COORD0_MIN, COORD1_MIN, COORD2_MIN};
   const auto posdirections = std::array<enum cell_boundary, 3>{COORD0_MAX, COORD1_MAX, COORD2_MAX};
 
-  // printout("checking inside cell boundary\n");
   for (int d = 0; d < ndim; d++) {
     // flip is either zero or one to indicate +ve and -ve boundaries along the selected axis
     for (int flip = 0; flip < 2; flip++) {
@@ -2526,14 +2519,6 @@ auto get_totmassradionuclide(const int z, const int a) -> double {
     }
   }
 
-  // printout("pkt.number %d\n", pkt.number);
-  // printout("delta1x %g delta2x %g\n",  (initpos[0] * globals::tmin/tstart)-grid::get_cellcoordmin(cellindex, 0),
-  // cellcoordmax[0] - (initpos[0] * globals::tmin/tstart)); printout("delta1y %g delta2y %g\n",  (initpos[1] *
-  // globals::tmin/tstart)-grid::get_cellcoordmin(cellindex, 1), cellcoordmax[1] - (initpos[1] *
-  // globals::tmin/tstart)); printout("delta1z %g delta2z %g\n",  (initpos[2] *
-  // globals::tmin/tstart)-grid::get_cellcoordmin(cellindex, 2), cellcoordmax[2] - (initpos[2] *
-  // globals::tmin/tstart)); printout("dir [%g, %g, %g]\n", dir[0],dir[1],dir[2]);
-
   // distance to reach the cell's upper boundary on each coordinate
   auto d_coordmaxboundary = std::array<double, 3>{-1};
 

input.cc

@@ -219,12 +219,9 @@ void read_phixs_data_table(std::fstream &phixsfile, const int nphixspoints_input
       // the photoionisation cross-sections in the database are given in Mbarn = 1e6 * 1e-28m^2
       // to convert to cgs units multiply by 1e-18
       levelphixstable[i] = phixs * 1e-18;
-      // fprintf(database_file,"%g %g\n", nutable[i], phixstable[i]);
     }
   }
 
-  // nbfcontinua++;
-  // printout("[debug] element %d, ion %d, level %d: phixs exists %g\n",element,lowerion,lowerlevel,phixs*1e-18);
   globals::nbfcontinua += get_nphixstargets(element, lowerion, lowerlevel);
   if (lowerlevel == 0 && get_nphixstargets(element, lowerion, lowerlevel) > 0) {
     globals::nbfcontinua_ground++;
@@ -446,8 +443,6 @@ void read_ion_transitions(std::fstream &ftransitiondata, const int tottransition
           if (tmplevel == prev_lower) {
             continue;
           }
-          // printout("+adding transition index %d Z=%02d ionstage %d lower %d upper %d\n", i, Z, ionstage, prev_lower,
-          // tmplevel);
           tottransitions++;
           assert_always(tmplevel >= 0);
           iontransitiontable.push_back(
@@ -457,9 +452,6 @@ void read_ion_transitions(std::fstream &ftransitiondata, const int tottransition
 
       iontransitiontable.push_back(
           {.lower = lower, .upper = upper, .A = A, .coll_str = coll_str, .forbidden = (intforbidden == 1)});
-      // printout("index %d, lower %d, upper %d, A %g\n",transitionindex,lower,upper,A);
-      //  printout("reading transition index %d lower %d upper %d\n", i, transitiontable[i].lower,
-      //  transitiontable[i].upper);
       prev_lower = lower;
       prev_upper = upper;
     }
@@ -1017,7 +1009,6 @@ void read_atomicdata_files() {
     double ionpot = 0.;
     for (int ion = 0; ion < nions; ion++) {
       int nlevelsmax = nlevelsmax_readin;
-      // printout("element %d ion %d\n", element, ion);
       // calculate the current levels ground level energy
       assert_always(ionpot >= 0);
       energyoffset += ionpot;
@@ -1646,8 +1637,6 @@ void input(int rank) {
 auto get_noncommentline(std::fstream &input, std::string &line) -> bool {
   while (true) {
     const bool linefound = !(!std::getline(input, line));
-    // printout("LINE: >%s<  linefound: %s commentonly: %s \n", line.c_str(), linefound ? "true" : "false",
-    // lineiscommentonly(line) ? "true" : "false");
     if (!linefound) {
       return false;
     }

ltepop.cc

@@ -185,7 +185,6 @@ auto calculate_levelpop_nominpop(const int modelgridindex, const int element, co
         // Case for when no NLTE level information is available yet
         nn = calculate_levelpop_lte(modelgridindex, element, ion, level);
       } else {
-        // printout("Using an nlte population!\n");
         nn = nltepop_over_rho * grid::get_rho(modelgridindex);
         if (!std::isfinite(nn)) {
           printout("[fatal] NLTE population failure.\n");
@@ -209,7 +208,6 @@ auto calculate_levelpop_nominpop(const int modelgridindex, const int element, co
         // Case for when no NLTE level information is available yet
         nn = calculate_levelpop_lte(modelgridindex, element, ion, level);
       } else {
-        // printout("Using a superlevel population!\n");
         nn = superlevelpop_over_rho * grid::get_rho(modelgridindex) *
              superlevel_boltzmann(modelgridindex, element, ion, level);
         if (!std::isfinite(nn)) {
@@ -544,8 +542,6 @@ __host__ __device__ auto calculate_sahafact(const int element, const int ion, co
   const double g_lower = stat_weight(element, ion, level);
   const double g_upper = stat_weight(element, ion + 1, upperionlevel);
   const double sf = SAHACONST * g_lower / g_upper * pow(T, -1.5) * exp(E_threshold / KB / T);
-  // printout("element %d, ion %d, level %d, T, %g, E %g has sf %g (g_l %g g_u %g)\n", element, ion, level, T,
-  // E_threshold, sf,stat_weight(element,ion,level),stat_weight(element,ion+1,0) );
   if (sf < 0) {
     printout(
         "[fatal] calculate_sahafact: Negative Saha factor. sfac %g element %d ion %d level %d upperionlevel %d "

macroatom.cc

@@ -682,8 +682,8 @@ auto rad_deexcitation_ratecoeff(const int modelgridindex, const int element, con
     }
 
     // printout("[debug] rad_rates_down: Z=%d, ionstage %d, upper %d, lower %d\n", get_atomicnumber(element),
-    // get_ionstage(element, ion), upper, lower); printout("[debug] rad_deexc: A_ul %g, tau_sobolev %g, n_u %g\n", A_ul,
-    // tau_sobolev, n_u);
+    // get_ionstage(element, ion), upper, lower); printout("[debug] rad_deexc: A_ul %g, tau_sobolev %g, n_u %g\n",
+    // A_ul, tau_sobolev, n_u);
     assert_testmodeonly(std::isfinite(R));
   }
 

nltepop.cc

@@ -1175,11 +1175,6 @@ void nltepop_write_to_file(const int modelgridindex, const int timestep) {
           }
 
           nnlevelnlte = slpopfactor * superlevel_partfunc;
-
-          // printout("nltepop_write_to_file: The Z=%d ionstage %d superlevel population is %g with rho %g and
-          // superlevel_partfunc %g Te %g scaled pop stored as %g\n", get_atomicnumber(element), get_ionstage(element,
-          // ion), nnlevelnlte, grid::modelgrid[modelgridindex].rho, superlevel_partfunc, grid::get_Te(modelgridindex),
-          // grid::modelgrid[modelgridindex].nlte_pops[ion_first_nlte + nlevels_nlte]);
         }
 
         const double ion_popfrac = nnlevelnlte / get_nnion(modelgridindex, element, ion);

nonthermal.cc

@@ -452,7 +452,6 @@ void read_auger_data() {
           assert_always(fabs(prob_sum - 1.0) < 0.001);
 
           printout("\n");
-          // printout("ionpot %g %g, g %d\n", colliondata[i].ionpot_ev, ionpot_ev, g);
           const bool found_existing_data = (collionrow.auger_g_accumulated > 0.);
 
           // keep existing data but update according to statistical weight represented by existing and new data
@@ -642,9 +641,6 @@ void read_collion_data() {
     collionrow.n_auger_elec_avg = 0.;
 
     colliondata.push_back(collionrow);
-
-    // printout("ci row: %2d %2d %1d %1d %lg %lg %lg %lg %lg\n", collionrow.Z, collionrow.nelec, collionrow.n,
-    //          collionrow.l, collionrow.ionpot_ev, collionrow.A, collionrow.B, collionrow.C, collionrow.D);
   }
   printout("Stored %zu of %d input shell cross sections\n", colliondata.size(), colliondatacount);
   for (int element = 0; element < get_nelements(); element++) {
@@ -1900,9 +1896,6 @@ void sfmatrix_add_ionization(std::vector<double> &sfmatrixuppertri, const int Z,
 
       assert_always(ionpot_ev >= SF_EMIN);
 
-      // printout("Z=%2d ionstage %d n %d l %d ionpot %g eV\n",
-      //          Z, ionstage, colliondata[n].n, colliondata[n].l, ionpot_ev);
-
       const int xsstartindex = get_xs_ionization_vector(vec_xs_ionization, collionrow);
       assert_always(xsstartindex >= 0);
       // Luke Shingles: the use of min and max on the epsilon limits keeps energies
@@ -2012,24 +2005,22 @@ auto sfmatrix_solve(const std::vector<double> &sfmatrix) -> std::array<double, S
   const auto &gsl_sfmatrix_LU = gsl_sfmatrix;
 
   // if the matrix is not upper triangular, then do a decomposition
-  // printout("Doing LU decomposition of SF matrix\n");
   // make a copy of the matrix for the LU decomp
   // std::array<double, SFPTS> sfmatrix_LU{};
   // auto gsl_sfmatrix_LU = gsl_matrix_view_array(sfmatrix_LU.data(), SFPTS, SFPTS).matrix;
   // gsl_matrix_memcpy(&gsl_sfmatrix_LU, &gsl_sfmatrix);
   // int s{};  // sign of the transformation
   // gsl_linalg_LU_decomp(&gsl_sfmatrix_LU, &p, &s);
 
-  // printout("Solving SF matrix equation\n");
-
   std::array<double, SFPTS> yvec_arr{};
   auto gsl_yvec = gsl_vector_view_array(yvec_arr.data(), SFPTS).vector;
 
   const auto gsl_rhsvec = gsl_vector_const_view_array(rhsvec.data(), SFPTS).vector;
 
   // solve matrix equation: sf_matrix * y_vec = rhsvec for yvec
   gsl_linalg_LU_solve(&gsl_sfmatrix_LU, &p, &gsl_rhsvec, &gsl_yvec);
-  // printout("Refining solution\n");
+
+  // refine the solution
 
   double error_best = -1.;
   std::array<double, SFPTS> yvec_best{};
@@ -2057,7 +2048,7 @@ auto sfmatrix_solve(const std::vector<double> &sfmatrix) -> std::array<double, S
       gsl_vector_memcpy(&gsl_yvec_best, &gsl_yvec);
       error_best = error;
     }
-    // printout("Linear algebra solver iteration %d has a maximum residual of %g\n",iteration,error);
+    // printout("Linear algebra solver iteration %d has a maximum residual of %g\n", iteration, error);
   }
   if (error_best >= 0.) {
     if (error_best > 1e-10) {
@@ -2370,9 +2361,6 @@ __host__ __device__ auto nt_ionization_ratecoeff(const int modelgridindex, const
       // probably because eff_ionpot = 0 because the solver hasn't been run yet, or no impact ionization cross sections
       // exist
       const double Y_nt_wfapprox = nt_ionization_ratecoeff_wfapprox(modelgridindex, element, ion);
-      // printout("Warning: Spencer-Fano solver gives non-finite ionization rate (%g) for element %d ionstage %d for
-      // cell %d. Using WF approx instead = %g\n",
-      //          Y_nt, get_atomicnumber(element), get_ionstage(element, ion), modelgridindex, Y_nt_wfapprox);
       return Y_nt_wfapprox;
     }
     if (Y_nt <= 0) {
@@ -2474,9 +2462,6 @@ __host__ __device__ void do_ntlepton_deposit(Packet &pkt) {
                                    pkt.e_cmf);
       }
 
-      // printout("NTLEPTON packet in cell %d selected ionization of Z=%d ionstage %d to %d\n",
-      //          modelgridindex, get_atomicnumber(element), get_ionstage(element, lowerion), get_ionstage(element,
-      //          upperion));
       do_macroatom(pkt, {.element = element, .ion = upperion, .level = 0, .activatingline = -99});
       return;
     }
@@ -2506,9 +2491,6 @@ __host__ __device__ void do_ntlepton_deposit(Packet &pkt) {
 
           stats::increment(stats::COUNTER_NT_STAT_TO_EXCITATION);
 
-          // printout("NTLEPTON packet selected in cell %d excitation of Z=%d ionstage %d level %d upperlevel %d\n",
-          //          modelgridindex, get_atomicnumber(element), get_ionstage(element, ion), lower, upper);
-
           do_macroatom(pkt, {.element = element, .ion = ion, .level = upper, .activatingline = -99});
           return;
         }
@@ -2664,10 +2646,8 @@ void solve_spencerfano(const int modelgridindex, const int timestep, const int i
   }
 
   // printout("SF matrix | RHS vector:\n");
-  // for (int row = 0; row < 10; row++)
-  // {
-  //   for (int column = 0; column < 10; column++)
-  //   {
+  // for (int row = 0; row < 10; row++) {
+  //   for (int column = 0; column < 10; column++) {
   //     char str[15];
   //     snprintf(str, 15, "%+.1e ", gsl_matrix_get(sfmatrix, row, column));
   //     printout(str);
@@ -2803,11 +2783,6 @@ void nt_MPI_Bcast(const int modelgridindex, const int root, const int root_node_
     return;
   }
 
-  // printout("nonthermal_MPI_Bcast cell %d before: ratecoeff(Z=%d ionstage %d): %g, eff_ionpot %g eV\n",
-  //          modelgridindex, logged_element_z, logged_ionstage,
-  //          nt_ionization_ratecoeff_sf(modelgridindex, logged_element_index, logged_ion_index),
-  //          get_eff_ionpot(modelgridindex, logged_element_index, logged_ion_index) / EV);
-
   MPI_Bcast(&deposition_rate_density_all_cells[modelgridindex], 1, MPI_DOUBLE, root, MPI_COMM_WORLD);
 
   if (NT_ON && NT_SOLVE_SPENCERFANO) {

radfield.cc

@@ -187,7 +187,6 @@ void add_detailed_line(const int lineindex) {
   }
   detailed_lineindicies[detailed_linecount] = lineindex;
   detailed_linecount++;
-  // printout("Added Jblue estimator for lineindex %d count %d\n", lineindex, detailed_linecount);
 }
 
 // get the normalised J_nu
@@ -688,13 +687,6 @@ auto get_Jblueindex(const int lineindex) -> int {
     }
   }
 
-  // const int element = linelist[lineindex].elementindex;
-  // const int ion = linelist[lineindex].ionindex;
-  // const int lower = linelist[lineindex].lowerlevelindex;
-  // const int upper = linelist[lineindex].upperlevelindex;
-  // printout("Could not find lineindex %d among %d items (Z=%02d ionstage %d lower %d upper %d)\n",
-  //          lineindex, detailed_linecount, get_atomicnumber(element), get_ionstage(element, ion), lower, upper);
-
   return -1;
 }
 

ratecoeff.cc

@@ -536,8 +536,6 @@ void read_recombrate_file() {
   int tablerows = 0;
 
   while (fscanf(recombrate_file, "%d %d %d\n", &atomicnumber, &upperionstage, &tablerows) > 0) {
-    // printout("%d %d %d\n", atomicnumber, upperionstage, tablerows);
-
     RRCRow T_highestbelow = {.log_Te = 0, .rrc_low_n = 0, .rrc_total = 0};
     RRCRow T_lowestabove = {.log_Te = 0, .rrc_low_n = 0, .rrc_total = 0};
     T_highestbelow.log_Te = -1;
@@ -984,8 +982,6 @@ __host__ __device__ auto get_spontrecombcoeff(int element, const int ion, const
 
     const double f_upper = spontrecombcoeffs[get_bflutindex(upperindex, element, ion, level, phixstargetindex)];
     const double f_lower = spontrecombcoeffs[get_bflutindex(lowerindex, element, ion, level, phixstargetindex)];
-    // printout("interpolate_spontrecombcoeff element %d, ion %d, level %d, upper %g, lower %g\n",
-    //          element,ion,level,f_upper,f_lower);
     Alpha_sp = (f_lower + (f_upper - f_lower) / (T_upper - T_lower) * (T_e - T_lower));
   } else {
     Alpha_sp = spontrecombcoeffs[get_bflutindex(TABLESIZE - 1, element, ion, level, phixstargetindex)];
@@ -1310,7 +1306,6 @@ auto iongamma_is_zero(const int nonemptymgi, const int element, const int ion) -
       }
 
       const double epsilon_trans = epsilon(element, ion + 1, upperlevel) - epsilon(element, ion, level);
-      // printout("%g %g %g\n", get_levelpop(n,element,ion,level),col_ionization(n,0,epsilon_trans),epsilon_trans);
 
       if (nnlevel * col_ionization_ratecoeff(T_e, nne, element, ion, level, phixstargetindex, epsilon_trans) > 0) {
         return false;
@@ -1348,8 +1343,6 @@ auto calculate_iongamma_per_gspop(const int modelgridindex, const int element, c
       Col_ion += nnlevel * col_ionization_ratecoeff(T_e, nne, element, ion, level, phixstargetindex, epsilon_trans);
     }
   }
-  // printout("element %d ion %d: col/gamma %g Te %g ne %g\n", element, ion, Col_ion/Gamma, grid::get_Te(n),
-  // grid::get_nne(n));
   Gamma += Col_ion;
   Gamma /= get_groundlevelpop(modelgridindex, element, ion);
   return Gamma;