Merge pull request sirocco-rt#1128 from sirocco-rt/dense

Handling of low frequency photons in macro mode

Author: jhmatthews

README.md

@@ -10,7 +10,7 @@
 The code is under active development, but we are looking for beta users to test the code, and potentially use it for their own research. If you are interested in using Sirocco please contact Knox Long via long[at]stsci[dot]edu 
 or James Matthews via james[dot]matthews[at]physics[dot]ox[dot]ac[dot]uk.
 
-Documentation is hosted on [ReadTheDocs]((https://sirocco-rt.readthedocs.io).
+Documentation is hosted on [ReadTheDocs](https://sirocco-rt.readthedocs.io).
 
 ## Installation
 

py_progs/plot_wind.py

@@ -140,6 +140,12 @@ def get_data(filename='fiducial_agn_master.txt', var='t_r',grid='ij',inwind='',s
     This routine reads and scales the data from a single variable that is read from an ascii table 
     representation of one or more of the parameters in a windsave file
 
+    The grid obptions:
+
+    ij
+    log
+    lin
+
     '''
 
     try:
@@ -182,8 +188,6 @@ def get_data(filename='fiducial_agn_master.txt', var='t_r',grid='ij',inwind='',s
         ylogmin=numpy.log10(xmin/10)
         x=numpy.select([x>1],[numpy.log10(x)],default=xlogmin)
         y=numpy.select([y>1],[numpy.log10(y)],default=ylogmin)
-        # x=numpy.log10(x)
-        # y=numpy.log10(y)
         xlabel='log(x)'
         ylabel='log(z)'
     else:
@@ -333,15 +337,51 @@ def doit(filename='fiducial_agn.master.txt', var='t_r',grid='ij',inwind='',scale
     return plotfile
 
 
+def steer(argv):
+    '''
+    Controlling routine for py_wind
+    '''
 
+    xfile=''
+    xvar=''
+    xgrid='ij'
+    xmin=-1e50
+    xmax=1e50
+    xscale='guess'
+
+    i=1
+    while i<len(argv):
+        if argv[i][0:2]=='-h':
+            print(__doc__)
+            return
+        elif argv[i]=='-log':
+            xgrid='log'
+        elif argv[i][0:4]=='-lin':
+            xgrid='linear'
+        elif argv[i][0:4]=='-min':
+            i+=1
+            xmin=eval(argv[i])
+        elif argv[i][0:4]=='-max':
+            i+=1
+            xmax=eval(argv[i])
+        elif argv[i][0]=='-':
+            print('Unknown options :',argv)
+            return
+        elif xfile=='':
+            xfile=argv[i]
+        elif xvar=='':
+            xvar=argv[i]
+        i+=1
 
 
+    doit(filename=xfile, var=xvar,grid=xgrid,inwind='',scale='log',zmin=xmin,zmax=xmax,
+        plot_dir='',root='')
 
 
 # Next lines permit one to run the routine from the command line
 if __name__ == "__main__":
     import sys
     if len(sys.argv)>2:
-        doit(sys.argv[1],sys.argv[2])
+         steer(sys.argv)
     else:
-        print('usage: plot_wind filename variable')
+        print(__doc__)

source/communicate_plasma.c

@@ -50,8 +50,13 @@ broadcast_plasma_grid (const int n_start, const int n_stop, const int n_cells_ra
 
   d_xsignal (files.root, "%-20s Begin communicating plasma grid\n", "NOK");
   const int n_cells_max = get_max_cells_per_rank (NPLASMA);
+//OLD  const int comm_buffer_size = calculate_comm_buffer_size (1 + n_cells_max * (1 + 20 + nphot_total + nions + NXBANDS + 2 * N_PHOT_PROC),
+//OLD                                                           n_cells_max * (71 + 11 * nions + nlte_levels + 2 * nphot_total + n_inner_tot +
+//OLD                                                                          11 * NXBANDS + NBINS_IN_CELL_SPEC + 6 * NFLUX_ANGLES +
+//OLD                                                                          N_DMO_DT_DIRECTIONS + 12 * NFORCE_DIRECTIONS));
+
   const int comm_buffer_size = calculate_comm_buffer_size (1 + n_cells_max * (1 + 20 + nphot_total + nions + NXBANDS + 2 * N_PHOT_PROC),
-                                                           n_cells_max * (71 + 11 * nions + nlte_levels + 2 * nphot_total + n_inner_tot +
+                                                           n_cells_max * (73 + 11 * nions + nlte_levels + 2 * nphot_total + n_inner_tot +
                                                                           11 * NXBANDS + NBINS_IN_CELL_SPEC + 6 * NFLUX_ANGLES +
                                                                           N_DMO_DT_DIRECTIONS + 12 * NFORCE_DIRECTIONS));
   char *comm_buffer = malloc (comm_buffer_size);
@@ -117,6 +122,8 @@ broadcast_plasma_grid (const int n_start, const int n_stop, const int n_cells_ra
         MPI_Pack (&cell->ntot_agn, 1, MPI_INT, comm_buffer, comm_buffer_size, &position, MPI_COMM_WORLD);
         MPI_Pack (&cell->nscat_es, 1, MPI_INT, comm_buffer, comm_buffer_size, &position, MPI_COMM_WORLD);
         MPI_Pack (&cell->nscat_res, 1, MPI_INT, comm_buffer, comm_buffer_size, &position, MPI_COMM_WORLD);
+        MPI_Pack (&cell->nscat_bf, 1, MPI_INT, comm_buffer, comm_buffer_size, &position, MPI_COMM_WORLD);
+        MPI_Pack (&cell->nscat_ff, 1, MPI_INT, comm_buffer, comm_buffer_size, &position, MPI_COMM_WORLD);
         MPI_Pack (&cell->mean_ds, 1, MPI_DOUBLE, comm_buffer, comm_buffer_size, &position, MPI_COMM_WORLD);
         MPI_Pack (&cell->n_ds, 1, MPI_INT, comm_buffer, comm_buffer_size, &position, MPI_COMM_WORLD);
         MPI_Pack (&cell->nrad, 1, MPI_INT, comm_buffer, comm_buffer_size, &position, MPI_COMM_WORLD);
@@ -272,6 +279,8 @@ broadcast_plasma_grid (const int n_start, const int n_stop, const int n_cells_ra
         MPI_Unpack (comm_buffer, comm_buffer_size, &position, &cell->ntot_agn, 1, MPI_INT, MPI_COMM_WORLD);
         MPI_Unpack (comm_buffer, comm_buffer_size, &position, &cell->nscat_es, 1, MPI_INT, MPI_COMM_WORLD);
         MPI_Unpack (comm_buffer, comm_buffer_size, &position, &cell->nscat_res, 1, MPI_INT, MPI_COMM_WORLD);
+        MPI_Unpack (comm_buffer, comm_buffer_size, &position, &cell->nscat_bf, 1, MPI_INT, MPI_COMM_WORLD);
+        MPI_Unpack (comm_buffer, comm_buffer_size, &position, &cell->nscat_ff, 1, MPI_INT, MPI_COMM_WORLD);
         MPI_Unpack (comm_buffer, comm_buffer_size, &position, &cell->mean_ds, 1, MPI_DOUBLE, MPI_COMM_WORLD);
         MPI_Unpack (comm_buffer, comm_buffer_size, &position, &cell->n_ds, 1, MPI_INT, MPI_COMM_WORLD);
         MPI_Unpack (comm_buffer, comm_buffer_size, &position, &cell->nrad, 1, MPI_INT, MPI_COMM_WORLD);
@@ -593,10 +602,12 @@ broadcast_updated_plasma_properties (const int n_start_rank, const int n_stop_ra
 
   d_xsignal (files.root, "%-20s Begin communicating updated plasma properties\n", "NOK");
   const int n_cells_max = get_max_cells_per_rank (NPLASMA);
-  const int num_ints = 1 + n_cells_max * (20 + nphot_total + 2 * NXBANDS + 2 * N_PHOT_PROC + nions);
+  //OLD  const int num_ints = 1 + n_cells_max * (20 + nphot_total + 2 * NXBANDS + 2 * N_PHOT_PROC + nions);
+  const int num_ints = 1 + n_cells_max * (22 + nphot_total + 2 * NXBANDS + 2 * N_PHOT_PROC + nions);
   const int num_doubles =
     n_cells_max * (71 + 1 * 3 + 9 * 4 + 6 * NFLUX_ANGLES + 3 * NFORCE_DIRECTIONS + 9 * nions + 1 * nlte_levels + 3 * nphot_total +
                    1 * n_inner_tot + 9 * NXBANDS + 1 * NBINS_IN_CELL_SPEC);
+
   const int size_of_comm_buffer = calculate_comm_buffer_size (num_ints, num_doubles);
   char *const comm_buffer = malloc (size_of_comm_buffer);
   if (comm_buffer == NULL)
@@ -662,6 +673,8 @@ broadcast_updated_plasma_properties (const int n_start_rank, const int n_stop_ra
         MPI_Pack (&plasmamain[n_plasma].ntot_wind, 1, MPI_INT, comm_buffer, size_of_comm_buffer, &position, MPI_COMM_WORLD);
         MPI_Pack (&plasmamain[n_plasma].ntot_agn, 1, MPI_INT, comm_buffer, size_of_comm_buffer, &position, MPI_COMM_WORLD);
         MPI_Pack (&plasmamain[n_plasma].nscat_es, 1, MPI_INT, comm_buffer, size_of_comm_buffer, &position, MPI_COMM_WORLD);
+        MPI_Pack (&plasmamain[n_plasma].nscat_bf, 1, MPI_INT, comm_buffer, size_of_comm_buffer, &position, MPI_COMM_WORLD);
+        MPI_Pack (&plasmamain[n_plasma].nscat_ff, 1, MPI_INT, comm_buffer, size_of_comm_buffer, &position, MPI_COMM_WORLD);
         MPI_Pack (&plasmamain[n_plasma].mean_ds, 1, MPI_DOUBLE, comm_buffer, size_of_comm_buffer, &position, MPI_COMM_WORLD);
         MPI_Pack (&plasmamain[n_plasma].n_ds, 1, MPI_INT, comm_buffer, size_of_comm_buffer, &position, MPI_COMM_WORLD);
         MPI_Pack (&plasmamain[n_plasma].nrad, 1, MPI_INT, comm_buffer, size_of_comm_buffer, &position, MPI_COMM_WORLD);
@@ -830,6 +843,8 @@ broadcast_updated_plasma_properties (const int n_start_rank, const int n_stop_ra
         MPI_Unpack (comm_buffer, size_of_comm_buffer, &position, &plasmamain[n_plasma].ntot_wind, 1, MPI_INT, MPI_COMM_WORLD);
         MPI_Unpack (comm_buffer, size_of_comm_buffer, &position, &plasmamain[n_plasma].ntot_agn, 1, MPI_INT, MPI_COMM_WORLD);
         MPI_Unpack (comm_buffer, size_of_comm_buffer, &position, &plasmamain[n_plasma].nscat_es, 1, MPI_INT, MPI_COMM_WORLD);
+        MPI_Unpack (comm_buffer, size_of_comm_buffer, &position, &plasmamain[n_plasma].nscat_bf, 1, MPI_INT, MPI_COMM_WORLD);
+        MPI_Unpack (comm_buffer, size_of_comm_buffer, &position, &plasmamain[n_plasma].nscat_ff, 1, MPI_INT, MPI_COMM_WORLD);
         MPI_Unpack (comm_buffer, size_of_comm_buffer, &position, &plasmamain[n_plasma].mean_ds, 1, MPI_DOUBLE, MPI_COMM_WORLD);
         MPI_Unpack (comm_buffer, size_of_comm_buffer, &position, &plasmamain[n_plasma].n_ds, 1, MPI_INT, MPI_COMM_WORLD);
         MPI_Unpack (comm_buffer, size_of_comm_buffer, &position, &plasmamain[n_plasma].nrad, 1, MPI_INT, MPI_COMM_WORLD);

source/macro_accelerate.c

@@ -808,7 +808,7 @@ f_kpkt_emit_accelerate (xplasma, freq_min, freq_max)
 
   /* JM 1511 -- Fix for issue 187. We need band limits for free free packet
      generation (see call to one_ff below) */
-  ff_freq_min = xband.f1[0];
+  ff_freq_min = 0.0;  // since #1128 the low frequency limit is always zero for free-free
   ff_freq_max = ALPHA_FF * xplasma->t_e / H_OVER_K;
 
   /* ksl This is a Bandaid for when the temperatures are very low */
@@ -887,30 +887,35 @@ f_kpkt_emit_accelerate (xplasma, freq_min, freq_max)
   /* consult issues #187, #492 regarding free-free */
   penorm += eprbs = mplasma->cooling_ff + mplasma->cooling_ff_lofreq;
 
-  total_ff_lofreq = total_free (xplasma, xplasma->t_e, 0, ff_freq_min);
+  //total_ff_lofreq = total_free (xplasma, xplasma->t_e, 0, ff_freq_min);
+  //total_ff = total_free (xplasma, xplasma->t_e, ff_freq_min, ff_freq_max);
   total_ff = total_free (xplasma, xplasma->t_e, ff_freq_min, ff_freq_max);
 
   /*
+   * Calculate the band-limited free-free luminosity by incrementing penorm_band. 
    * Do not increment penorm_band when the total free-free luminosity is zero
    */
-
-  if (freq_min > ff_freq_min)
-  {
-    if (total_ff > 0)
-      penorm_band += total_free (xplasma, xplasma->t_e, freq_min, freq_max) / total_ff * mplasma->cooling_ff;
-  }
-  else if (freq_max > ff_freq_min)
-  {
-    if (total_ff > 0)
-      penorm_band += total_free (xplasma, xplasma->t_e, ff_freq_min, freq_max) / total_ff * mplasma->cooling_ff;
-    if (total_ff_lofreq > 0)
-      penorm_band += total_free (xplasma, xplasma->t_e, freq_min, ff_freq_min) / total_ff_lofreq * mplasma->cooling_ff_lofreq;
-  }
-  else
-  {
-    if (total_ff_lofreq > 0)
-      penorm_band += total_free (xplasma, xplasma->t_e, freq_min, freq_max) / total_ff_lofreq * mplasma->cooling_ff_lofreq;
-  }
+  if (total_ff > 0)
+    penorm_band += total_free (xplasma, xplasma->t_e, freq_min, freq_max) / total_ff * mplasma->cooling_ff;
+
+  /* we used to have to do some extra stuff here due to low frequency free-free behaviour, see #1128 */
+  // if (freq_min > ff_freq_min)
+  // {
+  //   if (total_ff > 0)
+  //     penorm_band += total_free (xplasma, xplasma->t_e, freq_min, freq_max) / total_ff * mplasma->cooling_ff;
+  // }
+  // else if (freq_max > ff_freq_min)
+  // {
+  //   if (total_ff > 0)
+  //     penorm_band += total_free (xplasma, xplasma->t_e, ff_freq_min, freq_max) / total_ff * mplasma->cooling_ff;
+  //   if (total_ff_lofreq > 0)
+  //     penorm_band += total_free (xplasma, xplasma->t_e, freq_min, ff_freq_min) / total_ff_lofreq * mplasma->cooling_ff_lofreq;
+  // }
+  // else
+  // {
+  //   if (total_ff_lofreq > 0)
+  //     penorm_band += total_free (xplasma, xplasma->t_e, freq_min, freq_max) / total_ff_lofreq * mplasma->cooling_ff_lofreq;
+  // }
 
   penorm += eprbs = mplasma->cooling_adiabatic;
 

source/matom.c

@@ -1054,7 +1054,7 @@ kpkt (p, nres, escape, mode)
   }
   else if (destruction_choice < (mplasma->cooling_bftot + cooling_bbtot + mplasma->cooling_ff + mplasma->cooling_ff_lofreq))
   {
-    /*this is ff at a frequency that is so low frequency that it is not worth tracking further */
+/*this is ff at a frequency that is so low frequency that it is not worth tracking further */
     *escape = TRUE;
     *nres = NRES_FF;
     p->istat = P_LOFREQ_FF;
@@ -1065,7 +1065,7 @@ kpkt (p, nres, escape, mode)
   else if (destruction_choice <
            (mplasma->cooling_bftot + cooling_bbtot + mplasma->cooling_ff + mplasma->cooling_ff_lofreq + cooling_adiabatic))
   {
-    /* It is a k-packat that is destroyed by adiabatic cooling */
+/* It is a k-packat that is destroyed by adiabatic cooling */
 
     if (geo.adiabatic == 0 || mode == KPKT_MODE_CONTINUUM)
     {

source/photon2d.c

@@ -465,9 +465,23 @@ translate_in_wind (w, p, tau_scat, tau, nres)
     ds_current = calculate_ds (w, p, tau_scat, tau, nres, smax, &istat);
 
     if (p->nres == NRES_ES)
+    {
       xplasma->nscat_es++;
+    }
+    if (p->nres > NLINES)
+    {
+      xplasma->nscat_bf++;
+    }
+
     else if (p->nres > 0)
+    {
       xplasma->nscat_res++;
+    }
+    else if (p->nres == NRES_FF)
+    {
+      xplasma->nscat_ff++;
+    }
+
 
     /* We now increment the radiation field in the cell, translate the photon and wrap
      * things up.  For simple atoms, the routine radiation also reduces

source/resonate.c

@@ -425,14 +425,14 @@ calculate_ds (w, p, tau_scat, tau, nres, smax, istat)
  *
  *  * -1 implies electron scattering
  *  * -2 implies free-free
- *  * 0 or greater implies a specific photionization process was responsible (and 
- *  also that the program was operating in macro-atom mode.
+ *  * Greater that NLINES implies a specific photionization process was responsible  
  *
  * @details
  * This routine is called to determine which of several continuum proceseses
  * cause a photon to be scattered or absorbed.  In addition to electron scattering
  * and free-free absorption, the routine can identify which photoionization process
- * is implicated.
+ * is implicated assuming we are in macro-atom mode, when photoionization is considred
+ * as a scattering process..
  *
  * ### Notes ###
  *
@@ -443,10 +443,14 @@ calculate_ds (w, p, tau_scat, tau, nres, smax, istat)
  * which of the processes was responsible.  Data for the opacity due to
  * photonionization is passed remotely via the PlasmaPtr.
  *
- * In a program running in the two level approximation, only electron scattering
- * and ff and bf are treated as absorption processes.  In macro atom, ff and
- * photoionization are treated as a scattering 
- * processes.
+ * In a program running in the two level approximation, electron scattering
+ * and ff are treated as scattering processes; bf is treated as absorption processes
+ * and so in that case, the routine chooses between es and ff scatering.  
+ * In macro atom, ff and photoionization are treated as a scattering 
+ * processes, so any one of these can be returned
+ *
+ * In order to avoid confusion with bb processes, the bf option returns
+ * a value that is larger than NLINES.
  *
  **********************************************************/
 

source/sirocco.h

@@ -141,7 +141,7 @@ extern int NWAVE_NOW;         /**< Either NWAVE_IONIZ or NWAVE_EXTRACT depending
                                 */
 #define NWAVE_MIN   100       /**< The minimum number of wavelength bins in during spectral cycles
                                */
-#define MAXSCAT    2000
+#define MAXSCAT    20000
 
 /* Define the structures */
 #include "math_struc.h"
@@ -936,6 +936,8 @@ typedef struct plasma
 
   int nscat_es;                 /**<  The number of electrons scatters in the cell */
   int nscat_res;                /**<  The number of resonant line scatters in the cell */
+  int nscat_bf;                 /**< Number of bf scatters in the cell. (macro_only) */
+  int nscat_ff;                 /**< Number of ff scatters in the cell. (macro_only) */
 
   double mean_ds;               /**<  Mean photon path length in a cell. */
   int n_ds;                     /**<  Number of times a path lengyh was added; needed to compute mean_ds */

source/wind_updates2d.c

@@ -557,6 +557,8 @@ init_plasma_rad_properties (void)
     plasmamain[i].ntot_bl = 0;
     plasmamain[i].nscat_es = 0;
     plasmamain[i].nscat_res = 0;
+    plasmamain[i].nscat_bf = 0;
+    plasmamain[i].nscat_ff = 0;
     plasmamain[i].ntot_wind = 0;
     plasmamain[i].nrad = 0;
     plasmamain[i].nioniz = 0;

source/windsave2table_sub.c

@@ -154,8 +154,8 @@ create_master_table (ndom, rootname)
 {
   char filename[132];
   double *c[50], *converge;
-  char column_name[50][20];
-  char one_line[1024], start[1024], one_value[20];
+  char column_name[50][24];
+  char one_line[1024], start[1024], one_value[24];
   char name[132];               /* file name extension */
 
 
@@ -227,9 +227,21 @@ create_master_table (ndom, rootname)
   c[17] = get_one (ndom, "nioniz");
   strcpy (column_name[17], "nioniz");
 
+  c[18] = get_one (ndom, "nscat_es");
+  strcpy (column_name[18], "nscat_es");
+
+  c[19] = get_one (ndom, "nscat_res");
+  strcpy (column_name[19], "nscat_res");
+
+  c[20] = get_one (ndom, "nscat_ff");
+  strcpy (column_name[20], "nscat_ff");
+
+  c[21] = get_one (ndom, "nscat_bf");
+  strcpy (column_name[21], "nscat_bf");
+
 
   /* This should be the maxium number above +1 */
-  ncols = 18;
+  ncols = 22;
 
 
   converge = get_one (ndom, "converge");
@@ -1479,6 +1491,22 @@ get_one (ndom, variable_name)
       {
         x[n] = plasmamain[nplasma].nioniz;
       }
+      else if (strcmp (variable_name, "nscat_es") == 0)
+      {
+        x[n] = plasmamain[nplasma].nscat_es;
+      }
+      else if (strcmp (variable_name, "nscat_res") == 0)
+      {
+        x[n] = plasmamain[nplasma].nscat_res;
+      }
+      else if (strcmp (variable_name, "nscat_bf") == 0)
+      {
+        x[n] = plasmamain[nplasma].nscat_bf;
+      }
+      else if (strcmp (variable_name, "nscat_ff") == 0)
+      {
+        x[n] = plasmamain[nplasma].nscat_ff;
+      }
       else if (strcmp (variable_name, "heat_shock") == 0)
       {
         x[n] = plasmamain[nplasma].heat_shock;