Giter Site home page Giter Site logo

prorad's Introduction

PRORAD is a general-purpose proton radiography simulation tool developed at NIF/LLNL. It was written by summer intern Mason Black ([email protected]) and is maintained by Scott Wilks ([email protected]).

LLNL-CODE-739358

LICENSE INFO

Can be found in file https://github.com/LLNL/prorad/blob/master/LICENSE Additional information concerning the BSD License for this software can be found in file: https://github.com/LLNL/prorad/blob/master/AdditionalBSDNotice

USAGE DISCLAIMER

For purposes of speed, this code assumes that proton trajectories remain approximately paraxial. What this means for practical purposes is that the Larmor radius of your protons, at a given energy and field strength, should be significantly greater than the thickness of the region over which fields of that strength are found. If it is not, then large errors may be introduced.

COMPILATION

PRORAD is written primarily in Python 2 for ease of use, extensibility, and availability of libraries for reading the outputs of MHD and Rad-Hydro codes. However, much of the physics is handled by a compiled Python module written in C99, which improves performance drastically. To use this module you must first compile it, which is handled by a setup script. From the toplevel 'prorad' directory, run:

python setup.py build

or to compile to use OpenMP (if installed on your system), run:

python setup.py build -omp

This will generate a file named '_pmover.so'. This can be imported by Python and treated as any other module.

PRORAD can be used in two ways: interactively through a Python shell, or as a standalone script. In either case, the user needs to create an input file that defines the problem setup. The input file is itself written in Python, and must be saved in the 'inputs' directory. See existing input files for reference.

STANDALONE USE

To run as a standalone script, from the toplevel 'prorad' directory run:

python prorad.py inputfile

where inputfile is the name of an input file within the 'inputs' directory, omitting the '.py'. This will load a grid, propogate protons through it, and plot the resulting fluence at the film plane.

INTERACTIVE USE

To run PRORAD interactively, first start up an interactive Python session from the toplevel 'prorad' directory by typing 'python'.

The first step is to import prorad:

import prorad

Then, you will need to load an input file:

prorad.import_params('inputfile')

where 'inputfile' is whatever you named your file in the 'inputs' directory, omitting the '.py' (e.g. 'two_hohlraum' or 'scatter_test'). This is equivalent to running 'import inputs.inputfile as params' within the prorad script, in that all variables and functions defined in your input file will be accessible within the 'params' namespace (e.g. 'params.l_s2grid' or 'params.fformat').

Next, you will need to generate a grid through which to propogate the protons:

grid = prorad.load_grid()

Where the grid comes from depends on how you set 'params.fformat'. Current functional options are 'analytic_grid', 'FLASH', and 'HYDRA':

  • 'analytic_grid' uses the user-defined function 'params.fields' from the input file to populate a uniformly-spaced grid. 'fields' is a function of (x,y,z) and returns a 9-tuple of the form (Ex,Ey,Ez,Bx,By,Bz,nu,Z,N). E is electric field in CGS units (StatV/cm, or V/m multiplied by 0.3333e-4), B is magnetic field in CGS units (Gauss), nu is stopping power (currently unused), Z is atomic number (used for scattering only, so set to 0.0 unless in a solid), and N is ion number density (cm-3, used for scattering only).
  • 'FLASH' uses the yt library <http://yt-project.org/doc/installing.html> to read in the FLASH output file specified by 'params.fname'.
  • 'HYDRA' currently uses Yorick and Dave Strozzi's hyddjs tools to interpolate from a 2D HYDRA dump onto a 2D x-z uniform grid, and then rotates that slice around z to produce a 3D cylindrical grid. Contact David Strozzi ([email protected]) for access to the hyddjs tools.
  • 'LSP' uses the xdrlib tool to read in LSP grids from .p4 files. Contact Drew Higginson ([email protected]) for access.

Once you have a grid defined, you can run the simulation:

x,y,Ep,traces = prorad.push_protons(grid)

The output tuple contains four numpy arrays. The first three are length params.NP and contain the final x positions, y positions, and energies of the protons at the film plane. The fourth is an ndarray of shape (params.ntraces, params.nstep+3, 3) containing the x,y,z positions of a subset of the protons at each point along their trajectories, including the source and film.

Once you have run the simulation, you can analyze the final particle positions using whatever tools you like. To display a plot of proton fluence at the film plane, run:

prorad.plot_results(grid,x,y,Ep,traces)

INPUT FILES

Input files live in the 'inputs' directory. Several are included for reference. The easiest thing to do when making a new problem setup is to copy an existing input file and modify the values. Regardless of the type of grid being used, all input files must include some required parameters:

fformat : string
    An identifier string specifying the type of grid to be used.
fname : string
    If using a grid type that requires a data file, the absolute path to
    that file. If using an analytic grid, can leave as empty string.
r_source : float
    Effective radius of the proton source.
source_loc : array-like of size 3
    Location of the proton source (distances in cm).
prop_dir : array_like of size 3
    Rather than wastefully propogate protons in all directions, the code
    will initialize protons in a cone centered about this axis. Will be
    normalized automatically.
film_loc : array_like of size 3
    The center of the film on which final positions will be recorded.
    This, together with the two film axes, define a plane. Once the
    proton push finishes, proton trajectories are extrapolated to find
    the intersection with this plane, and the final 'x' and 'y' positions
    returned are locations on this plane relative to the two film axes.
film_axis1 : array_like of size 3
    Vector defining the first axis of the film. Length and direction
    determine the x axis of the fluence plot. 
film_axis2 : array_like of size 3
    Vector defining the second axis of the film. Length and direction
    determine the y axis of the fluence plot. Should be orthogonal to 
    film_axis1.
NP : int
    Number of protons.
ntraces : int
    Number of protons for which to record trajectories. Must be <= NP.
E0 : float
    Initial proton energy in MeV. Determines magnitude of initial velocity.
l_s2start : float
    Rather than being initialized at the source and propogated to the start
    of the grid, protons are instead initialized on a disc a distance
    l_s2start (in cm) along prop_dir from the source (hence the 'cone' 
    description). The distribution of positions on this disc and initial 
    velocities is set as though the protons actually did propogate here from 
    a source of the given size.
spread_angle : float
    Half-angle (in degrees) of the cone of protons. Determines the extent
    of the disc on which protons are initialized.
l_prop : float
    Distance (in cm) from starting plane (along prop_dir) through which to
    compute proton deflections. NOTE: this does not mean that protons only
    can move parallel to prop_dir. It only means that, during the push, 
    their net displacement projected onto prop_dir is of length l_prop. In
    addition, the displacement during each step projected onto prop_dir is 
    of equal length (l_prop/nsteps). This is forced to be the case by how
    the equations of motion are implemented in the pusher. This ensures
    that all protons start on one plane, and end on another parallel plane.
nsteps : float
    Number of steps each proton will take during the push. The displacement
    during each step projected onto prop_dir will be an equal distance. It
    is recommended that nsteps is greater than the expected number of grid
    cells traversed along prop_dir.

There may be other required/optional parameters depending on the grid type: * HYDRA: ngridx, ngridy, ngridz : int The number of radial cells, number of azimuthal (theta) cells, and number of z cells. HYDRA dump is interpolated onto R-z grid and then rotated about z axis. hyd_xrange : array-like of size 2, optional Min and max x values to read in from the HYDRA dump (currently the minimum is ignored because the cylindrical grid option requires the minimum radius to be 0.0). hyd_zrange : array-like of size 2, optional Min and max z values to read in from the HYDRA dump.

* analytic_grid: 
    ngridx, ngridy, ngridz : int
        The number of cells in x,y,z. 
    lx, ly, lz : float
        The spacial extend of your grid (in cm) in x,y,z.
    gridcorner : array_like of size 3
        Tuple containing the spatial position of the bottom/left/back
        corner of your grid, i.e. the position of index [0,0,0].
    fields(coord) : function, coord argument will (x,y,z) tuple.
        Returns a 9-tuple of fields at the given spatial coordinate. Will
        be called by load_grid function in prorad.py to populate a grid of
        the specified dimensions, position and extent. Try to make it fast,
        because loading large grids can take a while. The function can
        reference any other variables or functions you want to define in
        the input file. 
    grid_nthreads : int (optional)
        To speed things up, specifying a value for this will use Python's
        multiprocessing library to initialize the grid in parallel. If
        not defined the grid will be initialied in serial.
* LSP:
    The LSP reader library currently being used requires a directory name 
    and a step number, rather than an input filename. Just leave fname as
    an empty string, and define the strings 'lsp_dirname' (full path to
    directory containing LSP dump) and 'lsp_step' (int specifying step #).
    There is also the optional parameter 'lsp_ntile' that will tile the
    grid n times in x and z (used because 2D LSP grid I tested on was small 
    but periodic).

prorad's People

Contributors

mblack2012 avatar scottwilks avatar

Stargazers

 avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar

Watchers

 avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar

prorad's Issues

Suggestion: analytic_grid initialization speedup using NumPy vectorization

The section of prorad.py where it initializes the grid by looping over values in serial or even parallel is a choke point in speed:

            print("'params.grid_nthreads' not specified. Initializing grid in serial.")
            for i in field_xrange_idx:
                for j in field_yrange_idx:
                    for k in field_zrange_idx:
                        x = X_coords[i]
                        y = Y_coords[j]
                        z = Z_coords[k]
                        gridvals[i,j,k,:] = params.fields((x,y,z))   

There is a faster way to initialize the analytic grid, taking advantage of NumPy's vectorization capabilities. However, it requires more familiarity with the NumPy syntax by the user who writes the "fields()" function. I've written a local branch which leverages a new, vectorized "fields()" function in my input deck and then implements a new format in prorad.py called "analytic_grid2"; I can push it into this repo as an example if there is interest. I make a 300x300x300 analytic grid in a few seconds using the vectorized technique, but it a much smaller grid takes many times longer to generate with the original method.

Suggested addition to prorad.py:

    if fformat == "analytic_grid2":
        # Load user-defined analytic fields, Scott Feister's version using vectorized grid
        
        print("Generating "+str(ngridx)+"x"+str(ngridy)+"x"+str(ngridz)+" grid from analytic fields...")
        
        cyl_coords = False
        try: cyl_coords = params.cyl_coords
        except AttributeError: pass

        lx, ly, lz = params.lx, params.ly, params.lz
        xoffset, yoffset, zoffset = params.gridcorner

        X_coords = np.linspace(xoffset, xoffset+lx, ngridx)
        Y_coords = np.linspace(yoffset, yoffset+ly, ngridy)
        Z_coords = np.linspace(zoffset, zoffset+lz, ngridz)
        
        X, Y, Z = np.meshgrid(X_coords, Y_coords, Z_coords)

        print("Initializing grid in vectorized serial.")
        gridvals = params.fieldsvec((X, Y, Z), NUM_FIELDS)
        
        gridspacings = (lx/ngridx, ly/ngridy, lz/ngridz)
 
        grid = Grid(gridvals, gridspacings, (xoffset,yoffset,zoffset), (lx,ly,lz), cyl_coords=cyl_coords)

Corresponding example of definition of fields in input deck (a radial blob of E & B):

rad = 0.01 # Radius of blob

def fields(coord): # REQUIRED FOR ANALYTIC_GRID
    x,y,z = coord
    r = (x**2+y**2+z**2)
 
    if r < rad: # Create strongest field in the center
        Evec = (1 - r / rad) * np.array([1e4, 1e4, 1e4])
        Bvec = (1 - r / rad) * np.array([-4e4, -4e4, -4e4])
    else:
        Evec = np.array([0, 0, 0])
        Bvec = np.array([0, 0, 0])

    return (Evec[0], Evec[1], Evec[2], Bvec[0], Bvec[1], Bvec[2], 0.0,0.0,0.0)

def fieldsvec(coord, NUM_FIELDS): # REQUIRED FOR ANALYTIC_GRID2
    X, Y, Z = coord
    R = X**2 + Y**2 + Z**2
   
    Evec = np.zeros((X.shape[0],X.shape[1],X.shape[2],3))
    Bvec = np.zeros((X.shape[0],X.shape[1],X.shape[2],3))
    gridvals = np.zeros((X.shape[0],X.shape[1],X.shape[2],NUM_FIELDS))
    
    ct = R < rad
    Evec[ct,:] = np.expand_dims((1 - R[ct] / rad), axis=-1) * np.array([1e4, 1e4, 1e4])
    Bvec[ct,:] = np.expand_dims((1 - R[ct] / rad), axis=-1) * np.array([-4e4, -4e4, -4e4])
    
    gridvals[:,:,:,0:3] = Evec
    gridvals[:,:,:,3:6] = Bvec
    gridvals[:,:,:,6:9] = 0.0
    return gridvals

Headless plotting?

Headless plotting is not implemented; currently requires X11 forwarding if running on a server.

Traceback

 (most recent call last):
  File "prorad.py", line 734, in <module>
    main()
  File "prorad.py", line 78, in main
    plot_results(grid, film_x, film_y, Ep, traces, plot_fluence=True, plot_traces=True, plot_quiver=False, save_images=True)
  File "prorad.py", line 665, in plot_results
    fig = plt.figure(2,figsize=(8,8))
  File "/software/python-2.7.12-el7-x86_64/lib/python2.7/site-packages/matplotlib/pyplot.py", line 535, in figure
    **kwargs)
  File "/software/python-2.7.12-el7-x86_64/lib/python2.7/site-packages/matplotlib/backends/backend_tkagg.py", line 81, in new_figure_manager
    return new_figure_manager_given_figure(num, figure)
  File "/software/python-2.7.12-el7-x86_64/lib/python2.7/site-packages/matplotlib/backends/backend_tkagg.py", line 89, in new_figure_manager_given_figure
    window = Tk.Tk()
  File "/software/python-2.7.12-el7-x86_64/lib/python2.7/lib-tk/Tkinter.py", line 1815, in __init__
    self.tk = _tkinter.create(screenName, baseName, className, interactive, wantobjects, useTk, sync, use)
_tkinter.TclError: couldn't connect to display "localhost:16.0"

This could be corrected using the matplotlib 'Agg' plotter, which would also suppress GUI popups. I'll see about pushing a change.

Memory Error Regarding MultiProcessing When Loading Grid on Windows

Windows machines not using Cygwin run into a memory error which halts the program due to the following block of lines in prorad.py:

 try: 
            # If grid_nthreads is defined in input file, initialize the specified grid elements in parallel.
            # Using Python's multiprocessing library rather than mpi, so actual parallelism is currently limited 
            # to the number of processors on a single node.
            p = Pool(params.grid_nthreads)
            print("Using " + str(params.grid_nthreads) + " threads to initialize grid")
            coords = np.meshgrid(X_coords[field_xrange_idx],Y_coords[field_yrange_idx],Z_coords[field_zrange_idx],indexing='ij')
            coords = itertools.product(X_coords[field_xrange_idx],Y_coords[field_yrange_idx],Z_coords[field_zrange_idx])
            initialized_vals = np.array(p.map(params.fields,coords))
            initialized_vals = initialized_vals.reshape(len(field_xrange_idx),len(field_yrange_idx),len(field_zrange_idx),NUM_FIELDS)
            idx1,idx2,idx3 = np.meshgrid(field_xrange_idx,field_yrange_idx,field_zrange_idx,indexing='ij')
            gridvals[idx1,idx2,idx3] = initialized_vals
        except AttributeError:
            print("'params.grid_nthreads' not specified. Initializing grid in serial.")
            for i in field_xrange_idx:
                for j in field_yrange_idx:
                    for k in field_zrange_idx:
                        x = X_coords[i]
                        y = Y_coords[j]
                        z = Z_coords[k]
                        gridvals[i,j,k,:] = params.fields((x,y,z))

Specifically, the error exists in the line where the itertools.product (coords) are mapped to a function (params.fields), using Multiprocessing causes a memory error for Windows based machines most likely due to the fact that the regular C-compilers that are used (MinGW for example) do not have support for the map() function that is needed in the code above.

The specific line of the error is:

initialized_vals = np.array(p.map(params.fields,coords))

This can be fixed by changing the exception handling to catch the potential memory error that occurs before it crashes the program. With this fix, Linux based devices can still run MultiProcessing without error and Windows machines can run the program in serial with a small run-time difference.

The fixed code is shown below:

 try: 
            # If grid_nthreads is defined in input file, initialize the specified grid elements in parallel.
            # Using Python's multiprocessing library rather than mpi, so actual parallelism is currently limited 
            # to the number of processors on a single node.
            p = Pool(params.grid_nthreads)
            print("Using " + str(params.grid_nthreads) + " threads to initialize grid")
            coords = np.meshgrid(X_coords[field_xrange_idx],Y_coords[field_yrange_idx],Z_coords[field_zrange_idx],indexing='ij')
            coords = itertools.product(X_coords[field_xrange_idx],Y_coords[field_yrange_idx],Z_coords[field_zrange_idx])
            initialized_vals = np.array(p.map(params.fields,coords))
            initialized_vals = initialized_vals.reshape(len(field_xrange_idx),len(field_yrange_idx),len(field_zrange_idx),NUM_FIELDS)
            idx1,idx2,idx3 = np.meshgrid(field_xrange_idx,field_yrange_idx,field_zrange_idx,indexing='ij')
            gridvals[idx1,idx2,idx3] = initialized_vals

        # Deleted AttributeError catching to any generalized errors from Multiprocessing
        except:
            print("'params.grid_nthreads' not specified. Initializing grid in serial.")
            for i in field_xrange_idx:
                for j in field_yrange_idx:
                    for k in field_zrange_idx:
                        x = X_coords[i]
                        y = Y_coords[j]
                        z = Z_coords[k]
                        gridvals[i,j,k,:] = params.fields((x,y,z))

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    ๐Ÿ–– Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. ๐Ÿ“Š๐Ÿ“ˆ๐ŸŽ‰

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google โค๏ธ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.