README.txt
May 14, 2018 ยท View on GitHub
BEAM readme updated May 8, 2018
BEAM is a Monte Carlo (MC) based Radiative Transfer/Ray Tracing code written in Fortran 90. In addition, this package includes a few utility scripts written in python and R that we have found useful for diagnostic purposes. The MC code has been run and tested on MacOS and Linux systems (Pleiades at the NAS supercomputing facility).
The software consists of two parts; a set up code that generates a particle field and the "main" code that does the Monte Carlo based raytracing. The two parts are decoupled to allow for using a particle field obtained some other way.
prerequisites:
gnu fortran,
If running BEAM on a Mac:
If you don't have gcc installed. I recommend getting the binary from
You will also need to install Xtools from the Mac App store.
Open MPI,
if you want to run the parallel version locally (on your own machine) you will most likely need to install open MPI. Open MPI can be downloaded from,
https://www.open-mpi.org/software
When you download (get the most recent version), you will find a file called INSTALL that will guide you through getting it set up.
python,
Most systems (at least MacOS and Linux) have python pre-installed. NAS has a variety of python Versions available. The two scripts included here have been tested with python 2.7 but will probably work with newer versions, such as python 3. Some basic information about python is given in the document python/pythoninfo.txt included with this package.
R language,
In this package we include an R script as a utility for visualizing particle fields generated by BEAM. We recommend installing R studio. It has a nice matlab-like interface and can be downloaded for free from,
https://www.rstudio.com/products/rstudio/
for the program plotGranola.R you will need to add the RGL library this can be done easily from the tools menu in R studio.
Modules on NAS,
If running BEAM on Pleiades you will you will have to load the modules you need to compile and run the code. As of this writing the modules to load are,
>module load gcc/4.9.3
>module load mpi-sgi/mpt.2.15r20
>module load python
These seem to work but newer versions "should" work fine too.
Contents of package:
README.txt
-this file
in_progress.txt
- Document that details features of the software still being tested
run_mpi4_S0tau05D001Phase60sha45Ntarg300.pbs
-this a sample pbs script for runing BEAM on NAS it
executes the MPI version of the fortran and python
code mentioned below to average over multiple runs.
runTau10D01S0_Phase20_50targs.sh
-Sample BASH shell script for averaging over multiple runs
that can be run on any unix system. This is meant to duplicate
what the pbs script does, but averages over far fewer runs
and requires a much longer runtime.
source
Fortran source code
make_particles_wakes.F90
-original particle set up code. Contains only
a moderate capability to ensure particles do not overlap.
Works well with low D and tau (density less
then ~ 0.1, tau less then about 1.0) and narrow size
distributions (r=5 to 500 cm, for example).
mpw_iquad_xy.f90
-particle set up code that does a better job
of ensuring that no particles overlap by seperating
the coordinate grid into quadrants. Better to use this
code for high density cases with a wide size distribution.
mc.F90
-main ray tracing code
mc_mpi.F90
-main code madified to use message massing interface (MPI) on NAS
planarflux_simpsons.f90
-program used to generate the Apf tables found in support files
directory.
misc
plotGranola.R
-R program that plots a particle field in 3d. Used for visualizing
output from make_particles_wakes.F90 or mpw_iquad_xy.f90. For
different particle files you just need to change line 4.
supportfiles
-Tables of planar flux for a k=1,1.05, and 1.15 for S=0 to 2.0
used by montecarlo code. The files are,
apf_table_k115.dat
apf_table_k105.dat
apf_table_k1.dat
these were generated by planarflux_simpsons.f90 assuming a
hapke-like shadowing function (f(a) = exp(-shad*sqrt(tan(a/2)))).
If you want to use a simple exponential (f(a)=exp(-shad*a))
in the main code, you will need to generate new tables.
python
python scripts,
measureTau.py
-Can be used to check the particle file that it accurately
represents the desired tau and D values. Because of the
option in model to have particles gaussian distributed
vertically, the D values output from this script may
be less accurate the the tau values.
SumFlux.py
-Used by pbs and shell script to average flux and I/F over
multiple simulations.
sphericalAlbedo.py
-script that plots the spherical albedo based/inspired by
the paper,
"Rough surfaces: Is the dark stuff just shadow?"(2016)
by Jeffrey N. Cuzzi, Lindsey B. Chambers, Amanda R. Hendrix
Used to understand how shadowing and mineart parameters
effect reflectivity and not required for running BEAM.
pythoninfo.txt
-a document that describes where to obtain a python interpreter
and installing modules that are used in our utility programs.
Also, includes where to get more information on python online.
bin
-compiled binaries of source code go in this directory
mcx - monte carlo program
wakex - particle field set up program
doc
-supporting documentation
input
-sample input files for BEAM
-input files for monte carlo code
Lambert surface element scattering
mc-beamS0Phase20sha100tau10D01LitGeom.in
mc-beamS0Phase100sha100tau10D01LitGeom.in
mc-beamS0Phase150sha100tau10D01LitGeom.in
Callisto phase function scattering
mc-beamS0Phase20sha100tau10D01Callisto.in
-input files for setup code
part-no_wakesD01tau10_nbins40.in
part-wakesDw01taug01tauw10_phi0_nbins40.in
part-wakesDw01taug01tauw10_phi30_nbins40.in
-template files used by python script SumFlux.py
if_sum.out
flux_sum.out
output
this directory contains sample output.Details follow
Ouput from Monte Carlo code,
All of these use obslat = 41.7 degrees, sunlat = 26.7,
distance from saturn = 88000 km
Lambert_D01tau10Phase20_mcx - output files from run with
lambert surface element scattering with phase angle of
20 degrees
Lambert_D01tau10Phase100_mcx - output files from run with
lambert surface element scattering with phase angle of
100 degrees
Callisto_D01tau10Phase20_mcx - output files from run with
Callisto phase function scattering with phase angle of
20 degrees
Lambert_D01tau10Phase20_mc_mpi - output files from run with
lambert surface element scattering with phase angle of
20 degrees using MPI version of the code.
Ntarg300D01tau10_S0_phase20sha100 - results generated by
the pbs script run_mpi4_S0tau10D01Phase20sha45Ntarg300.pbs
on Pleides.
BEAM_tau10D01S0Phase20SHA100Rev3_50targs - results generated
the shell script runTau10D01S0_Phase20_50targs.sh
Output from mpw_iquad_xy.f90 or make_particles_wakes.F90
nowakes_D01_tau10_nbins40.partdat - particle field with
with D=0.1, tau=1.0 (nowakes case), rmin = 5 cm, and
rmax = 500 cm.
wakes_Dw01taug01tauw10_phi0_nbins40.partdat - particle field
with wake/gap regions. In the gap tau is 0.1, in the wake
D=0.1, tau = 1.0, rmin = 5 cm, and
rmax = 500 cm. The wake orientation is at 0.0 degrees.
wakes_Dw01taug01tauw10_phi30_nbins40.partdat - particle field
with wake/gap regions. In the gap tau is 0.1, in the wake
D=0.1, tau = 1.0, rmin = 5 cm, and rmax = 500 cm.
The wake orientation is at 30.0 degrees.
Compiling code:
cd into source directory
particle code should be compiled as,
>gfortran -O3 make_particles_wakes.F90 -o ../bin/wakex
or the "quad" set up code as,
>gfortran -O3 -mcmodel=medium mpw_iquad_xy.f90 -o ../bin/wakex
*theoretically, -mcmodel=large can be used, but this is not fully
implemented in gfortran and will not work in a MacOS environment.
monte carlo code (serial version) should be compiled as,
> gfortran -O3 mc.F90 -o ../bin/mcx
the MPI (parallel) version can be compiled on MacOS as,
> mpifort -O3 mc_mpi.F90 -o ../bin/mc_mpix
On NAS you will need to load a couple of modules first. As of
this writing they are,
>module load gcc/4.9.3
>module load mpi-sgi/mpt.2.15r20
then compile MPI code as,
>gfortran -O3 -I/nasa/sgi/mpt/2.15r20/include -o ../bin/mc_mpix mc_mpi.F90 -lmpi
if you use a different version of the sgi/mpt library you will
have to change the compile statement accordingly.
Running the model
First run the setup code, the steps below assume you are have cd'd
into the output directory,
type,
>../bin/wakex < ../input/part-no_wakesD01tau10_nbins40.in > output_particles
the file output_particles will contain messages normally output to the screen.
I use the included file ../input/part-no_wakesD01tau10_nbins40.in as an example.
This generates a non-wake particle field with the smallest particle at r=5cm and largest
at r=500 cm. The particle sizes start with the largest and are put in
geometrically spaced bins (this example has 40 bins). The physical parameters are
tau - 1.0 and D = 0.1. After running you should see a file called:
nowakes_D01_tau10_nbins40.partdat. This will be used by mcx
**** Actual tau and D ****************************************************
A bit of caution here the actual particle field generated will generally
have a tau then the input tau. How much depends on the number of size bins
used. The same is true of D. I find with 40 bins tau is off by 2% and D is
off by 7-8%. If you want to validate the results of BEAM with classical
results or compare with other models you will want to measure the tau and D
in the particle file. You may find the script measureTau.py useful for this.
However, be aware that because, in rare cases, there are particles whose z
coordinates are outliers the calculated D values are misleading.
**************************************************************************
Running mcx (the serial version),
type > ../bin/mcx < ../input/mc-beamS0Phase20sha100tau10D01LitGeom.in > output
The file output contains messages printed while the program runs.
In this example the sun latitude is 26.7 degrees observer latitude is
41.7 degrees, phase angle is 20 degrees, and the solar hour angle
is 100 degrees. The distance to the ring patch is 88000 km with
10000 solar and satshine photons are used. We use the the minneart
scatering law as well. There is also a particle roughness parameter
but this is set to zero here.
One side note, for the viewing geometry to make sense, the phase angle
can not be larger then 180 - (obslat+sunlat) and can not be smaller
then abs(obslat - sunlat). If you try using the input file
mc-beamS0Phase150sha100tau10D01LitGeom.in you will get an error message
and the program will halt execution for this reason.
Four files are generated,
flux.out - a table of flux due to orders of scattering from solar photons
and satshine
if.out - tabulation of I/F as a function of single scattering albedo
nphotpass.out - a file that is useful for understanding output in
some cases. It reports 6 numbers;
* Total number of photons that pass through without hitting any particles
* Satshine photons that that did not interact
* Solar photons that did not interact
* Satshine photons used to compute flux
* Solar photons used to compute flux
* Number of scattering events that result in flux not seen by the observer
due to the scattering particle being in the way.
* Number of scattering events that result in flux not seen by the observer
because of another particle being in the way.
output - contains messages related to progress of program run
Running mc_mpix (the parallel version),
Most higher end Macs have 2 to 8 cores so the parallel version can
be run locally by first linking the input file to a hardcoded filename,
>ln -s ../input/mc-beamS0Phase20sha100tau10D01LitGeom.in input_file.in
and then (if you have four cores) execute the code as,
>mpiexec -np 4 ../bin/mc_mpix
generates same files as serial version
running the model on Pleidies
As aluded to above, to speed up execution time we converted the serial
version to a MPI version that splits up a large number of
photons among multiple processors. Also we found averaging runs
among unique particle fields (or targets) with the same physical
parameters eliminates bias. So we developed a pbs script that
loops through 100 (for high density cases) or 300 (for low density)
targets. The PBS script assumes we are running it in a CSHELL environment.
The beginning of the code has some shell variable declarations used
form input filenames and output directories.
So after compiling both codes on NAS one can execute the sample pbs
script in the root directory as,
>qsub -q normal run_mpi4_S0tau10D01Phase20sha45Ntarg300.pbs
this puts the "job" into the normal queue. This is for jobs that take
less then 8 hours to run.
as it stands output files are written to /nobackupp8/dmolson/
on lines 37 and 38 you will need to modify those lines to work in your
account.
When executed the script will put results in a folder such as the one
you find in /output/Ntarg300D01tau10_S0_phase20sha100/
with seperate results for each target as well as averaged results in flux_sum.out
and if_sum.out.
Shell script
We also provide a shell script that does the same thing as the PBS script
but uses the serial version of the code rather then the mpi version.
While in the BASH shell (the script is written using BASH commands)
and in the root directory type,
> ./runTau10D01S0_Phase20_50targs.sh
While running, the script echos the progress as it goes through
each run. As with the PBS script, this script first creates a
parent folder for the results with a subfolder for the results
from individual runs. An example of output can be found in
/output/BEAM_tau10D01S0Phase20SHA100Rev3_50targs.