Merge branch 'master' into cmfd_parallel

Conflicts:
	src/Makefile
This commit is contained in:
Paul Romano 2012-10-16 15:01:21 -04:00
commit a414f9295e
16 changed files with 594 additions and 14 deletions

View file

@ -13,5 +13,7 @@ as debugging.
:maxdepth: 2
structures
styleguide
workflow
xml-fortran
statepoint

View file

@ -4,3 +4,86 @@
Data Structures
===============
The purpose of this section is to give you an overview of the major data
structures in OpenMC and how they are logically related. A majority of variables
in OpenMC are `derived types`_ (similar to a struct in C). These derived types
are defined in the various header modules, e.g. src/geometry_header.F90. Most
important variables are found in the `global module`_. Have a look through that
module to get a feel for what variables you'll often come across when looking at
OpenMC code.
--------
Particle
--------
Perhaps the variable that you will see most often is simply called ``p`` and is
of type(Particle). This variable stores information about a particle's physical
characteristics (coordinates, direction, energy), what cell and material it's
currently in, how many collisions it has undergone, etc. In practice, only one
particle is followed at a time so there is no array of type(Particle). The
Particle type is defined in the `particle_header module`_.
You will notice that the direction and angle of the particle is stored in a
linked list of type(LocalCoord). In geometries with multiple :ref:`universes`,
the coordinates in each universe are stored in this linked list. If universes or
lattices are not used in a geometry, only one LocalCoord is present in the
linked list.
The LocalCoord type has a component called cell which gives the index in the
``cells`` array in the `global module`_. The ``cells`` array is of type(Cell)
and stored information about each region defined by the user.
----
Cell
----
The Cell type is defined in the `geometry_header module`_ along with other
geometry-related derived types. Each cell in the problem is described in terms
of its bounding surfaces, which are listed on the ``surfaces`` component. The
absolute value of each item in the ``surfaces`` component contains the index of
the corresponding surface in the ``surfaces`` array defined in the `global
module`_. The sign on each item in the ``surfaces`` component indicates whether
the cell exists on the positive or negative side of the surface (see
:ref:`methods_geometry`).
Each cell can either be filled with another universe/lattice or with a
material. If it is filled with a material, the ``material`` component gives the
index of the material in the ``materials`` array defined in the `global
module`_.
-------
Surface
-------
The Surface type is defined in the `geometry_header module`_. A surface is
defined by a type (sphere, cylinder, etc.) and a list of coefficients for that
surface type. The simplest example would be a plane perpendicular to the xy, yz,
or xz plane which needs only one parameter. The ``type`` component indicates the
type through integer parameters such as SURF_SPHERE or SURF_CYL_Y (these are
defined in the `constants module`_). The ``coeffs`` component gives the
necessary coefficients to parameterize the surface type (see
:ref:`surface_element`).
--------
Material
--------
The Material type is defined in the `material_header module`_. Each material
contains a number of nuclides at a given atom density. Each item in the
``nuclide`` component corresponds to the index in the global ``nuclides`` array
(as usual, found in the `global module`_). The ``atom_density`` component is the
same length as the ``nuclides`` component and lists the corresponding atom
density in atom/barn-cm for each nuclide in the ``nuclides`` component.
If the material contains nuclides for which binding effects are important in
low-energy scattering, a :math:`S(\alpha,\beta)` can be associated with that
material through the ``sab_table`` component. Again, this component contains the
index in the ``sab_tables`` array from the `global module`_.
.. _derived types: http://nf.nci.org.au/training/FortranAdvanced/slides/slides.025.html
.. _global module: https://github.com/mit-crpg/openmc/blob/master/src/global.F90
.. _particle_header module: https://github.com/mit-crpg/openmc/blob/master/src/particle_header.F90
.. _geometry_header module: https://github.com/mit-crpg/openmc/blob/master/src/geometry_header.F90
.. _constants module: https://github.com/mit-crpg/openmc/blob/master/src/constants.F90
.. _material_header module: https://github.com/mit-crpg/openmc/blob/master/src/material_header.F90

View file

@ -0,0 +1,142 @@
.. _devguide_styleguide:
======================
Style Guide for OpenMC
======================
In order to keep the OpenMC code base consistent in style, this guide specifies
a number of rules which should be adhered to when modified existing code or
adding new code in OpenMC.
-------------
General Rules
-------------
Conform to the Fortran 2008 standard.
Make sure code can be compiled with most common compilers, especially gfortran
and the Intel Fortran compiler. This supercedes the previous rule --- if a
Fortran 2003/2008 feature is not implemented in a common compiler, do not use
it.
Do not use special extensions that can be only be used from certain compilers.
In general, write your code in lower-case. Having code in all caps does not
enhance code readability or otherwise.
Always include comments to describe what your code is doing. Do not be afraid of
using copious amounts of comments.
Use <, >, <=, >=, ==, and /= rather than .lt., .gt., .le., .ge., .eq., and .ne.
Try to keep code within 80 columns when possible.
Don't use ``print *`` or ``write(*,*)``. If writing to a file, use a specific
unit. Writing to standard output or standard error should be handled by the
``write_message`` subroutine or functionality in the error module.
-------------------------
Subroutines and Functions
-------------------------
Above each subroutine/function, include a comment block giving a brief
description of what the subroutine or function does.
Arguments to subroutines/functions should be explicitly specified as intent(in),
intent(out), or intent(inout).
Include a comment describing what each argument to a subroutine/function is.
---------
Variables
---------
Never, under any circumstances, should implicit variables be used! Always
include ``implicit none`` and define all your variables.
Variable names should be all lower-case and descriptive, i.e. not a random
assortment of letters that doesn't give any information to someone seeing it for
the first time. Variables consisting of multiple words should be separated by
underscores, not hyphens or in camel case.
Constant (parameter) variables should be in ALL CAPITAL LETTERS and defined in
in the constants.F90 module.
32-bit reals (real(4)) should never be used. Always use 64-bit reals (real(8)).
For arbitrary length character variables, use the pre-defined lengths
MAX_LINE_LEN, MAX_WORD_LEN, and MAX_FILE_LEN if possible.
Integer values being used to indicate a certain state should be defined as named
constants (see the constants.F90 module for many examples).
Yes:
.. code-block:: fortran
if (boundary_condition == BC_VACUUM) then
No:
.. code-block:: fortran
if (boundary_condition == -10) then
Never create arrays with a pre-defined maximum length. Always use dynamic memory
allocation. Use allocatable variables instead of pointer variables when
possible.
Shared/Module Variables
-----------------------
Always put shared variables in modules. Access module variables through a
``use`` statement. Always use the ``only`` specifier on the ``use`` statement
except for variables from the global, constants, and various header modules.
Never use ``equivalence`` statements, ``common`` blocks, or ``data`` statements.
-------------------------
Derived Types and Classes
-------------------------
Derived types and classes should have CamelCase names with words not separated
by underscores or hyphens.
-----------
Indentation
-----------
Never use tab characters. Indentation should always be applied using
spaces. Emacs users should include the following line in their .emacs file:
.. code-block:: common-lisp
(setq-default indent-tabs-mode nil)
vim users should include the following line in their .vimrc file::
set expandtab
Use 2 spaces per indentation level. This applies to all constructs such as
program, subroutine, function, if, associate, etc. Emacs users should set the
variables f90-if-indent, f90-do-indent, f90-continuation-indent,
f90-type-indent, f90-associate-indent, and f90-program indent to 2.
Continuation lines should be indented by an extra 5 spaces. This is the default
value of f90-continuation-indent in Emacs.
-------------------------
Whitespace in Expressions
-------------------------
Avoid extraneous whitespace in the following situations:
- In subroutine/function calls::
Yes: call somesub(x, y(2), z)
No: call somesub( x, y( 2 ), z )
- In logical expressions, use one space around operators but nowhere else::
Yes: if (variable == 2) then
No: if ( variable==2 ) then

View file

@ -0,0 +1,45 @@
.. _devguide_workflow:
====================
Development Workflow
====================
Anyone wishing to make contributions to OpenMC should be fully acquianted and
comfortable working with git_ and GitHub_. The primary means of modifying and
making contributions to OpenMC is through GitHub `pull requests`_. This is
what's known as a fork and pull development model. The steps for this are as
follows:
1. Fork the main openmc repository from `mit-crpg/openmc`_. This will create a
repository with the same name under your personal account. As such, you can
commit to it as you please without disrupting other developers.
2. Create a branch that you want merged back to `mit-crpg/openmc`_ and make
commits that you intend to go back. If you have made other changes that should
not be merged back, those changes should be on another branch.
3. Issue a pull request from GitHub and select the branch you want merged.
4. The OpenMC integration manager will review your pull request and make sure it
conforms to the :ref:`devguide_styleguide`, compiles correctly, runs in parallel
correctly, does not break other features in the code, etc. Any issues with the
pull request can be discussed directly on the pull request page itself.
5. After the pull request has been thoroughly vetted, it is merged back into
`mit-crpg/openmc`_.
While the process above depends on the fork of the OpenMC repository being
publicly available on GitHub, you may also wish to do development on a private
repository for research or commercial purposes. The proper way to do this is to
create a complete copy of the OpenMC repository (not a fork from GitHub). The
private repository can then either be stored just locally or in conjunction with
a private repository on Github (this requires a `paid plan`_). If you want to
merge some changes you've made in your private repository back to
`mit-crpg/openmc`_ repository, simply follow the steps above with an extra step
of pulling a branch from your private repository into your public fork.
.. _git: http://git-scm.com/
.. _GitHub: https://github.com/
.. _pull requests: https://help.github.com/articles/using-pull-requests
.. _mit-crpg/openmc: https://github.com/mit-crpg/openmc
.. _paid plan: https://github.com/plans

View file

@ -4,3 +4,37 @@
xml-fortran Input Parsing
=========================
OpenMC relies on the xml-fortran package for reading and intrepreting the XML
input files for geometry, materials, settings, tallies, etc. The use of an XML
format makes writing input files considerably more flexible than would otherwise
be possible.
With the xml-fortran package, extending the user input files to include new tags
is fairly straightforward. A "template" file exists for each diferent type of
input file that tells xml-fortran what to expect in a file. These template files
can be found in the src/templates directory. The steps for modifying/adding
input are as follows:
1. Add a ``<variable>``` tag to the desired template file,
e.g. src/templates/geometry_t.xml. See the `xml-fortran documentation`_ for a
description of the acceptable fields.
2. In the input_xml module, any input given in your new tag will be read
automatically through a call to, e.g. read_xml_file_geometry_t. Whatever
variable name you specified should have the data available.
3. Add code in the appropriate subroutine to check the variable for any possible
errors.
4. Add a variable in OpenMC to copy the temporary variable into if there are no
errors.
A set of `RELAX NG`_ schemata exists that enables real-time validation of input
files when using the GNU Emacs text editor. You should also modify the RELAX NG
schema for the template you changed (e.g. src/templates/geometry.rnc) so that
those who use Emacs can confirm whether their input is valid before they
run. You will need to be familiar with RELAX NG `compact syntax`_.
.. _xml-fortran documentation: http://xml-fortran.sourceforge.net/documentation.html
.. _RELAX NG: http://relaxng.org/
.. _compact syntax: http://relaxng.org/compact-tutorial-20030326.html

View file

@ -64,6 +64,8 @@ In OpenMC, any second-order surface of the form
can be modeled in OpenMC. For example, the equation for a sphere centered at
:math:`(\bar{x},\bar{y},\bar{z})` and of radius :math:`R` can be written as
.. _universes:
Universes
---------

View file

@ -5,18 +5,18 @@ Publications
============
- Paul K. Romano and Benoit Forget, "Reducing Parallel Communication in Monte
Carlo Simulations via Batch Statistics," *Trans. Am. Nucl. Soc.*, Accepted
(2012).
Carlo Simulations via Batch Statistics," *Trans. Am. Nucl. Soc.*, **107**,
xxx--xxx (2012).
- Paul K. Romano and Benoit Forget, "The OpenMC Monte Carlo Particle Transport
Code," *Annals of Nuclear Energy*, Accepted (2012).
Code," *Ann. Nucl. Energy*, **51**, 274--281
(2012). `<http://dx.doi.org/10.1016/j.anucene.2012.06.040>`_
- Andrew R. Siegel, Kord Smith, Paul K. Romano, Benoit Forget, and Kyle Felker,
"The effect of load imbalances on the performance of Monte Carlo codes in LWR
analysis", *Journal of Computational Physics*,
analysis", *J. Comput. Phys.*,
`<http://dx.doi.org/10.1016/j.jcp.2012.06.012>`_ (2012).
- Paul K. Romano and Benoit Forget, "Parallel Fission Bank Algorithms in Monte
Carlo Criticality Calculations," *Nuclear Science and Engineering*, **170**,
pp. 125--135 (2012). [`PDF
<http://web.mit.edu/romano7/www/nse_v170_n2_pp125-135.pdf>`_]
Carlo Criticality Calculations," *Nucl. Sci. Eng.*, **170**, 125--135
(2012). `<http://hdl.handle.net/1721.1/73569>`_

View file

@ -454,6 +454,8 @@ could be written as:
</geometry>
.. _surface_element:
``<surface>`` Element
---------------------
@ -519,6 +521,21 @@ The following quadratic surfaces can be modeled:
A sphere of the form :math:`(x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2 =
R^2`. The coefficients specified are ":math:`x_0 \: y_0 \: z_0 \: R`".
:x-cone:
A cone parallel to the x-axis of the form :math:`(y - y_0)^2 + (z - z_0)^2 =
R^2 (x - x_0)^2`. The coefficients specified are ":math:`x_0 \: y_0 \: z_0
\: R^2`".
:y-cone:
A cone parallel to the y-axis of the form :math:`(x - x_0)^2 + (z - z_0)^2 =
R^2 (y - y_0)^2`. The coefficients specified are ":math:`x_0 \: y_0 \: z_0
\: R^2`".
:z-cone:
A cone parallel to the x-axis of the form :math:`(x - x_0)^2 + (y - y_0)^2 =
R^2 (z - z_0)^2`. The coefficients specified are ":math:`x_0 \: y_0 \: z_0
\: R^2`".
``<cell>`` Element
------------------

View file

@ -107,11 +107,14 @@ module constants
SURF_CYL_Y = 6, & ! Cylinder along y-axis
SURF_CYL_Z = 7, & ! Cylinder along z-axis
SURF_SPHERE = 8, & ! Sphere
SURF_BOX_X = 9, & ! Box extending infinitely in x-direction
SURF_BOX_Y = 10, & ! Box extending infinitely in y-direction
SURF_BOX_Z = 11, & ! Box extending infinitely in z-direction
SURF_BOX = 12, & ! Rectangular prism
SURF_GQ = 13 ! General quadratic surface
SURF_CONE_X = 9, & ! Cone parallel to x-axis
SURF_CONE_Y = 10, & ! Cone parallel to y-axis
SURF_CONE_Z = 11, & ! Cone parallel to z-axis
SURF_BOX_X = 12, & ! Box extending infinitely in x-direction
SURF_BOX_Y = 13, & ! Box extending infinitely in y-direction
SURF_BOX_Z = 14, & ! Box extending infinitely in z-direction
SURF_BOX = 15, & ! Rectangular prism
SURF_GQ = 16 ! General quadratic surface
! ============================================================================
! CROSS SECTION RELATED CONSTANTS

View file

@ -376,6 +376,48 @@ contains
v = v - 2*dot_prod*y/(R*R)
w = w - 2*dot_prod*z/(R*R)
case (SURF_CONE_X)
! Find x-x0, y-y0, z-z0 and dot product of direction and surface
! normal
x = p % coord0 % xyz(1) - surf % coeffs(1)
y = p % coord0 % xyz(2) - surf % coeffs(2)
z = p % coord0 % xyz(3) - surf % coeffs(3)
R = surf % coeffs(4)
dot_prod = (v*y + w*z - R*u*x)/((R + ONE)*R*x*x)
! Reflect direction according to normal
u = u + 2*dot_prod*R*x
v = v - 2*dot_prod*y
w = w - 2*dot_prod*z
case (SURF_CONE_Y)
! Find x-x0, y-y0, z-z0 and dot product of direction and surface
! normal
x = p % coord0 % xyz(1) - surf % coeffs(1)
y = p % coord0 % xyz(2) - surf % coeffs(2)
z = p % coord0 % xyz(3) - surf % coeffs(3)
R = surf % coeffs(4)
dot_prod = (u*x + w*z - R*v*y)/((R + ONE)*R*y*y)
! Reflect direction according to normal
u = u - 2*dot_prod*x
v = v + 2*dot_prod*R*y
w = w - 2*dot_prod*z
case (SURF_CONE_Z)
! Find x-x0, y-y0, z-z0 and dot product of direction and surface
! normal
x = p % coord0 % xyz(1) - surf % coeffs(1)
y = p % coord0 % xyz(2) - surf % coeffs(2)
z = p % coord0 % xyz(3) - surf % coeffs(3)
R = surf % coeffs(4)
dot_prod = (u*x + v*y - R*w*z)/((R + ONE)*R*z*z)
! Reflect direction according to normal
u = u - 2*dot_prod*x
v = v - 2*dot_prod*y
w = w + 2*dot_prod*R*z
case default
message = "Reflection not supported for surface " // &
trim(to_str(surf % id))
@ -853,6 +895,147 @@ contains
end if
case (SURF_CONE_X)
x0 = surf % coeffs(1)
y0 = surf % coeffs(2)
z0 = surf % coeffs(3)
r = surf % coeffs(4)
x = x - x0
y = y - y0
z = z - z0
a = v*v + w*w - r*u*u
k = y*v + z*w + r*x*u
c = y*y + z*z - r*x*x
quad = k*k - c
if (quad < ZERO .or. a == ZERO) then
! no intersection with cone
d = INFINITY
elseif (on_surface) then
! particle is on the cone, thus one distance is
! positive/negative and the other is zero. The sign of k
! determines if we are facing in or out
if (k >= ZERO) then
d = INFINITY
else
d = (-k + sqrt(quad))/a
end if
elseif (c < ZERO) then
! particle is inside the cone, thus one distance must be
! negative and one must be positive. The positive distance will
! be the one with negative sign on sqrt(quad)
d = (-k + sqrt(quad))/a
else
! particle is outside the cone, thus both distances are either
! positive or negative. If positive, the smaller distance is the
! one with positive sign on sqrt(quad)
d = (-k - sqrt(quad))/a
if (d <= ZERO) d = INFINITY
end if
case (SURF_CONE_Y)
x0 = surf % coeffs(1)
y0 = surf % coeffs(2)
z0 = surf % coeffs(3)
r = surf % coeffs(4)
x = x - x0
y = y - y0
z = z - z0
a = u*u + w*w - r*v*v
k = x*u + z*w + r*y*v
c = x*x + z*z - r*y*y
quad = k*k - c
if (quad < ZERO .or. a == ZERO) then
! no intersection with cone
d = INFINITY
elseif (on_surface) then
! particle is on the cone, thus one distance is
! positive/negative and the other is zero. The sign of k
! determines if we are facing in or out
if (k >= ZERO) then
d = INFINITY
else
d = (-k + sqrt(quad))/a
end if
elseif (c < ZERO) then
! particle is inside the cone, thus one distance must be
! negative and one must be positive. The positive distance will
! be the one with negative sign on sqrt(quad)
d = (-k + sqrt(quad))/a
else
! particle is outside the cone, thus both distances are either
! positive or negative. If positive, the smaller distance is the
! one with positive sign on sqrt(quad)
d = (-k - sqrt(quad))/a
if (d <= ZERO) d = INFINITY
end if
case (SURF_CONE_Z)
x0 = surf % coeffs(1)
y0 = surf % coeffs(2)
z0 = surf % coeffs(3)
r = surf % coeffs(4)
x = x - x0
y = y - y0
z = z - z0
a = u*u + v*v - r*w*w
k = x*u + y*v + r*z*w
c = x*x + y*y - r*z*z
quad = k*k - c
if (quad < ZERO .or. a == ZERO) then
! no intersection with cone
d = INFINITY
elseif (on_surface) then
! particle is on the cone, thus one distance is
! positive/negative and the other is zero. The sign of k
! determines if we are facing in or out
if (k >= ZERO) then
d = INFINITY
else
d = (-k + sqrt(quad))/a
end if
elseif (c < ZERO) then
! particle is inside the cone, thus one distance must be
! negative and one must be positive. The positive distance will
! be the one with negative sign on sqrt(quad)
d = (-k + sqrt(quad))/a
else
! particle is outside the cone, thus both distances are either
! positive or negative. If positive, the smaller distance is the
! one with positive sign on sqrt(quad)
d = (-k - sqrt(quad))/a
if (d <= ZERO) d = INFINITY
end if
case (SURF_GQ)
message = "Surface distance not yet implement for general quadratic."
call fatal_error()
@ -1035,6 +1218,36 @@ contains
z = z - z0
func = x*x + y*y + z*z - r*r
case (SURF_CONE_X)
x0 = surf % coeffs(1)
y0 = surf % coeffs(2)
z0 = surf % coeffs(3)
r = surf % coeffs(4)
x = x - x0
y = y - y0
z = z - z0
func = y*y + z*z - r*x*x
case (SURF_CONE_Y)
x0 = surf % coeffs(1)
y0 = surf % coeffs(2)
z0 = surf % coeffs(3)
r = surf % coeffs(4)
x = x - x0
y = y - y0
z = z - z0
func = x*x + z*z - r*y*y
case (SURF_CONE_Z)
x0 = surf % coeffs(1)
y0 = surf % coeffs(2)
z0 = surf % coeffs(3)
r = surf % coeffs(4)
x = x - x0
y = y - y0
z = z - z0
func = x*x + y*y - r*z*z
case (SURF_BOX_X)
y0 = surf % coeffs(1)
z0 = surf % coeffs(2)

View file

@ -247,6 +247,12 @@ contains
call h5ltmake_dataset_string_f(temp_group, "type", "Z Cylinder", hdf5_err)
case (SURF_SPHERE)
call h5ltmake_dataset_string_f(temp_group, "type", "Sphere", hdf5_err)
case (SURF_CONE_X)
call h5ltmake_dataset_string_f(temp_group, "type", "X Cone", hdf5_err)
case (SURF_CONE_Y)
call h5ltmake_dataset_string_f(temp_group, "type", "Y Cone", hdf5_err)
case (SURF_CONE_Z)
call h5ltmake_dataset_string_f(temp_group, "type", "Z Cone", hdf5_err)
case (SURF_BOX_X)
case (SURF_BOX_Y)
case (SURF_BOX_Z)

View file

@ -766,6 +766,15 @@ contains
case ('sphere')
s % type = SURF_SPHERE
coeffs_reqd = 4
case ('x-cone')
s % type = SURF_CONE_X
coeffs_reqd = 4
case ('y-cone')
s % type = SURF_CONE_Y
coeffs_reqd = 4
case ('z-cone')
s % type = SURF_CONE_Z
coeffs_reqd = 4
case ('box-x')
s % type = SURF_BOX_X
coeffs_reqd = 4

View file

@ -3,7 +3,7 @@ module material_header
implicit none
!===============================================================================
! MATERIAL describes a material by its constituent isotopes
! MATERIAL describes a material by its constituent nuclides
!===============================================================================
type Material

View file

@ -486,6 +486,12 @@ contains
string = "Z Cylinder"
case (SURF_SPHERE)
string = "Sphere"
case (SURF_CONE_X)
string = "X Cone"
case (SURF_CONE_Y)
string = "Y Cone"
case (SURF_CONE_Z)
string = "Z Cone"
case (SURF_BOX_X)
case (SURF_BOX_Y)
case (SURF_BOX_Z)

View file

@ -8,7 +8,7 @@ module source_header
!===============================================================================
type ExtSource
integer :: type_space ! spacial distributione, e.g. 'box' or 'point'
integer :: type_space ! spacial distribution, e.g. 'box' or 'point'
integer :: type_angle ! angle distribution, e.g. 'isotropic'
integer :: type_energy ! energy distribution, e.g. 'Watt'
real(8), allocatable :: params_space(:) ! parameters for spatial distribution

View file

@ -1,5 +1,23 @@
module state_point
!===============================================================================
! STATE_POINT -- This module handles writing and reading binary state point
! files. State points are contain complete tally results, source sites, and
! various other data. They can be used to restart a run or to reconstruct
! confidence intervals for tallies (this requires post-processing via Python
! scripts).
!
! Modifications to this module should be made with care. There are essentially
! three different ways to write or read state points: 1) normal Fortran file
! I/O, 2) MPI-IO, and 3) HDF5. The HDF5 functionality is contained in the
! hdf5_interface module. If you plan to change the state point, you will need to
! change all methods. You should also increment REVISION_STATEPOINT in the
! constants module.
!
! State points can be written at any batch during a simulation, or at specified
! intervals, using the <state_point ... /> tag.
!===============================================================================
use error, only: warning, fatal_error
use global
use math, only: t_percentile