discretize SOPHIA

[MHoerbe]

currently discretizing the SOPHIA code, aiming to fix the saturation of the PhotoPionProduction beyond 8 cores.

[adundovi]

I once had started to do something similar, but never finished. Maybe the following will help.

This is Makefile that generates a python interface for SOPHIA (you need f2py):

F2PY=f2py

all: sophia.so

sophia.so: sophia.pyf sophia_interface.f
        $(F2PY) -c sophia_interface.f sophia.pyf

clean:
        rm sophia.so

And this is my naive approach to generate tables of events

from sophia import *
import numpy as np
import sys

import itertools as it
from multiprocessing import Pool

# parallelize the work
max_cores = 4

# input
nature = [0,1] # proton, neutron
Ein_div = 10000
Ein = {'CMB': 3.72*np.logspace(9, 13, num=Ein_div),
       'IRB': 5.83*np.logspace(6, 12, num=Ein_div) }  # Gev

redshift = np.array([0, 0.01, 0.03, 0.05, 0.1, 0.2, 0.3, 0.4,
                     0.5, 0.6, 0.8, 1.0, 1.25, 1.5, 2.0, 2.5,
                     3.0, 4.0, 5.0, 6.0]) # usual redshifts
backgrounds = { 'CMB': 1,'IRB': 2 }
maxRedshift = 100 # hardcoded in CRPropa's source

def sevent( n, E, z, bg ):
    """Wrapper for ``sophiaevent``.
       Parameters
       ----------
       nature : input int
       ein : input float
       outpart : in/output rank-2 array('d') with bounds (2000,5)
       outparttype : in/output rank-1 array('i') with bounds (2000)
       nboutpart : in/output rank-0 array(int,'i')
       z_particle : input float
       bgflag : input int
       zmax_irb : input float
       en_data : input rank-1 array('d') with bounds (idatamax)
       flux_data : input rank-1 array('d') with bounds (idatamax)
    """

    # output variables
    particleList = np.zeros( 2000, dtype=np.int32 )
    momentaList = np.zeros( (2000,5), dtype=np.float_, order='F' )
    nParticles = np.array(0) # number of particles, scalar field

    # dummy variables
    dummy1 = np.int_()
    dummy2 = np.zeros(10, dtype=np.float_)
    dummy3 = np.zeros(10, dtype=np.float_)

    sophiaevent( n, E, momentaList, particleList, nParticles, z, bg, maxRedshift, dummy2, dummy3 )

    return '%d\t%.6e\t%.2f\t%d\t%s\t%s\t' % (n, E, z, nParticles, np.array2string(momentaList[0:nParticles,3], separator=','), np.array2string(particleList[0:nParticles], separator=',') )

def swrapper( args ):
  """ wrapper function because Pool.map() doesn't accept multiple arguments
      see http://stackoverflow.com/a/5443941
  """
  return sevent( *args )

if __name__ == '__main__':
  """The entry point for the script"""

  pool = Pool( processes=max_cores )

  np.set_printoptions(threshold=np.inf, linewidth=np.inf) # for saving np arrays without summaries
  header = 'Photo-pion production data by SOPHIA\nnature\tinput_energy\tredshift\tsecondaries_number\tsecondaries_energies\tsecondaries_type\t'
  fmt = '%d\t%.6e\t%s\t%s\t%.0f\t%.2f'

  counter = 0
  tot = float( len(backgrounds)*len(nature)*len(redshift) )/100

  for bg_name, bg in backgrounds.items():
    fname = 'sophia_%s.txt' % bg_name

    data = []

    for n in nature:
      for z in redshift:
	counter += 1
	sys.stdout.write("\r%.2f%%" % (counter/tot) )
	sys.stdout.flush()
	
	data += pool.map(swrapper, it.izip( it.repeat(n), Ein[bg_name], it.repeat(z), it.repeat(bg) ) )

    np.savetxt(fname, data, fmt='%s', header=header)

  print " Completed."

If this can help great, if not - well, you proceed as planned.

[MHoerbe]

I came up with an event generator that generates SOPHIA events from histogram data. These histograms have been created from raw SOPHIA output data. The event generator may be operated in "preserve" mode which preserves quantum numbers at the expense of run time or a "stochastic" mode which is faster but does not preserve quantum numbers.

Subsequently, you can find an example of this generator's output taken for a primary proton at 10^12 GeV interacting with a photon of 1 meV (roughly the CMB).

If you have questions or would like to see the files which have generated the data, feel free to get in touch!

[TobiasWinchen]

Looks great - I am looking forward to your pull request. Probably it would be best to have this as a pair of flags (e.g. useTabulatedData=True/False, preserveQuantumNumbers=True/False) in the existing module module.

It would be nice to also have the performance improvements of the modes compared to standard as function of the number of cores?

[MHoerbe]

I have implemented the discretized SOPHIA event generator into CRPropa's PhotoPionProduction. This module now has an additional bool useTabulatedData. The option preserveQuantumNumbers has been discarded for now as this option appears to crash CRPropa.

In the figure below you find preliminary simulation run time measurements of the two options in dependence of the number of cores used. This was made on my local machine as due to the holiday shutdown I have no access to more powerful computers. What we can learn though is a rough two-fold run time improvement of the discretized (quantum-number-preserving) event generator over regular SOPHIA. If the saturation problem at 8 cores is solved with this generator has to be investigated yet.

Simualtion details:
The simulation comprised 400 protons with energies of 1e11 GeV that were passed to the PhotoPionProduction module with a blackbody photon background of T=10 eV. The SimplePropagation was used over a total distance of 10 kpc in steps of 1e-1 pc. The photon field's density decreased according to 1/r² from the source and was constant in time.