Merge remote-tracking branch 'upstream/develop' into squash-distribcells

This commit is contained in:
Will Boyd 2015-05-11 11:20:45 -04:00
commit 24d7502af8
37 changed files with 152 additions and 2935 deletions

View file

@ -10,7 +10,6 @@ before_install:
- export MPI_DIR=$PWD/mpich_install
- export PHDF5_DIR=$PWD/phdf5_install
- export HDF5_DIR=$PWD/hdf5_install
- export PETSC_DIR=$PWD/petsc_install
before_script:
- cd data

View file

@ -0,0 +1,64 @@
.. _devguide_docbuild:
=============================
Building Sphinx Documentation
=============================
In order to build the documentation in the ``docs`` directory, you will need to
have the Sphinx_ third-party Python package. The easiest way to install Sphinx
is via pip:
.. code-block:: sh
sudo pip install sphinx
Additionally, you will also need two Sphinx extensions for TikZ support and
numbering figures. The sphinxcontrib-tikz_ package should be installed directly
from the git repository as such:
.. code-block:: sh
git clone https://bitbucket.org/philexander/tikz.git
cd tikz
sudo python setup.py install
The Numfig_ package can be installed directly with pip:
.. code-block:: sh
sudo pip install sphinx-numfig
-----------------------------------
Building Documentation as a Webpage
-----------------------------------
To build the documentation as a webpage (what appears at
http://mit-crpg.github.io/openmc), simply go to the ``docs`` directory and run:
.. code-block:: sh
make html
-------------------------------
Building Documentation as a PDF
-------------------------------
To build PDF documentation, you will need to have a LaTeX distribution installed
on your computer as well as Inkscape_, which is used to convert .svg files to
.pdf files. Inkscape can be installed in a Debian-derivative with:
.. code-block:: sh
sudo apt-get install inkscape
One the pre-requisites are installed, simply go to the ``docs`` directory and
run:
.. code-block:: sh
make latexpdf
.. _Sphinx: http://sphinx-doc.org
.. _sphinxcontrib-tikz: https://bitbucket.org/philexander/tikz
.. _Numfig: https://pypi.python.org/pypi/sphinx_numfig
.. _Inkscape: https://inkscape.org

View file

@ -18,3 +18,4 @@ as debugging.
xml-parsing
statepoint
voxel
docbuild

View file

@ -153,10 +153,6 @@ variables should be set if the default paths are incorrect:
* Default - */opt/phdf5/1.8.14-gnu*
* **PETSC_DIR** - The path to the PETSc directory.
* Default - */opt/petsc/3.5.2-gnu*
To run the full test suite, the following command can be executed in the
tests directory:

View file

@ -22,14 +22,13 @@ package manager`_. Simply enter the following commands into the terminal:
sudo apt-get update
sudo apt-get install openmc
Currently, the binary package does not allow for parallel simulations, HDF5_, or
CMFD acceleration through PETSc_. Users who need such capabilities should build
OpenMC from source as is described in :ref:`usersguide_install`.
Currently, the binary package does not allow for parallel simulations or use of
HDF5_. Users who need such capabilities should build OpenMC from source as is
described in :ref:`usersguide_install`.
.. _Personal Package Archive: https://launchpad.net/~paulromano/+archive/staging
.. _APT package manager: https://help.ubuntu.com/community/AptGet/Howto
.. _HDF5: http://www.hdfgroup.org/HDF5/
.. _PETSc: http://www.mcs.anl.gov/petsc/
-------------------------------------------
Installing from Source on Linux or Mac OS X

View file

@ -1061,7 +1061,7 @@ Here is an example of a properly defined 2d hexagonal lattice:
<pitch> 1.0 </pitch>
<universes>
202
202 202
202 202
202 202 202
202 202
202 101 202
@ -1718,20 +1718,10 @@ It can be turned on with "true" and off with "false".
The ``<gauss_seidel_tolerance>`` element specifies two parameters. The first is
the absolute inner tolerance for Gauss-Seidel iterations when performing CMFD
and the second is the relative inner tolerance for Gauss-Seidel iterations
for CMFD calculations. It is only used in the standalone CMFD power iteration
solver and not when PETSc is active.
for CMFD calculations.
*Default*: 1.e-10 1.e-5
``<ksp_monitor>`` Element
-------------------------
The ``<ksp_monitor>`` element is used to view the convergence of linear GMRES
iterations in PETSc. This option can be turned on with "true" and turned off
with "false".
*Default*: false
``<ktol>`` Element
--------------------
@ -1820,26 +1810,16 @@ iteration. This option can be turned on with "true" and turned off with "false".
-------------------------
The ``<run_adjoint>`` element can be turned on with "true" to have an adjoint
calculation be performed on the last batch when CMFD is active. OpenMC should be
compiled with PETSc when using this option.
calculation be performed on the last batch when CMFD is active.
*Default*: false
``<solver>`` Element
--------------------
The ``<solver>`` element controls whether the CMFD eigenproblem is solved with
standard power iteration or nonlinear Jacobian-free Newton Krylov (JFNK).
By setting "power", power iteration is used and by setting "jfnk", JFNK is used.
*Default*: power
``<shift>`` Element
--------------------
The ``<shfit>`` element specifies an optional Wielandt shift parameter for
accelerating power iterations. It can only be used when PETSc is not active.
It is by default very large so the impact of the shift is effectively zero.
The ``<shift>`` element specifies an optional Wielandt shift parameter for
accelerating power iterations. It is by default very large so the impact of the
shift is effectively zero.
*Default*: 1e6
@ -1848,10 +1828,9 @@ It is by default very large so the impact of the shift is effectively zero.
The ``<spectral>`` element specifies an optional spectral radius that can be set to
accelerate the convergence of Gauss-Seidel iterations during CMFD power iteration
solve. Note this is only used in the standalone CMFD solver and does not affect
the calculation when PETSc is active.
solve.
*Default*: power
*Default*: 0.0
``<stol>`` Element
------------------
@ -1872,10 +1851,9 @@ should be reset.
``<write_matrices>`` Element
----------------------------
The ``<write_matrices>`` element is used to view the PETSc sparse matrices
created when solving CMFD equations. These binary output files can be imported
into MATLAB using PETSc-MATLAB utilities. This option can be
turned on with "true" and off with "false".
The ``<write_matrices>`` element is used to write the sparse matrices created
when solving CMFD equations. This option can be turned on with "true" and off
with "false".
*Default*: false

View file

@ -86,21 +86,6 @@ Prerequisites
You may omit ``--enable-parallel`` if you want to compile HDF5_ in serial.
* PETSc_ for CMFD acceleration
To enable some features of CMFD acceleration, you will need to have
PETSc_ (3.4.2 or higher) installed on your computer. The installed version
will need to have been compiled with the same compiler you intend to
compile OpenMC with. OpenMC requires PETSc_ to be configured with Fortran
datatypes. An example of configuring PETSc_ is listed below::
./configure --prefix=/opt/petsc/3.4.4 --download-f-blas-lapack \
--with-mpi-dir=/opt/mpich/3.1 --with-shared-libraries \
--with-fortran-datatypes
The BLAS/LAPACK library is not required to be downloaded and can be linked
explicitly (e.g., Intel MKL library).
* git_ version control software for obtaining source code
.. _gfortran: http://gcc.gnu.org/wiki/GFortran
@ -108,7 +93,6 @@ Prerequisites
.. _OpenMPI: http://www.open-mpi.org
.. _MPICH: http://www.mpich.org
.. _HDF5: http://www.hdfgroup.org/HDF5/
.. _PETSc: http://www.mcs.anl.gov/petsc/
Obtaining the Source
--------------------
@ -139,10 +123,10 @@ Build Configuration
-------------------
Compiling OpenMC with CMake is carried out in two steps. First, ``cmake`` is run
to determine the compiler, whether optional packages (MPI, HDF5, PETSc) are
available, to generate a list of dependencies between source files so that they
may be compiled in the correct order, and to generate a normal Makefile. The
Makefile is then used by ``make`` to actually carry out the compile and linking
to determine the compiler, whether optional packages (MPI, HDF5) are available,
to generate a list of dependencies between source files so that they may be
compiled in the correct order, and to generate a normal Makefile. The Makefile
is then used by ``make`` to actually carry out the compile and linking
commands. A typical out-of-source build would thus look something like the
following
@ -177,10 +161,6 @@ openmp
Enables shared-memory parallelism using the OpenMP API. The Fortran compiler
being used must support OpenMP.
petsc
Enables PETSc for use in CMFD acceleration. The PETSC_DIR variable should be
set to the base directory of the PETSc installation.
To set any of these options (e.g. turning on debug mode), the following form
should be used:

View file

@ -21,7 +21,6 @@ endif()
option(openmp "Enable shared-memory parallelism with OpenMP" OFF)
option(profile "Compile with profiling flags" OFF)
option(petsc "Enable PETSC for use in CMFD acceleration" OFF)
option(debug "Compile with debug flags" OFF)
option(optimize "Turn on all compiler optimization flags" OFF)
option(verbose "Create verbose Makefiles" OFF)
@ -144,78 +143,6 @@ elseif(CMAKE_Fortran_COMPILER_ID STREQUAL "Cray")
endif()
#===============================================================================
# PETSc for CMFD functionality
#===============================================================================
set (PETSC_ENABLED FALSE)
if(petsc)
set(PETSC_ENABLED TRUE)
find_package(PETSc REQUIRED HINTS $ENV{PETSC_DIR}/conf)
find_library(libpetsc petsc $ENV{PETSC_DIR}/lib)
# If libfblas wasn't found, search the PETSc lib directory
if(PETSC_FBLAS_LIB STREQUAL "PETSC_FBLAS_LIB-NOTFOUND")
find_library(PETSC_FBLAS_LIB fblas $ENV{PETSC_DIR}/lib)
list(REMOVE_ITEM PETSC_PACKAGE_LIBS PETSC_FBLAS_LIB-NOTFOUND)
list(INSERT PETSC_PACKAGE_LIBS 0 ${PETSC_FBLAS_LIB})
endif()
# If libflapack wasn't found, search the PETSc lib directory
if(PETSC_FLAPACK_LIB STREQUAL "PETSC_FLAPACK_LIB-NOTFOUND")
find_library(PETSC_FLAPACK_LIB flapack $ENV{PETSC_DIR}/lib)
list(REMOVE_ITEM PETSC_PACKAGE_LIBS PETSC_FLAPACK_LIB-NOTFOUND)
list(INSERT PETSC_PACKAGE_LIBS 0 ${PETSC_FLAPACK_LIB})
endif()
# If libdl wasn't found, search /usr/lib64
if(PETSC_DL_LIB STREQUAL "PETSC_DL_LIB-NOTFOUND")
find_library(PETSC_DL_LIB libdl.so /usr/lib64)
list(REMOVE_ITEM PETSC_PACKAGE_LIBS PETSC_DL_LIB-NOTFOUND)
list(INSERT PETSC_PACKAGE_LIBS 0 ${PETSC_DL_LIB})
endif()
# If libm wasn't found, search /usr/lib64
if(PETSC_M_LIB STREQUAL "PETSC_M_LIB-NOTFOUND")
find_library(PETSC_M_LIB libm.so /usr/lib64)
list(REMOVE_ITEM PETSC_PACKAGE_LIBS PETSC_M_LIB-NOTFOUND)
list(INSERT PETSC_PACKAGE_LIBS 0 ${PETSC_M_LIB})
endif()
# If libpthread wasn't found, search /usr/lib64
if(PETSC_PTHREAD_LIB STREQUAL "PETSC_PTHREAD_LIB-NOTFOUND")
find_library(PETSC_PTHREAD_LIB libpthread.so /usr/lib64)
list(REMOVE_ITEM PETSC_PACKAGE_LIBS PETSC_PTHREAD_LIB-NOTFOUND)
list(INSERT PETSC_PACKAGE_LIBS 0 ${PETSC_PTHREAD_LIB})
endif()
# If librt wasn't found, search /usr/lib64
if(PETSC_RT_LIB STREQUAL "PETSC_RT_LIB-NOTFOUND")
find_library(PETSC_RT_LIB librt.so /usr/lib64)
list(REMOVE_ITEM PETSC_PACKAGE_LIBS PETSC_RT_LIB-NOTFOUND)
list(INSERT PETSC_PACKAGE_LIBS 0 ${PETSC_RT_LIB})
endif()
# If libssl wasn't found, search /usr/lib64
if(PETSC_SSL_LIB STREQUAL "PETSC_SSL_LIB-NOTFOUND")
find_library(PETSC_SSL_LIB libssl.so /usr/lib64)
list(REMOVE_ITEM PETSC_PACKAGE_LIBS PETSC_SSL_LIB-NOTFOUND)
list(INSERT PETSC_PACKAGE_LIBS 0 ${PETSC_SSL_LIB})
endif()
# If libcrypto wasn't found, search /usr/lib64
if(PETSC_CRYPTO_LIB STREQUAL "PETSC_CRYPTO_LIB-NOTFOUND")
find_library(PETSC_CRYPTO_LIB libcrypto.so /usr/lib64)
list(REMOVE_ITEM PETSC_PACKAGE_LIBS PETSC_CRYPTO_LIB-NOTFOUND)
list(INSERT PETSC_PACKAGE_LIBS 0 ${PETSC_CRYPTO_LIB})
endif()
message("-- Using PETSC: ${libpetsc}")
add_definitions(-DPETSC)
include_directories($ENV{PETSC_DIR}/include)
set(libraries "${libpetsc};${PETSC_PACKAGE_LIBS};${libraries}")
endif()
#===============================================================================
# git SHA1 hash
#===============================================================================
@ -295,14 +222,6 @@ include(CTest)
# Get a list of all the tests to run
file(GLOB_RECURSE TESTS ${CMAKE_CURRENT_SOURCE_DIR}/../tests/test_*.py)
# Check to see if PETSC is compiled for CMFD tests
if (NOT ${PETSC_ENABLED})
file(GLOB_RECURSE CMFD_TESTS ${CMAKE_CURRENT_SOURCE_DIR}/../tests/test_cmfd_jfnk.py)
foreach(cmfd_test in ${CMFD_TESTS})
list(REMOVE_ITEM TESTS ${cmfd_test})
endforeach(cmfd_test)
endif(NOT ${PETSC_ENABLED})
# Check for MEM_CHECK and COVERAGE variables
if (DEFINED ENV{MEM_CHECK})
set(MEM_CHECK $ENV{MEM_CHECK})

View file

@ -1,4 +1,4 @@
module cmfd_execute
module cmfd_execute
!==============================================================================
! CMFD_EXECUTE -- This module is the highest level cmfd module that controls the
@ -20,10 +20,8 @@ contains
subroutine execute_cmfd()
use cmfd_data, only: set_up_cmfd
use cmfd_power_solver, only: cmfd_power_execute
use cmfd_jfnk_solver, only: cmfd_jfnk_execute
use cmfd_solver, only: cmfd_solver_execute
use error, only: warning, fatal_error
use error, only: warning, fatal_error
! CMFD single processor on master
if (master) then
@ -31,35 +29,18 @@ contains
! Start cmfd timer
call time_cmfd % start()
! Create cmfd data from OpenMC tallies
! Create cmfd data from OpenMC tallies
call set_up_cmfd()
! Process solver options
call process_cmfd_options()
! Call solver
#ifdef PETSC
if (trim(cmfd_solver_type) == 'power') then
call cmfd_power_execute()
elseif (trim(cmfd_solver_type) == 'jfnk') then
call cmfd_jfnk_execute()
else
call fatal_error('solver type became invalid after input processing')
end if
#else
call cmfd_solver_execute()
#endif
! Save k-effective
cmfd % k_cmfd(current_batch) = cmfd % keff
! check to perform adjoint on last batch
if (current_batch == n_batches .and. cmfd_run_adjoint) then
if (trim(cmfd_solver_type) == 'power') then
call cmfd_power_execute(adjoint = .true.)
elseif (trim(cmfd_solver_type) == 'jfnk') then
call cmfd_jfnk_execute(adjoint = .true.)
end if
call cmfd_solver_execute(adjoint=.true.)
end if
end if
@ -102,26 +83,6 @@ contains
end subroutine cmfd_init_batch
!==============================================================================
! PROCESS_CMFD_OPTIONS processes user options that interface with PETSc
!==============================================================================
subroutine process_cmfd_options()
#ifdef PETSC
use global, only: cmfd_snes_monitor, cmfd_ksp_monitor, mpi_err
! Check for snes monitor
if (cmfd_snes_monitor) call PetscOptionsSetValue("-snes_monitor", &
"stdout", mpi_err)
! Check for ksp monitor
if (cmfd_ksp_monitor) call PetscOptionsSetValue("-ksp_monitor", &
"stdout", mpi_err)
#endif
end subroutine process_cmfd_options
!===============================================================================
! CALC_FISSION_SOURCE calculates the cmfd fission source
!===============================================================================
@ -130,7 +91,7 @@ contains
use constants, only: CMFD_NOACCEL, ZERO, TWO
use global, only: cmfd, cmfd_coremap, master, entropy_on, current_batch
use string, only: to_str
use string, only: to_str
#ifdef MPI
use global, only: mpi_err
@ -291,11 +252,11 @@ contains
ng = cmfd % indices(4)
! allocate arrays in cmfd object (can take out later extend to multigroup)
if (.not.allocated(cmfd%sourcecounts)) then
if (.not.allocated(cmfd%sourcecounts)) then
allocate(cmfd%sourcecounts(ng,nx,ny,nz))
cmfd % sourcecounts = 0
end if
if (.not.allocated(cmfd % weightfactors)) then
if (.not.allocated(cmfd % weightfactors)) then
allocate(cmfd % weightfactors(ng,nx,ny,nz))
cmfd % weightfactors = ONE
end if
@ -303,7 +264,7 @@ contains
! Compute new weight factors
if (new_weights) then
! Set weight factors to a default 1.0
! Set weight factors to a default 1.0
cmfd%weightfactors = ONE
! Count bank sites in mesh and reverse due to egrid structure

View file

@ -2,10 +2,6 @@ module cmfd_input
use global
#ifdef PETSC
use petscsys
#endif
implicit none
private
public :: configure_cmfd
@ -13,7 +9,7 @@ module cmfd_input
contains
!===============================================================================
! CONFIGURE_CMFD initializes PETSc and CMFD parameters
! CONFIGURE_CMFD initializes CMFD parameters
!===============================================================================
subroutine configure_cmfd()
@ -32,17 +28,6 @@ contains
color = 2
end if
! Split up procs
#ifdef PETSC
call MPI_COMM_SPLIT(MPI_COMM_WORLD, color, 0, cmfd_comm, mpi_err)
! assign to PETSc
PETSC_COMM_WORLD = cmfd_comm
! Initialize PETSc on all procs
call PetscInitialize(PETSC_NULL_CHARACTER, mpi_err)
#endif
! Initialize timers
call time_cmfd % reset()
call time_cmfdbuild % reset()
@ -165,7 +150,7 @@ contains
cmfd_downscatter = .true.
end if
! Reset dhat parameters
! Reset dhat parameters
if (check_for_node(doc, "dhat_reset")) then
call get_node_value(doc, "dhat_reset", temp_str)
temp_str = to_lower(temp_str)
@ -173,23 +158,7 @@ contains
dhat_reset = .true.
end if
! Set the solver type
if (check_for_node(doc, "solver")) &
call get_node_value(doc, "solver", cmfd_solver_type)
! Set monitoring
if (check_for_node(doc, "snes_monitor")) then
call get_node_value(doc, "snes_monitor", temp_str)
temp_str = to_lower(temp_str)
if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') &
cmfd_snes_monitor = .true.
end if
if (check_for_node(doc, "ksp_monitor")) then
call get_node_value(doc, "ksp_monitor", temp_str)
temp_str = to_lower(temp_str)
if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') &
cmfd_ksp_monitor = .true.
end if
if (check_for_node(doc, "power_monitor")) then
call get_node_value(doc, "power_monitor", temp_str)
temp_str = to_lower(temp_str)
@ -210,9 +179,6 @@ contains
call get_node_value(doc, "run_adjoint", temp_str)
temp_str = to_lower(temp_str)
if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') &
#ifndef PETSC
call fatal_error('Must use PETSc when running adjoint option.')
#endif
cmfd_run_adjoint = .true.
end if
@ -238,12 +204,6 @@ contains
! Get display
if (check_for_node(doc, "display")) &
call get_node_value(doc, "display", cmfd_display)
if (trim(cmfd_display) == 'dominance' .and. &
trim(cmfd_solver_type) /= 'power') then
if (master) call warning('Dominance Ratio only aviable with power &
&iteration solver')
cmfd_display = ''
end if
! Read in spectral radius estimate and tolerances
if (check_for_node(doc, "spectral")) &

View file

@ -1,415 +0,0 @@
module cmfd_jfnk_solver
use cmfd_loss_operator, only: init_loss_matrix, build_loss_matrix
use cmfd_power_solver, only: cmfd_power_execute
use cmfd_prod_operator, only: init_prod_matrix, build_prod_matrix
use matrix_header, only: Matrix
use solver_interface, only: JFNKSolver, Jfnk_ctx
use vector_header, only: Vector
implicit none
private
public :: cmfd_jfnk_execute
logical :: adjoint_calc ! true if an adjoint is to be calculated
type(Jfnk_ctx) :: jfnk_data ! object that holds pointers to user routines
type(Matrix) :: jac_prec ! Jacobian preconditioner object
type(Matrix) :: loss ! CMFD loss matrix
type(Matrix) :: prod ! CMFD production matrix
#ifdef PETSC
type(JFNKSolver) :: jfnk ! JFNK solver object
#endif
type(Vector) :: resvec ! JFNK residual vector
type(Vector) :: xvec ! JFNK solution vector
contains
!===============================================================================
! CMFD_JFNK_EXECUTE main routine for JFNK solver
!===============================================================================
subroutine cmfd_jfnk_execute(adjoint)
use global, only: time_cmfdbuild, time_cmfdsolve
logical, intent(in), optional :: adjoint ! adjoint calculation
! Check for adjoint
adjoint_calc = .false.
if (present(adjoint)) adjoint_calc = adjoint
! Seed with power iteration to help converge to fundamental mode
call cmfd_power_execute(k_tol=1.E-3_8, s_tol=1.E-3_8, adjoint=adjoint_calc)
! Start timer for build
call time_cmfdbuild % start()
! Initialize data
call init_data()
! Initialize solver
#ifdef PETSC
call jfnk % create()
#endif
! Set up residual and jacobian routines
#ifdef PETSC
jfnk_data % res_proc_ptr => compute_nonlinear_residual
jfnk_data % jac_proc_ptr => build_jacobian_matrix
call jfnk % set_functions(jfnk_data, resvec, jac_prec)
#endif
! Stop timer for build
call time_cmfdbuild % stop()
! Solve the system
call time_cmfdsolve % start()
#ifdef PETSC
call jfnk % solve(xvec)
#endif
call time_cmfdsolve % stop()
! Extracts results to cmfd object
call extract_results()
! Deallocate all slepc data
call finalize()
end subroutine cmfd_jfnk_execute
!===============================================================================
! INIT_DATA allocates matrices vectors for CMFD solution
!===============================================================================
subroutine init_data()
use constants, only: ZERO, ONE
use global, only: cmfd, cmfd_adjoint_type, current_batch
logical :: physical_adjoint ! physical adjoing calculation logical
integer :: n ! size of matrices
! Set up all matrices
call init_loss_matrix(loss)
call init_prod_matrix(prod)
call init_jacobian_matrix()
! Check for physical adjoint
physical_adjoint = .false.
if (adjoint_calc .and. trim(cmfd_adjoint_type) == 'physical') &
physical_adjoint = .true.
! Create matrix operators
call build_loss_matrix(loss, adjoint = physical_adjoint)
call build_prod_matrix(prod, adjoint = physical_adjoint)
! Assemble matrices and use PETSc
call loss % assemble()
call prod % assemble()
#ifdef PETSC
call loss % setup_petsc()
call prod % setup_petsc()
#endif
! Check for mathematical adjoint
if (adjoint_calc .and. trim(cmfd_adjoint_type) == 'math') &
call compute_adjoint()
! Get problem size
n = loss % n
! Create problem vectors
call resvec % create(n + 1)
call xvec % create(n + 1)
! Set flux in guess from rough power iteration
if (adjoint_calc) then
xvec % val(1:n) = cmfd % adj_phi
else
xvec % val(1:n) = cmfd % phi
end if
! Set keff in guess from rough power iteration
if (adjoint_calc) then
xvec % val(n + 1) = ONE/cmfd % adj_keff
else
xvec % val(n + 1) = ONE/cmfd % keff
end if
! Set up vectors for PETSc
#ifdef PETSC
call resvec % setup_petsc()
call xvec % setup_petsc()
#endif
! Build jacobian from initial guess
call build_jacobian_matrix(xvec)
! Set up Jacobian for PETSc
#ifdef PETSC
call jac_prec % setup_petsc()
#endif
! Set dominance ratio to 0
cmfd % dom(current_batch) = ZERO
end subroutine init_data
!==============================================================================
! INIT_JACOBIAN_MATRIX preallocates jacobian matrix and initializes it
!==============================================================================
subroutine init_jacobian_matrix()
integer :: nnz ! number of nonzeros in matrix
integer :: n ! dimension of matrix
! Get length of matrix and number of nonzeros total in loss matrix
nnz = loss % nnz
n = loss % n
! There is one more nonzero for each row and and additional row of nonzeros
nnz = nnz + 2*n + 1
! We need 1 more row for the jacobian
n = n + 1
! Configure Jacobian matrix
call jac_prec % create(n, nnz)
end subroutine init_jacobian_matrix
!===============================================================================
! BUILD_JACOBIAN_MATRIX creates the Jacobian of nonlinear eigenvalue problem
!===============================================================================
subroutine build_jacobian_matrix(x)
use constants, only: ONE
type(Vector), intent(in) :: x ! solution vector
integer :: i ! loop counter for jacobian rows
integer :: jjac ! loop counter for jacobian cols
integer :: jloss ! loop counter for loss matrix cols
integer :: jprod ! loop counter for prod matrix cols
integer :: n ! problem size
real(8) :: lambda ! eigenvalue
real(8) :: val ! temporary real scalar
type(Vector) :: flux ! flux vector
type(Vector) :: fsrc ! fission source vector
! Get the problem size
n = loss % n
! Get flux and eigenvalue
flux % val => x % val(1:n)
lambda = x % val(n + 1)
! Create fission source vector
call fsrc % create(n)
call prod % vector_multiply(flux, fsrc)
! Reset counters in Jacobian matrix
jac_prec % n_count = 1
jac_prec % nz_count = 1
! Begin loop around rows of Jacobian matrix
ROWS: do i = 1, n
! Add a new row in Jacobian
call jac_prec % new_row()
! Begin loop around columns of loss matrix
COLS_LOSS: do jloss = loss % get_row(i), loss % get_row(i + 1) - 1
! Start with the value in the loss matrix
val = loss % val(jloss)
! Loop around columns of prod matrix
COLS_PROD: do jprod = prod % get_row(i), prod % get_row(i + 1) - 1
! See if columns agree with loss matrix
if (prod % get_col(jprod) == loss % get_col(jloss)) then
val = val - lambda*prod % val(jprod)
exit
end if
end do COLS_PROD
! Record value in jacobian
call jac_prec % add_value(loss % get_col(jloss), val)
end do COLS_LOSS
! Add fission source value
val = -fsrc % val(i)
call jac_prec % add_value(n + 1, val)
end do ROWS
! Need to add negative transpose of flux vector on last row
call jac_prec % new_row()
do jjac = 1, n
val = -flux % val(jjac)
call jac_prec % add_value(jjac, val)
end do
! Add unity on bottom right corner of matrix
call jac_prec % add_value(n + 1, ONE)
! CRS requires a final value in row
call jac_prec % new_row()
! Free all allocated memory
flux % val => null()
call fsrc % destroy()
end subroutine build_jacobian_matrix
!===============================================================================
! COMPUTE_NONLINEAR_RESIDUAL computes the residual of the nonlinear equations
!===============================================================================
#ifdef PETSC
subroutine compute_nonlinear_residual(x, res)
use global, only: cmfd_write_matrices
type(Vector), intent(in) :: x ! solution vector
type(Vector), intent(inout) :: res ! residual vector
character(len=25) :: filename
integer :: n
real(8) :: lambda
type(Vector) :: res_loss
type(Vector) :: res_prod
type(Vector) :: flux
! Get problem size
n = loss % n
! Set up temporary vectors
call res_loss % create(n)
call res_prod % create(n)
! Extract flux
flux % n = n
flux % val => x % val(1:n)
! Extract eigenvalue
lambda = x % val(n + 1)
! Calculate M*flux then F*flux
call loss % vector_multiply(flux, res_loss)
call prod % vector_multiply(flux, res_prod)
! Put flux component values in residual vector
res % val(1:n) = res_loss % val - lambda*res_prod % val
! Put eigenvalue component in residual vector
res % val(n+1) = 0.5_8 - 0.5_8*sum(flux % val * flux % val)
! Write out data in binary files (debugging)
if (cmfd_write_matrices) then
! Write out residual vector
if (adjoint_calc) then
filename = 'adj_res.bin'
else
filename = 'res.bin'
end if
#ifdef PETSC
call res % write_petsc_binary(filename)
#endif
! Write out solution vector
if (adjoint_calc) then
filename = 'adj_x.bin'
else
filename = 'x.bin'
end if
#ifdef PETSC
call x % write_petsc_binary(filename)
#endif
end if
end subroutine compute_nonlinear_residual
#endif
!===============================================================================
! COMPUTE_ADJOINT calculates mathematical adjoint by taking transposes
!===============================================================================
subroutine compute_adjoint()
use global, only: cmfd_write_matrices
! Transpose matrices
#ifdef PETSC
call loss % transpose()
call prod % transpose()
#endif
! Write out matrix in binary file (debugging)
if (cmfd_write_matrices) then
#ifdef PETSC
call loss % write_petsc_binary('adj_lossmat.bin')
call loss % write_petsc_binary('adj_prodmat.bin')
#endif
end if
end subroutine compute_adjoint
!===============================================================================
! EXTRACT_RESULTS gets the results and puts them in global CMFD object
!===============================================================================
subroutine extract_results()
use constants, only: ZERO, ONE
use global, only: cmfd
integer :: n ! problem size
! Get problem size
n = loss % n
! Also allocate in cmfd object
if (adjoint_calc) then
if (.not. allocated(cmfd % adj_phi)) allocate(cmfd % adj_phi(n))
else
if (.not. allocated(cmfd % phi)) allocate(cmfd % phi(n))
end if
! Get flux and eigenvalue
if (adjoint_calc) then
cmfd % adj_phi = xvec % val(1:n)
cmfd % adj_keff = ONE / xvec % val(n+1)
else
cmfd % phi = xvec % val(1:n)
cmfd % keff = ONE / xvec % val(n+1)
end if
end subroutine extract_results
!===============================================================================
! FINALIZE frees all memory from a JFNK calculation
!===============================================================================
subroutine finalize()
call loss % destroy()
call prod % destroy()
call jac_prec % destroy()
call xvec % destroy()
call resvec % destroy()
#ifdef PETSC
call jfnk % destroy()
#endif
nullify(jfnk_data % res_proc_ptr)
nullify(jfnk_data % jac_proc_ptr)
end subroutine finalize
end module cmfd_jfnk_solver

View file

@ -1,356 +0,0 @@
module cmfd_power_solver
! This module contains routines to execute the power iteration solver
use cmfd_loss_operator, only: init_loss_matrix, build_loss_matrix
use cmfd_prod_operator, only: init_prod_matrix, build_prod_matrix
use matrix_header, only: Matrix
use solver_interface, only: GMRESSolver
use vector_header, only: Vector
implicit none
private
public :: cmfd_power_execute
logical :: iconv ! did the problem converged
real(8) :: k_n ! new k-eigenvalue
real(8) :: k_o ! old k-eigenvalue
real(8) :: ktol = 1.e-8_8 ! tolerance on keff
real(8) :: stol = 1.e-8_8 ! tolerance on source
real(8) :: norm_n ! current norm of source vector
real(8) :: norm_o ! old norm of source vector
real(8) :: kerr ! error in keff
real(8) :: serr ! error in source
logical :: adjoint_calc ! run an adjoint calculation
type(Matrix) :: loss ! cmfd loss matrix
type(Matrix) :: prod ! cmfd prod matrix
type(Vector) :: phi_n ! new flux vector
type(Vector) :: phi_o ! old flux vector
type(Vector) :: s_n ! new source vector
type(Vector) :: s_o ! old flux vector
type(Vector) :: serr_v ! error in source
#ifdef PETSC
type(GMRESSolver) :: gmres ! gmres solver
#endif
contains
!===============================================================================
! CMFD_POWER_EXECUTE sets up and runs power iteration solver for CMFD
!===============================================================================
subroutine cmfd_power_execute(k_tol, s_tol, adjoint)
use global, only: cmfd_adjoint_type, time_cmfdbuild, time_cmfdsolve
real(8), intent(in), optional :: k_tol ! tolerance on keff
real(8), intent(in), optional :: s_tol ! tolerance on source
logical, intent(in), optional :: adjoint ! adjoint calc
logical :: physical_adjoint = .false. ! physical adjoint default false
! Set tolerances if present
if (present(k_tol)) ktol = k_tol
if (present(s_tol)) stol = s_tol
! Check for adjoint execution
adjoint_calc = .false.
if (present(adjoint)) adjoint_calc = adjoint
! Check for physical adjoint
if (adjoint_calc .and. trim(cmfd_adjoint_type) == 'physical') &
physical_adjoint = .true.
! Start timer for build
call time_cmfdbuild % start()
! Initialize solver
#ifdef PETSC
call gmres % create()
#endif
! Initialize matrices and vectors
call init_data(physical_adjoint)
! Check for mathematical adjoint calculation
if (adjoint_calc .and. trim(cmfd_adjoint_type) == 'math') &
call compute_adjoint()
! Set up krylov info
#ifdef PETSC
call gmres % set_oper(loss, loss)
#endif
! Stop timer for build
call time_cmfdbuild % stop()
! Begin power iteration
call time_cmfdsolve % start()
call execute_power_iter()
call time_cmfdsolve % stop()
! Extract results
call extract_results()
! Deallocate data
call finalize()
end subroutine cmfd_power_execute
!===============================================================================
! INIT_DATA allocates matrices and vectors for CMFD solution
!===============================================================================
subroutine init_data(adjoint)
use constants, only: ONE, ZERO
#ifdef PETSC
use global, only: cmfd_write_matrices
#endif
logical, intent(in) :: adjoint ! adjoint calcualtion
integer :: n ! problem size
real(8) :: guess ! initial guess
! Set up matrices
call init_loss_matrix(loss)
call init_prod_matrix(prod)
! Get problem size
n = loss % n
! Set up flux vectors
call phi_n % create(n)
call phi_o % create(n)
! Set up source vectors
call s_n % create(n)
call s_o % create(n)
call serr_v % create(n)
! Set initial guess
guess = ONE
phi_n % val = guess
phi_o % val = guess
k_n = guess
k_o = guess
! Fill in loss matrix
call build_loss_matrix(loss, adjoint=adjoint)
! Fill in production matrix
call build_prod_matrix(prod, adjoint=adjoint)
! Setup petsc for everything
call loss % assemble()
call prod % assemble()
#ifdef PETSC
call loss % setup_petsc()
call prod % setup_petsc()
call phi_n % setup_petsc()
call phi_o % setup_petsc()
call s_o % setup_petsc()
call s_n % setup_petsc()
if (cmfd_write_matrices) call loss % write_petsc_binary('loss.bin')
if (cmfd_write_matrices) call prod % write_petsc_binary('prod.bin')
#endif
! Set norms to 0
norm_n = ZERO
norm_o = ZERO
end subroutine init_data
!===============================================================================
! COMPUTE_ADJOINT computes a mathematical adjoint of CMFD problem
!===============================================================================
subroutine compute_adjoint()
use global, only: cmfd_write_matrices
! Transpose matrices
#ifdef PETSC
call loss % transpose()
call prod % transpose()
#endif
! Write out matrix in binary file (debugging)
if (cmfd_write_matrices) then
#ifdef PETSC
call loss % write_petsc_binary('adj_lossmat.bin')
call prod % write_petsc_binary('adj_prodmat.bin')
#endif
end if
end subroutine compute_adjoint
!===============================================================================
! EXECUTE_POWER_ITER is the main power iteration routine
! for the cmfd calculation
!===============================================================================
subroutine execute_power_iter()
integer :: i ! iteration counter
! Reset convergence flag
iconv = .false.
! Begin power iteration
do i = 1, 10000
! Compute source vector
call prod % vector_multiply(phi_o, s_o)
! Normalize source vector
s_o % val = s_o % val / k_o
! Compute new flux vector
#ifdef PETSC
call gmres % solve(s_o, phi_n)
#endif
! Compute new source vector
call prod % vector_multiply(phi_n, s_n)
! Compute new k-eigenvalue
k_n = sum(s_n % val) / sum(s_o % val)
! Renormalize the old source
s_o % val = s_o % val * k_o
! Check convergence
call convergence(i)
! Break loop if converged
if (iconv) exit
! Record old values
phi_o % val = phi_n % val
k_o = k_n
norm_o = norm_n
end do
end subroutine execute_power_iter
!===============================================================================
! CONVERGENCE checks the convergence of the CMFD problem
!===============================================================================
subroutine convergence(iter)
use constants, only: ONE, ZERO
use global, only: cmfd_power_monitor, master
use, intrinsic :: ISO_FORTRAN_ENV
integer, intent(in) :: iter ! iteration number
! Reset convergence flag
iconv = .false.
! Calculate error in keff
kerr = abs(k_o - k_n)/k_n
! Calculate max error in source
where (s_n % val > ZERO)
serr_v % val = ((s_n % val - s_o % val)/s_n % val)**2
end where
serr = sqrt(ONE/dble(s_n % n) * sum(serr_v % val))
! Check for convergence
if(kerr < ktol .and. serr < stol) iconv = .true.
! Save the L2 norm of the source
norm_n = serr
! Print out to user
if (cmfd_power_monitor .and. master) then
write(OUTPUT_UNIT,FMT='(I0,":",T10,"k-eff: ",F0.8,T30,"k-error: ", &
&1PE12.5,T55, "src-error: ",1PE12.5)') iter, k_n, kerr, serr
end if
end subroutine convergence
!===============================================================================
! EXTRACT_RESULTS takes results and puts them in CMFD global data object
!===============================================================================
subroutine extract_results()
use global, only: cmfd, cmfd_write_matrices, current_batch
character(len=25) :: filename ! name of file to write data
integer :: n ! problem size
! Get problem size
n = loss % n
! Allocate in cmfd object if not already allocated
if (adjoint_calc) then
if (.not. allocated(cmfd%adj_phi)) allocate(cmfd%adj_phi(n))
else
if (.not. allocated(cmfd%phi)) allocate(cmfd%phi(n))
end if
! Save values
if (adjoint_calc) then
cmfd % adj_phi = phi_n % val
else
cmfd % phi = phi_n % val
end if
! Save eigenvalue
if(adjoint_calc) then
cmfd%adj_keff = k_n
else
cmfd%keff = k_n
end if
! Normalize phi to 1
if (adjoint_calc) then
cmfd%adj_phi = cmfd%adj_phi/sqrt(sum(cmfd%adj_phi*cmfd%adj_phi))
else
cmfd%phi = cmfd%phi/sqrt(sum(cmfd%phi*cmfd%phi))
end if
! Save dominance ratio
cmfd % dom(current_batch) = norm_n/norm_o
! Write out results
if (cmfd_write_matrices) then
if (adjoint_calc) then
filename = 'adj_fluxvec.bin'
else
filename = 'fluxvec.bin'
end if
#ifdef PETSC
call phi_n % write_petsc_binary(filename)
#endif
end if
end subroutine extract_results
!===============================================================================
! FINALIZE frees all memory associated with power iteration
!===============================================================================
subroutine finalize()
! Destroy all objects
#ifdef PETSC
call gmres % destroy()
#endif
call loss % destroy()
call prod % destroy()
call phi_n % destroy()
call phi_o % destroy()
call s_n % destroy()
call s_o % destroy()
call serr_v % destroy
end subroutine finalize
end module cmfd_power_solver

View file

@ -1,214 +0,0 @@
module cmfd_slepc_solver
#ifdef SLEPC
use cmfd_loss_operator, only: loss_operator, init_M_operator, &
build_loss_matrix, destroy_M_operator
use cmfd_prod_operator, only: prod_operator, init_F_operator, &
build_prod_matrix, destroy_F_operator
implicit none
#include <finclude/petsc.h90>
#include <finclude/slepcsys.h>
#include <finclude/slepceps.h>
type(loss_operator) :: loss
type(prod_operator) :: prod
Vec :: phi ! eigenvector
EPS :: eps ! slepc eigenvalue object
ST :: st ! slepc spectral trans object
KSP :: ksp ! linear solver object
PC :: pc ! preconditioner object
integer :: ierr ! error flag
real(8) :: keff ! the converged eigenvalue
contains
!===============================================================================
! CMFD_SLEPC_EXECUTE
!===============================================================================
subroutine cmfd_slepc_execute()
use global, only: time_cmfd, master
call time_cmfd % start()
! initialize data
call init_data()
! initialize solver
call init_solver()
!all timer_stop(time_cmfd)
!f(master) print *,'Init Time:',time_cmfd%elapsed
!all timer_reset(time_cmfd)
!all timer_start(time_cmfd)
! build operators
call build_loss_matrix(loss)
call build_prod_matrix(prod)
!all timer_stop(time_cmfd)
!f(master) print *,'Build Time',time_cmfd%elapsed
!all timer_reset(time_cmfd)
!all timer_start(time_cmfd)
! set operators to EPS object
call EPSSetOperators(eps, prod%F, loss%M, ierr)
! set EPS options
call EPSSetFromOptions(eps, ierr)
! solve the system
call EPSSolve(eps, ierr)
!all timer_stop(time_cmfd)
!f(master) print *,'Solve Time:',time_cmfd%elapsed
!all timer_reset(time_cmfd)
!all timer_start(time_cmfd)
! extracts results to cmfd object
call extract_results()
!all timer_stop(time_cmfd)
!f(master) print *,'Extraction Time:',time_cmfd%elapsed
!all timer_reset(time_cmfd)
! deallocate all slepc data
call finalize()
end subroutine cmfd_slepc_execute
!===============================================================================
! INIT_DATA allocates matrices vectors for CMFD solution
!===============================================================================
subroutine init_data()
integer :: n ! problem size
! set up matrices
call init_M_operator(loss)
call init_F_operator(prod)
! get problem size
n = loss%n
! set up eigenvector
call VecCreateMPI(PETSC_COMM_WORLD, PETSC_DECIDE, n, phi, ierr)
end subroutine init_data
!===============================================================================
! INIT_SOLVER
!===============================================================================
subroutine init_solver()
character(LEN=20) :: epstype, sttype, ksptype, pctype
! create EPS Object
call EPSCreate(PETSC_COMM_WORLD, eps, ierr)
call EPSSetProblemType(eps, EPS_GNHEP, ierr)
call EPSSetType(eps, EPSPOWER, ierr)
call EPSSetWhichEigenpairs(eps, EPS_LARGEST_MAGNITUDE, ierr)
! get ST, KSP and PC objects
call EPSGetST(eps, st, ierr)
call STGetKSP(st, ksp, ierr)
call KSPGetPC(ksp, pc, ierr)
! set GMRES default
call KSPSetType(ksp, KSPGMRES, ierr)
! set precursor type
call PCSetType(pc, PCHYPRE, ierr)
call PCSetFromOptions(pc, ierr)
! get all types and print
call EPSGetType(eps, epstype, ierr)
call STGetType(st, sttype, ierr)
call KSPGetType(ksp, ksptype, ierr)
call PCGetType(pc, pctype, ierr)
! display information to user
! write(*,'(/,A)') 'SLEPC SOLVER OPTIONS:'
! write(*,*) '---------------------'
! write(*,*) 'EPS TYPE IS: ',epstype
! write(*,*) 'ST TYPE IS: ',sttype
! write(*,*) 'KSP TYPE IS: ',ksptype
! write(*,*) 'PC TYPE IS: ',pctype
end subroutine init_solver
!==============================================================================
! EXTRACT_RESULTS
!==============================================================================
subroutine extract_results()
use constants, only: ZERO
use global, only: cmfd, path_input, master
integer :: n ! problem size
integer :: i_eig = 0 ! eigenvalue to extract
integer :: row_start ! starting row
integer :: row_end ! ending row
real(8),allocatable :: mybuf(:) ! temp buffer
! PetscViewer :: viewer ! petsc output object
PetscScalar, pointer :: phi_v(:) ! pointer to eigenvector info
! get problem size
n = loss%n
! also allocate in cmfd object
if (.not. allocated(cmfd%phi)) allocate(cmfd%phi(n))
if (.not. allocated(mybuf)) allocate(mybuf(n))
! zero out cmfd object
cmfd%phi = ZERO
! extract run information
call EPSGetEigenpair(eps, i_eig, keff, PETSC_NULL_DOUBLE, phi, &
PETSC_NULL_OBJECT, ierr)
! get ownership range
call VecGetOwnershipRange(phi, row_start, row_end, ierr)
! convert petsc phi_object to cmfd_obj
call VecGetArrayF90(phi, phi_v, ierr)
cmfd%phi(row_start+1:row_end) = phi_v
call VecRestoreArrayF90(phi, phi_v, ierr)
! save eigenvalue
cmfd%keff = keff
! reduce result to master
mybuf = ZERO
call MPI_ALLREDUCE(cmfd%phi, mybuf, n, MPI_REAL8, MPI_SUM, &
PETSC_COMM_WORLD, ierr)
! move buffer to object and deallocate
cmfd%phi = mybuf
if(allocated(mybuf)) deallocate(mybuf)
! write out results
! call PetscViewerBinaryOpen(PETSC_COMM_SELF,trim(path_input)//'fluxvec.bin' &
! & ,FILE_MODE_WRITE,viewer,ierr)
! call VecView(phi,viewer,ierr)
! call PetscViewerDestroy(viewer,ierr)
end subroutine extract_results
!==============================================================================
! FINALIZE
!==============================================================================
subroutine finalize()
! finalize data objects
call destroy_M_operator(loss)
call destroy_F_operator(prod)
call VecDestroy(phi, ierr)
! finalize solver objects
call EPSDestroy(eps, ierr)
end subroutine finalize
#endif
end module cmfd_slepc_solver

View file

@ -33,7 +33,6 @@ module cmfd_solver
type(Vector) :: serr_v ! Error in source
! CMFD linear solver interface
procedure(linsolve), pointer :: cmfd_linsolver => null()
abstract interface
subroutine linsolve(A, b, x, tol, i)
import :: Matrix
@ -101,7 +100,7 @@ contains
subroutine init_data(adjoint)
use constants, only: ONE, ZERO
use global, only: cmfd, cmfd_shift, keff, cmfd_ktol, cmfd_stol, &
use global, only: cmfd_shift, keff, cmfd_ktol, cmfd_stol, &
cmfd_write_matrices
logical, intent(in) :: adjoint
@ -155,16 +154,6 @@ contains
norm_n = ZERO
norm_o = ZERO
! Set up solver
select case(cmfd % indices(4))
case(1)
cmfd_linsolver => cmfd_linsolver_1g
case(2)
cmfd_linsolver => cmfd_linsolver_2g
case default
cmfd_linsolver => cmfd_linsolver_ng
end select
! Set tolerances
ktol = cmfd_ktol
stol = cmfd_stol
@ -178,23 +167,17 @@ contains
subroutine compute_adjoint()
use error, only: fatal_error
#ifdef PETSC
use global, only: cmfd_write_matrices
#endif
#ifdef PETSC
! Transpose matrices
call loss % transpose()
call prod % transpose()
loss = loss % transpose()
prod = prod % transpose()
! Write out matrix in binary file (debugging)
if (cmfd_write_matrices) then
call loss % write_petsc_binary('adj_lossmat.bin')
call prod % write_petsc_binary('adj_prodmat.bin')
call loss % write('adj_loss.dat')
call prod % write('adj_prod.dat')
end if
#else
call fatal_error('Adjoint calculations only allowed with PETSc')
#endif
end subroutine compute_adjoint
@ -207,7 +190,7 @@ contains
use constants, only: ONE
use error, only: fatal_error
use global, only: cmfd_atoli, cmfd_rtoli
use global, only: cmfd, cmfd_atoli, cmfd_rtoli
integer :: i ! iteration counter
integer :: innerits ! # of inner iterations
@ -245,7 +228,14 @@ contains
s_o % val = s_o % val / k_lo
! Compute new flux vector
call cmfd_linsolver(loss, s_o, phi_n, toli, innerits)
select case(cmfd % indices(4))
case(1)
call cmfd_linsolver_1g(loss, s_o, phi_n, toli, innerits)
case(2)
call cmfd_linsolver_2g(loss, s_o, phi_n, toli, innerits)
case default
call cmfd_linsolver_ng(loss, s_o, phi_n, toli, innerits)
end select
! Compute new source vector
call prod % vector_multiply(phi_n, s_n)
@ -743,13 +733,11 @@ contains
! Write out results
if (cmfd_write_matrices) then
if (adjoint_calc) then
filename = 'adj_fluxvec.bin'
filename = 'adj_fluxvec.dat'
else
filename = 'fluxvec.bin'
filename = 'fluxvec.dat'
end if
#ifdef PETSC
call phi_n % write_petsc_binary(filename)
#endif
! TODO: call phi_n % write(filename)
end if
end subroutine extract_results

View file

@ -38,14 +38,6 @@ contains
if (check_overlaps) call reduce_overlap_count()
end if
#ifdef PETSC
! Finalize PETSc
if (cmfd_run) then
call PetscFinalize(mpi_err)
call MPI_COMM_FREE(cmfd_comm, mpi_err)
end if
#endif
! Stop timers and show timing statistics
call time_finalize % stop()
call time_total % stop()

View file

@ -317,9 +317,6 @@ module global
! Is CMFD active
logical :: cmfd_run = .false.
! CMFD communicator
integer :: cmfd_comm
! Timing objects
type(Timer) :: time_cmfd ! timer for whole cmfd calculation
type(Timer) :: time_cmfdbuild ! timer for matrix build
@ -338,9 +335,6 @@ module global
integer :: n_cmfd_meshes = 1 ! # of structured meshes
integer :: n_cmfd_tallies = 3 ! # of user-defined tallies
! Eigenvalue solver type
character(len=10) :: cmfd_solver_type = 'power'
! Adjoint method type
character(len=10) :: cmfd_adjoint_type = 'physical'
@ -358,8 +352,6 @@ module global
logical :: cmfd_downscatter = .false.
! Convergence monitoring
logical :: cmfd_snes_monitor = .false.
logical :: cmfd_ksp_monitor = .false.
logical :: cmfd_power_monitor = .false.
! Cmfd output

View file

@ -438,9 +438,6 @@ contains
case ('-v', '-version', '--version')
call print_version()
stop
case ('-eps_tol', '-ksp_gmres_restart')
! Handle options that would be based to PETSC
i = i + 1
case ('-t', '-track', '--track')
write_all_tracks = .true.
case default
@ -561,7 +558,7 @@ contains
integer :: k ! loop index for lattices
integer :: m ! loop index for lattices
integer :: lid ! lattice IDs
integer :: i_array ! index in surfaces/materials array
integer :: i_array ! index in surfaces/materials array
integer :: id ! user-specified id
type(Cell), pointer :: c => null()
class(Lattice), pointer :: lat => null()

View file

@ -2934,7 +2934,7 @@ contains
end select
! Append the score to the list of possible trigger scores
call trigger_scores % add_key(trim(score_name), l)
if (trigger_on) call trigger_scores % add_key(trim(score_name), l)
end do
t % n_score_bins = n_scores

View file

@ -1,9 +1,5 @@
module matrix_header
#ifdef PETSC
use petscmat
#endif
implicit none
private
@ -15,10 +11,6 @@ module matrix_header
integer, allocatable :: row(:) ! csr row vector
integer, allocatable :: col(:) ! column vector
real(8), allocatable :: val(:) ! matrix value vector
# ifdef PETSC
type(mat) :: petsc_mat
# endif
logical :: petsc_active
contains
procedure :: create => matrix_create
procedure :: destroy => matrix_destroy
@ -31,17 +23,9 @@ module matrix_header
procedure :: search_indices => matrix_search_indices
procedure :: write => matrix_write
procedure :: copy => matrix_copy
# ifdef PETSC
procedure :: setup_petsc => matrix_setup_petsc
procedure :: write_petsc_binary => matrix_write_petsc_binary
procedure :: transpose => matrix_transpose
# endif
end type matrix
#ifdef PETSC
integer :: petsc_err
#endif
contains
!===============================================================================
@ -67,9 +51,6 @@ contains
self % n = n
self % nnz = nnz
! Set PETSc active by default to false
self % petsc_active = .false.
end subroutine matrix_create
!===============================================================================
@ -80,10 +61,6 @@ contains
class(Matrix), intent(inout) :: self ! matrix instance
#ifdef PETSC
if (self % petsc_active) call MatDestroy(self % petsc_mat, petsc_err)
#endif
if (allocated(self % row)) deallocate(self % row)
if (allocated(self % col)) deallocate(self % col)
if (allocated(self % val)) deallocate(self % val)
@ -104,10 +81,6 @@ contains
self % col(self % nz_count) = col
self % val(self % nz_count) = val
! Need to adjust column indices if PETSc is active
if (self % petsc_active) self % col(self % nz_count) = &
self % col(self % nz_count) - 1
! Increment the number of nonzeros currently stored
self % nz_count = self % nz_count + 1
@ -124,10 +97,6 @@ contains
! Record the current number of nonzeros
self % row(self % n_count) = self % nz_count
! If PETSc is active, we have to reference indices off 0
if (self % petsc_active) self % row(self % n_count) = &
self % row(self % n_count) - 1
! Increment the current row that we are on
self % n_count = self % n_count + 1
@ -239,7 +208,7 @@ contains
end subroutine split
!===============================================================================
! MATRIX_GET_ROW is a method to get row and checks for PETSc C indexing use
! MATRIX_GET_ROW is a method to get row
!===============================================================================
function matrix_get_row(self, i) result(row)
@ -250,12 +219,10 @@ contains
row = self % row(i)
if (self % petsc_active) row = row + 1
end function matrix_get_row
!===============================================================================
! MATRIX_GET_COL is a method to get column and checks for PETSc C index use
! MATRIX_GET_COL is a method to get column
!===============================================================================
function matrix_get_col(self, i) result(col)
@ -266,67 +233,45 @@ contains
col = self % col(i)
if (self % petsc_active) col = col + 1
end function matrix_get_col
!===============================================================================
! MATRIX_SETUP_PETSC configures the row/col vectors and links to a PETSc object
! MATRIX_TRANSPOSE transposes a sparse matrix
!===============================================================================
#ifdef PETSC
subroutine matrix_setup_petsc(self)
function matrix_transpose(self) result(mat)
class(Matrix), intent(inout) :: self ! matrix instance
class(Matrix), intent(in) :: self ! matrix instance
type(Matrix) :: mat ! transposed matrix
! change indices to c notation
self % row = self % row - 1
self % col = self % col - 1
integer :: i, j ! loop indices for row/column
integer :: i_transpose ! loop index for row in transposed matrix
integer :: first, last ! indices in col() array
! Link to petsc
call MatCreateSeqAIJWithArrays(PETSC_COMM_WORLD, self % n, self % n, &
self % row, self % col, self % val, self % petsc_mat, petsc_err)
call mat % create(self % n, self % nnz)
do i_transpose = 1, mat % n
! Set up row in transposed matrix
call mat % new_row()
! Petsc is now active
self % petsc_active = .true.
ROW: do i = 1, self % n
! Get range of columns for row i
first = self % row(i)
last = self % row(i+1) - 1
end subroutine matrix_setup_petsc
#endif
COL: do j = first, last
if (self % col(j) == i_transpose) then
! If column in original matrix matches row in transposed, add value
call mat % add_value(i, self % val(j))
elseif (self % col(j) > i_transpose) then
exit COL
end if
end do COL
end do ROW
end do
!===============================================================================
! MATRIX_WRITE_PETSC_BINARY writes a PETSc matrix binary file
!===============================================================================
call mat % new_row()
#ifdef PETSC
subroutine matrix_write_petsc_binary(self, filename)
character(*), intent(in) :: filename ! file name to write to
class(Matrix), intent(in) :: self ! matrix instance
type(PetscViewer) :: viewer ! a petsc viewer instance
call PetscViewerBinaryOpen(PETSC_COMM_WORLD, trim(filename), &
FILE_MODE_WRITE, viewer, petsc_err)
call MatView(self % petsc_mat, viewer, petsc_err)
call PetscViewerDestroy(viewer, petsc_err)
end subroutine matrix_write_petsc_binary
#endif
!===============================================================================
! MATRIX_TRANSPOSE uses PETSc to transpose a matrix
!===============================================================================
#ifdef PETSC
subroutine matrix_transpose(self)
class(Matrix), intent(inout) :: self ! matrix instance
call MatTranspose(self % petsc_mat, MAT_REUSE_MATRIX, self % petsc_mat, &
petsc_err)
end subroutine matrix_transpose
#endif
end function matrix_transpose
!===============================================================================
! MATRIX_VECTOR_MULTIPLY allow a vector to multiply the matrix
@ -363,7 +308,7 @@ contains
end subroutine matrix_vector_multiply
!===============================================================================
! MATRIX_SEARCH_INDICES searches for an index in column corresponding to a row
! MATRIX_SEARCH_INDICES searches for an index in column corresponding to a row
!===============================================================================
subroutine matrix_search_indices(self, row, col, idx, found)
@ -435,9 +380,6 @@ contains
if (.not.allocated(self % col)) allocate(self % col(self % nnz))
if (.not.allocated(self % val)) allocate(self % val(self % nnz))
! Set PETSc active to false
self % petsc_active = .false.
! Copy over data
self % row = mattocopy % row
self % col = mattocopy % col

View file

@ -26,12 +26,6 @@ element cmfd {
element dhat_reset { xsd:boolean }? &
element solver { xsd:string }? &
element snes_monitor { xsd:boolean }? &
element ksp_monitor { xsd:boolean }? &
element power_monitor { xsd:boolean }? &
element write_matrices { xsd:boolean }? &

View file

@ -143,21 +143,6 @@
<data type="boolean"/>
</element>
</optional>
<optional>
<element name="solver">
<data type="string"/>
</element>
</optional>
<optional>
<element name="snes_monitor">
<data type="boolean"/>
</element>
</optional>
<optional>
<element name="ksp_monitor">
<data type="boolean"/>
</element>
</optional>
<optional>
<element name="power_monitor">
<data type="boolean"/>

View file

@ -1,317 +0,0 @@
module solver_interface
use error, only: fatal_error
use matrix_header, only: Matrix
use vector_header, only: Vector
#ifdef PETSC
use petscksp
use petscsnes
# include <petscversion.h>
#endif
implicit none
private
! GMRES solver type
type, public :: GMRESSolver
#ifdef PETSC
type(ksp) :: ksp_ ! Krylov linear solver instance
type(pc) :: pc_ ! Preconditioner instance
#endif
contains
#ifdef PETSC
procedure :: create => petsc_gmres_create
procedure :: set_oper => petsc_gmres_set_oper
procedure :: destroy => petsc_gmres_destroy
procedure :: solve => petsc_gmres_solve
#endif
end type GMRESSolver
! Derived type to contain list of data needed to be passed to procedures
type, public :: Jfnk_ctx
procedure (res_interface), pointer, nopass :: res_proc_ptr
procedure (jac_interface), pointer, nopass :: jac_proc_ptr
end type Jfnk_ctx
! JFNK solver type
type, public :: JFNKSolver
#ifdef PETSC
type(ksp) :: ksp_ ! Krylov linear solver instance
type(pc) :: pc_ ! Preconditioner instance
type(snes) :: snes_ ! Nonlinear solver instance
type(mat) :: jac_mf ! Matrix free jacobian instance
integer :: ls ! Line search instance
#endif
contains
#ifdef PETSC
procedure :: create => petsc_jfnk_create
procedure :: destroy => petsc_jfnk_destroy
procedure :: set_functions => petsc_jfnk_set_functions
procedure :: solve => petsc_jfnk_solve
#endif
end type JFNKSolver
! Abstract interface stating how jacobian and residual routines look
abstract interface
subroutine res_interface(x, res)
import :: Vector
type(Vector), intent(in) :: x ! solution vector
type(Vector), intent(inout) :: res ! residual vector
end subroutine res_interface
subroutine jac_interface(x)
import :: Vector
type(Vector), intent(in) :: x ! solution vector
end subroutine jac_interface
end interface
#ifdef PETSC
integer :: petsc_err ! petsc error code
! Checks for PETSc version to handle 3.5 changes
# if (PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR > 4)
# define PETSC_GREATER_34
# endif
#endif
contains
#ifdef PETSC
!===============================================================================
! PETSC_GMRES_CREATE sets up a PETSc GMRES solver
!===============================================================================
subroutine petsc_gmres_create(self)
class(GMRESSolver), intent(inout) :: self ! GMRES solver instance
integer :: ilu_levels = 5
real(8) :: rtol = 1.0e-10_8
real(8) :: atol = 1.0e-10_8
call KSPCreate(PETSC_COMM_WORLD, self % ksp_, petsc_err)
#ifdef PETSC_GREATER_34
call KSPSetTolerances(self % ksp_, rtol, atol, &
PETSC_DEFAULT_REAL, PETSC_DEFAULT_INTEGER, petsc_err)
#else
call KSPSetTolerances(self % ksp_, rtol, atol, &
PETSC_DEFAULT_DOUBLE_PRECISION, PETSC_DEFAULT_INTEGER, petsc_err)
#endif
call KSPSetType(self % ksp_, 'gmres', petsc_err)
call KSPSetInitialGuessNonzero(self % ksp_, PETSC_TRUE, petsc_err)
call KSPGetPC(self % ksp_, self % pc_, petsc_err)
call PCFactorSetLevels(self % pc_, ilu_levels, petsc_err)
call KSPSetFromOptions(self % ksp_, petsc_err)
end subroutine petsc_gmres_create
!===============================================================================
! PETSC_GMRES_SET_OPER sets the matrix opetors for the GMRES solver
!===============================================================================
subroutine petsc_gmres_set_oper(self, prec_mat, mat_in)
class(GMRESSolver), intent(inout) :: self ! GMRES solver instanace
type(Matrix), intent(inout) :: prec_mat ! preconditioner matrix
type(Matrix), intent(inout) :: mat_in ! coefficient matrix
#ifdef PETSC_GREATER_34
call KSPSetOperators(self % ksp_, mat_in % petsc_mat, prec_mat % petsc_mat, &
petsc_err)
#else
call KSPSetOperators(self % ksp_, mat_in % petsc_mat, prec_mat % petsc_mat, &
SAME_NONZERO_PATTERN, petsc_err)
#endif
call KSPSetUp(self % ksp_, petsc_err)
end subroutine petsc_gmres_set_oper
!===============================================================================
! PETSC_GMRES_DESTROY frees all memory associated with the GMRES solver
!===============================================================================
subroutine petsc_gmres_destroy(self)
class(GMRESSolver), intent(inout) :: self ! GMRES solver instance
call KSPDestroy(self % ksp_, petsc_err)
end subroutine petsc_gmres_destroy
!===============================================================================
! PETSC_GMRES_SOLVE solves the linear system
!===============================================================================
subroutine petsc_gmres_solve(self, b, x)
class(GMRESSolver), intent(inout) :: self ! GMRES solver instance
type(Vector), intent(inout) :: b ! right hand side vector
type(Vector), intent(inout) :: x ! solution vector
call KSPSolve(self % ksp_, b % petsc_vec, x % petsc_vec, petsc_err)
end subroutine petsc_gmres_solve
!===============================================================================
! PETSC_JFNK_CREATE sets up a JFNK solver using PETSc SNES
!===============================================================================
subroutine petsc_jfnk_create(self)
class(JFNKSolver), intent(inout) :: self ! JFNK solver instance
! Turn on mf_operator option for matrix free jacobian
call PetscOptionsSetValue("-snes_mf_operator", "TRUE", petsc_err)
! Create the SNES context
call SNESCreate(PETSC_COMM_WORLD, self % snes_, petsc_err)
! Set up the GMRES solver
call SNESGetKSP(self % snes_, self % ksp_, petsc_err)
call KSPSetType(self % ksp_, 'gmres', petsc_err)
! Apply options
call SNESGetLineSearch(self % snes_, self % ls, petsc_err)
call SNESLineSearchSetType(self % ls, 'basic', petsc_err)
call SNESSetFromOptions(self % snes_, petsc_err)
! Set up preconditioner
call KSPGetPC(self % ksp_, self % pc_, petsc_err)
call PCSetType(self % pc_, 'ilu', petsc_err)
call PCFactorSetLevels(self % pc_, 5, petsc_err)
call PCSetFromOptions(self % pc_, petsc_err)
call KSPSetFromOptions(self % ksp_, petsc_err)
end subroutine petsc_jfnk_create
!===============================================================================
! PETSC_JFNK_DESTROY frees all memory associated with JFNK solver
!===============================================================================
subroutine petsc_jfnk_destroy(self)
class(JFNKSolver), intent(inout) :: self ! JFNK solver instance
call MatDestroy(self % jac_mf, petsc_err)
call SNESDestroy(self % snes_, petsc_err)
end subroutine petsc_jfnk_destroy
!===============================================================================
! PETSC_JFNK_SET_FUNCTIONS sets user functions and matrix free objects
!===============================================================================
subroutine petsc_jfnk_set_functions(self, ctx, res, jac_prec)
class(JFNKSolver), intent(inout) :: self ! JFNK solver instance
type(Jfnk_ctx), intent(inout) :: ctx ! JFNK context instance
type(Vector), intent(inout) :: res ! residual vector
type(Matrix), intent(inout) :: jac_prec ! preconditioner matrix
! Set residual procedure
call SNESSetFunction(self % snes_, res % petsc_vec, &
petsc_jfnk_compute_residual, ctx, petsc_err)
! Create the matrix free jacobian
call MatCreateSNESMF(self % snes_, self % jac_mf, petsc_err)
! Set Jacobian procedure
call SNESSetJacobian(self % snes_, self % jac_mf, jac_prec % petsc_mat, &
petsc_jfnk_compute_jacobian, ctx, petsc_err)
! Set up Jacobian Lags
call SNESSetLagJacobian(self % snes_, -2, petsc_err)
call SNESSetLagPreconditioner(self % snes_, -1, petsc_err)
end subroutine petsc_jfnk_set_functions
!===============================================================================
! PETSC_JFNK_SOLVE solves the nonlinear system
!===============================================================================
subroutine petsc_jfnk_solve(self, xvec)
class(JFNKSolver), intent(inout) :: self ! JFNK instance
type(Vector), intent(inout) :: xvec ! solution vector
call SNESSolve(self % snes_, PETSC_NULL_DOUBLE, xvec % petsc_vec, petsc_err)
end subroutine petsc_jfnk_solve
!===============================================================================
! PETSC_JFNK_COMPUTE_RESIDUAL buffer routine to user specifed residual routine
!===============================================================================
subroutine petsc_jfnk_compute_residual(snes_, x, res, ctx, ierr)
type(snes), intent(inout) :: snes_ ! PETSc SNES object
type(vec), intent(inout) :: x ! PETSc solution vector
type(vec), intent(inout) :: res ! PETSc residual vector
integer, intent(inout) :: ierr ! error code
type(Jfnk_ctx), intent(inout) :: ctx ! JFNK context instance
type(Vector) :: xvec ! solution vector
type(Vector) :: resvec ! residual vector
! We need to use the x and res that come from PETSc because the pointer
! location changes and therefore we can not use module variables
! for the residual and x vectors here
! Need to point an OpenMC vector to PETSc vector
call VecGetArrayF90(x, xvec % val, ierr)
call VecGetArrayF90(res, resvec % val, ierr)
! Call user residual routine
call ctx % res_proc_ptr(xvec, resvec)
! Need to restore the PETSc vector
call VecRestoreArrayF90(x, xvec % val, ierr)
call VecRestoreArrayF90(res, resvec % val, ierr)
end subroutine petsc_jfnk_compute_residual
!===============================================================================
! PETSC_JFNK_COMPUTE_JACOBIAN buffer routine to user specified jacobian routine
!===============================================================================
#ifdef PETSC_GREATER_34
subroutine petsc_jfnk_compute_jacobian(snes_, x, jac_mf, jac_prec, &
ctx, ierr)
#else
subroutine petsc_jfnk_compute_jacobian(snes_, x, jac_mf, jac_prec, flag, &
ctx, ierr)
#endif
type(snes), intent(inout) :: snes_ ! PETSc snes instance
type(vec), intent(inout) :: x ! PETSc solution vector
type(mat), intent(inout) :: jac_mf ! PETSc matrix free jacobian
type(mat), intent(inout) :: jac_prec ! PETSc matrix jacobian precond.
type(Jfnk_ctx), intent(inout) :: ctx ! JFNK context instance
integer, intent(inout) :: ierr ! error code
#ifndef PETSC_GREATER_34
integer, intent(inout) :: flag ! unused madatory flag
#endif
type(Vector) :: xvec ! solution vector
! Again, we use the vector that comes from Petsc to build the Jacobian
! matrix
! Need to point OpenMC vector to PETSc Vector
call VecGetArrayF90(x, xvec % val, ierr)
! Evaluate user jacobian routine
call ctx % jac_proc_ptr(xvec)
! Restore the PETSc vector
call VecRestoreArrayF90(x, xvec % val, ierr)
call MatAssemblyBegin(jac_mf, MAT_FINAL_ASSEMBLY, petsc_err)
call MatAssemblyEnd(jac_mf, MAT_FINAL_ASSEMBLY, petsc_err)
end subroutine petsc_jfnk_compute_jacobian
#endif
end module solver_interface

View file

@ -306,14 +306,11 @@ class CMFDFile(object):
self._feedback = None
self._inactive = None
self._inactive_flush = None
self._ksp_monitor = None
self._cmfd_mesh = None
self._norm = None
self._num_flushes = None
self._power_monitor = None
self._run_adjoint = None
self._snes_monitor = None
self._solver = None
self._write_matrices = None
self._cmfd_file = ET.Element("cmfd")
@ -350,11 +347,6 @@ class CMFDFile(object):
return self._inactive_flush
@property
def ksp_monitor(self):
return self._ksp_monitor
@property
def cmfd_mesh(self):
return self._cmfd_mesh
@ -380,11 +372,6 @@ class CMFDFile(object):
return self._run_adjoint
@property
def snes_monitor(self):
return self._snes_monitor
@property
def solver(self):
return self._solver
@ -481,17 +468,6 @@ class CMFDFile(object):
self._inactive_flush = inactive_flush
@ksp_monitor.setter
def ksp_monitor(self, ksp_monitor):
if not isinstance(ksp_monitor, bool):
msg = 'Unable to set CMFD ksp monitor to {0} which is a ' \
'non-boolean value'.format(ksp_monitor)
raise ValueError(msg)
self._ksp_monitor = ksp_monitor
@cmfd_mesh.setter
def cmfd_mesh(self, mesh):
@ -552,28 +528,6 @@ class CMFDFile(object):
self._run_adjoint = run_adjoint
@snes_monitor.setter
def snes_monitor(self, snes_monitor):
if not isinstance(snes_monitor, bool):
msg = 'Unable to set CMFD snes monitor to {0} which is a ' \
'non-boolean value'.format(snes_monitor)
raise ValueError(msg)
self._snes_monitor = snes_monitor
@solver.setter
def solver(self, solver):
if not solver in ['power', 'jfnk']:
msg = 'Unable to set CMFD solver to {0} which is not ' \
'"power" or "jfnk"'.format(solver)
raise ValueError(msg)
self._solver = solver
@write_matrices.setter
def write_matrices(self, write_matrices):
@ -627,13 +581,6 @@ class CMFDFile(object):
element.text = '{0}'.format(str(self._inactive_flush))
def create_ksp_monitor_subelement(self):
if not self._ksp_monitor is None:
element = ET.SubElement(self._cmfd_file, "ksp_monitor")
element.text = '{0}'.format(str(self._ksp_monitor).lower())
def create_mesh_subelement(self):
if not self._mesh is None:
@ -669,20 +616,6 @@ class CMFDFile(object):
element.text = '{0}'.format(str(self._run_adjoint).lower())
def create_snes_monitor_subelement(self):
if not self._snes_monitor is None:
element = ET.SubElement(self._cmfd_file, "snes_monitor")
element.text = '{0}'.format(str(self._snes_monitor).lower())
def create_solver_subelement(self):
if not self._solver is None:
element = ET.SubElement(self._cmfd_file, "solver")
element.text = '{0}'.format(str(self._solver))
def create_write_matrices_subelement(self):
if not self._write_matrices is None:
@ -698,14 +631,11 @@ class CMFDFile(object):
self.create_feedback_subelement()
self.create_inactive_subelement()
self.create_inactive_flush_subelement()
self.create_ksp_monitor_subelement()
self.create_mesh_subelement()
self.create_norm_subelement()
self.create_num_flushes_subelement()
self.create_power_monitor_subelement()
self.create_run_adjoint_subelement()
self.create_snes_monitor_subelement()
self.create_solver_subelement()
self.create_write_matrices_subelement()
# Clean the indentation in the file to be user-readable

View file

@ -2,40 +2,25 @@ module vector_header
use constants, only: ZERO
#ifdef PETSC
use petscvec
#endif
implicit none
private
type, public :: Vector
type, public :: Vector
integer :: n ! number of rows/cols in matrix
real(8), allocatable :: data(:) ! where vector data is stored
real(8), pointer :: val(:) ! pointer to vector data
# ifdef PETSC
type(vec) :: petsc_vec ! PETSc vector
# endif
logical :: petsc_active ! Logical if PETSc is being used
contains
procedure :: create => vector_create
procedure :: destroy => vector_destroy
procedure :: add_value => vector_add_value
procedure :: copy => vector_copy
# ifdef PETSC
procedure :: setup_petsc => vector_setup_petsc
procedure :: write_petsc_binary => vector_write_petsc_binary
# endif
! TODO: procedure :: write => vector_write
end type Vector
#ifdef PETSC
integer :: petsc_err ! petsc error code
#endif
contains
!===============================================================================
! VECTOR_CREATE allocates and initializes a vector
! VECTOR_CREATE allocates and initializes a vector
!===============================================================================
subroutine vector_create(self, n)
@ -53,9 +38,6 @@ contains
! Initialize to zero
self % val = ZERO
! Petsc is default not active
self % petsc_active = .false.
end subroutine vector_create
!===============================================================================
@ -66,10 +48,6 @@ contains
class(Vector), intent(inout) :: self ! vector instance
#ifdef PETSC
if (self % petsc_active) call VecDestroy(self % petsc_vec, petsc_err)
#endif
if (associated(self % val)) nullify(self % val)
if (allocated(self % data)) deallocate(self % data)
@ -89,45 +67,6 @@ contains
end subroutine vector_add_value
!===============================================================================
! VECTOR_SETUP_PETSC links the data to a PETSc vector
!===============================================================================
#ifdef PETSC
subroutine vector_setup_petsc(self)
class(Vector), intent(inout) :: self ! vector instance
! Link to PETSc
call VecCreateSeqWithArray(PETSC_COMM_WORLD, 1, self % n, self % val, &
self % petsc_vec, petsc_err)
! Set that PETSc is now active
self % petsc_active = .true.
end subroutine vector_setup_petsc
#endif
!===============================================================================
! VECTOR_WRITE_PETSC_BINARY writes the PETSc vector to a binary file
!===============================================================================
#ifdef PETSC
subroutine vector_write_petsc_binary(self, filename)
character(*), intent(in) :: filename ! name of file to write to
class(Vector), intent(in) :: self ! vector instance
type(PetscViewer) :: viewer ! PETSc viewer instance
call PetscViewerBinaryOpen(PETSC_COMM_WORLD, trim(filename), &
FILE_MODE_WRITE, viewer, petsc_err)
call VecView(self % petsc_vec, viewer, petsc_err)
call PetscViewerDestroy(viewer, petsc_err)
end subroutine vector_write_petsc_binary
#endif
!===============================================================================
! VECTOR_COPY allocates a separate vector and copies
!===============================================================================
@ -144,12 +83,9 @@ contains
! Set n
self % n = vectocopy % n
! Copy values
! Copy values
self % val = vectocopy % val
! Petsc is default not active
self % petsc_active = .false.
end subroutine vector_copy
end module vector_header

View file

@ -45,7 +45,6 @@ FC='gfortran'
MPI_DIR='/opt/mpich/3.1.3-gnu'
HDF5_DIR='/opt/hdf5/1.8.14-gnu'
PHDF5_DIR='/opt/phdf5/1.8.14-gnu'
PETSC_DIR='/opt/petsc/3.5.2-gnu'
# Script mode for extra capability
script_mode = False
@ -59,8 +58,6 @@ if 'HDF5_DIR' in os.environ:
HDF5_DIR = os.environ['HDF5_DIR']
if 'PHDF5_DIR' in os.environ:
PHDF5_DIR = os.environ['PHDF5_DIR']
if 'PETSC_DIR' in os.environ:
PETSC_DIR = os.environ['PETSC_DIR']
# CTest script template
ctest_str = """set (CTEST_SOURCE_DIRECTORY "{source_dir}")
@ -110,14 +107,13 @@ tests = OrderedDict()
class Test(object):
def __init__(self, name, debug=False, optimize=False, mpi=False, openmp=False,
hdf5=False, petsc=False, valgrind=False, coverage=False):
hdf5=False, valgrind=False, coverage=False):
self.name = name
self.debug = debug
self.optimize = optimize
self.mpi = mpi
self.openmp = openmp
self.hdf5 = hdf5
self.petsc = petsc
self.valgrind = valgrind
self.coverage = coverage
self.success = True
@ -152,8 +148,6 @@ class Test(object):
build_str += "-Doptimize=ON "
if self.openmp:
build_str += "-Dopenmp=ON "
if self.petsc:
build_str += "-Dpetsc=ON "
if self.coverage:
build_str += "-Dcoverage=ON "
self.build_opts = build_str
@ -167,8 +161,6 @@ class Test(object):
# Runs the ctest script which performs all the cmake/ctest/cdash
def run_ctest_script(self):
os.environ['FC'] = self.fc
if self.petsc:
os.environ['PETSC_DIR'] = PETSC_DIR
if self.mpi:
os.environ['MPI_DIR'] = MPI_DIR
rc = call(['ctest', '-S', 'ctestscript.run','-V'])
@ -179,8 +171,6 @@ class Test(object):
# Runs cmake when in non-script mode
def run_cmake(self):
os.environ['FC'] = self.fc
if self.petsc:
os.environ['PETSC_DIR'] = PETSC_DIR
if self.mpi:
os.environ['MPI_DIR'] = MPI_DIR
build_opts = self.build_opts.split()
@ -272,9 +262,9 @@ class Test(object):
# Simple function to add a test to the global tests dictionary
def add_test(name, debug=False, optimize=False, mpi=False, openmp=False,\
hdf5=False, petsc=False, valgrind=False, coverage=False):
tests.update({name:Test(name, debug, optimize, mpi, openmp, hdf5, petsc,
valgrind, coverage)})
hdf5=False, valgrind=False, coverage=False):
tests.update({name: Test(name, debug, optimize, mpi, openmp, hdf5,
valgrind, coverage)})
# List of all tests that may be run. User can add -C to command line to specify
# a subset of these configurations
@ -302,22 +292,11 @@ add_test('phdf5-optimize', mpi=True, hdf5=True, optimize=True)
add_test('phdf5-omp-normal', mpi=True, hdf5=True, openmp=True)
add_test('phdf5-omp-debug', mpi=True, hdf5=True, openmp=True, debug=True)
add_test('phdf5-omp-optimize', mpi=True, hdf5=True, openmp=True, optimize=True)
add_test('petsc-normal', petsc=True, mpi=True)
add_test('petsc-debug', petsc=True, mpi=True, debug=True)
add_test('petsc-optimize', petsc=True, mpi=True, optimize=True)
add_test('phdf5-petsc-normal', mpi=True, hdf5=True, petsc=True)
add_test('phdf5-petsc-debug', mpi=True, hdf5=True, petsc=True, debug=True)
add_test('phdf5-petsc-optimize', mpi=True, hdf5=True, petsc=True, optimize=True)
add_test('omp-phdf5-petsc-normal', openmp=True, mpi=True, hdf5=True, petsc=True)
add_test('omp-phdf5-petsc-debug', openmp=True, mpi=True, hdf5=True, petsc=True, debug=True)
add_test('omp-phdf5-petsc-optimize', openmp=True, mpi=True, hdf5=True, petsc=True, optimize=True)
add_test('basic-debug_valgrind', debug=True, valgrind=True)
add_test('hdf5-debug_valgrind', hdf5=True, debug=True, valgrind=True)
add_test('petsc-debug_valgrind', petsc=True, mpi=True, debug=True, valgrind=True)
add_test('basic-debug_coverage', debug=True, coverage=True)
add_test('hdf5-debug_coverage', debug=True, hdf5=True, coverage=True)
add_test('mpi-debug_coverage', debug=True, mpi=True, coverage=True)
add_test('petsc-debug_coverage', debug=True, petsc=True, mpi=True, coverage=True)
# Check to see if we should just print build configuration information to user
if options.list_build_configs:
@ -328,7 +307,6 @@ if options.list_build_configs:
print(' HDF5 Active:..........{0}'.format(tests[key].hdf5))
print(' MPI Active:...........{0}'.format(tests[key].mpi))
print(' OpenMP Active:........{0}'.format(tests[key].openmp))
print(' PETSc Active:.........{0}'.format(tests[key].petsc))
print(' Valgrind Test:........{0}'.format(tests[key].valgrind))
print(' Coverage Test:........{0}\n'.format(tests[key].coverage))
exit()

View file

@ -10,7 +10,6 @@
<begin> 5 </begin>
<display> dominance </display>
<solver> power </solver>
<feedback> true </feedback>
<gauss_seidel_tolerance> 1.e-15 1.e-20 </gauss_seidel_tolerance>
</cmfd>

View file

@ -1,16 +0,0 @@
<?xml version="1.0"?>
<cmfd>
<mesh>
<lower_left> -10 -1 -1 </lower_left>
<upper_right> 10 1 1 </upper_right>
<dimension> 10 1 1 </dimension>
<albedo> 0.0 0.0 1.0 1.0 1.0 1.0 </albedo>
</mesh>
<begin> 5 </begin>
<display> balance </display>
<solver> jfnk </solver>
<feedback> true </feedback>
</cmfd>

View file

@ -1,43 +0,0 @@
<?xml version="1.0" encoding="UTF-8"?>
<geometry>
<!-- Definition of Cells -->
<cell id="1">
<universe>0</universe>
<surfaces>-1 2 -3 4 -5 6</surfaces>
<material>1</material>
</cell>
<!-- Defition of Surfaces -->
<surface id="1">
<type>x-plane</type>
<coeffs>10</coeffs>
<boundary> vacuum </boundary>
</surface>
<surface id="2">
<type>x-plane</type>
<coeffs>-10</coeffs>
<boundary> vacuum </boundary>
</surface>
<surface id="3">
<type>y-plane</type>
<coeffs>1</coeffs>
<boundary>reflective</boundary>
</surface>
<surface id="4">
<type>y-plane</type>
<coeffs>-1</coeffs>
<boundary>reflective</boundary>
</surface>
<surface id="5">
<type>z-plane</type>
<coeffs>1</coeffs>
<boundary>reflective</boundary>
</surface>
<surface id="6">
<type>z-plane</type>
<coeffs>-1</coeffs>
<boundary>reflective</boundary>
</surface>
</geometry>

View file

@ -1,12 +0,0 @@
<?xml version="1.0"?>
<materials>
<!-- Definition of materials -->
<material id="1">
<density value="19" units="g/cc" />
<nuclide name="U-235" xs="71c" wo="0.21" />
<nuclide name="U-238" xs="71c" wo="0.68" />
<nuclide name="O-16" xs="71c" wo="0.11" />
</material>
</materials>

View file

@ -1,94 +0,0 @@
#!/usr/bin/env python
import sys
import numpy as np
sys.path.insert(0, '../../src/utils')
# import statepoint
import openmc
from openmc.statepoint import StatePoint
# read in statepoint file
if len(sys.argv) > 1:
sp = StatePoint(sys.argv[1])
else:
sp = StatePoint('statepoint.20.binary')
sp.read_results()
# extract tally results and convert to vector
tally1 = sp.get_tally('flux', [openmc.Filter('mesh', [1])], \
[-1], estimator='tracklength')
tally2 = sp.get_tally('flux', [openmc.Filter('mesh', [2])], \
[-1], estimator='analog')
tally3 = sp.get_tally('nu-fission', [openmc.Filter('mesh', [2])], \
[-1], estimator='analog')
tally4 = sp.get_tally('current', [openmc.Filter('mesh', [2]), \
openmc.Filter('surface', [1,2,3,4,5,6])], \
[-1], estimator='analog')
results1 = np.zeros((tally1.sum.size + tally1.sum.size, ))
results1[0::2] = tally1.sum.ravel()
results1[1::2] = tally1.sum_sq.ravel()
results2 = np.zeros((tally2.sum.size + tally2.sum.size, ))
results2[0::2] = tally2.sum.ravel()
results2[1::2] = tally2.sum_sq.ravel()
results3 = np.zeros((tally3.sum.size + tally3.sum.size, ))
results3[0::2] = tally3.sum.ravel()
results3[1::2] = tally3.sum_sq.ravel()
results4 = np.zeros((tally4.sum.size + tally4.sum.size, ))
results4[0::2] = tally4.sum.ravel()
results4[1::2] = tally4.sum_sq.ravel()
# set up output string
outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp._k_combined[0], sp._k_combined[1])
# write out tally results
outstr += 'tally 1:\n'
for item in results1:
outstr += "{0:12.6E}\n".format(item)
outstr += 'tally 2:\n'
for item in results2:
outstr += "{0:12.6E}\n".format(item)
outstr += 'tally 3:\n'
for item in results3:
outstr += "{0:12.6E}\n".format(item)
outstr += 'tally 4:\n'
for item in results4:
outstr += "{0:12.6E}\n".format(item)
# write out cmfd answers
outstr += 'cmfd indices\n'
for item in sp._cmfd_indices:
outstr += "{0:12.6E}\n".format(item)
outstr += 'k cmfd\n'
for item in sp._k_cmfd:
outstr += "{0:12.6E}\n".format(item)
outstr += 'cmfd entropy\n'
for item in sp._cmfd_entropy:
outstr += "{0:12.6E}\n".format(item)
outstr += 'cmfd balance\n'
for item in sp._cmfd_balance:
outstr += "{0:12.6E}\n".format(item)
outstr += 'cmfd dominance ratio\n'
for item in sp._cmfd_dominance:
outstr += "{0:10.3E}\n".format(item)
outstr += 'cmfd openmc source comparison\n'
for item in sp._cmfd_srccmp:
outstr += "{0:12.6E}\n".format(item)
outstr += 'cmfd source\n'
cmfdsrc = np.reshape(sp._cmfd_src, np.product(sp._cmfd_indices), order='F')
for item in cmfdsrc:
outstr += "{0:12.6E}\n".format(item)
# write results to file
with open('results_test.dat', 'w') as fh:
fh.write(outstr)

View file

@ -1,775 +0,0 @@
k-combined:
1.177396E+00 4.883437E-03
tally 1:
1.107335E+01
1.235221E+01
2.002676E+01
4.017045E+01
2.844184E+01
8.111731E+01
3.428492E+01
1.177848E+02
3.735839E+01
1.397523E+02
3.894180E+01
1.517824E+02
3.528006E+01
1.246309E+02
2.863448E+01
8.222321E+01
2.125384E+01
4.535192E+01
1.112124E+01
1.241089E+01
tally 2:
2.243113E+01
2.535057E+01
1.557550E+01
1.222813E+01
2.164704E+00
2.389690E-01
4.049814E+01
8.218116E+01
2.854669E+01
4.085264E+01
3.934452E+00
7.833208E-01
5.706024E+01
1.633739E+02
4.049737E+01
8.234850E+01
5.231368E+00
1.386110E+00
6.798682E+01
2.320807E+02
4.819178E+01
1.166820E+02
6.247056E+00
1.968990E+00
7.449317E+01
2.782374E+02
5.315156E+01
1.417287E+02
6.667584E+00
2.250879E+00
7.530107E+01
2.849330E+02
5.378916E+01
1.454619E+02
7.071012E+00
2.522919E+00
6.820303E+01
2.335083E+02
4.850330E+01
1.181013E+02
6.241747E+00
1.973368E+00
5.831373E+01
1.704779E+02
4.149124E+01
8.633656E+01
5.385636E+00
1.468642E+00
4.184350E+01
8.792430E+01
2.971699E+01
4.435896E+01
3.784806E+00
7.343249E-01
2.202673E+01
2.444732E+01
1.531988E+01
1.183160E+01
2.073873E+00
2.272561E-01
tally 3:
1.500439E+01
1.135466E+01
9.942586E-01
5.051831E-02
2.744453E+01
3.776958E+01
1.781992E+00
1.612021E-01
3.901013E+01
7.642401E+01
2.472924E+00
3.077607E-01
4.642527E+01
1.082730E+02
2.924369E+00
4.335819E-01
5.109083E+01
1.309660E+02
3.325865E+00
5.592642E-01
5.182512E+01
1.350762E+02
3.259912E+00
5.350517E-01
4.672480E+01
1.095974E+02
3.097309E+00
4.854729E-01
3.997405E+01
8.014326E+01
2.589176E+00
3.374615E-01
2.861624E+01
4.114210E+01
1.806237E+00
1.656944E-01
1.474531E+01
1.096355E+01
9.807860E-01
4.934364E-02
tally 4:
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
3.050887E+00
4.698799E-01
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
5.388951E+00
1.456464E+00
2.664528E+00
3.577345E-01
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
7.165280E+00
2.573230E+00
4.930263E+00
1.218988E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
8.550278E+00
3.668687E+00
6.985934E+00
2.450395E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
9.236725E+00
4.281659E+00
8.429932E+00
3.565264E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
9.396613E+00
4.431004E+00
9.307625E+00
4.351276E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
8.721885E+00
3.821463E+00
9.438995E+00
4.479948E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
7.096946E+00
2.525113E+00
8.605114E+00
3.710134E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
5.222832E+00
1.373166E+00
7.480684E+00
2.804650E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
2.726690E+00
3.773397E-01
5.470375E+00
1.504111E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
3.026290E+00
4.597502E-01
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
cmfd indices
1.000000E+01
1.000000E+00
1.000000E+00
1.000000E+00
k cmfd
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
1.150583E+00
1.159578E+00
1.162798E+00
1.171003E+00
1.162732E+00
1.165591E+00
1.166032E+00
1.165387E+00
1.165451E+00
1.170726E+00
1.170820E+00
1.167945E+00
1.166050E+00
1.165757E+00
1.167631E+00
1.168366E+00
cmfd entropy
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
3.215201E+00
3.212850E+00
3.214094E+00
3.208054E+00
3.212891E+00
3.213349E+00
3.208290E+00
3.206895E+00
3.208786E+00
3.209630E+00
3.212321E+00
3.211750E+00
3.212453E+00
3.213370E+00
3.211791E+00
3.212685E+00
cmfd balance
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
6.348026E-03
5.102747E-03
4.043947E-03
3.952742E-03
2.533934E-03
2.215030E-03
2.945032E-03
2.597862E-03
2.300873E-03
1.991974E-03
1.679430E-03
1.636923E-03
1.581336E-03
1.431472E-03
1.419168E-03
1.256242E-03
cmfd dominance ratio
0.000E+00
0.000E+00
0.000E+00
0.000E+00
0.000E+00
0.000E+00
0.000E+00
0.000E+00
0.000E+00
0.000E+00
0.000E+00
0.000E+00
0.000E+00
0.000E+00
0.000E+00
0.000E+00
0.000E+00
0.000E+00
0.000E+00
0.000E+00
cmfd openmc source comparison
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
8.388554E-03
7.868057E-03
6.387210E-03
6.979018E-03
5.634468E-03
5.579467E-03
5.647780E-03
5.289856E-03
4.550548E-03
4.373716E-03
4.042349E-03
3.813510E-03
3.749150E-03
3.358125E-03
3.562360E-03
3.904943E-03
cmfd source
4.142294E-02
7.484381E-02
1.050784E-01
1.256771E-01
1.444717E-01
1.420230E-01
1.350735E-01
1.118744E-01
7.755343E-02
4.198177E-02

View file

@ -1,32 +0,0 @@
<?xml version="1.0" encoding="UTF-8"?>
<settings>
<!-- Parameters for criticality calculation -->
<eigenvalue>
<batches>20</batches>
<inactive>10</inactive>
<particles>1000</particles>
</eigenvalue>
<!-- How verbose output should be -->
<verbosity value="7" />
<!-- Starting source -->
<source>
<space>
<type>box</type>
<parameters>-10 -1 -1 10 1 1</parameters>
</space>
</source>
<!-- Shannon Entropy -->
<entropy>
<dimension> 10 1 1 </dimension>
<lower_left> -10.0 -1.0 -1.0 </lower_left>
<upper_right> 10.0 1.0 1.0 </upper_right>
</entropy>
<!-- Run CMFD -->
<run_cmfd> true </run_cmfd>
</settings>

View file

@ -1,16 +0,0 @@
<?xml version="1.0"?>
<tallies>
<mesh id="1">
<type>rectangular</type>
<lower_left>-10 -1 -1 </lower_left>
<upper_right>10 1 1</upper_right>
<dimension>10 1 1</dimension>
</mesh>
<tally id="1">
<filter type="mesh" bins="1" />
<scores>flux</scores>
</tally>
</tallies>

View file

@ -1,68 +0,0 @@
#!/usr/bin/env python
import os
from subprocess import Popen, STDOUT, PIPE, call
import filecmp
import glob
from optparse import OptionParser
parser = OptionParser()
parser.add_option('--mpi_exec', dest='mpi_exec', default='')
parser.add_option('--mpi_np', dest='mpi_np', default='3')
parser.add_option('--exe', dest='exe')
(opts, args) = parser.parse_args()
cwd = os.getcwd()
def test_run():
if opts.mpi_exec != '':
proc = Popen([opts.mpi_exec, '-np', opts.mpi_np, opts.exe, cwd],
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([opts.exe, cwd], stderr=STDOUT, stdout=PIPE)
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0, 'OpenMC did not exit successfully.'
def test_created_statepoint():
statepoint = glob.glob(os.path.join(cwd, 'statepoint.20.*'))
assert len(statepoint) == 1, 'Either multiple or no statepoint files exist.'
assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5'),\
'Statepoint file is not a binary or hdf5 file.'
def test_output_exists():
assert os.path.exists(os.path.join(cwd, 'tallies.out')), 'Tally output file does not exist.'
def test_results():
statepoint = glob.glob(os.path.join(cwd, 'statepoint.20.*'))
call(['python', 'results.py', statepoint[0]])
compare = filecmp.cmp('results_test.dat', 'results_true.dat')
if not compare:
os.rename('results_test.dat', 'results_error.dat')
assert compare, 'Results do not agree.'
def teardown():
output = glob.glob(os.path.join(cwd, 'statepoint.20.*'))
output.append(os.path.join(cwd, 'tallies.out'))
output.append(os.path.join(cwd, 'results_test.dat'))
for f in output:
if os.path.exists(f):
os.remove(f)
if __name__ == '__main__':
# test for openmc executable
if opts.exe is None:
raise Exception('Must specify OpenMC executable from command line with --exe.')
# run tests
try:
test_run()
test_created_statepoint()
test_output_exists()
test_results()
finally:
teardown()

View file

@ -10,7 +10,6 @@
<begin> 5 </begin>
<display> dominance </display>
<solver> power </solver>
<feedback> false </feedback>
<gauss_seidel_tolerance> 1.e-15 1.e-20 </gauss_seidel_tolerance>

View file

@ -4,7 +4,7 @@ set -ev
# Run all debug tests
if [ "$TRAVIS_PULL_REQUEST" != "false" ]; then
./run_tests.py -C "^basic-debug$|^hdf5-debug$|^mpi-omp-debug$|^phdf5-omp-debug$|^omp-phdf5-petsc-debug$" -j 4 -s
./run_tests.py -C "^basic-debug$|^hdf5-debug$|^mpi-omp-debug$|^phdf5-omp-debug$" -j 4 -s
else
./run_tests.py -C "^basic-debug$" -j 4
fi

View file

@ -2,7 +2,7 @@
set -ev
# Build HDF5 and PETSc for rest of debug tests
# Build MPICH and HDF5 for rest of debug tests
if [ "$TRAVIS_PULL_REQUEST" != "false" ]; then
# Build MPICH
@ -35,18 +35,4 @@ if [ "$TRAVIS_PULL_REQUEST" != "false" ]; then
make install >/dev/null 2>&1
cd ..
# Build PETSc
wget -q http://ftp.mcs.anl.gov/pub/petsc/release-snapshots/petsc-lite-3.5.3.tar.gz
tar -xzvf petsc-lite-3.5.3.tar.gz >/dev/null 2>&1
cd petsc-3.5.3
./configure --prefix=$PWD/../petsc_install -q --download-fblaslapack \
--with-mpi-dir=$PWD/../mpich_install --with-share-libraries \
--with-fortran-datatypes
make PETSC_DIR=$PWD \
PETSC_ARCH=arch-linux2-c-debug all >/dev/null 2>&1
make PETSC_DIR=$PWD \
PETSC_ARCH=arch-linux2-c-debug install >/dev/null 2>&1
make install >/dev/null 2>&1
cd ..
fi