Merge pull request #199 from mit-crpg/cmfd_clean

New CMFD Routines
This commit is contained in:
Paul Romano 2013-10-16 20:55:51 -07:00
commit 14b3db3cd7
111 changed files with 5083 additions and 3098 deletions

View file

@ -4,6 +4,544 @@
State Point Binary File Specifications
======================================
-----------
Revision 10
-----------
**integer(4) FILETYPE_STATEPOINT**
Flags whether this file is a statepoint file or a particle restart file.
**integer(4) REVISION_STATEPOINT**
Revision of the binary state point file. Any time a change is made in the
format of the state-point file, this integer is incremented.
**integer(4) VERSION_MAJOR**
Major version number for OpenMC
**integer(4) VERSION_MINOR**
Minor version number for OpenMC
**integer(4) VERSION_RELEASE**
Release version number for OpenMC
**character(19) time_stamp**
Date and time the state point was written.
**character(255) path**
Absolute path to directory containing input files.
**integer(8) seed**
Pseudo-random number generator seed.
**integer(4) run_mode**
run mode used. The modes are described in constants.F90.
**integer(8) n_particles**
Number of particles used per generation.
**integer(4) n_batches**
Total number of batches (active + inactive).
**integer(4) current_batch**
The number of batches already simulated.
if (run_mode == MODE_EIGENVALUE)
**integer(4) n_inactive**
Number of inactive batches
**integer(4) gen_per_batch**
Number of generations per batch for criticality calculations
*do i = 1, current_batch \* gen_per_batch*
**real(8) k_generation(i)**
k-effective for the i-th total generation
*do i = 1, current_batch \* gen_per_batch*
**real(8) entropy(i)**
Shannon entropy for the i-th total generation
**real(8) k_col_abs**
Sum of product of collision/absorption estimates of k-effective
**real(8) k_col_tra**
Sum of product of collision/track-length estimates of k-effective
**real(8) k_abs_tra**
Sum of product of absorption/track-length estimates of k-effective
**real(8) k_combined(2)**
Mean and standard deviation of a combined estimate of k-effective
**integer(4) cmfd_on**
Flag that cmfd is on
if (cmfd_on)
**integer(4) cmfd % indices**
Indices for cmfd mesh (i,j,k,g)
**real(8) cmfd % k_cmfd(1:current_batch)**
CMFD eigenvalues
**real(8) cmfd % src(1:I,1:J,1:K,1:G)**
CMFD fission source
**real(8) cmfd % entropy(1:current_batch)**
CMFD estimate of Shannon entropy
**real(8) cmfd % balance(1:current_batch)**
RMS of the residual neutron balance equation on CMFD mesh
**real(8) cmfd % dom(1:current_batch)**
CMFD estimate of dominance ratio
**real(8) cmfd % scr_cmp(1:current_batch)**
RMS comparison of difference between OpenMC and CMFD fission source
**integer(4) n_meshes**
Number of meshes in tallies.xml file
*do i = 1, n_meshes*
**integer(4) meshes(i) % id**
Unique ID of mesh.
**integer(4) meshes(i) % type**
Type of mesh.
**integer(4) meshes(i) % n_dimension**
Number of dimensions for mesh (2 or 3).
**integer(4) meshes(i) % dimension(:)**
Number of mesh cells in each dimension.
**real(8) meshes(i) % lower_left(:)**
Coordinates of lower-left corner of mesh.
**real(8) meshes(i) % upper_right(:)**
Coordinates of upper-right corner of mesh.
**real(8) meshes(i) % width(:)**
Width of each mesh cell in each dimension.
**integer(4) n_tallies**
*do i = 1, n_tallies*
**integer(4) tallies(i) % id**
Unique ID of tally.
**integer(4) tallies(i) % n_realizations**
Number of realizations for the i-th tally.
**integer(4) size(tallies(i) % scores, 1)**
Total number of score bins for the i-th tally
**integer(4) size(tallies(i) % scores, 2)**
Total number of filter bins for the i-th tally
**integer(4) tallies(i) % n_filters**
*do j = 1, tallies(i) % n_filters*
**integer(4) tallies(i) % filter(j) % type**
Type of tally filter.
**integer(4) tallies(i) % filter(j) % n_bins**
Number of bins for filter.
**integer(4)/real(8) tallies(i) % filter(j) % bins(:)**
Value for each filter bin of this type.
**integer(4) tallies(i) % n_nuclide_bins**
Number of nuclide bins. If none are specified, this is just one.
*do j = 1, tallies(i) % n_nuclide_bins*
**integer(4) tallies(i) % nuclide_bins(j)**
Values of specified nuclide bins
**integer(4) tallies(i) % n_score_bins**
Number of scoring bins.
*do j = 1, tallies(i) % n_score_bins*
**integer(4) tallies(i) % score_bins(j)**
Values of specified scoring bins (e.g. SCORE_FLUX).
*do j = 1, tallies(i) % n_score_bins*
**integer(4) tallies(i) % scatt_order(j)**
Scattering Order specified scoring bins.
**integer(4) tallies(i) % n_score_bins**
Number of scoring bins without accounting for those added by
the scatter-pn command.
**integer(4) n_realizations**
Number of realizations for global tallies.
**integer(4) N_GLOBAL_TALLIES**
Number of global tally scores
*do i = 1, N_GLOBAL_TALLIES*
**real(8) global_tallies(i) % sum**
Accumulated sum for the i-th global tally
**real(8) global_tallies(i) % sum_sq**
Accumulated sum of squares for the i-th global tally
**integer(4) tallies_on**
Flag indicated if tallies are present in the file.
if (tallies_on > 0)
*do i = 1, n_tallies*
*do k = 1, size(tallies(i) % scores, 2)*
*do j = 1, size(tallies(i) % scores, 1)*
**real(8) tallies(i) % scores(j,k) % sum**
Accumulated sum for the j-th score and k-th filter of the
i-th tally
**real(8) tallies(i) % scores(j,k) % sum_sq**
Accumulated sum of squares for the j-th score and k-th
filter of the i-th tally
if (run_mode == MODE_EIGENVALUE)
*do i = 1, n_particles*
**real(8) source_bank(i) % wgt**
Weight of the i-th source particle
**real(8) source_bank(i) % xyz(1:3)**
Coordinates of the i-th source particle.
**real(8) source_bank(i) % uvw(1:3)**
Direction of the i-th source particle
**real(8) source_bank(i) % E**
Energy of the i-th source particle.
----------
Revision 9
----------
**integer(4) FILETYPE_STATEPOINT**
Flags whether this file is a statepoint file or a particle restart file.
**integer(4) REVISION_STATEPOINT**
Revision of the binary state point file. Any time a change is made in the
format of the state-point file, this integer is incremented.
**integer(4) VERSION_MAJOR**
Major version number for OpenMC
**integer(4) VERSION_MINOR**
Minor version number for OpenMC
**integer(4) VERSION_RELEASE**
Release version number for OpenMC
**character(19) time_stamp**
Date and time the state point was written.
**character(255) path**
Absolute path to directory containing input files.
**integer(8) seed**
Pseudo-random number generator seed.
**integer(4) run_mode**
run mode used. The modes are described in constants.F90.
**integer(8) n_particles**
Number of particles used per generation.
**integer(4) n_batches**
Total number of batches (active + inactive).
**integer(4) current_batch**
The number of batches already simulated.
if (run_mode == MODE_EIGENVALUE)
**integer(4) n_inactive**
Number of inactive batches
**integer(4) gen_per_batch**
Number of generations per batch for criticality calculations
*do i = 1, current_batch \* gen_per_batch*
**real(8) k_generation(i)**
k-effective for the i-th total generation
*do i = 1, current_batch \* gen_per_batch*
**real(8) entropy(i)**
Shannon entropy for the i-th total generation
**real(8) k_col_abs**
Sum of product of collision/absorption estimates of k-effective
**real(8) k_col_tra**
Sum of product of collision/track-length estimates of k-effective
**real(8) k_abs_tra**
Sum of product of absorption/track-length estimates of k-effective
**real(8) k_combined(2)**
Mean and standard deviation of a combined estimate of k-effective
**integer(4) n_meshes**
Number of meshes in tallies.xml file
*do i = 1, n_meshes*
**integer(4) meshes(i) % id**
Unique ID of mesh.
**integer(4) meshes(i) % type**
Type of mesh.
**integer(4) meshes(i) % n_dimension**
Number of dimensions for mesh (2 or 3).
**integer(4) meshes(i) % dimension(:)**
Number of mesh cells in each dimension.
**real(8) meshes(i) % lower_left(:)**
Coordinates of lower-left corner of mesh.
**real(8) meshes(i) % upper_right(:)**
Coordinates of upper-right corner of mesh.
**real(8) meshes(i) % width(:)**
Width of each mesh cell in each dimension.
**integer(4) n_tallies**
*do i = 1, n_tallies*
**integer(4) tallies(i) % id**
Unique ID of tally.
**integer(4) tallies(i) % n_realizations**
Number of realizations for the i-th tally.
**integer(4) size(tallies(i) % scores, 1)**
Total number of score bins for the i-th tally
**integer(4) size(tallies(i) % scores, 2)**
Total number of filter bins for the i-th tally
**integer(4) tallies(i) % n_filters**
*do j = 1, tallies(i) % n_filters*
**integer(4) tallies(i) % filter(j) % type**
Type of tally filter.
**integer(4) tallies(i) % filter(j) % n_bins**
Number of bins for filter.
**integer(4)/real(8) tallies(i) % filter(j) % bins(:)**
Value for each filter bin of this type.
**integer(4) tallies(i) % n_nuclide_bins**
Number of nuclide bins. If none are specified, this is just one.
*do j = 1, tallies(i) % n_nuclide_bins*
**integer(4) tallies(i) % nuclide_bins(j)**
Values of specified nuclide bins
**integer(4) tallies(i) % n_score_bins**
Number of scoring bins.
*do j = 1, tallies(i) % n_score_bins*
**integer(4) tallies(i) % score_bins(j)**
Values of specified scoring bins (e.g. SCORE_FLUX).
*do j = 1, tallies(i) % n_score_bins*
**integer(4) tallies(i) % scatt_order(j)**
Scattering Order specified scoring bins.
**integer(4) tallies(i) % n_score_bins**
Number of scoring bins without accounting for those added by
the scatter-pn command.
**integer(4) n_realizations**
Number of realizations for global tallies.
**integer(4) N_GLOBAL_TALLIES**
Number of global tally scores
*do i = 1, N_GLOBAL_TALLIES*
**real(8) global_tallies(i) % sum**
Accumulated sum for the i-th global tally
**real(8) global_tallies(i) % sum_sq**
Accumulated sum of squares for the i-th global tally
**integer(4) tallies_on**
Flag indicated if tallies are present in the file.
if (tallies_on > 0)
*do i = 1, n_tallies*
*do k = 1, size(tallies(i) % scores, 2)*
*do j = 1, size(tallies(i) % scores, 1)*
**real(8) tallies(i) % scores(j,k) % sum**
Accumulated sum for the j-th score and k-th filter of the
i-th tally
**real(8) tallies(i) % scores(j,k) % sum_sq**
Accumulated sum of squares for the j-th score and k-th
filter of the i-th tally
if (run_mode == MODE_EIGENVALUE)
*do i = 1, n_particles*
**real(8) source_bank(i) % wgt**
Weight of the i-th source particle
**real(8) source_bank(i) % xyz(1:3)**
Coordinates of the i-th source particle.
**real(8) source_bank(i) % uvw(1:3)**
Direction of the i-th source particle
**real(8) source_bank(i) % E**
Energy of the i-th source particle.
----------
Revision 8
----------

View file

@ -1139,6 +1139,22 @@ The ``<begin>`` element controls what batch CMFD calculations should begin.
*Default*: 1
``<display>`` Element
---------------------
The ``<display>`` element sets one additional CMFD output column. Options are:
* "balance" - prints the RMS [%] of the resdiual from the neutron balance equation
on CMFD tallies.
* "dominance" - prints the estimated dominance ratio from the CMFD iterations.
**This will only work for power iteration eigensolver**.
* "entropy" - prints the *entropy* of the CMFD predicted fission source.
**Can only be used if OpenMC entropy is active as well**.
* "source" - prints the RMS [%] between the OpenMC fission source and CMFD
fission source.
*Default*: None
``<feedback>`` Element
----------------------
@ -1158,14 +1174,14 @@ with "true" and off with "false"
*Default*: true
``<keff_tol>`` Element
----------------------
``<inactive_flush>`` Element
----------------------------
The ``<keff_tol>`` element specifies acceptance criteria of a CMFD eigenvalue.
If the CMFD eigenvalue and OpenMC batch eigenvalue are within this tolerance,
CMFD is allowed to modify source neutron weights.
The ``<inactive_flush>`` element controls when CMFD tallies are reset during
inactive batches. The integer set here is the interval at which this reset
occurs. The amout of resets is controlled with the ``<num_flushes>`` element.
*Default*: 0.005
*Defualt*: 9999
``<ksp_monitor>`` Element
-------------------------
@ -1249,14 +1265,13 @@ not impact the calculation.
*Default*: 1.0
``<n_procs_cmfd>`` Element
--------------------------
``<num_flushes>`` Element
-------------------------
The ``<n_procs_cmfd>`` element is used to set the number of processors used
for CMFD calculation. It should be less than or equal to the number of
processors used during OpenMC.
The ``<num_flushes>`` element controls the number of CMFD tally resets that
occur during inactive CMFD batches.
*Default*: 1
*Default*: 9999
``<power_monitor>`` Element
---------------------------
@ -1264,26 +1279,33 @@ processors used during OpenMC.
The ``<power_monitor>`` element is used to view the convergence of power iteration.
This option can be turned on with "true" and turned off with "false".
*Default*: false
``<run_adjoint>`` Element
-------------------------
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.
*Default*: false
``<write_balance>`` Element
---------------------------
``<snes_monitor>`` Element
--------------------------
The ``<write_balance>`` element is used to view the balance of OpenMC tally
residuals for every coarse mesh region and energy group. This option can be
turned on with "true" and off with "false".
The ``<snes_monitor>`` element is used to view the convergence of the nonlinear SNES
function in PETSc. This option can be turned on with "true" and turned off with "false".
*Default*: false
``<write_hdf5>`` Element
------------------------
``<solver>`` Element
--------------------
The ``<write_hdf5>`` element can be turned on with "true" to get an
HDF5 output file of CMFD results.
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*: false
*Default*: power
``<write_matrices>`` Element
----------------------------

View file

@ -71,13 +71,29 @@ Prerequisites
To compile with support for HDF5_ output (highly recommended), you will
need to have HDF5 installed on your computer. The installed version will
need to have been compiled with the same compiler you intend to compile
OpenMC with.
OpenMC with. HDF5_ must be built with parallel I/O features if you intend
to use HDF5_ with MPI. An example of configuring HDF5_ is listed below::
FC=/opt/mpich/3.0.4-gnu/bin/mpif90 CC=/opt/mpich/3.0.4-gnu/bin/mpicc \
./configure --prefix=/opt/hdf5/1.8.11-gnu --enable-fortran \
--enable-fortran2003 --enable-parallel
You may omit '--enable-parallel' if you want to compile HDF5_ in serial.
* PETSc_ for CMFD acceleration
To enable CMFD acceleration, you will need to have PETSc_ installed on
your computer. The installed version will need to have been compiled with
the same compiler you intend to compile OpenMC with.
To enable 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.2-gnu --download-f-blas-lapack \
--with-mpi-dir=/opt/mpich/3.0.4-gnu/ --with-shared-libraries=0 \
--with-fortran-datatypes
The BLAS/LAPACK library is not required to be downloaded and can be linked
explicitly (e.g., Intel MLK library).
* git_ version control software for obtaining source code

View file

@ -1,37 +1,29 @@
set_header.o: constants.o
set_header.o: list_header.o
ace.o: ace_header.o
ace.o: constants.o
ace.o: endf.o
ace.o: error.o
ace.o: fission.o
ace.o: global.o
ace.o: material_header.o
ace.o: output.o
ace.o: set_header.o
ace.o: string.o
energy_grid.o: constants.o
energy_grid.o: global.o
energy_grid.o: list_header.o
energy_grid.o: output.o
ace_header.o: constants.o
ace_header.o: endf_header.o
list_header.o: constants.o
cmfd_slepc_solver.o: cmfd_loss_operator.o
cmfd_slepc_solver.o: cmfd_prod_operator.o
cmfd_slepc_solver.o: constants.o
cmfd_slepc_solver.o: global.o
cmfd_loss_operator.o: constants.o
cmfd_loss_operator.o: global.o
particle_restart.o: bank_header.o
particle_restart.o: constants.o
particle_restart.o: geometry_header.o
particle_restart.o: global.o
particle_restart.o: output.o
particle_restart.o: output_interface.o
particle_restart.o: particle_header.o
particle_restart.o: random_lcg.o
particle_restart.o: tracking.o
doppler.o: constants.o
cmfd_data.o: cmfd_header.o
cmfd_data.o: constants.o
cmfd_data.o: error.o
cmfd_data.o: global.o
cmfd_data.o: mesh.o
cmfd_data.o: mesh_header.o
cmfd_data.o: string.o
cmfd_data.o: tally_header.o
cmfd_execute.o: cmfd_data.o
cmfd_execute.o: cmfd_message_passing.o
cmfd_execute.o: cmfd_jfnk_solver.o
cmfd_execute.o: cmfd_power_solver.o
cmfd_execute.o: cmfd_snes_solver.o
cmfd_execute.o: constants.o
cmfd_execute.o: error.o
cmfd_execute.o: global.o
@ -41,27 +33,99 @@ cmfd_execute.o: output.o
cmfd_execute.o: search.o
cmfd_execute.o: tally.o
cmfd_header.o: constants.o
cmfd_input.o: cmfd_header.o
cmfd_input.o: error.o
cmfd_input.o: global.o
cmfd_input.o: mesh_header.o
cmfd_input.o: output.o
cmfd_input.o: string.o
cmfd_input.o: tally.o
cmfd_input.o: tally_header.o
cmfd_input.o: tally_initialize.o
cmfd_input.o: templates/cmfd_t.o
cmfd_jfnk_solver.o: cmfd_loss_operator.o
cmfd_jfnk_solver.o: cmfd_power_solver.o
cmfd_jfnk_solver.o: cmfd_prod_operator.o
cmfd_jfnk_solver.o: constants.o
cmfd_jfnk_solver.o: global.o
cmfd_jfnk_solver.o: matrix_header.o
cmfd_jfnk_solver.o: solver_interface.o
cmfd_jfnk_solver.o: vector_header.o
cmfd_loss_operator.o: constants.o
cmfd_loss_operator.o: global.o
cmfd_loss_operator.o: matrix_header.o
cmfd_power_solver.o: cmfd_loss_operator.o
cmfd_power_solver.o: cmfd_prod_operator.o
cmfd_power_solver.o: constants.o
cmfd_power_solver.o: global.o
cmfd_power_solver.o: string.o
cmfd_power_solver.o: matrix_header.o
cmfd_power_solver.o: solver_interface.o
cmfd_power_solver.o: vector_header.o
cmfd_message_passing.o: cmfd_header.o
cmfd_message_passing.o: global.o
cmfd_prod_operator.o: constants.o
cmfd_prod_operator.o: global.o
cmfd_prod_operator.o: matrix_header.o
random_lcg.o: global.o
cmfd_slepc_solver.o: cmfd_loss_operator.o
cmfd_slepc_solver.o: cmfd_prod_operator.o
cmfd_slepc_solver.o: constants.o
cmfd_slepc_solver.o: global.o
plot.o: constants.o
plot.o: error.o
plot.o: geometry.o
plot.o: geometry_header.o
plot.o: global.o
plot.o: output.o
plot.o: particle_header.o
plot.o: plot_header.o
plot.o: ppmlib.o
plot.o: string.o
cross_section.o: ace_header.o
cross_section.o: constants.o
cross_section.o: error.o
cross_section.o: fission.o
cross_section.o: global.o
cross_section.o: material_header.o
cross_section.o: particle_header.o
cross_section.o: random_lcg.o
cross_section.o: search.o
doppler.o: constants.o
eigenvalue.o: cmfd_execute.o
eigenvalue.o: constants.o
eigenvalue.o: error.o
eigenvalue.o: global.o
eigenvalue.o: math.o
eigenvalue.o: mesh.o
eigenvalue.o: mesh_header.o
eigenvalue.o: output.o
eigenvalue.o: particle_header.o
eigenvalue.o: random_lcg.o
eigenvalue.o: search.o
eigenvalue.o: source.o
eigenvalue.o: state_point.o
eigenvalue.o: string.o
eigenvalue.o: tally.o
eigenvalue.o: tracking.o
endf.o: constants.o
endf.o: string.o
energy_grid.o: constants.o
energy_grid.o: global.o
energy_grid.o: list_header.o
energy_grid.o: output.o
error.o: global.o
finalize.o: global.o
finalize.o: hdf5_interface.o
finalize.o: output.o
finalize.o: tally.o
fission.o: ace_header.o
fission.o: constants.o
fission.o: error.o
fission.o: global.o
fission.o: interpolation.o
fission.o: search.o
fixed_source.o: constants.o
fixed_source.o: global.o
@ -74,114 +138,15 @@ fixed_source.o: string.o
fixed_source.o: tally.o
fixed_source.o: tracking.o
source.o: bank_header.o
source.o: constants.o
source.o: error.o
source.o: geometry_header.o
source.o: global.o
source.o: math.o
source.o: output.o
source.o: particle_header.o
source.o: random_lcg.o
source.o: string.o
cmfd_prod_operator.o: constants.o
cmfd_prod_operator.o: global.o
ace_header.o: constants.o
ace_header.o: endf_header.o
fission.o: ace_header.o
fission.o: constants.o
fission.o: error.o
fission.o: global.o
fission.o: interpolation.o
fission.o: search.o
cmfd_jacobian_operator.o: cmfd_loss_operator.o
cmfd_jacobian_operator.o: cmfd_prod_operator.o
cmfd_jacobian_operator.o: constants.o
cmfd_jacobian_operator.o: global.o
cmfd_snes_solver.o: cmfd_jacobian_operator.o
cmfd_snes_solver.o: cmfd_loss_operator.o
cmfd_snes_solver.o: cmfd_power_solver.o
cmfd_snes_solver.o: cmfd_prod_operator.o
cmfd_snes_solver.o: constants.o
cmfd_snes_solver.o: global.o
cmfd_snes_solver.o: string.o
cmfd_input.o: cmfd_message_passing.o
cmfd_input.o: error.o
cmfd_input.o: global.o
cmfd_input.o: mesh_header.o
cmfd_input.o: output.o
cmfd_input.o: string.o
cmfd_input.o: tally.o
cmfd_input.o: tally_header.o
cmfd_input.o: tally_initialize.o
cmfd_input.o: templates/cmfd_t.o
input_xml.o: cmfd_input.o
input_xml.o: constants.o
input_xml.o: dict_header.o
input_xml.o: error.o
input_xml.o: geometry_header.o
input_xml.o: global.o
input_xml.o: list_header.o
input_xml.o: mesh_header.o
input_xml.o: output.o
input_xml.o: plot_header.o
input_xml.o: random_lcg.o
input_xml.o: string.o
input_xml.o: tally_header.o
input_xml.o: tally_initialize.o
input_xml.o: templates/cross_sections_t.o
input_xml.o: templates/geometry_t.o
input_xml.o: templates/materials_t.o
input_xml.o: templates/plots_t.o
input_xml.o: templates/settings_t.o
input_xml.o: templates/tallies_t.o
main.o: constants.o
main.o: eigenvalue.o
main.o: finalize.o
main.o: fixed_source.o
main.o: global.o
main.o: initialize.o
main.o: particle_restart.o
main.o: plot.o
particle_restart_write.o: bank_header.o
particle_restart_write.o: global.o
particle_restart_write.o: output_interface.o
particle_restart_write.o: particle_header.o
particle_restart_write.o: string.o
timer_header.o: constants.o
math.o: constants.o
math.o: random_lcg.o
interpolation.o: constants.o
interpolation.o: endf_header.o
interpolation.o: error.o
interpolation.o: global.o
interpolation.o: search.o
interpolation.o: string.o
cmfd_data.o: cmfd_header.o
cmfd_data.o: constants.o
cmfd_data.o: error.o
cmfd_data.o: global.o
cmfd_data.o: mesh.o
cmfd_data.o: mesh_header.o
cmfd_data.o: tally_header.o
cmfd_output.o: cmfd_data.o
cmfd_output.o: cmfd_header.o
cmfd_output.o: constants.o
cmfd_output.o: global.o
geometry.o: constants.o
geometry.o: error.o
geometry.o: geometry_header.o
geometry.o: global.o
geometry.o: output.o
geometry.o: particle_header.o
geometry.o: particle_restart_write.o
geometry.o: string.o
geometry.o: tally.o
global.o: ace_header.o
global.o: bank_header.o
@ -198,51 +163,17 @@ global.o: source_header.o
global.o: tally_header.o
global.o: timer_header.o
string.o: constants.o
string.o: error.o
string.o: global.o
finalize.o: cmfd_output.o
finalize.o: global.o
finalize.o: hdf5_interface.o
finalize.o: output.o
finalize.o: tally.o
tally_initialize.o: constants.o
tally_initialize.o: global.o
tally_initialize.o: tally_header.o
tally.o: ace_header.o
tally.o: constants.o
tally.o: error.o
tally.o: global.o
tally.o: math.o
tally.o: mesh.o
tally.o: mesh_header.o
tally.o: output.o
tally.o: particle_header.o
tally.o: search.o
tally.o: string.o
tally.o: tally_header.o
particle_header.o: constants.o
particle_header.o: geometry_header.o
output_interface.o: constants.o
output_interface.o: error.o
output_interface.o: global.o
output_interface.o: hdf5_interface.o
output_interface.o: mpiio_interface.o
output_interface.o: tally_header.o
mesh.o: constants.o
mesh.o: global.o
mesh.o: mesh_header.o
mesh.o: particle_header.o
mesh.o: search.o
endf.o: constants.o
endf.o: string.o
hdf5_summary.o: ace_header.o
hdf5_summary.o: constants.o
hdf5_summary.o: endf.o
hdf5_summary.o: geometry_header.o
hdf5_summary.o: global.o
hdf5_summary.o: material_header.o
hdf5_summary.o: mesh_header.o
hdf5_summary.o: output.o
hdf5_summary.o: output_interface.o
hdf5_summary.o: string.o
hdf5_summary.o: tally_header.o
initialize.o: ace.o
initialize.o: bank_header.o
@ -265,97 +196,56 @@ initialize.o: string.o
initialize.o: tally_header.o
initialize.o: tally_initialize.o
cross_section.o: ace_header.o
cross_section.o: constants.o
cross_section.o: error.o
cross_section.o: fission.o
cross_section.o: global.o
cross_section.o: material_header.o
cross_section.o: particle_header.o
cross_section.o: random_lcg.o
cross_section.o: search.o
input_xml.o: cmfd_input.o
input_xml.o: constants.o
input_xml.o: dict_header.o
input_xml.o: error.o
input_xml.o: geometry_header.o
input_xml.o: global.o
input_xml.o: list_header.o
input_xml.o: mesh_header.o
input_xml.o: output.o
input_xml.o: plot_header.o
input_xml.o: random_lcg.o
input_xml.o: string.o
input_xml.o: tally_header.o
input_xml.o: tally_initialize.o
input_xml.o: templates/cross_sections_t.o
input_xml.o: templates/geometry_t.o
input_xml.o: templates/materials_t.o
input_xml.o: templates/plots_t.o
input_xml.o: templates/settings_t.o
input_xml.o: templates/tallies_t.o
state_point.o: constants.o
state_point.o: error.o
state_point.o: global.o
state_point.o: output.o
state_point.o: output_interface.o
state_point.o: string.o
state_point.o: tally_header.o
interpolation.o: constants.o
interpolation.o: endf_header.o
interpolation.o: error.o
interpolation.o: global.o
interpolation.o: search.o
interpolation.o: string.o
eigenvalue.o: cmfd_execute.o
eigenvalue.o: constants.o
eigenvalue.o: error.o
eigenvalue.o: global.o
eigenvalue.o: math.o
eigenvalue.o: mesh.o
eigenvalue.o: mesh_header.o
eigenvalue.o: output.o
eigenvalue.o: particle_header.o
eigenvalue.o: random_lcg.o
eigenvalue.o: search.o
eigenvalue.o: source.o
eigenvalue.o: state_point.o
eigenvalue.o: string.o
eigenvalue.o: tally.o
eigenvalue.o: tracking.o
list_header.o: constants.o
search.o: error.o
search.o: global.o
main.o: constants.o
main.o: eigenvalue.o
main.o: finalize.o
main.o: fixed_source.o
main.o: global.o
main.o: initialize.o
main.o: particle_restart.o
main.o: plot.o
tracking.o: cross_section.o
tracking.o: error.o
tracking.o: geometry.o
tracking.o: geometry_header.o
tracking.o: global.o
tracking.o: output.o
tracking.o: particle_header.o
tracking.o: physics.o
tracking.o: random_lcg.o
tracking.o: string.o
tracking.o: tally.o
math.o: constants.o
math.o: random_lcg.o
ace.o: ace_header.o
ace.o: constants.o
ace.o: endf.o
ace.o: error.o
ace.o: fission.o
ace.o: global.o
ace.o: material_header.o
ace.o: output.o
ace.o: set_header.o
ace.o: string.o
matrix_header.o: constants.o
matrix_header.o: vector_header.o
geometry.o: constants.o
geometry.o: error.o
geometry.o: geometry_header.o
geometry.o: global.o
geometry.o: output.o
geometry.o: particle_header.o
geometry.o: particle_restart_write.o
geometry.o: string.o
geometry.o: tally.o
plot_header.o: constants.o
cmfd_header.o: constants.o
tally_header.o: constants.o
hdf5_summary.o: ace_header.o
hdf5_summary.o: constants.o
hdf5_summary.o: endf.o
hdf5_summary.o: geometry_header.o
hdf5_summary.o: global.o
hdf5_summary.o: hdf5_interface.o
hdf5_summary.o: material_header.o
hdf5_summary.o: mesh_header.o
hdf5_summary.o: output.o
hdf5_summary.o: output_interface.o
hdf5_summary.o: string.o
hdf5_summary.o: tally_header.o
error.o: global.o
mesh.o: constants.o
mesh.o: global.o
mesh.o: mesh_header.o
mesh.o: particle_header.o
mesh.o: search.o
output.o: ace_header.o
output.o: constants.o
@ -371,6 +261,32 @@ output.o: plot_header.o
output.o: string.o
output.o: tally_header.o
output_interface.o: constants.o
output_interface.o: error.o
output_interface.o: global.o
output_interface.o: hdf5_interface.o
output_interface.o: mpiio_interface.o
output_interface.o: tally_header.o
particle_header.o: constants.o
particle_header.o: geometry_header.o
particle_restart.o: bank_header.o
particle_restart.o: constants.o
particle_restart.o: geometry_header.o
particle_restart.o: global.o
particle_restart.o: output.o
particle_restart.o: output_interface.o
particle_restart.o: particle_header.o
particle_restart.o: random_lcg.o
particle_restart.o: tracking.o
particle_restart_write.o: bank_header.o
particle_restart_write.o: global.o
particle_restart_write.o: output_interface.o
particle_restart_write.o: particle_header.o
particle_restart_write.o: string.o
physics.o: ace_header.o
physics.o: constants.o
physics.o: endf.o
@ -388,3 +304,87 @@ physics.o: random_lcg.o
physics.o: search.o
physics.o: string.o
plot.o: constants.o
plot.o: error.o
plot.o: geometry.o
plot.o: geometry_header.o
plot.o: global.o
plot.o: output.o
plot.o: particle_header.o
plot.o: plot_header.o
plot.o: ppmlib.o
plot.o: string.o
plot_header.o: constants.o
random_lcg.o: global.o
search.o: error.o
search.o: global.o
set_header.o: constants.o
set_header.o: list_header.o
solver_interface.o: error.o
solver_interface.o: global.o
solver_interface.o: matrix_header.o
solver_interface.o: vector_header.o
source.o: bank_header.o
source.o: constants.o
source.o: error.o
source.o: geometry_header.o
source.o: global.o
source.o: math.o
source.o: output.o
source.o: particle_header.o
source.o: random_lcg.o
source.o: string.o
state_point.o: constants.o
state_point.o: error.o
state_point.o: global.o
state_point.o: output.o
state_point.o: output_interface.o
state_point.o: string.o
state_point.o: tally_header.o
string.o: constants.o
string.o: error.o
string.o: global.o
tally.o: ace_header.o
tally.o: constants.o
tally.o: error.o
tally.o: global.o
tally.o: math.o
tally.o: mesh.o
tally.o: mesh_header.o
tally.o: output.o
tally.o: particle_header.o
tally.o: search.o
tally.o: string.o
tally.o: tally_header.o
tally_header.o: constants.o
tally_initialize.o: constants.o
tally_initialize.o: global.o
tally_initialize.o: tally_header.o
timer_header.o: constants.o
tracking.o: cross_section.o
tracking.o: error.o
tracking.o: geometry.o
tracking.o: geometry_header.o
tracking.o: global.o
tracking.o: output.o
tracking.o: particle_header.o
tracking.o: physics.o
tracking.o: random_lcg.o
tracking.o: string.o
tracking.o: tally.o
vector_header.o: constants.o

View file

@ -42,12 +42,12 @@ GIT_SHA1 = $(shell git log -1 2>/dev/null | head -n 1 | awk '{print $$2}')
ifeq ($(COMPILER),gnu)
F90 = gfortran
F90FLAGS := -cpp -fbacktrace
F90FLAGS := -cpp -std=f2008 -fbacktrace
LDFLAGS =
# Debugging
ifeq ($(DEBUG),yes)
F90FLAGS += -g -Wall -pedantic -std=f2008 -fbounds-check \
F90FLAGS += -g -Wall -pedantic -fbounds-check \
-ffpe-trap=invalid,overflow,underflow
LDFLAGS += -g
endif
@ -70,12 +70,12 @@ endif
ifeq ($(COMPILER),intel)
F90 = ifort
F90FLAGS := -fpp -warn -assume byterecl -traceback
F90FLAGS := -fpp -std08 -assume byterecl -traceback
LDFLAGS =
# Debugging
ifeq ($(DEBUG),yes)
F90FLAGS += -g -ftrapuv -fp-stack-check -check all -fpe0
F90FLAGS += -g -warn -ftrapuv -fp-stack-check -check all -fpe0
LDFLAGS += -g
endif

File diff suppressed because it is too large Load diff

View file

@ -5,81 +5,57 @@ module cmfd_execute
! cross section generation, diffusion calculation, and source re-weighting
!==============================================================================
use global
implicit none
private
public :: execute_cmfd, cmfd_init_batch
# ifdef PETSC
# include <finclude/petsc.h90>
# endif
contains
!==============================================================================
! EXECUTE_CMFD
! EXECUTE_CMFD runs the CMFD calculation
!==============================================================================
subroutine execute_cmfd()
# ifdef PETSC
use cmfd_data, only: set_up_cmfd
use cmfd_message_passing, only: petsc_init_mpi, cmfd_bcast
use cmfd_power_solver, only: cmfd_power_execute
use cmfd_snes_solver, only: cmfd_snes_execute
use cmfd_jfnk_solver, only: cmfd_jfnk_execute
use error, only: warning, fatal_error
use global, only: n_procs_cmfd, cmfd, &
cmfd_solver_type, time_cmfd, &
cmfd_run_adjoint, cmfd_write_hdf5, &
cmfd_feedback,cmfd_hold_weights, &
cmfd_inact_flush, cmfd_keff_tol, &
cmfd_act_flush, current_batch, keff, &
n_batches, message, master, mpi_err, rank
! stop cmfd timer
! CMFD single processor on master
if (master) then
! Start cmfd timer
call time_cmfd % start()
end if
! filter processors (lowest PETSc group)
if (rank < n_procs_cmfd) then
! Create cmfd data from OpenMC tallies
call set_up_cmfd()
! set up cmfd data (master only)
if (master) call set_up_cmfd()
! broadcast cmfd to all petsc procs
call cmfd_bcast()
! process solver options
! Process solver options
call process_cmfd_options()
end if
! filter processors (lowest PETSc group)
if (rank < n_procs_cmfd) then
! call solver
! Call solver
if (trim(cmfd_solver_type) == 'power') then
call cmfd_power_execute()
elseif (trim(cmfd_solver_type) == 'jfnk') then
call cmfd_snes_execute()
call cmfd_jfnk_execute()
else
message = 'solver type became invalid after input processing'
call fatal_error()
end if
! perform any last batch tasks
if (current_batch == n_batches) then
! Save k-effective
cmfd % k_cmfd(current_batch) = cmfd % keff
! check for adjoint run
if (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_snes_execute(adjoint = .true.)
end if
! 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
end if
end if
@ -91,19 +67,12 @@ contains
if (cmfd_feedback) call cmfd_reweight(.true.)
! stop cmfd timer
if (master) then
call time_cmfd % stop()
end if
! wait here for all procs
call MPI_Barrier(MPI_COMM_WORLD, mpi_err)
# endif
if (master) call time_cmfd % stop()
end subroutine execute_cmfd
!==============================================================================
! CMFD_INIT_BATCH
! CMFD_INIT_BATCH handles cmfd options at the start of every batch
!==============================================================================
subroutine cmfd_init_batch()
@ -112,7 +81,7 @@ contains
cmfd_inact_flush, cmfd_act_flush, cmfd_run, &
current_batch, cmfd_hold_weights
! check to activate CMFD diffusion and possible feedback
! Check to activate CMFD diffusion and possible feedback
! this guarantees that when cmfd begins at least one batch of tallies are
! accumulated
if (cmfd_run .and. cmfd_begin == current_batch) then
@ -120,18 +89,19 @@ contains
cmfd_tally_on = .true.
end if
! check to flush cmfd tallies for active batches, no more inactive flush
! If this is a restart run and we are just replaying batches leave
if (restart_run .and. current_batch <= restart_batch) return
! Check to flush cmfd tallies for active batches, no more inactive flush
if (cmfd_run .and. cmfd_act_flush == current_batch) then
call cmfd_tally_reset()
cmfd_tally_on = .true.
cmfd_inact_flush(2) = -1
end if
! check to flush cmfd tallies during inactive batches (>= on number of
! Check to flush cmfd tallies during inactive batches (>= on number of
! flushes important as the code will flush on the first batch which we
! dont want to count)
! if (cmfd_run .and. current_batch < n_inactive .and. mod(current_batch-1,cmfd_inact_flush(1)) &
! == 0 .and. cmfd_inact_flush(2) >= 0) then
if (cmfd_run .and. mod(current_batch,cmfd_inact_flush(1)) &
== 0 .and. cmfd_inact_flush(2) > 0 .and. cmfd_begin < current_batch) then
cmfd_hold_weights = .true.
@ -141,23 +111,23 @@ contains
end subroutine cmfd_init_batch
# ifdef PETSC
!==============================================================================
! PROCESS_CMFD_OPTIONS
! PROCESS_CMFD_OPTIONS processes user options that interface with PETSc
!==============================================================================
subroutine process_cmfd_options()
use global, only: cmfd_snes_monitor, cmfd_ksp_monitor, mpi_err
! check for snes monitor
#ifdef PETSC
! Check for snes monitor
if (cmfd_snes_monitor) call PetscOptionsSetValue("-snes_monitor", &
"stdout", mpi_err)
! check for ksp monitor
! Check for ksp monitor
if (cmfd_ksp_monitor) call PetscOptionsSetValue("-ksp_monitor", &
"stdout", mpi_err)
#endif
end subroutine process_cmfd_options
@ -168,39 +138,44 @@ contains
subroutine calc_fission_source()
use constants, only: CMFD_NOACCEL, ZERO, TWO
use global, only: cmfd, cmfd_coremap, master, mpi_err, entropy_on
use global, only: cmfd, cmfd_coremap, master, mpi_err, entropy_on, &
current_batch
#ifdef MPI
use mpi
#endif
integer :: nx ! maximum number of cells in x direction
integer :: ny ! maximum number of cells in y direction
integer :: nz ! maximum number of cells in z direction
integer :: ng ! maximum number of energy groups
integer :: n ! total size
integer :: i ! iteration counter for x
integer :: j ! iteration counter for y
integer :: k ! iteration counter for z
integer :: g ! iteration counter for groups
integer :: idx ! index in vector
integer :: nx ! maximum number of cells in x direction
integer :: ny ! maximum number of cells in y direction
integer :: nz ! maximum number of cells in z direction
integer :: ng ! maximum number of energy groups
integer :: n ! total size
integer :: i ! iteration counter for x
integer :: j ! iteration counter for y
integer :: k ! iteration counter for z
integer :: g ! iteration counter for groups
integer :: idx ! index in vector
real(8) :: hxyz(3) ! cell dimensions of current ijk cell
real(8) :: vol ! volume of cell
real(8),allocatable :: source(:,:,:,:) ! tmp source array for entropy
! get maximum of spatial and group indices
nx = cmfd%indices(1)
ny = cmfd%indices(2)
nz = cmfd%indices(3)
ng = cmfd%indices(4)
! Get maximum of spatial and group indices
nx = cmfd % indices(1)
ny = cmfd % indices(2)
nz = cmfd % indices(3)
ng = cmfd % indices(4)
n = ng*nx*ny*nz
! allocate cmfd source if not already allocated and allocate buffer
if (.not. allocated(cmfd%cmfd_src)) allocate(cmfd%cmfd_src(ng,nx,ny,nz))
! Allocate cmfd source if not already allocated and allocate buffer
if (.not. allocated(cmfd % cmfd_src)) &
allocate(cmfd % cmfd_src(ng,nx,ny,nz))
! reset cmfd source to 0
cmfd%cmfd_src = ZERO
! Reset cmfd source to 0
cmfd % cmfd_src = ZERO
! only perform for master
! Only perform for master
if (master) then
! loop around indices to map to cmfd object
! Loop around indices to map to cmfd object
ZLOOP: do k = 1, nz
YLOOP: do j = 1, ny
@ -209,25 +184,25 @@ contains
GROUP: do g = 1, ng
! check for core map
! Check for core map
if (cmfd_coremap) then
if (cmfd%coremap(i,j,k) == CMFD_NOACCEL) then
if (cmfd % coremap(i,j,k) == CMFD_NOACCEL) then
cycle
end if
end if
! get dimensions of cell
hxyz = cmfd%hxyz(:,i,j,k)
! Get dimensions of cell
hxyz = cmfd % hxyz(:,i,j,k)
! calculate volume
! Calculate volume
vol = hxyz(1)*hxyz(2)*hxyz(3)
! get first index
! Get first index
idx = get_matrix_idx(1,i,j,k,ng,nx,ny)
! compute fission source
cmfd%cmfd_src(g,i,j,k) = sum(cmfd%nfissxs(:,g,i,j,k) * &
cmfd%phi(idx:idx+(ng-1)))*vol
! Compute fission source
cmfd % cmfd_src(g,i,j,k) = sum(cmfd % nfissxs(:,g,i,j,k) * &
cmfd % phi(idx:idx + (ng - 1)))*vol
end do GROUP
@ -237,43 +212,49 @@ contains
end do ZLOOP
! normalize source such that it sums to 1.0
cmfd%cmfd_src = cmfd%cmfd_src/sum(cmfd%cmfd_src)
! Normalize source such that it sums to 1.0
cmfd % cmfd_src = cmfd % cmfd_src/sum(cmfd % cmfd_src)
! compute entropy
! Compute entropy
if (entropy_on) then
! allocate tmp array
! Allocate tmp array
if (.not.allocated(source)) allocate(source(ng,nx,ny,nz))
! initialize the source
! Initialize the source
source = ZERO
! compute log
where (cmfd%cmfd_src > ZERO)
source = cmfd%cmfd_src*log(cmfd%cmfd_src)/log(TWO)
! Compute log
where (cmfd % cmfd_src > ZERO)
source = cmfd % cmfd_src*log(cmfd % cmfd_src)/log(TWO)
end where
! sum that source
cmfd%entropy = -sum(source)
! Sum that source
cmfd % entropy(current_batch) = -sum(source)
! deallocate tmp array
! Deallocate tmp array
if (allocated(source)) deallocate(source)
end if
! normalize source so average is 1.0
cmfd%cmfd_src = cmfd%cmfd_src/sum(cmfd%cmfd_src)*cmfd%norm
! Normalize source so average is 1.0
cmfd % cmfd_src = cmfd % cmfd_src/sum(cmfd % cmfd_src)*cmfd % norm
! Calculate differences between normalized sources
cmfd % src_cmp(current_batch) = sqrt(ONE/cmfd % norm * &
sum((cmfd % cmfd_src - cmfd % openmc_src)**2))
end if
! broadcast full source to all procs
call MPI_BCAST(cmfd%cmfd_src, n, MPI_REAL8, 0, MPI_COMM_WORLD, mpi_err)
#ifdef MPI
! Broadcast full source to all procs
call MPI_BCAST(cmfd % cmfd_src, n, MPI_REAL8, 0, MPI_COMM_WORLD, mpi_err)
#endif
end subroutine calc_fission_source
!===============================================================================
! CMFD_REWEIGHT
! CMFD_REWEIGHT performs weighting of particles in the source bank
!===============================================================================
subroutine cmfd_reweight(new_weights)
@ -286,110 +267,116 @@ contains
use mesh, only: count_bank_sites, get_mesh_indices
use search, only: binary_search
! local variables
integer :: nx ! maximum number of cells in x direction
integer :: ny ! maximum number of cells in y direction
integer :: nz ! maximum number of cells in z direction
integer :: ng ! maximum number of energy groups
integer :: i ! iteration counter
integer :: ijk(3) ! spatial bin location
integer :: e_bin ! energy bin of source particle
integer :: n_groups ! number of energy groups
logical :: outside ! any source sites outside mesh
logical :: in_mesh ! source site is inside mesh
logical :: new_weights ! calcualte new weights
type(StructuredMesh), pointer :: m ! point to mesh
real(8), allocatable :: egrid(:)
#ifdef MPI
use mpi
#endif
! associate pointer
logical, intent(in) :: new_weights ! calcualte new weights
integer :: nx ! maximum number of cells in x direction
integer :: ny ! maximum number of cells in y direction
integer :: nz ! maximum number of cells in z direction
integer :: ng ! maximum number of energy groups
integer :: i ! iteration counter
integer :: ijk(3) ! spatial bin location
integer :: e_bin ! energy bin of source particle
integer :: n_groups ! number of energy groups
logical :: outside ! any source sites outside mesh
logical :: in_mesh ! source site is inside mesh
type(StructuredMesh), pointer :: m ! point to mesh
real(8), allocatable :: egrid(:) ! energy grid
! Associate pointer
m => meshes(n_user_meshes + 1)
! get maximum of spatial and group indices
nx = cmfd%indices(1)
ny = cmfd%indices(2)
nz = cmfd%indices(3)
ng = cmfd%indices(4)
! Get maximum of spatial and group indices
nx = cmfd % indices(1)
ny = cmfd % indices(2)
nz = cmfd % indices(3)
ng = cmfd % indices(4)
! allocate arrays in cmfd object (can take out later extend to multigroup)
if (.not.allocated(cmfd%sourcecounts)) then
allocate(cmfd%sourcecounts(ng,nx,ny,nz))
cmfd % sourcecounts = 0
end if
if (.not.allocated(cmfd%weightfactors)) then
allocate(cmfd%weightfactors(ng,nx,ny,nz))
if (.not.allocated(cmfd % weightfactors)) then
allocate(cmfd % weightfactors(ng,nx,ny,nz))
cmfd % weightfactors = ONE
end if
! allocate energy grid and reverse cmfd energy grid
if (.not. allocated(egrid)) allocate(egrid(ng+1))
egrid = (/(cmfd%egrid(ng-i+2),i = 1,ng+1)/)
! Allocate energy grid and reverse cmfd energy grid
if (.not. allocated(egrid)) allocate(egrid(ng + 1))
egrid = (/(cmfd % egrid(ng - i + 2), i = 1, ng + 1)/)
! compute new weight factors
! Compute new weight factors
if (new_weights) then
! zero out weights
! Zero out weights
cmfd%weightfactors = ZERO
! count bank sites in mesh
! Count bank sites in mesh
call count_bank_sites(m, source_bank, cmfd%sourcecounts, egrid, &
sites_outside=outside, size_bank=work)
! check for sites outside of the mesh
! Check for sites outside of the mesh
if (master .and. outside) then
message = "Source sites outside of the CMFD mesh!"
call fatal_error()
end if
! have master compute weight factors
! Have master compute weight factors (watch for 0s)
if (master) then
where(cmfd%cmfd_src > ZERO .and. cmfd%sourcecounts > ZERO)
cmfd%weightfactors = cmfd%cmfd_src/sum(cmfd%cmfd_src)* &
sum(cmfd%sourcecounts) / cmfd%sourcecounts
where(cmfd % cmfd_src > ZERO .and. cmfd % sourcecounts > ZERO)
cmfd % weightfactors = cmfd % cmfd_src/sum(cmfd % cmfd_src)* &
sum(cmfd % sourcecounts) / cmfd % sourcecounts
end where
end if
! broadcast weight factors to all procs
call MPI_BCAST(cmfd%weightfactors, ng*nx*ny*nz, MPI_REAL8, 0, &
! Broadcast weight factors to all procs
#ifdef MPI
call MPI_BCAST(cmfd % weightfactors, ng*nx*ny*nz, MPI_REAL8, 0, &
MPI_COMM_WORLD, mpi_err)
#endif
end if
! begin loop over source bank
do i = 1, int(work,4)
! determine spatial bin
call get_mesh_indices(m, source_bank(i)%xyz, ijk, in_mesh)
! Determine spatial bin
call get_mesh_indices(m, source_bank(i) % xyz, ijk, in_mesh)
! determine energy bin
n_groups = size(cmfd%egrid) - 1
if (source_bank(i) % E < cmfd%egrid(1)) then
! Determine energy bin
n_groups = size(cmfd % egrid) - 1
if (source_bank(i) % E < cmfd % egrid(1)) then
e_bin = 1
message = 'source pt below energy grid'
message = 'Source pt below energy grid'
call warning()
elseif (source_bank(i) % E > cmfd%egrid(n_groups+1)) then
elseif (source_bank(i) % E > cmfd % egrid(n_groups + 1)) then
e_bin = n_groups
message = 'source pt above energy grid'
message = 'Source pt above energy grid'
call warning()
else
e_bin = binary_search(cmfd%egrid, n_groups + 1, source_bank(i) % E)
e_bin = binary_search(cmfd % egrid, n_groups + 1, source_bank(i) % E)
end if
! reverese energy bin (lowest grp is highest energy bin)
! Reverese energy bin (lowest grp is highest energy bin)
e_bin = n_groups - e_bin + 1
! check for outside of mesh
! Check for outside of mesh
if (.not. in_mesh) then
message = 'Source site found outside of CMFD mesh!'
message = 'Source site found outside of CMFD mesh'
call fatal_error()
end if
! reweight particle
source_bank(i)%wgt = source_bank(i)%wgt * &
cmfd%weightfactors(e_bin,ijk(1),ijk(2),ijk(3))
! Reweight particle
source_bank(i) % wgt = source_bank(i) % wgt * &
cmfd % weightfactors(e_bin, ijk(1), ijk(2), ijk(3))
end do
! deallocate
! Deallocate all
if (allocated(egrid)) deallocate(egrid)
end subroutine cmfd_reweight
@ -403,33 +390,31 @@ contains
use global, only: cmfd, cmfd_coremap
integer :: matidx ! the index location in matrix
integer :: i ! current x index
integer :: j ! current y index
integer :: k ! current z index
integer :: g ! current group index
integer :: nx ! maximum number of cells in x direction
integer :: ny ! maximum number of cells in y direction
integer :: ng ! maximum number of energy groups
integer, intent(in) :: i ! current x index
integer, intent(in) :: j ! current y index
integer, intent(in) :: k ! current z index
integer, intent(in) :: g ! current group index
integer, intent(in) :: nx ! maximum number of cells in x direction
integer, intent(in) :: ny ! maximum number of cells in y direction
integer, intent(in) :: ng ! maximum number of energy groups
! check if coremap is used
! Check if coremap is used
if (cmfd_coremap) then
! get idx from core map
! Get idx from core map
matidx = ng*(cmfd % coremap(i,j,k)) - (ng - g)
else
! compute index
! Compute index
matidx = g + ng*(i - 1) + ng*nx*(j - 1) + ng*nx*ny*(k - 1)
end if
end function get_matrix_idx
# endif
!===============================================================================
! CMFD_TALLY_RESET
! CMFD_TALLY_RESET resets all cmfd tallies
!===============================================================================
subroutine cmfd_tally_reset()
@ -440,14 +425,14 @@ contains
integer :: i ! loop counter
! print message
! Print message
message = "CMFD tallies reset"
call write_message(7)
! begin loop around CMFD tallies
! Begin loop around CMFD tallies
do i = 1, n_cmfd_tallies
! reset that tally
! Reset that tally
cmfd_tallies(i) % n_realizations = 0
call reset_result(cmfd_tallies(i) % results)

View file

@ -8,93 +8,106 @@ module cmfd_header
type, public :: cmfd_type
! indices for problem
! Indices for problem
integer :: indices(4)
! albedo boundary condition
! Albedo boundary condition
real(8) :: albedo(6)
! core overlay map
! Core overlay map
integer, allocatable :: coremap(:,:,:)
integer, allocatable :: indexmap(:,:)
integer :: mat_dim = CMFD_NOACCEL
! energy grid
! Energy grid
real(8), allocatable :: egrid(:)
! cross sections
! Cross sections
real(8), allocatable :: totalxs(:,:,:,:)
real(8), allocatable :: p1scattxs(:,:,:,:)
real(8), allocatable :: scattxs(:,:,:,:,:)
real(8), allocatable :: nfissxs(:,:,:,:,:)
! diffusion coefficient
! Diffusion coefficient
real(8), allocatable :: diffcof(:,:,:,:)
! current
! Current
real(8), allocatable :: current(:,:,:,:,:)
! flux
! Flux
real(8), allocatable :: flux(:,:,:,:)
! coupling coefficients and equivalence parameters
! Coupling coefficients and equivalence parameters
real(8), allocatable :: dtilde(:,:,:,:,:)
real(8), allocatable :: dhat(:,:,:,:,:)
! dimensions of mesh cells ([hu,hv,hw],xloc,yloc,zloc)
! Dimensions of mesh cells ([hu,hv,hw],xloc,yloc,zloc)
real(8), allocatable :: hxyz(:,:,:,:)
! source distributions
! Source distributions
real(8), allocatable :: cmfd_src(:,:,:,:)
real(8), allocatable :: openmc_src(:,:,:,:)
! source sites in each mesh box
! Source sites in each mesh box
real(8), allocatable :: sourcecounts(:,:,:,:)
! weight adjustment factors
! Weight adjustment factors
real(8), allocatable :: weightfactors(:,:,:,:)
! eigenvector/eigenvalue from cmfd run
! Eigenvector/eigenvalue from cmfd run
real(8), allocatable :: phi(:)
real(8) :: keff = ZERO
! eigenvector/eigenvalue from adjoint run
! Eigenvector/eigenvalue from adjoint run
real(8), allocatable :: adj_phi(:)
real(8) :: adj_keff = ZERO
! residual for neutron balance
! Residual for neutron balance
real(8), allocatable :: resnb(:,:,:,:)
! openmc source normalization factor
! Openmc source normalization factor
real(8) :: norm = ONE
! Shannon entropy from cmfd fission source
real(8) :: entropy
! "Shannon entropy" from cmfd fission source
real(8), allocatable :: entropy(:)
! RMS of neutron balance equations
real(8), allocatable :: balance(:)
! RMS deviation of OpenMC and CMFD normalized source
real(8), allocatable :: src_cmp(:)
! Dominance ratio
real(8), allocatable :: dom(:)
! List of CMFD k
real(8), allocatable :: k_cmfd(:)
end type cmfd_type
contains
!==============================================================================
! ALLOCATE_CMFD
! ALLOCATE_CMFD allocates all data in of cmfd type
!==============================================================================
subroutine allocate_cmfd(this)
subroutine allocate_cmfd(this, n_batches)
type(cmfd_type) :: this
integer, intent(in) :: n_batches ! number of batches in calc
type(cmfd_type), intent(inout) :: this ! cmfd instance
integer :: nx ! number of mesh cells in x direction
integer :: ny ! number of mesh cells in y direction
integer :: nz ! number of mesh cells in z direction
integer :: ng ! number of energy groups
! extract spatial and energy indices from object
! Extract spatial and energy indices from object
nx = this % indices(1)
ny = this % indices(2)
nz = this % indices(3)
ng = this % indices(4)
! allocate flux, cross sections and diffusion coefficient
! Allocate flux, cross sections and diffusion coefficient
if (.not. allocated(this % flux)) allocate(this % flux(ng,nx,ny,nz))
if (.not. allocated(this % totalxs)) allocate(this % totalxs(ng,nx,ny,nz))
if (.not. allocated(this % p1scattxs)) allocate(this % p1scattxs(ng,nx,ny,nz))
@ -102,25 +115,32 @@ contains
if (.not. allocated(this % nfissxs)) allocate(this % nfissxs(ng,ng,nx,ny,nz))
if (.not. allocated(this % diffcof)) allocate(this % diffcof(ng,nx,ny,nz))
! allocate dtilde and dhat
! Allocate dtilde and dhat
if (.not. allocated(this % dtilde)) allocate(this % dtilde(6,ng,nx,ny,nz))
if (.not. allocated(this % dhat)) allocate(this % dhat(6,ng,nx,ny,nz))
! allocate dimensions for each box (here for general case)
! Allocate dimensions for each box (here for general case)
if (.not. allocated(this % hxyz)) allocate(this % hxyz(3,nx,ny,nz))
! allocate surface currents
! Allocate surface currents
if (.not. allocated(this % current)) allocate(this % current(12,ng,nx,ny,nz))
! allocate source distributions
! Allocate source distributions
if (.not. allocated(this % cmfd_src)) allocate(this % cmfd_src(ng,nx,ny,nz))
if (.not. allocated(this % openmc_src)) allocate(this % openmc_src(ng,nx,ny,nz))
! allocate source weight modification vars
! Allocate source weight modification vars
if (.not. allocated(this % sourcecounts)) allocate(this % sourcecounts(ng,nx,ny,nz))
if (.not. allocated(this % weightfactors)) allocate(this % weightfactors(ng,nx,ny,nz))
! set everthing to 0 except weight multiply factors if feedback isnt on
! Allocate batchwise parameters
if (.not. allocated(this % entropy)) allocate(this % entropy(n_batches))
if (.not. allocated(this % balance)) allocate(this % balance(n_batches))
if (.not. allocated(this % src_cmp)) allocate(this % src_cmp(n_batches))
if (.not. allocated(this % dom)) allocate(this % dom(n_batches))
if (.not. allocated(this % k_cmfd)) allocate(this % k_cmfd(n_batches))
! Set everthing to 0 except weight multiply factors if feedback isnt on
this % flux = ZERO
this % totalxs = ZERO
this % p1scattxs = ZERO
@ -135,16 +155,21 @@ contains
this % openmc_src = ZERO
this % sourcecounts = ZERO
this % weightfactors = ONE
this % balance = ZERO
this % src_cmp = ZERO
this % dom = ZERO
this % k_cmfd = ZERO
this % entropy = ZERO
end subroutine allocate_cmfd
!===============================================================================
! DEALLOCATE_CMFD
! DEALLOCATE_CMFD frees all memory of cmfd type
!===============================================================================
subroutine deallocate_cmfd(this)
type(cmfd_type) :: this
type(cmfd_type), intent(inout) :: this ! cmfd instance
if (allocated(this % egrid)) deallocate(this % egrid)
if (allocated(this % totalxs)) deallocate(this % totalxs)
@ -164,6 +189,11 @@ contains
if (allocated(this % weightfactors)) deallocate(this % weightfactors)
if (allocated(this % cmfd_src)) deallocate(this % cmfd_src)
if (allocated(this % openmc_src)) deallocate(this % openmc_src)
if (allocated(this % balance)) deallocate(this % balance)
if (allocated(this % src_cmp)) deallocate(this % src_cmp)
if (allocated(this % dom)) deallocate(this % dom)
if (allocated(this % k_cmfd)) deallocate(this % k_cmfd)
if (allocated(this % entropy)) deallocate(this % entropy)
end subroutine deallocate_cmfd

View file

@ -1,5 +1,11 @@
module cmfd_input
use global
#ifdef PETSC
use petscsys
#endif
implicit none
private
public :: configure_cmfd
@ -7,23 +13,47 @@ module cmfd_input
contains
!===============================================================================
! CONFIGURE_CMFD
! CONFIGURE_CMFD initializes PETSc and CMFD parameters
!===============================================================================
subroutine configure_cmfd()
# ifdef PETSC
use cmfd_message_passing, only: petsc_init_mpi
# endif
use cmfd_header, only: allocate_cmfd
! read in cmfd input file
integer :: new_comm ! new mpi communicator
integer :: color ! color group of processor
! Read in cmfd input file
call read_cmfd_xml()
! initialize petsc on mpi
# ifdef PETSC
call petsc_init_mpi()
! Assign color
if (master) then
color = 1
else
color = 2
end if
! Split up procs
# ifdef PETSC
call MPI_COMM_SPLIT(MPI_COMM_WORLD, color, 0, new_comm, mpi_err)
# endif
! assign to PETSc
# ifdef PETSC
PETSC_COMM_WORLD = new_comm
! Initialize PETSc on all procs
call PetscInitialize(PETSC_NULL_CHARACTER, mpi_err)
# endif
! Initialize timers
call time_cmfd % reset()
call time_cmfdbuild % reset()
call time_cmfdsolve % reset()
! Allocate cmfd object
call allocate_cmfd(cmfd, n_batches)
end subroutine configure_cmfd
!===============================================================================
@ -32,18 +62,17 @@ contains
subroutine read_cmfd_xml()
use error, only: fatal_error
use global
use error, only: fatal_error, warning
use output, only: write_message
use string, only: lower_case
use xml_data_cmfd_t
use, intrinsic :: ISO_FORTRAN_ENV
integer :: ng
logical :: file_exists ! does cmfd.xml exist?
character(MAX_LINE_LEN) :: filename
integer :: ng ! number of energy groups
logical :: file_exists ! does cmfd.xml exist?
character(MAX_LINE_LEN) :: filename ! name of input file
! read cmfd infput file
! Read cmfd input file
filename = trim(path_input) // "cmfd.xml"
inquire(FILE=filename, EXIST=file_exists)
if (.not. file_exists) then
@ -55,79 +84,78 @@ contains
return
else
! tell user
! Tell user
message = "Reading CMFD XML file..."
call write_message(5)
end if
! parse cmfd.xml file
! Parse cmfd.xml file
call read_xml_file_cmfd_t(filename)
! set spatial dimensions in cmfd object
cmfd % indices(1:3) = mesh_ % dimension(1:3) ! sets spatial dimensions
! Set spatial dimensions in cmfd object (structed Cartesian mesh)
cmfd % indices(1:3) = mesh_ % dimension(1:3)
! get number of energy groups
! Get number of energy groups or set to 1 group default
if (associated(mesh_ % energy)) then
ng = size(mesh_ % energy)
if(.not.allocated(cmfd%egrid)) allocate(cmfd%egrid(ng))
cmfd%egrid = mesh_ % energy
if(.not.allocated(cmfd % egrid)) allocate(cmfd % egrid(ng))
cmfd % egrid = mesh_ % energy
cmfd % indices(4) = ng - 1 ! sets energy group dimension
else
if(.not.allocated(cmfd%egrid)) allocate(cmfd%egrid(2))
cmfd%egrid = (/0.0_8,20.0_8/)
if(.not.allocated(cmfd % egrid)) allocate(cmfd % egrid(2))
cmfd % egrid = (/0.0_8,20.0_8/)
cmfd % indices(4) = 1 ! one energy group
end if
! set global albedo
! Set global albedo, these can be overwritten by coremap
if (associated(mesh_ % albedo)) then
cmfd % albedo = mesh_ % albedo
else
cmfd % albedo = (/1.0, 1.0, 1.0, 1.0, 1.0, 1.0/)
end if
! get acceleration map
! Get core map overlay for a subset mesh for CMFD
if (associated(mesh_ % map)) then
! Allocate a core map with appropriate dimensions
allocate(cmfd % coremap(cmfd % indices(1), cmfd % indices(2), &
cmfd % indices(3)))
! Check and make sure it is of correct size
if (size(mesh_ % map) /= product(cmfd % indices(1:3))) then
message = 'FATAL==>CMFD coremap not to correct dimensions'
message = 'CMFD coremap not to correct dimensions'
call fatal_error()
end if
cmfd % coremap = reshape(mesh_ % map,(cmfd % indices(1:3)))
cmfd_coremap = .true.
end if
! check for core map activation by printing note
if (cmfd_coremap .and. master) then
! Reshape core map vector into (x,y,z) array
cmfd % coremap = reshape(mesh_ % map,(cmfd % indices(1:3)))
! Indicate to cmfd that a core map overlay is active
cmfd_coremap = .true.
! Write to stdout that a core map overlay is active
message = "Core Map Overlay Activated"
call write_message()
call write_message(10)
end if
! check for normalization constant
! Get normalization constant for source (default is 1.0 from XML)
cmfd % norm = norm_
! set feedback logical
! Set feedback logical
call lower_case(feedback_)
if (feedback_ == 'true' .or. feedback_ == '1') cmfd_feedback = .true.
! set balance logical
! call lower_case(balance_)
! if (balance_ == 'true' .or. balance == '1') cmfd_balance = .true.
! Set downscatter logical
call lower_case(downscatter_)
if (downscatter_ == 'true' .or. downscatter_ == '1') &
cmfd_downscatter = .true.
! set downscatter logical
! call lower_case(downscatter_)
! if (downscatter_ == 'true' .or. downscatter == '1') &
! cmfd_downscatter = downscatter_
! set 2 group fix
call lower_case(run_2grp_)
if (run_2grp_ == 'true' .or. run_2grp_ == '1') cmfd_run_2grp = .true.
! set the solver type
! Set the solver type (default power from XML)
cmfd_solver_type = solver_(1:10)
! set monitoring
! Set convergence monitoring
call lower_case(snes_monitor_)
call lower_case(ksp_monitor_)
call lower_case(power_monitor_)
@ -138,46 +166,41 @@ contains
if (power_monitor_ == 'true' .or. power_monitor_ == '1') &
cmfd_power_monitor = .true.
! output logicals
call lower_case(write_balance_)
! Output logicals
call lower_case(write_matrices_)
! call lower_case(write_hdf5_)
if (write_balance_ == 'true' .or. write_balance_ == '1') &
cmfd_write_balance = .true.
if (write_matrices_ == 'true' .or. write_matrices_ == '1') &
cmfd_write_matrices = .true.
! if (write_hdf5_ == 'true' .or. write_hdf5_ == '1') &
! cmfd_write_hdf5 = .true.
! run an adjoint calc
! Run an adjoint calc
call lower_case(run_adjoint_)
if (run_adjoint_ == 'true' .or. run_adjoint_ == '1') &
cmfd_run_adjoint = .true.
! batch to begin cmfd
! Batch to begin cmfd (default is 1 from XML)
cmfd_begin = begin_
! tally during inactive batches
! Tally during inactive batches (by default we will always tally from 1)
call lower_case(inactive_)
if (inactive_ == 'false' .or. inactive_ == '0') cmfd_tally_on = .false.
! inactive batch flush window
cmfd_inact_flush(1) = inactive_flush_
cmfd_inact_flush(2) = num_flushes_
! Inactive batch flush window
cmfd_inact_flush(1) = inactive_flush_ ! the interval of batches
cmfd_inact_flush(2) = num_flushes_ ! number of times to do this
! last flush before active batches
! Last flush before active batches
cmfd_act_flush = active_flush_
! tolerance on keff
cmfd_keff_tol = keff_tol_
! create tally objects
call create_cmfd_tally()
! Get display
cmfd_display = display_
if (trim(cmfd_display) == 'dominance' .and. &
trim(cmfd_solver_type) /= 'power') then
message = 'Dominance Ratio only aviable with power iteration solver'
call warning()
cmfd_display = ''
end if
! set number of CMFD processors and report to user
n_procs_cmfd = n_cmfd_procs_
if (master) write(OUTPUT_UNIT,'(A,1X,I0,1X,A)') "CMFD Running on", &
n_procs_cmfd," processors."
! Create tally objects
call create_cmfd_tally()
end subroutine read_cmfd_xml
@ -193,7 +216,6 @@ contains
subroutine create_cmfd_tally()
use error, only: fatal_error, warning
use global
use mesh_header, only: StructuredMesh
use string
use tally, only: setup_active_cmfdtallies
@ -206,27 +228,27 @@ contains
integer :: ng ! number of energy groups (default 1)
integer :: n_filters ! number of filters
integer :: i_filter_mesh ! index for mesh filter
character(MAX_LINE_LEN) :: filename
character(MAX_LINE_LEN) :: filename ! name of cmfd file
type(TallyObject), pointer :: t => null()
type(StructuredMesh), pointer :: m => null()
type(TallyFilter) :: filters(N_FILTER_TYPES) ! temporary filters
! parse cmfd.xml file
! Parse cmfd.xml file
filename = trim(path_input) // "cmfd.xml"
call read_xml_file_cmfd_t(filename)
! set global variables if they are 0 (this can happen if there is no tally
! Set global variables if they are 0 (this can happen if there is no tally
! file)
if (n_meshes == 0) n_meshes = n_user_meshes + n_cmfd_meshes
! allocate mesh
! Allocate mesh
if (.not. allocated(meshes)) allocate(meshes(n_meshes))
m => meshes(n_user_meshes+1)
! set mesh id
! Set mesh id
m % id = n_user_meshes + 1
! set mesh type to rectangular
! Set mesh type to rectangular
m % type = LATTICE_RECT
! Determine number of dimensions for mesh
@ -321,20 +343,20 @@ contains
! Add mesh to dictionary
call mesh_dict % add_key(m % id, n_user_meshes + 1)
! allocate tallies
! Allocate tallies
call add_tallies("cmfd", n_cmfd_tallies)
! begin loop around tallies
! Begin loop around tallies
do i = 1, n_cmfd_tallies
! point t to tally variable
! Point t to tally variable
t => cmfd_tallies(i)
! set reset property
! Set reset property
call lower_case(reset_)
if (reset_ == 'true' .or. reset_ == '1') t % reset = .true.
! set up mesh filter
! Set up mesh filter
n_filters = 1
filters(n_filters) % type = FILTER_MESH
filters(n_filters) % n_bins = product(m % dimension)
@ -342,7 +364,7 @@ contains
filters(n_filters) % int_bins(1) = n_user_meshes + 1
t % find_filter(FILTER_MESH) = n_filters
! read and set incoming energy mesh filter
! Read and set incoming energy mesh filter
if (associated(mesh_ % energy)) then
n_filters = n_filters + 1
filters(n_filters) % type = FILTER_ENERGYIN
@ -353,40 +375,40 @@ contains
t % find_filter(FILTER_ENERGYIN) = n_filters
end if
! set number of nucilde bins
! Set number of nucilde bins
allocate(t % nuclide_bins(1))
t % nuclide_bins(1) = -1
t % n_nuclide_bins = 1
! record tally id which is equivalent to loop number
! Record tally id which is equivalent to loop number
t % id = i_cmfd_tallies + i
if (i == 1) then
! set label
! Set label
t % label = "CMFD flux, total, scatter-1"
! set tally estimator to analog
! Set tally estimator to analog
t % estimator = ESTIMATOR_ANALOG
! set tally type to volume
! Set tally type to volume
t % type = TALLY_VOLUME
! allocate and set filters
! Allocate and set filters
t % n_filters = n_filters
allocate(t % filters(n_filters))
t % filters = filters(1:n_filters)
! allocate scoring bins
! Allocate scoring bins
allocate(t % score_bins(3))
t % n_score_bins = 3
t % n_user_score_bins = 3
! allocate scattering order data
! Allocate scattering order data
allocate(t % scatt_order(3))
t % scatt_order = 0
! set macro_bins
! Set macro_bins
t % score_bins(1) = SCORE_FLUX
t % score_bins(2) = SCORE_TOTAL
t % score_bins(3) = SCORE_SCATTER_N
@ -394,16 +416,16 @@ contains
else if (i == 2) then
! set label
! Set label
t % label = "CMFD neutron production"
! set tally estimator to analog
! Set tally estimator to analog
t % estimator = ESTIMATOR_ANALOG
! set tally type to volume
! Set tally type to volume
t % type = TALLY_VOLUME
! read and set outgoing energy mesh filter
! Read and set outgoing energy mesh filter
if (associated(mesh_ % energy)) then
n_filters = n_filters + 1
filters(n_filters) % type = FILTER_ENERGYOUT
@ -414,34 +436,34 @@ contains
t % find_filter(FILTER_ENERGYOUT) = n_filters
end if
! allocate and set filters
! Allocate and set filters
t % n_filters = n_filters
allocate(t % filters(n_filters))
t % filters = filters(1:n_filters)
! deallocate filters bins array
! Deallocate filters bins array
if (associated(mesh_ % energy)) &
deallocate(filters(n_filters) % real_bins)
! allocate macro reactions
! Allocate macro reactions
allocate(t % score_bins(2))
t % n_score_bins = 2
t % n_user_score_bins = 2
! allocate scattering order data
! Allocate scattering order data
allocate(t % scatt_order(2))
t % scatt_order = 0
! set macro_bins
! Set macro_bins
t % score_bins(1) = SCORE_NU_SCATTER
t % score_bins(2) = SCORE_NU_FISSION
else if (i == 3) then
! set label
! Set label
t % label = "CMFD surface currents"
! set tally estimator to analog
! Set tally estimator to analog
t % estimator = ESTIMATOR_ANALOG
! Add extra filter for surface
@ -458,40 +480,41 @@ contains
end if
t % find_filter(FILTER_SURFACE) = n_filters
! allocate and set filters
! Allocate and set filters
t % n_filters = n_filters
allocate(t % filters(n_filters))
t % filters = filters(1:n_filters)
! deallocate filters bins array
! Deallocate filters bins array
deallocate(filters(n_filters) % int_bins)
! allocate macro reactions
! Allocate macro reactions
allocate(t % score_bins(1))
t % n_score_bins = 1
t % n_user_score_bins = 1
! allocate scattering order data
! Allocate scattering order data
allocate(t % scatt_order(1))
t % scatt_order = 0
! set macro bins
! Set macro bins
t % score_bins(1) = SCORE_CURRENT
t % type = TALLY_SURFACE_CURRENT
! we need to increase the dimension by one since we also need
! We need to increase the dimension by one since we also need
! currents coming into and out of the boundary mesh cells.
i_filter_mesh = t % find_filter(FILTER_MESH)
t % filters(i_filter_mesh) % n_bins = product(m % dimension + 1)
end if
! deallocate filter bins
! Deallocate filter bins
deallocate(filters(1) % int_bins)
if (associated(mesh_ % energy)) deallocate(filters(2) % real_bins)
end do
! Put cmfd tallies into active tally array and turn tallies on
call setup_active_cmfdtallies()
tallies_on = .true.

View file

@ -1,357 +0,0 @@
module cmfd_jacobian_operator
# ifdef PETSC
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
private
public :: init_J_operator, build_jacobian_matrix, destroy_J_operator
# include <finclude/petsc.h90>
integer :: nx ! maximum number of x cells
integer :: ny ! maximum number of y cells
integer :: nz ! maximum number of z cells
integer :: ng ! maximum number of groups
integer :: ierr ! petsc error code
type, public :: jacobian_operator
Mat :: J ! petsc matrix for neutronic prod operator
integer :: n ! dimensions of matrix
integer :: nnz ! max number of nonzeros
integer :: localn ! local size on proc
integer, allocatable :: d_nnz(:) ! vector of diagonal preallocation
integer, allocatable :: o_nnz(:) ! vector of off-diagonal preallocation
end type jacobian_operator
type, public :: operators
type(loss_operator) :: loss
type(prod_operator) :: prod
end type operators
contains
!==============================================================================
! INIT_J_OPERATOR
!==============================================================================
subroutine init_J_operator(this,ctx)
type(jacobian_operator) :: this
type(operators) :: ctx
! get indices
call get_J_indices(this)
! get preallocation
call preallocate_jacobian_matrix(this,ctx)
! set up M operator
call MatCreateAIJ(PETSC_COMM_WORLD, this%localn, this%localn, PETSC_DECIDE,&
PETSC_DECIDE, PETSC_NULL_INTEGER, this%d_nnz, PETSC_NULL_INTEGER, &
this%o_nnz, this%J,ierr)
call MatSetOption(this%J, MAT_NEW_NONZERO_LOCATIONS, PETSC_TRUE, ierr)
call MatSetOption(this%J, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE, ierr)
end subroutine init_J_operator
!==============================================================================
! GET_J_INDICES
!==============================================================================
subroutine get_J_indices(this)
use global, only: cmfd, cmfd_coremap
type(jacobian_operator) :: this
! get maximum number of cells in each direction
nx = cmfd%indices(1)
ny = cmfd%indices(2)
nz = cmfd%indices(3)
ng = cmfd%indices(4)
! get number of nonzeros
this%nnz = 7 + ng - 1
! calculate dimensions of matrix
this%n = nx*ny*nz*ng
! calculate dimensions of matrix
if (cmfd_coremap) then
this%n = cmfd % mat_dim * ng
else
this%n = nx*ny*nz*ng
end if
! add 1 for eigenvalue row
! this%n = this%n + 1
end subroutine get_J_indices
!===============================================================================
! PREALLOCATE_JACOBIAN_MATRIX
!===============================================================================
subroutine preallocate_jacobian_matrix(this,ctx)
use global, only: cmfd, n_procs_cmfd, rank
type(jacobian_operator) :: this
type(operators) :: ctx
integer :: n ! the extent of the matrix
integer :: row_start ! index of local starting row
integer :: row_end ! index of local final row
! get local problem size
n = this%n
! determine local size, divide evenly between all other procs
this%localn = n/(n_procs_cmfd)
! add 1 more if less proc id is less than mod
if (rank < mod(n,n_procs_cmfd)) this%localn = this%localn + 1
! add another 1 on last proc
if (rank == n_procs_cmfd - 1) this%localn = this%localn + 1
! determine local starting row
row_start = 0
if (rank < mod(n,n_procs_cmfd)) then
row_start = rank*(n/n_procs_cmfd+1)
else
row_start = min(mod(n,n_procs_cmfd)*(n/n_procs_cmfd+1) + (rank - &
mod(n,n_procs_cmfd))*(n/n_procs_cmfd),n)
end if
! determine local final row
row_end = row_start + this%localn - 1
! allocate counters
if (.not. allocated(this%d_nnz)) allocate(this%d_nnz(row_start:row_end))
if (.not. allocated(this%o_nnz)) allocate(this%o_nnz(row_start:row_end))
this % d_nnz = 0
this % o_nnz = 0
! start with pattern from loss matrix
if (rank == n_procs_cmfd - 1) then
this % d_nnz(row_start:row_end-1) = ctx%loss%d_nnz
this % o_nnz(row_start:row_end-1) = ctx%loss%o_nnz
else
this % d_nnz = ctx%loss%d_nnz
this % o_nnz = ctx%loss%o_nnz
end if
! append -F*phi term for last processor will take care of 1 for lambda
if (rank == n_procs_cmfd - 1) then
this%d_nnz = this%d_nnz + 1
else
this%o_nnz = this%o_nnz + 1
end if
! do last row which has all filled (already did lower left corner above)
if (rank == n_procs_cmfd - 1) then
this % d_nnz(row_end) = this % d_nnz(row_end) + (row_end - row_start)
this % o_nnz(row_end) = this % o_nnz(row_end) + (this%n - (row_end - &
row_start))
end if
end subroutine preallocate_jacobian_matrix
!===============================================================================
! BUILD_JACOBIAN_MATRIX creates the matrix representing loss of neutrons
!===============================================================================
subroutine build_jacobian_matrix(snes,x,jac,jac_prec,flag,ctx,ierr)
use constants, only: ZERO, ONE
use global, only: n_procs_cmfd, cmfd_write_matrices, rank
SNES :: snes ! the snes context
Vec :: x ! the solution vector
Mat :: jac ! the jacobian matrix
Mat :: jac_prec ! the jacobian preconditioner
MatStructure :: flag ! not used
type(operators) :: ctx ! not used
integer :: ierr ! petsc error flag
Vec :: phi ! flux vector
Vec :: source ! source vector
integer :: n ! problem size
integer :: k ! implied do loop counter
integer :: ncols ! number of nonzeros in cols
integer :: irow ! row counter
integer :: row_start! starting local row on process
integer :: row_end ! ending local row on process
integer, allocatable :: dims(:) ! vec of starting and ending rows
integer, allocatable :: dims1(:) ! vec of sizes on each proc
integer, allocatable :: cols(:) ! vector of column numbers
real(8) :: lambda ! eigenvalue
real(8), pointer :: xptr(:) ! pointer to solution vector
real(8), pointer :: sptr(:) ! pointer to source vector
real(8), allocatable :: vals(:) ! vector of row values
real(8), allocatable :: phi_tmp(:) ! temp buffer for flux
! create operators
call build_loss_matrix(ctx%loss)
call build_prod_matrix(ctx%prod)
! get problem size
n = ctx%loss%n
! get local size on each processor
call MatGetOwnershipRange(jac_prec, row_start, row_end, ierr)
! allocate cols and initialize to zero
if (.not. allocated(cols)) allocate(cols(&
maxval(ctx%loss%d_nnz + ctx%loss%o_nnz)))
if (.not. allocated(vals)) allocate(vals(&
maxval(ctx%loss%d_nnz + ctx%loss%o_nnz)))
cols = 0
vals = ZERO
! get pointers to residual vector
call VecGetArrayF90(x, xptr, ierr)
! create petsc vector for flux
call VecCreateMPI(PETSC_COMM_WORLD, ctx%loss%localn, PETSC_DECIDE, phi, ierr)
! extract flux and eigenvalue
call VecPlaceArray(phi, xptr, ierr)
if (rank == n_procs_cmfd - 1) lambda = xptr(size(xptr))
call MPI_BCAST(lambda, 1, MPI_REAL8, n_procs_cmfd-1, PETSC_COMM_WORLD, ierr)
! compute math (M-lambda*F) M is overwritten here
call MatAXPY(ctx%loss%M, -lambda, ctx%prod%F, &
DIFFERENT_NONZERO_PATTERN, ierr)
! create tmp petsc vector for source
call VecCreateMPI(PETSC_COMM_WORLD, ctx%loss%localn, PETSC_DECIDE, &
source, ierr)
! perform math (-F*phi --> source)
call MatMult(ctx%prod%F, phi, source, ierr)
call VecScale(source, -ONE, ierr)
! get pointer to source
call VecGetArrayF90(source, sptr, ierr)
! begin loop to insert things into matrix
do irow = row_start, row_end - 1
! don't do last row
if (irow == n) cycle
! get row of matrix
call MatGetRow(ctx%loss%M, irow, ncols, cols, vals, ierr)
! set that row to Jacobian matrix
call MatSetValues(jac_prec, 1, (/irow/), ncols, cols(1:ncols), vals, &
INSERT_VALUES, ierr)
! restore the row
call MatRestoreRow(ctx%loss%M, irow, ncols, cols, vals, ierr)
! insert source value
call MatSetValue(jac_prec, irow, n, sptr(irow-row_start+1), &
INSERT_VALUES, ierr)
end do
! allocate space for flux vector buffer
if (rank == n_procs_cmfd - 1) then
if (.not. allocated(phi_tmp)) allocate(phi_tmp(0:n-1))
end if
! get size on each proc
if (.not. allocated(dims)) allocate(dims(0:n_procs_cmfd))
if (.not. allocated(dims1)) allocate(dims1(0:n_procs_cmfd-1))
call VecGetOwnershipRanges(phi, dims, ierr)
do k = 0, n_procs_cmfd-1
dims1(k) = dims(k+1) - dims(k)
end do
! gather data on all procs (will truncate xptr if needed for last proc)
call MPI_GATHERV(xptr, dims1(rank), MPI_REAL8, phi_tmp, dims1, &
dims(0:n_procs_cmfd-1), MPI_REAL8, n_procs_cmfd-1, &
PETSC_COMM_WORLD, ierr)
! set values in last row of matrix
if (rank == n_procs_cmfd - 1) then
phi_tmp = -phi_tmp ! negate the transpose
call MatSetValues(jac_prec, 1, (/n/), n, (/(k,k=0,n-1)/), phi_tmp, &
INSERT_VALUES, ierr)
call MatSetValue(jac_prec, n, n, ONE, INSERT_VALUES, ierr)
end if
! assemble matrix
call MatAssemblyBegin(jac_prec, MAT_FINAL_ASSEMBLY, ierr)
call MatAssemblyEnd(jac_prec, MAT_FINAL_ASSEMBLY, ierr)
call MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY, ierr)
call MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY, ierr)
! reset all vectors
call VecResetArray(phi,ierr)
! restore all vectors
call VecRestoreArrayF90(x, xptr, ierr)
call VecRestoreArrayF90(source, sptr, ierr)
! destroy all temporary objects
call VecDestroy(phi, ierr)
call VecDestroy(source, ierr)
! deallocate all temporary space
if (allocated(cols)) deallocate(cols)
if (allocated(vals)) deallocate(vals)
if (allocated(phi_tmp)) deallocate(phi_tmp)
if (allocated(dims)) deallocate(dims)
if (allocated(dims1)) deallocate(dims1)
! print jacobian out
if (cmfd_write_matrices) call print_J_operator(jac_prec)
end subroutine build_jacobian_matrix
!===============================================================================
! PRINT_J_OPERATOR
!===============================================================================
subroutine print_J_operator(jac)
Mat :: jac
PetscViewer :: viewer
! write out matrix in binary file (debugging)
call PetscViewerBinaryOpen(PETSC_COMM_WORLD, 'jacobian.bin', &
FILE_MODE_WRITE, viewer, ierr)
call MatView(jac, viewer, ierr)
call PetscViewerDestroy(viewer, ierr)
end subroutine print_J_operator
!==============================================================================
! DESTROY_J_OPERATOR
!==============================================================================
subroutine destroy_J_operator(this)
type(jacobian_operator) :: this
! deallocate matrix
call MatDestroy(this%J,ierr)
! deallocate other parameters
if (allocated(this%d_nnz)) deallocate(this%d_nnz)
if (allocated(this%o_nnz)) deallocate(this%o_nnz)
end subroutine destroy_J_operator
# endif
end module cmfd_jacobian_operator

398
src/cmfd_jfnk_solver.F90 Normal file
View file

@ -0,0 +1,398 @@
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
type(JFNKSolver) :: jfnk ! JFNK solver object
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
jfnk_data % res_proc_ptr => compute_nonlinear_residual
jfnk_data % jac_proc_ptr => build_jacobian_matrix
#ifdef PETSC
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()
call loss % setup_petsc()
call prod % setup_petsc()
! 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
call resvec % setup_petsc()
call xvec % setup_petsc()
! Build jacobian from initial guess
call build_jacobian_matrix(xvec)
! Set up Jacobian for PETSc
call jac_prec % setup_petsc()
! 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
!===============================================================================
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
call res % write_petsc_binary(filename)
! Write out solution vector
if (adjoint_calc) then
filename = 'adj_x.bin'
else
filename = 'x.bin'
end if
call x % write_petsc_binary(filename)
end if
end subroutine compute_nonlinear_residual
!===============================================================================
! COMPUTE_ADJOINT calculates mathematical adjoint by taking transposes
!===============================================================================
subroutine compute_adjoint()
use global, only: cmfd_write_matrices
! Transpose matrices
call loss % transpose()
call prod % transpose()
! Write out matrix in binary file (debugging)
if (cmfd_write_matrices) then
call loss % write_petsc_binary('adj_lossmat.bin')
call loss % write_petsc_binary('adj_prodmat.bin')
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,103 +1,86 @@
module cmfd_loss_operator
# ifdef PETSC
use constants, only: CMFD_NOACCEL, ZERO
use global, only: cmfd, cmfd_coremap
use matrix_header, only: Matrix
implicit none
private
public :: init_M_operator, build_loss_matrix, destroy_M_operator
# include <finclude/petsc.h90>
integer :: nx ! maximum number of x cells
integer :: ny ! maximum number of y cells
integer :: nz ! maximum number of z cells
integer :: ng ! maximum number of groups
integer :: ierr ! petsc error code
type, public :: loss_operator
Mat :: M ! petsc matrix for neutronic loss operator
integer :: n ! dimensions of matrix
integer :: nnz ! max number of nonzeros
integer :: localn ! local size on proc
integer, allocatable :: d_nnz(:) ! vector of diagonal preallocation
integer, allocatable :: o_nnz(:) ! vector of off-diagonal preallocation
end type loss_operator
logical :: adjoint_calc = .false. ! adjoint calculation
public :: init_loss_matrix, build_loss_matrix
contains
!===============================================================================
! INIT_M_OPERATOR
! INIT_LOSS_MATRIX preallocates loss matrix and initializes it
!===============================================================================
subroutine init_M_operator(this)
subroutine init_loss_matrix(loss_matrix)
type(loss_operator) :: this
type(Matrix), intent(inout) :: loss_matrix ! cmfd loss matrix
! get indices
call get_M_indices(this)
integer :: nx ! maximum number of x cells
integer :: ny ! maximum number of y cells
integer :: nz ! maximum number of z cells
integer :: ng ! maximum number of groups
integer :: n ! total length of matrix
integer :: nnz ! number of nonzeros in matrix
integer :: n_i ! number of interior cells
integer :: n_c ! number of corner cells
integer :: n_s ! number side cells
integer :: nz_c ! number of non-zero corner cells
integer :: nz_s ! number of non-zero side cells
integer :: nz_i ! number of non-zero interior cells
! get preallocation
call preallocate_loss_matrix(this)
! set up M operator
call MatCreateAIJ(PETSC_COMM_WORLD, this%localn, this%localn, PETSC_DECIDE,&
PETSC_DECIDE, PETSC_NULL_INTEGER, this%d_nnz, PETSC_NULL_INTEGER, &
this%o_nnz, this%M,ierr)
call MatSetOption(this%M, MAT_NEW_NONZERO_LOCATIONS, PETSC_TRUE, ierr)
call MatSetOption(this%M, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE, ierr)
end subroutine init_M_operator
!===============================================================================
! GET_M_INDICES
!===============================================================================
subroutine get_M_indices(this)
use global, only: cmfd, cmfd_coremap
type(loss_operator) :: this
! get maximum number of cells in each direction
! Get maximum number of cells in each direction
nx = cmfd%indices(1)
ny = cmfd%indices(2)
nz = cmfd%indices(3)
ng = cmfd%indices(4)
! get number of nonzeros
this%nnz = 7 + ng - 1
! calculate dimensions of matrix
! Calculate dimensions of matrix
if (cmfd_coremap) then
this%n = cmfd % mat_dim * ng
n = cmfd % mat_dim * ng
else
this%n = nx*ny*nz*ng
n = nx*ny*nz*ng
end if
end subroutine get_M_indices
! Calculate number of nonzeros, if core map -> need to determine manually
if (cmfd_coremap) then
nnz = preallocate_loss_matrix(nx, ny, nz, ng, n)
else ! structured Cartesian grid
n_c = 4 ! define # of corners
n_s = 2*(nx + ny) - 8 ! define # of sides
n_i = nx*ny - (n_c + n_S) ! define # of interiors
nz_c = ng*n_c*(3 + ng - 1) ! define # nonzero corners
nz_s = ng*n_s*(4 + ng - 1) ! define # nonzero sides
nz_i = ng*n_i*(5 + ng - 1) ! define # nonzero interiors
nnz = nz_c + nz_s + nz_i
end if
! Configure loss matrix
call loss_matrix % create(n, nnz)
end subroutine init_loss_matrix
!===============================================================================
! PREALLOCATE_LOSS_MATRIX
! PREALLOCATE_LOSS_MATRIX manually preallocates the loss matrix
!===============================================================================
subroutine preallocate_loss_matrix(this)
function preallocate_loss_matrix(nx, ny, nz, ng, n) result(nnz)
use constants, only: CMFD_NOACCEL
use global, only: cmfd, cmfd_coremap
integer, intent(in) :: nx ! maximum number of x cells
integer, intent(in) :: ny ! maximum number of y cells
integer, intent(in) :: nz ! maximum number of z cells
integer, intent(in) :: ng ! maximum number of groups
integer, intent(in) :: n ! total length of matrix
integer :: nnz ! number of nonzeros
type(loss_operator) :: this
integer :: rank ! rank of processor
integer :: sizen ! number of procs
integer :: i ! iteration counter for x
integer :: j ! iteration counter for y
integer :: k ! iteration counter for z
integer :: g ! iteration counter for groups
integer :: l ! iteration counter for leakages
integer :: h ! energy group when doing scattering
integer :: n ! the extent of the matrix
integer :: irow ! row counter
integer :: bound(6) ! vector for comparing when looking for bound
integer :: xyz_idx ! index for determining if x,y or z leakage
@ -105,113 +88,68 @@ contains
integer :: neig_idx(3) ! spatial indices of neighbour
integer :: nxyz(3,2) ! single vector containing bound. locations
integer :: shift_idx ! parameter to shift index by +1 or -1
integer :: row_start ! index of local starting row
integer :: row_end ! index of local final row
integer :: neig_mat_idx ! matrix index of neighbor cell
integer :: scatt_mat_idx ! matrix index for h-->g scattering terms
! initialize size and rank
sizen = 0
rank = 0
! Reset number of nonzeros to 0
nnz = 0
! get rank and max rank of procs
call MPI_COMM_RANK(PETSC_COMM_WORLD, rank, ierr)
call MPI_COMM_SIZE(PETSC_COMM_WORLD, sizen, ierr)
! get local problem size
n = this%n
! determine local size, divide evenly between all other procs
this%localn = n/(sizen)
! add 1 more if less proc id is less than mod
if (rank < mod(n,sizen)) this%localn = this%localn + 1
! determine local starting row
row_start = 0
if (rank < mod(n,sizen)) then
row_start = rank*(n/sizen+1)
else
row_start = min(mod(n,sizen)*(n/sizen+1) + &
(rank - mod(n,sizen))*(n/sizen),n)
end if
! determine local final row
row_end = row_start + this%localn - 1
! allocate counters
if (.not. allocated(this%d_nnz)) allocate(this%d_nnz(row_start:row_end))
if (.not. allocated(this%o_nnz)) allocate(this%o_nnz(row_start:row_end))
this % d_nnz = 0
this % o_nnz = 0
! create single vector of these indices for boundary calculation
! Create single vector of these indices for boundary calculation
nxyz(1,:) = (/1,nx/)
nxyz(2,:) = (/1,ny/)
nxyz(3,:) = (/1,nz/)
! begin loop around local rows
ROWS: do irow = row_start,row_end
! Begin loop around local rows
ROWS: do irow = 1, n
! initialize counters
this%d_nnz(irow) = 1 ! already add in matrix diagonal
this%o_nnz(irow) = 0
! Set a nonzero for diagonal
nnz = nnz + 1
! get location indices
call matrix_to_indices(irow, g, i, j, k)
! Get location indices
call matrix_to_indices(irow, g, i, j, k, ng, nx, ny, nz)
! create boundary vector
! Create boundary vector
bound = (/i,i,j,j,k,k/)
! begin loop over leakages
! Begin loop over leakages
LEAK: do l = 1,6
! define (x,y,z) and (-,+) indices
! Define (x,y,z) and (-,+) indices
xyz_idx = int(ceiling(real(l)/real(2))) ! x=1, y=2, z=3
dir_idx = 2 - mod(l,2) ! -=1, +=2
! calculate spatial indices of neighbor
! Calculate spatial indices of neighbor
neig_idx = (/i,j,k/) ! begin with i,j,k
shift_idx = -2*mod(l,2) +1 ! shift neig by -1 or +1
neig_idx(xyz_idx) = shift_idx + neig_idx(xyz_idx)
! check for global boundary
! Check for global boundary
if (bound(l) /= nxyz(xyz_idx,dir_idx)) then
! check for coremap
! Check for coremap
if (cmfd_coremap) then
! check for neighbor that is non-acceleartred
! Check for neighbor that is non-acceleartred
if (cmfd % coremap(neig_idx(1),neig_idx(2),neig_idx(3)) /= &
CMFD_NOACCEL) then
! get neighbor matrix index
! Get neighbor matrix index
call indices_to_matrix(g,neig_idx(1), neig_idx(2), &
neig_idx(3), neig_mat_idx)
neig_idx(3), neig_mat_idx, ng, nx, ny)
! record nonzero
if (((neig_mat_idx-1) >= row_start) .and. &
((neig_mat_idx-1) <= row_end)) then
this%d_nnz(irow) = this%d_nnz(irow) + 1
else
this%o_nnz(irow) = this%o_nnz(irow) + 1
end if
! Record nonzero
nnz = nnz + 1
end if
else
! get neighbor matrix index
! Get neighbor matrix index
call indices_to_matrix(g, neig_idx(1), neig_idx(2), neig_idx(3), &
neig_mat_idx)
neig_mat_idx, ng, nx, ny)
! record nonzero
if (((neig_mat_idx-1) >= row_start) .and. &
((neig_mat_idx-1) <= row_end)) then
this%d_nnz(irow) = this%d_nnz(irow) + 1
else
this%o_nnz(irow) = this%o_nnz(irow) + 1
end if
! Record nonzero
nnz = nnz + 1
end if
@ -219,214 +157,206 @@ contains
end do LEAK
! begin loop over off diagonal in-scattering
! Begin loop over off diagonal in-scattering
SCATTR: do h = 1, ng
! cycle though if h=g
! Cycle though if h=g, it was already banked in removal xs
if (h == g) cycle
! get neighbor matrix index
call indices_to_matrix(h, i, j, k, scatt_mat_idx)
! Get neighbor matrix index
call indices_to_matrix(h, i, j, k, scatt_mat_idx, ng, nx, ny)
! record nonzero
if (((scatt_mat_idx-1) >= row_start) .and. &
((scatt_mat_idx-1) <= row_end)) then
this%d_nnz(irow) = this%d_nnz(irow) + 1
else
this%o_nnz(irow) = this%o_nnz(irow) + 1
end if
! Record nonzero
nnz = nnz + 1
end do SCATTR
end do ROWS
end subroutine preallocate_loss_matrix
end function preallocate_loss_matrix
!===============================================================================
! BUILD_LOSS_MATRIX creates the matrix representing loss of neutrons
!===============================================================================
subroutine build_loss_matrix(this, adjoint)
subroutine build_loss_matrix(loss_matrix, adjoint)
use constants, only: CMFD_NOACCEL, ZERO
use global, only: cmfd, cmfd_coremap, cmfd_write_matrices
type(Matrix), intent(inout) :: loss_matrix ! cmfd loss matrix
logical, intent(in), optional :: adjoint ! set up the adjoint
type(loss_operator) :: this
logical, optional :: adjoint ! set up the adjoint
integer :: nxyz(3,2) ! single vector containing bound. locations
integer :: i ! iteration counter for x
integer :: j ! iteration counter for y
integer :: k ! iteration counter for z
integer :: g ! iteration counter for groups
integer :: l ! iteration counter for leakages
integer :: h ! energy group when doing scattering
integer :: nx ! maximum number of x cells
integer :: ny ! maximum number of y cells
integer :: nz ! maximum number of z cells
integer :: ng ! maximum number of groups
integer :: neig_mat_idx ! matrix index of neighbor cell
integer :: scatt_mat_idx ! matrix index for h-->g scattering terms
integer :: bound(6) ! vector for comparing when looking for bound
integer :: xyz_idx ! index for determining if x,y or z leakage
integer :: dir_idx ! index for determining - or + face of cell
integer :: neig_idx(3) ! spatial indices of neighbour
integer :: shift_idx ! parameter to shift index by +1 or -1
integer :: irow ! iteration counter over row
logical :: adjoint_calc ! is this a physical adjoint calculation?
real(8) :: totxs ! total macro cross section
real(8) :: scattxsgg ! scattering macro cross section g-->g
real(8) :: scattxshg ! scattering macro cross section h-->g
real(8) :: dtilde(6) ! finite difference coupling parameter
real(8) :: dhat(6) ! nonlinear coupling parameter
real(8) :: hxyz(3) ! cell lengths in each direction
real(8) :: jn ! direction dependent leakage coeff to neig
real(8) :: jo(6) ! leakage coeff in front of cell flux
real(8) :: jnet ! net leakage from jo
real(8) :: val ! temporary variable before saving to
integer :: nxyz(3,2) ! single vector containing bound. locations
integer :: i ! iteration counter for x
integer :: j ! iteration counter for y
integer :: k ! iteration counter for z
integer :: g ! iteration counter for groups
integer :: l ! iteration counter for leakages
integer :: h ! energy group when doing scattering
integer :: neig_mat_idx ! matrix index of neighbor cell
integer :: scatt_mat_idx ! matrix index for h-->g scattering terms
integer :: bound(6) ! vector for comparing when looking for bound
integer :: xyz_idx ! index for determining if x,y or z leakage
integer :: dir_idx ! index for determining - or + face of cell
integer :: neig_idx(3) ! spatial indices of neighbour
integer :: shift_idx ! parameter to shift index by +1 or -1
integer :: row_start ! the first local row on the processor
integer :: row_finish ! the last local row on the processor
integer :: irow ! iteration counter over row
real(8) :: totxs ! total macro cross section
real(8) :: scattxsgg ! scattering macro cross section g-->g
real(8) :: scattxshg ! scattering macro cross section h-->g
real(8) :: dtilde(6) ! finite difference coupling parameter
real(8) :: dhat(6) ! nonlinear coupling parameter
real(8) :: hxyz(3) ! cell lengths in each direction
real(8) :: jn ! direction dependent leakage coeff to neig
real(8) :: jo(6) ! leakage coeff in front of cell flux
real(8) :: jnet ! net leakage from jo
real(8) :: val ! temporary variable before saving to
! check for adjoint
! Check for adjoint
adjoint_calc = .false.
if (present(adjoint)) adjoint_calc = adjoint
! create single vector of these indices for boundary calculation
! Get maximum number of cells in each direction
nx = cmfd%indices(1)
ny = cmfd%indices(2)
nz = cmfd%indices(3)
ng = cmfd%indices(4)
! Create single vector of these indices for boundary calculation
nxyz(1,:) = (/1,nx/)
nxyz(2,:) = (/1,ny/)
nxyz(3,:) = (/1,nz/)
! initialize row start and finish
row_start = 0
row_finish = 0
! Begin iteration loops
ROWS: do irow = 1, loss_matrix % n
! get row bounds for this processor
call MatGetOwnershipRange(this%M, row_start, row_finish, ierr)
! Set up a new row in matrix
call loss_matrix % new_row()
! begin iteration loops
ROWS: do irow = row_start, row_finish-1
! Get indices for that row
call matrix_to_indices(irow, g, i, j, k, ng, nx, ny, nz)
! get indices for that row
call matrix_to_indices(irow, g, i, j, k)
! retrieve cell data
! Retrieve cell data
totxs = cmfd%totalxs(g,i,j,k)
scattxsgg = cmfd%scattxs(g,g,i,j,k)
dtilde = cmfd%dtilde(:,g,i,j,k)
hxyz = cmfd%hxyz(:,i,j,k)
! check and get dhat
! Check and get dhat
if (allocated(cmfd%dhat)) then
dhat = cmfd%dhat(:,g,i,j,k)
else
dhat = ZERO
end if
! create boundary vector
! Create boundary vector
bound = (/i,i,j,j,k,k/)
! begin loop over leakages
! Begin loop over leakages
! 1=-x, 2=+x, 3=-y, 4=+y, 5=-z, 6=+z
LEAK: do l = 1,6
! define (x,y,z) and (-,+) indices
! Define (x,y,z) and (-,+) indices
xyz_idx = int(ceiling(real(l)/real(2))) ! x=1, y=2, z=3
dir_idx = 2 - mod(l,2) ! -=1, +=2
! calculate spatial indices of neighbor
! Calculate spatial indices of neighbor
neig_idx = (/i,j,k/) ! begin with i,j,k
shift_idx = -2*mod(l,2) +1 ! shift neig by -1 or +1
neig_idx(xyz_idx) = shift_idx + neig_idx(xyz_idx)
! check for global boundary
! Check for global boundary
if (bound(l) /= nxyz(xyz_idx,dir_idx)) then
! check for core map
! Check for core map
if (cmfd_coremap) then
! check that neighbor is not reflector
! Check that neighbor is not reflector
if (cmfd % coremap(neig_idx(1),neig_idx(2),neig_idx(3)) /= &
CMFD_NOACCEL) then
! compute leakage coefficient for neighbor
! Compute leakage coefficient for neighbor
jn = -dtilde(l) + shift_idx*dhat(l)
! get neighbor matrix index
! Get neighbor matrix index
call indices_to_matrix(g, neig_idx(1), neig_idx(2), neig_idx(3), &
neig_mat_idx)
neig_mat_idx, ng, nx, ny)
! compute value and record to bank
! Compute value and record to bank
val = jn/hxyz(xyz_idx)
! record value in matrix
call MatSetValue(this%M, irow, neig_mat_idx-1, val, &
INSERT_VALUES, ierr)
! Record value in matrix
call loss_matrix % add_value(neig_mat_idx, val)
end if
else
! compute leakage coefficient for neighbor
! Compute leakage coefficient for neighbor
jn = -dtilde(l) + shift_idx*dhat(l)
! get neighbor matrix index
! Get neighbor matrix index
call indices_to_matrix(g, neig_idx(1), neig_idx(2), neig_idx(3), &
neig_mat_idx)
neig_mat_idx, ng, nx, ny)
! compute value and record to bank
! Compute value and record to bank
val = jn/hxyz(xyz_idx)
! record value in matrix
call MatSetValue(this%M, irow, neig_mat_idx-1, val, &
INSERT_VALUES, ierr)
! Record value in matrix
call loss_matrix % add_value(neig_mat_idx, val)
end if
end if
! compute leakage coefficient for target
! Compute leakage coefficient for target
jo(l) = shift_idx*dtilde(l) + dhat(l)
end do LEAK
! calate net leakage coefficient for target
! Calculate net leakage coefficient for target
jnet = (jo(2) - jo(1))/hxyz(1) + (jo(4) - jo(3))/hxyz(2) + &
(jo(6) - jo(5))/hxyz(3)
! calculate loss of neutrons
! Calculate loss of neutrons
val = jnet + totxs - scattxsgg
! record diagonal term
call MatSetValue(this%M, irow, irow, val, INSERT_VALUES, ierr)
! Record diagonal term
call loss_matrix % add_value(irow, val)
! begin loop over off diagonal in-scattering
! Begin loop over off diagonal in-scattering
SCATTR: do h = 1, ng
! cycle though if h=g
! Cycle though if h=g, value already banked in removal xs
if (h == g) cycle
! get neighbor matrix index
call indices_to_matrix(h, i, j, k, scatt_mat_idx)
! Get neighbor matrix index
call indices_to_matrix(h, i, j, k, scatt_mat_idx, ng, nx, ny)
! get scattering macro xs
scattxshg = cmfd%scattxs(h, g, i, j, k)
! Check for adjoint
if (adjoint_calc) then
! Get scattering macro xs, transposed!
scattxshg = cmfd%scattxs(g, h, i, j, k)
else
! Get scattering macro xs
scattxshg = cmfd%scattxs(h, g, i, j, k)
end if
! record value in matrix (negate it)
! Negate the scattering xs
val = -scattxshg
! check for adjoint and bank value
if (adjoint_calc) then
call MatSetValue(this%M, scatt_mat_idx-1, irow, val, &
INSERT_VALUES, ierr)
else
call MatSetValue(this%M, irow, scatt_mat_idx-1, val, &
INSERT_VALUES, ierr)
end if
! Record value in matrix
call loss_matrix % add_value(scatt_mat_idx, val)
end do SCATTR
end do ROWS
! assemble matrix
call MatAssemblyBegin(this%M, MAT_FLUSH_ASSEMBLY, ierr)
call MatAssemblyEnd(this%M, MAT_FINAL_ASSEMBLY, ierr)
! print out operator to file
if (cmfd_write_matrices) call print_M_operator(this)
! CSR requires n+1 row
call loss_matrix % new_row()
end subroutine build_loss_matrix
@ -434,25 +364,26 @@ contains
! INDICES_TO_MATRIX takes (x,y,z,g) indices and computes location in matrix
!===============================================================================
subroutine indices_to_matrix(g, i, j, k, matidx)
subroutine indices_to_matrix(g, i, j, k, matidx, ng, nx, ny)
use global, only: cmfd, cmfd_coremap
integer, intent(out) :: matidx ! the index location in matrix
integer, intent(in) :: i ! current x index
integer, intent(in) :: j ! current y index
integer, intent(in) :: k ! current z index
integer, intent(in) :: g ! current group index
integer, intent(in) :: nx ! maximum number of x cells
integer, intent(in) :: ny ! maximum number of y cells
integer, intent(in) :: ng ! maximum number of groups
integer :: matidx ! the index location in matrix
integer :: i ! current x index
integer :: j ! current y index
integer :: k ! current z index
integer :: g ! current group index
! check if coremap is used
! Check if coremap is used
if (cmfd_coremap) then
! get idx from core map
! Get idx from core map
matidx = ng*(cmfd % coremap(i,j,k)) - (ng - g)
else
! compute index
! Compute index
matidx = g + ng*(i - 1) + ng*nx*(j - 1) + ng*nx*ny*(k - 1)
end if
@ -460,80 +391,40 @@ contains
end subroutine indices_to_matrix
!===============================================================================
! MATRIX_TO_INDICES
! MATRIX_TO_INDICES converts a matrix index to spatial and group indicies
!===============================================================================
subroutine matrix_to_indices(irow, g, i, j, k)
subroutine matrix_to_indices(irow, g, i, j, k, ng, nx, ny, nz)
use global, only: cmfd, cmfd_coremap
integer, intent(out) :: i ! iteration counter for x
integer, intent(out) :: j ! iteration counter for y
integer, intent(out) :: k ! iteration counter for z
integer, intent(out) :: g ! iteration counter for groups
integer, intent(in) :: irow ! iteration counter over row (0 reference)
integer, intent(in) :: nx ! maximum number of x cells
integer, intent(in) :: ny ! maximum number of y cells
integer, intent(in) :: nz ! maximum number of z cells
integer, intent(in) :: ng ! maximum number of groups
integer :: i ! iteration counter for x
integer :: j ! iteration counter for y
integer :: k ! iteration counter for z
integer :: g ! iteration counter for groups
integer :: irow ! iteration counter over row (0 reference)
! check for core map
! Check for core map
if (cmfd_coremap) then
! get indices from indexmap
g = mod(irow, ng) + 1
i = cmfd % indexmap(irow/ng+1,1)
j = cmfd % indexmap(irow/ng+1,2)
k = cmfd % indexmap(irow/ng+1,3)
! Get indices from indexmap
g = mod(irow-1, ng) + 1
i = cmfd % indexmap((irow-1)/ng+1,1)
j = cmfd % indexmap((irow-1)/ng+1,2)
k = cmfd % indexmap((irow-1)/ng+1,3)
else
! compute indices
g = mod(irow, ng) + 1
i = mod(irow, ng*nx)/ng + 1
j = mod(irow, ng*nx*ny)/(ng*nx)+ 1
k = mod(irow, ng*nx*ny*nz)/(ng*nx*ny) + 1
! Compute indices
g = mod(irow-1, ng) + 1
i = mod(irow-1, ng*nx)/ng + 1
j = mod(irow-1, ng*nx*ny)/(ng*nx)+ 1
k = mod(irow-1, ng*nx*ny*nz)/(ng*nx*ny) + 1
end if
end subroutine matrix_to_indices
!===============================================================================
! PRINT_M_OPERATOR
!===============================================================================
subroutine print_M_operator(this)
type(loss_operator) :: this
PetscViewer :: viewer
! write out matrix in binary file (debugging)
if (adjoint_calc) then
call PetscViewerBinaryOpen(PETSC_COMM_WORLD, 'adj_lossmat.bin', &
FILE_MODE_WRITE, viewer, ierr)
else
call PetscViewerBinaryOpen(PETSC_COMM_WORLD, 'lossmat.bin', &
FILE_MODE_WRITE, viewer, ierr)
end if
call MatView(this%M, viewer, ierr)
call PetscViewerDestroy(viewer, ierr)
end subroutine print_M_operator
!==============================================================================
! DESTROY_M_OPERATOR
!==============================================================================
subroutine destroy_M_operator(this)
type(loss_operator) :: this
! deallocate matrix
call MatDestroy(this%M, ierr)
! deallocate other parameters
if (allocated(this%d_nnz)) deallocate(this%d_nnz)
if (allocated(this%o_nnz)) deallocate(this%o_nnz)
end subroutine destroy_M_operator
# endif
end module cmfd_loss_operator

View file

@ -1,98 +0,0 @@
module cmfd_message_passing
# ifdef PETSC
! This module contains MPI routines related to the CMFD calculation
implicit none
private
public :: petsc_init_mpi, cmfd_bcast
# include <finclude/petsc.h90>
contains
!===============================================================================
! PETSC_INIT_MPI
!===============================================================================
subroutine petsc_init_mpi()
use global, only: n_procs_cmfd, time_cmfd, rank, mpi_err
integer :: new_comm ! new communicator
integer :: color ! communicator color
! assign color
if (rank < n_procs_cmfd) then
color = 1
else
color = 2
end if
! split up procs
call MPI_COMM_SPLIT(MPI_COMM_WORLD,color,0,new_comm,mpi_err)
! assign to PETSc
PETSC_COMM_WORLD = new_comm
! initialize petsc on all procs
call PetscInitialize(PETSC_NULL_CHARACTER,mpi_err)
! initialize timer
call time_cmfd % reset()
end subroutine petsc_init_mpi
!===============================================================================
! CMFD_BCAST
!===============================================================================
subroutine cmfd_bcast()
use global, only: cmfd_coremap, cmfd, mpi_err
use cmfd_header, only: allocate_cmfd
integer :: nx ! number of mesh cells in x direction
integer :: ny ! number of mesh cells in y direction
integer :: nz ! number of mesh cells in z direction
integer :: ng ! number of energy groups
! extract spatial and energy indices from object
nx = cmfd % indices(1)
ny = cmfd % indices(2)
nz = cmfd % indices(3)
ng = cmfd % indices(4)
! initialize data
if (.not. allocated(cmfd%flux)) call allocate_cmfd(cmfd)
! sync up procs
call MPI_Barrier(PETSC_COMM_WORLD, mpi_err)
! broadcast all data
call MPI_BCAST(cmfd%flux, ng*nx*ny*nz, MPI_REAL8, 0, PETSC_COMM_WORLD, mpi_err)
call MPI_BCAST(cmfd%totalxs, ng*nx*ny*nz, MPI_REAL8, 0, PETSC_COMM_WORLD, mpi_err)
call MPI_BCAST(cmfd%p1scattxs, ng*nx*ny*nz, MPI_REAL8, 0, PETSC_COMM_WORLD, mpi_err)
call MPI_BCAST(cmfd%scattxs, ng*ng*nx*ny*nz, MPI_REAL8, 0, PETSC_COMM_WORLD, mpi_err)
call MPI_BCAST(cmfd%nfissxs, ng*ng*nx*ny*nz, MPI_REAL8, 0, PETSC_COMM_WORLD, mpi_err)
call MPI_BCAST(cmfd%diffcof, ng*nx*ny*nz, MPI_REAL8, 0, PETSC_COMM_WORLD, mpi_err)
call MPI_BCAST(cmfd%dtilde, 6*ng*nx*ny*nz, MPI_REAL8, 0, PETSC_COMM_WORLD, mpi_err)
call MPI_BCAST(cmfd%dhat, 6*ng*nx*ny*nz, MPI_REAL8, 0, PETSC_COMM_WORLD, mpi_err)
call MPI_BCAST(cmfd%hxyz, 3*nx*ny*nz, MPI_REAL8, 0, PETSC_COMM_WORLD, mpi_err)
call MPI_BCAST(cmfd%current, 12*ng*nx*ny*nz, MPI_REAL8, 0, PETSC_COMM_WORLD, mpi_err)
! broadcast coremap info
if (cmfd_coremap) then
call MPI_BCAST(cmfd%coremap, nx*ny*nz, MPI_INTEGER8, 0, PETSC_COMM_WORLD, mpi_err)
call MPI_BCAST(cmfd%mat_dim, 1, MPI_INTEGER8, 0, PETSC_COMM_WORLD, mpi_err)
if (.not. allocated(cmfd % indexmap)) &
allocate(cmfd % indexmap(cmfd % mat_dim,3))
call MPI_BCAST(cmfd%indexmap, cmfd%mat_dim*3, MPI_INTEGER8, 0, PETSC_COMM_WORLD, mpi_err)
end if
end subroutine cmfd_bcast
# endif
end module cmfd_message_passing

View file

@ -1,61 +0,0 @@
module cmfd_output
#ifdef PETSC
! This modules cleans up cmfd objects and echos the results
implicit none
private
public :: finalize_cmfd
contains
!===============================================================================
! FINALIZE_CMFD
!===============================================================================
subroutine finalize_cmfd()
use global, only: cmfd, cmfd_write_balance, cmfd_write_hdf5, &
master, mpi_err
use cmfd_header, only: deallocate_cmfd
! finalize petsc
call PetscFinalize(mpi_err)
! write out final neutron balance
if (master .and. cmfd_write_balance) call write_neutron_balance()
! deallocate cmfd object
call deallocate_cmfd(cmfd)
end subroutine finalize_cmfd
!===============================================================================
! WRITE_NEUTRON_BALANCE
!===============================================================================
subroutine write_neutron_balance()
use cmfd_data, only: neutron_balance
use constants, only: CMFD_BALANCE, MAX_FILE_LEN
use global, only: path_output
character(MAX_FILE_LEN) :: filename
filename = trim(path_output) // 'neutron_balance.out'
! open file for output
open(UNIT=CMFD_BALANCE, FILE=filename, ACTION='write')
! write out the tally
call neutron_balance(CMFD_BALANCE)
! close file
close(CMFD_BALANCE)
end subroutine write_neutron_balance
#endif
end module cmfd_output

View file

@ -1,445 +1,346 @@
module cmfd_power_solver
# ifdef PETSC
! This module contains routines to execute the power iteration solver
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
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
# include <finclude/petsc.h90>
type(loss_operator) :: loss ! M loss matrix
type(prod_operator) :: prod ! F production matrix
Vec :: phi_n ! new flux eigenvector
Vec :: phi_o ! old flux eigenvector
Vec :: S_n ! new source vector
Vec :: S_o ! old source vector
KSP :: krylov ! krylov solver
KSP :: sub_krylov ! sub-ksp for bjacobi
PC :: prec ! preconditioner for krylov
PC :: sub_prec ! sub-prec for bjacobi
integer :: ierr ! error flag
integer :: nlocal
integer :: first
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
logical :: iconv ! did the problem converged
logical :: adjoint_calc = .false. ! run an adjoint calculation
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
type(GMRESSolver) :: gmres ! gmres solver
contains
!===============================================================================
! CMFD_POWER_EXECUTE
! 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, n_procs_cmfd
use global, only: cmfd_adjoint_type, time_cmfdbuild, time_cmfdsolve
real(8), optional :: k_tol ! tolerance on keff
real(8), optional :: s_tol ! tolerance on source
logical, optional :: adjoint ! adjoint calc
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.
logical :: physical_adjoint = .false. ! physical adjoint default false
! set tolerances if present
! Set tolerances if present
if (present(k_tol)) ktol = k_tol
if (present(s_tol)) stol = s_tol
! check for adjoint execution
! Check for adjoint execution
adjoint_calc = .false.
if (present(adjoint)) adjoint_calc = adjoint
! check for physical adjoint
! Check for physical adjoint
if (adjoint_calc .and. trim(cmfd_adjoint_type) == 'physical') &
physical_adjoint = .true.
! initialize solver
call init_solver()
! Start timer for build
call time_cmfdbuild % start()
! initialize matrices and vectors
call init_data()
! Initialize solver
#ifdef PETSC
call gmres % create()
#endif
! set up M loss matrix
call build_loss_matrix(loss, adjoint=physical_adjoint)
! Initialize matrices and vectors
call init_data(physical_adjoint)
! set up F production matrix
call build_prod_matrix(prod, adjoint=physical_adjoint)
! check for adjoint calculation
! Check for mathematical adjoint calculation
if (adjoint_calc .and. trim(cmfd_adjoint_type) == 'math') &
call compute_adjoint()
! set up krylov info
call KSPSetOperators(krylov, loss%M, loss%M, SAME_NONZERO_PATTERN, ierr)
! Set up krylov info
#ifdef PETSC
call gmres % set_oper(loss, loss)
#endif
! precondition matrix
call precondition_matrix()
! Stop timer for build
call time_cmfdbuild % stop()
! begin power iteration
! Begin power iteration
call time_cmfdsolve % start()
call execute_power_iter()
call time_cmfdsolve % stop()
! extract results
! Extract results
call extract_results()
! deallocate petsc objects
! Deallocate data
call finalize()
end subroutine cmfd_power_execute
!===============================================================================
! INIT_DATA allocates matrices vectors for CMFD solution
! INIT_DATA allocates matrices and vectors for CMFD solution
!===============================================================================
subroutine init_data()
subroutine init_data(adjoint)
use constants, only: ONE
use constants, only: ONE, ZERO
use global, only: cmfd_write_matrices
logical, intent(in) :: adjoint ! adjoint calcualtion
integer :: n ! problem size
real(8) :: guess ! initial guess
! set up matrices
call init_M_operator(loss)
call init_F_operator(prod)
! Set up matrices
call init_loss_matrix(loss)
call init_prod_matrix(prod)
! get problem size
n = loss%localn
! Get problem size
n = loss % n
! set up flux vectors
call VecCreateMPI(PETSC_COMM_WORLD, n, PETSC_DECIDE, phi_n, ierr)
call VecCreateMPI(PETSC_COMM_WORLD, n, PETSC_DECIDE, phi_o, ierr)
! Set up flux vectors
call phi_n % create(n)
call phi_o % create(n)
! set up source vectors
call VecCreateMPI(PETSC_COMM_WORLD, n, PETSC_DECIDE, S_n, ierr)
call VecCreateMPI(PETSC_COMM_WORLD, n, PETSC_DECIDE, S_o, ierr)
! Set up source vectors
call s_n % create(n)
call s_o % create(n)
call serr_v % create(n)
! set initial guess
! Set initial guess
guess = ONE
call VecSet(phi_n, guess, ierr)
call VecSet(phi_o, guess, ierr)
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
!===============================================================================
! INIT_SOLVER
!===============================================================================
subroutine init_solver()
real(8) :: rtol ! relative tolerance
real(8) :: atol ! absolute tolerance
! set tolerance
rtol = 1.0e-10_8
atol = 1.0e-10_8
! set up krylov solver
call KSPCreate(PETSC_COMM_WORLD, krylov, ierr)
call KSPSetTolerances(krylov, rtol, atol, PETSC_DEFAULT_DOUBLE_PRECISION, &
PETSC_DEFAULT_INTEGER, ierr)
call KSPSetType(krylov, KSPGMRES, ierr)
call KSPSetInitialGuessNonzero(krylov, PETSC_TRUE, ierr)
end subroutine init_solver
!===============================================================================
! PRECONDITION_MATRIX
!===============================================================================
subroutine precondition_matrix()
use global, only: cmfd_power_monitor, cmfd_solver_type, &
n_procs_cmfd, cmfd_ilu_levels, master
use string, only: to_str
use, intrinsic :: ISO_FORTRAN_ENV
character(len=20) :: ksptype,pctype
! set up preconditioner
call KSPGetPC(krylov, prec, ierr)
if (n_procs_cmfd == 1) then
call PCSetType(prec, PCILU, ierr)
call PCFactorSetLevels(prec, cmfd_ilu_levels, ierr)
call KSPSetUp(krylov, ierr)
call PCFactorGetMatrix(prec, loss%M, ierr)
else
call PetscOptionsSetValue("-pc_type", "bjacobi", ierr)
call PetscOptionsSetValue("-sub_pc_type", "ilu", ierr)
call PetscOptionsSetValue("-sub_pc_factor_levels", &
trim(to_str(cmfd_ilu_levels)), ierr)
call PCSetFromOptions(prec, ierr)
call KSPSetUp(krylov, ierr)
end if
! get options
if (trim(cmfd_solver_type) == 'power') call KSPSetFromOptions(krylov, ierr)
! get all types and print
call KSPGetType(krylov, ksptype, ierr)
call PCGetType(prec, pctype, ierr)
! print heading information
if (cmfd_power_monitor .and. master) then
write(OUTPUT_UNIT,'(A)') ''
write(OUTPUT_UNIT,'(A)') '########################################################'
write(OUTPUT_UNIT,'(A)') '################ Power Iteration Solver ################'
write(OUTPUT_UNIT,'(A)') '########################################################'
write(OUTPUT_UNIT,'(A)')
write(OUTPUT_UNIT,100) 'Eigenvalue Tolerance:', ktol
write(OUTPUT_UNIT,100) 'Source Tolerance: ', stol
write(OUTPUT_UNIT,'(A)') ''
write(OUTPUT_UNIT,102) 'Linear Solver Type: ', ksptype
write(OUTPUT_UNIT,102) 'Preconditioner Type: ', pctype
write(OUTPUT_UNIT,101) 'ILU Fill Levels:', cmfd_ilu_levels
write(OUTPUT_UNIT,'(A)') ''
write(OUTPUT_UNIT,'(A)') '---------------------------------------------'
write(OUTPUT_UNIT,'(A)') ''
end if
100 FORMAT(A,1X,1PE11.4)
101 FORMAT(A,1X,I0)
102 FORMAT(A,1X,A)
end subroutine precondition_matrix
!===============================================================================
! COMPUTE_ADJOINT
! COMPUTE_ADJOINT computes a mathematical adjoint of CMFD problem
!===============================================================================
subroutine compute_adjoint()
use global, only: cmfd_write_matrices
PetscViewer :: viewer
! Transpose matrices
call loss % transpose()
call prod % transpose()
! transpose matrices
call MatTranspose(loss%M, MAT_REUSE_MATRIX, loss%M, ierr)
call MatTranspose(prod%F, MAT_REUSE_MATRIX, prod%F, ierr)
! write out matrix in binary file (debugging)
! Write out matrix in binary file (debugging)
if (cmfd_write_matrices) then
call PetscViewerBinaryOpen(PETSC_COMM_WORLD, 'adj_lossmat.bin', &
FILE_MODE_WRITE, viewer, ierr)
call MatView(loss%M, viewer, ierr)
call PetscViewerDestroy(viewer, ierr)
call PetscViewerBinaryOpen(PETSC_COMM_WORLD, 'adj_prodmat.bin', &
FILE_MODE_WRITE, viewer, ierr)
call MatView(prod%F, viewer, ierr)
call PetscViewerDestroy(viewer, ierr)
call loss % write_petsc_binary('adj_lossmat.bin')
call prod % write_petsc_binary('adj_prodmat.bin')
end if
end subroutine compute_adjoint
!===============================================================================
! EXECUTE_POWER_ITER in the main power iteration routine
! EXECUTE_POWER_ITER is the main power iteration routine
! for the cmfd calculation
!===============================================================================
subroutine execute_power_iter()
use constants, only: ONE
integer :: i ! iteration counter
real(8) :: num ! numerator for eigenvalue update
real(8) :: den ! denominator for eigenvalue update
integer :: i ! iteration counter
! reset convergence flag
! Reset convergence flag
iconv = .false.
! begin power iteration
! Begin power iteration
do i = 1, 10000
! compute source vector
call MatMult(prod%F, phi_o, S_o, ierr)
! Compute source vector
call prod % vector_multiply(phi_o, s_o)
! normalize source vector
call VecScale(S_o, ONE/k_o, ierr)
! Normalize source vector
s_o % val = s_o % val / k_o
! compute new flux vector
call KSPSolve(krylov, S_o, phi_n, ierr)
! Compute new flux vector
#ifdef PETSC
call gmres % solve(s_o, phi_n)
#endif
! compute new source vector
call MatMult(prod%F, phi_n, S_n, ierr)
! Compute new source vector
call prod % vector_multiply(phi_n, s_n)
! compute new k-eigenvalue
call VecSum(S_n, num, ierr)
call VecSum(S_o, den, ierr)
k_n = num/den
! Compute new k-eigenvalue
k_n = sum(s_n % val) / sum(s_o % val)
! renormalize the old source
call VecScale(S_o, k_o, ierr)
! Renormalize the old source
s_o % val = s_o % val * k_o
! check convergence
! Check convergence
call convergence(i)
! to break or not to break
! Break loop if converged
if (iconv) exit
! record old values
call VecCopy(phi_n, phi_o, ierr)
! 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
! CONVERGENCE checks the convergence of the CMFD problem
!===============================================================================
subroutine convergence(iter)
use global, only: cmfd_power_monitor, master
use constants, only: ONE, TINY_BIT
use global, only: cmfd_power_monitor, master
use, intrinsic :: ISO_FORTRAN_ENV
integer :: iter ! iteration number
integer, intent(in) :: iter ! iteration number
real(8) :: kerr ! error in keff
real(8) :: serr ! error in source
real(8) :: norm_n ! L2 norm of new source
real(8) :: norm_o ! L2 norm of old source
integer :: ierr ! petsc error code
! reset convergence flag
! Reset convergence flag
iconv = .false.
! calculate error in keff
! Calculate error in keff
kerr = abs(k_o - k_n)/k_n
! calculate max error in source
call VecNorm(S_n, NORM_2, norm_n, ierr)
call VecNorm(S_o, NORM_2, norm_o, ierr)
serr = abs(norm_n-norm_o)/norm_n
! Calculate max error in source
where (s_n % val > TINY_BIT)
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
! Check for convergence
if(kerr < ktol .and. serr < stol) iconv = .true.
! print out to user
! 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
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
! EXTRACT_RESULTS takes results and puts them in CMFD global data object
!===============================================================================
subroutine extract_results()
use constants, only: ZERO
use global, only: cmfd, n_procs_cmfd, cmfd_write_matrices
use global, only: cmfd, cmfd_write_matrices, current_batch
character(len=25) :: filename ! name of file to write data
integer :: n ! problem size
integer :: row_start ! local row start
integer :: row_end ! local row end
real(8),allocatable :: mybuf(:) ! temp buffer
PetscScalar, pointer :: phi_v(:) ! pointer to eigenvector info
PetscViewer :: viewer ! petsc viewer for binary write
! get problem size
n = loss%n
! Get problem size
n = loss % n
! also allocate in cmfd object
! 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
if (.not. allocated(mybuf)) allocate(mybuf(n))
! get ownership range
call VecGetOwnershipRange(phi_n, row_start, row_end, ierr)
! convert petsc phi_object to cmfd_obj
call VecGetArrayF90(phi_n, phi_v, ierr)
! Save values
if (adjoint_calc) then
cmfd%adj_phi(row_start+1:row_end) = phi_v
cmfd % adj_phi = phi_n % val
else
cmfd%phi(row_start+1:row_end) = phi_v
cmfd % phi = phi_n % val
end if
call VecRestoreArrayF90(phi_n, phi_v, ierr)
! save eigenvalue
! Save eigenvalue
if(adjoint_calc) then
cmfd%adj_keff = k_n
else
cmfd%keff = k_n
end if
! reduce result to all
mybuf = ZERO
if (adjoint_calc) then
call MPI_ALLREDUCE(cmfd%adj_phi, mybuf, n, MPI_REAL8, MPI_SUM, &
PETSC_COMM_WORLD, ierr)
else
call MPI_ALLREDUCE(cmfd%phi, mybuf, n, MPI_REAL8, MPI_SUM, &
PETSC_COMM_WORLD, ierr)
end if
! normalize phi to 1
! 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
! write out results
! Save dominance ratio
cmfd % dom(current_batch) = norm_n/norm_o
! Write out results
if (cmfd_write_matrices) then
if (adjoint_calc) then
call PetscViewerBinaryOpen(PETSC_COMM_WORLD, 'adj_fluxvec.bin', &
FILE_MODE_WRITE, viewer, ierr)
filename = 'adj_fluxvec.bin'
else
call PetscViewerBinaryOpen(PETSC_COMM_WORLD, 'fluxvec.bin', &
FILE_MODE_WRITE, viewer, ierr)
filename = 'fluxvec.bin'
end if
call VecView(phi_n, viewer, ierr)
call PetscViewerDestroy(viewer, ierr)
call phi_n % write_petsc_binary(filename)
end if
! nullify pointer and deallocate local vars
if (associated(phi_v)) nullify(phi_v)
if (allocated(mybuf)) deallocate(mybuf)
end subroutine extract_results
!===============================================================================
! FINALIZE
! FINALIZE frees all memory associated with power iteration
!===============================================================================
subroutine finalize()
use global, only: n_procs_cmfd
! finalize solver objects
call KSPDestroy(krylov, ierr)
call KSPDestroy(sub_krylov, ierr)
! finalize data objects
if (n_procs_cmfd > 1) call destroy_M_operator(loss) ! only destroy for jacobi
call destroy_F_operator(prod)
call VecDestroy(phi_n, ierr)
call VecDestroy(phi_o, ierr)
call VecDestroy(S_n, ierr)
call VecDestroy(S_o, ierr)
! 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
# endif
end module cmfd_power_solver

View file

@ -1,66 +1,72 @@
module cmfd_prod_operator
# ifdef PETSC
use constants, only: CMFD_NOACCEL
use global, only: cmfd, cmfd_coremap
use matrix_header, only: Matrix
implicit none
private
public :: init_F_operator,build_prod_matrix,destroy_F_operator
# include <finclude/petsc.h90>
integer :: nx ! maximum number of x cells
integer :: ny ! maximum number of y cells
integer :: nz ! maximum number of z cells
integer :: ng ! maximum number of groups
integer :: ierr ! petsc error code
type, public :: prod_operator
Mat :: F ! petsc matrix for neutronic prod operator
integer :: n ! dimensions of matrix
integer :: nnz ! max number of nonzeros
integer :: localn ! local size on proc
integer, allocatable :: d_nnz(:) ! vector of diagonal preallocation
integer, allocatable :: o_nnz(:) ! vector of off-diagonal preallocation
end type prod_operator
logical :: adjoint_calc = .false. ! adjoint calculation
public :: init_prod_matrix, build_prod_matrix
contains
!==============================================================================
! INIT_F_OPERATOR
! INIT_PROD_MATRIX preallocates prod matrix and initializes it
!==============================================================================
subroutine init_F_operator(this)
subroutine init_prod_matrix(prod_matrix)
type(prod_operator) :: this
type(Matrix), intent(inout) :: prod_matrix ! production matrix
! get indices
call get_F_indices(this)
integer :: nx ! maximum number of x cells
integer :: ny ! maximum number of y cells
integer :: nz ! maximum number of z cells
integer :: ng ! maximum number of groups
integer :: n ! total length of matrix
integer :: nnz ! number of nonzeros in matrix
! get preallocation
call preallocate_prod_matrix(this)
! Get maximum number of cells in each direction
nx = cmfd%indices(1)
ny = cmfd%indices(2)
nz = cmfd%indices(3)
ng = cmfd%indices(4)
! set up M operator
call MatCreateAIJ(PETSC_COMM_WORLD, this%localn, this%localn, PETSC_DECIDE, &
PETSC_DECIDE, PETSC_NULL_INTEGER, this%d_nnz, PETSC_NULL_INTEGER, &
this%o_nnz, this%F, ierr)
call MatSetOption(this%F, MAT_NEW_NONZERO_LOCATIONS, PETSC_TRUE, ierr)
call MatSetOption(this%F, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE, ierr)
! Calculate dimensions and number of nonzeros in matrix
if (cmfd_coremap) then
n = cmfd % mat_dim * ng
else
n = nx*ny*nz*ng
end if
nnz = n * ng
end subroutine init_F_operator
! Configure prod matrix
call prod_matrix % create(n, nnz)
!==============================================================================
! GET_F_INDICES
!==============================================================================
end subroutine init_prod_matrix
subroutine get_F_indices(this)
!===============================================================================
! BUILD_PROD_MATRIX creates the matrix representing production of neutrons
!===============================================================================
use global, only: cmfd, cmfd_coremap
subroutine build_prod_matrix(prod_matrix, adjoint)
type(prod_operator) :: this
type(Matrix), intent(inout) :: prod_matrix ! production matrix
logical, intent(in), optional :: adjoint ! adjoint calculation logical
integer :: i ! iteration counter for x
integer :: j ! iteration counter for y
integer :: k ! iteration counter for z
integer :: g ! iteration counter for groups
integer :: h ! energy group when doing scattering
integer :: nx ! maximum number of x cells
integer :: ny ! maximum number of y cells
integer :: nz ! maximum number of z cells
integer :: ng ! maximum number of groups
integer :: hmat_idx ! index in matrix for energy group h
integer :: irow ! iteration counter over row
logical :: adjoint_calc ! is this a physical adjoint?
real(8) :: nfissxs ! nufission cross section h-->g
real(8) :: val ! temporary variable for nfissxs
! get maximum number of cells in each direction
nx = cmfd%indices(1)
@ -68,162 +74,18 @@ contains
nz = cmfd%indices(3)
ng = cmfd%indices(4)
! get number of nonzeros
this%nnz = 7 + ng - 1
! calculate dimensions of matrix
if (cmfd_coremap) then
this%n = cmfd % mat_dim * ng
else
this%n = nx*ny*nz*ng
end if
end subroutine get_F_indices
!===============================================================================
! PREALLOCATE_PROD_MATRIX
!===============================================================================
subroutine preallocate_prod_matrix(this)
use constants, only: CMFD_NOACCEL
use global, only: cmfd, cmfd_coremap
type(prod_operator) :: this
integer :: rank ! rank of processor
integer :: sizen ! number of procs
integer :: i ! iteration counter for x
integer :: j ! iteration counter for y
integer :: k ! iteration counter for z
integer :: g ! iteration counter for groups
integer :: h ! energy group when doing scattering
integer :: n ! the extent of the matrix
integer :: irow ! row counter
integer :: row_start ! index of local starting row
integer :: row_end ! index of local final row
integer :: hmat_idx ! index in matrix for energy group h
! initialize size and rank
rank = 0
sizen = 0
! get rank and max rank of procs
call MPI_COMM_RANK(PETSC_COMM_WORLD, rank, ierr)
call MPI_COMM_SIZE(PETSC_COMM_WORLD, sizen, ierr)
! get local problem size
n = this%n
! determine local size, divide evenly between all other procs
this%localn = n/(sizen)
! add 1 more if less proc id is less than mod
if (rank < mod(n, sizen)) this%localn = this%localn + 1
! determine local starting row
row_start = 0
if (rank < mod(n, sizen)) then
row_start = rank*(n/sizen+1)
else
row_start = min(mod(n, sizen)*(n/sizen+1) + &
(rank - mod(n, sizen))*(n/sizen), n)
end if
! determine local final row
row_end = row_start + this%localn - 1
! allocate counters
if (.not. allocated(this%d_nnz)) allocate(this%d_nnz(row_start:row_end))
if (.not. allocated(this%o_nnz)) allocate(this%o_nnz(row_start:row_end))
this % d_nnz = 0
this % o_nnz = 0
! begin loop around local rows
ROWS: do irow = row_start, row_end
! initialize counters
this%d_nnz(irow) = 1 ! already add in matrix diagonal
this%o_nnz(irow) = 0
! get location indices
call matrix_to_indices(irow, g, i, j, k)
! check if not including reflector
if (cmfd_coremap) then
! check if at a reflector
if (cmfd % coremap(i,j,k) == CMFD_NOACCEL) then
cycle
end if
end if
! begin loop over off diagonal in-scattering
NFISS: do h = 1, ng
! cycle though if h=g
if (h == g) then
cycle
end if
! get neighbor matrix index
call indices_to_matrix(h, i, j, k, hmat_idx)
! record nonzero
if (((hmat_idx-1) >= row_start) .and. &
((hmat_idx-1) <= row_end)) then
this%d_nnz(irow) = this%d_nnz(irow) + 1
else
this%o_nnz(irow) = this%o_nnz(irow) + 1
end if
end do NFISS
end do ROWS
end subroutine preallocate_prod_matrix
!===============================================================================
! BUILD_PROD_MATRIX creates the matrix representing loss of neutrons
!===============================================================================
subroutine build_prod_matrix(this, adjoint)
use constants, only: CMFD_NOACCEL
use global, only: cmfd, cmfd_coremap, cmfd_write_matrices
type(prod_operator) :: this
logical, optional :: adjoint
integer :: i ! iteration counter for x
integer :: j ! iteration counter for y
integer :: k ! iteration counter for z
integer :: g ! iteration counter for groups
integer :: h ! energy group when doing scattering
integer :: hmat_idx ! index in matrix for energy group h
integer :: ierr ! Petsc error code
integer :: row_start ! the first local row on the processor
integer :: row_finish ! the last local row on the processor
integer :: irow ! iteration counter over row
real(8) :: nfissxs ! nufission cross section h-->g
real(8) :: val ! temporary variable for nfissxs
! check for adjoint
adjoint_calc = .false.
if (present(adjoint)) adjoint_calc = adjoint
! initialize row start and finish
row_start = 0
row_finish = 0
! get row bounds for this processor
call MatGetOwnershipRange(this%F, row_start, row_finish, ierr)
! begin iteration loops
ROWS: do irow = row_start, row_finish-1
ROWS: do irow = 1, prod_matrix % n
! add a new row to matrix
call prod_matrix % new_row()
! get indices for that row
call matrix_to_indices(irow, g, i, j, k)
call matrix_to_indices(irow, g, i, j, k, ng, nx, ny, nz)
! check if not including reflector
if (cmfd_coremap) then
@ -236,34 +98,34 @@ contains
end if
! loop around all other groups
NFISS: do h = 1, ng
! get cell data
nfissxs = cmfd%nfissxs(h,g,i,j,k)
! get matrix column location
call indices_to_matrix(h, i, j, k, hmat_idx)
! reocrd value in matrix
val = nfissxs
call indices_to_matrix(h, i, j, k, hmat_idx, ng, nx, ny)
! check for adjoint and bank val
if (adjoint_calc) then
call MatSetValue(this%F, hmat_idx-1, irow, val, INSERT_VALUES, ierr)
! get nu-fission cross section from cell
nfissxs = cmfd%nfissxs(g,h,i,j,k)
else
call MatSetValue(this%F, irow, hmat_idx-1, val, INSERT_VALUES, ierr)
! get nu-fission cross section from cell
nfissxs = cmfd%nfissxs(h,g,i,j,k)
end if
! set as value to be recorded
val = nfissxs
! record value in matrix
call prod_matrix % add_value(hmat_idx, val)
end do NFISS
end do ROWS
! assemble matrix
call MatAssemblyBegin(this%F, MAT_FLUSH_ASSEMBLY, ierr)
call MatAssemblyEnd(this%F, MAT_FINAL_ASSEMBLY, ierr)
! print out operator to file
if (cmfd_write_matrices) call print_F_operator(this)
! CSR requires n+1 row
call prod_matrix % new_row()
end subroutine build_prod_matrix
@ -271,19 +133,20 @@ contains
! INDICES_TO_MATRIX takes (x,y,z,g) indices and computes location in matrix
!===============================================================================
subroutine indices_to_matrix(g, i, j, k, matidx)
subroutine indices_to_matrix(g, i, j, k, matidx, ng, nx, ny)
use global, only: cmfd, cmfd_coremap
integer, intent(out) :: matidx ! the index location in matrix
integer, intent(in) :: i ! current x index
integer, intent(in) :: j ! current y index
integer, intent(in) :: k ! current z index
integer, intent(in) :: g ! current group index
integer, intent(in) :: nx ! maximum number of x cells
integer, intent(in) :: ny ! maximum number of y cells
integer, intent(in) :: ng ! maximum number of groups
integer :: matidx ! the index location in matrix
integer :: i ! current x index
integer :: j ! current y index
integer :: k ! current z index
integer :: g ! current group index
! check if coremap is used
if (cmfd_coremap) then
! get idx from core map
matidx = ng*(cmfd % coremap(i,j,k)) - (ng - g)
@ -297,80 +160,40 @@ contains
end subroutine indices_to_matrix
!===============================================================================
! MATRIX_TO_INDICES
! MATRIX_TO_INDICES converts matrix index to spatial and group indices
!===============================================================================
subroutine matrix_to_indices(irow, g, i, j, k)
subroutine matrix_to_indices(irow, g, i, j, k, ng, nx, ny, nz)
use global, only: cmfd, cmfd_coremap
integer :: i ! iteration counter for x
integer :: j ! iteration counter for y
integer :: k ! iteration counter for z
integer :: g ! iteration counter for groups
integer :: irow ! iteration counter over row (0 reference)
integer, intent(out) :: i ! iteration counter for x
integer, intent(out) :: j ! iteration counter for y
integer, intent(out) :: k ! iteration counter for z
integer, intent(out) :: g ! iteration counter for groups
integer, intent(in) :: irow ! iteration counter over row (0 reference)
integer, intent(in) :: nx ! maximum number of x cells
integer, intent(in) :: ny ! maximum number of y cells
integer, intent(in) :: nz ! maximum number of z cells
integer, intent(in) :: ng ! maximum number of groups
! check for core map
if (cmfd_coremap) then
! get indices from indexmap
g = mod(irow,ng) + 1
i = cmfd % indexmap(irow/ng+1,1)
j = cmfd % indexmap(irow/ng+1,2)
k = cmfd % indexmap(irow/ng+1,3)
g = mod(irow-1, ng) + 1
i = cmfd % indexmap((irow-1)/ng+1,1)
j = cmfd % indexmap((irow-1)/ng+1,2)
k = cmfd % indexmap((irow-1)/ng+1,3)
else
! compute indices
g = mod(irow, ng) + 1
i = mod(irow, ng*nx)/ng + 1
j = mod(irow, ng*nx*ny)/(ng*nx)+ 1
k = mod(irow, ng*nx*ny*nz)/(ng*nx*ny) + 1
g = mod(irow-1, ng) + 1
i = mod(irow-1, ng*nx)/ng + 1
j = mod(irow-1, ng*nx*ny)/(ng*nx)+ 1
k = mod(irow-1, ng*nx*ny*nz)/(ng*nx*ny) + 1
end if
end subroutine matrix_to_indices
!===============================================================================
! PRINT_F_OPERATOR
!===============================================================================
subroutine print_F_operator(this)
type(prod_operator) :: this
PetscViewer :: viewer
! write out matrix in binary file (debugging)
if (adjoint_calc) then
call PetscViewerBinaryOpen(PETSC_COMM_WORLD, 'adj_prodmat.bin', &
FILE_MODE_WRITE, viewer, ierr)
else
call PetscViewerBinaryOpen(PETSC_COMM_WORLD, 'prodmat.bin', &
FILE_MODE_WRITE, viewer, ierr)
end if
call MatView(this%F, viewer, ierr)
call PetscViewerDestroy(viewer, ierr)
end subroutine print_F_operator
!==============================================================================
! DESTROY_F_OPERATOR
!==============================================================================
subroutine destroy_F_operator(this)
type(prod_operator) :: this
! deallocate matrix
call MatDestroy(this%F, ierr)
! deallocate other parameters
if (allocated(this%d_nnz)) deallocate(this%d_nnz)
if (allocated(this%o_nnz)) deallocate(this%o_nnz)
end subroutine destroy_F_operator
# endif
end module cmfd_prod_operator

View file

@ -1,473 +0,0 @@
module cmfd_snes_solver
# ifdef PETSC
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
use cmfd_jacobian_operator, only: jacobian_operator,init_J_operator, &
build_jacobian_matrix,destroy_J_operator,&
operators
use cmfd_power_solver, only: cmfd_power_execute
implicit none
private
public :: cmfd_snes_execute
# include <finclude/petsc.h90>
type(jacobian_operator) :: jac_prec
type(operators) :: ctx
Mat :: jac ! jacobian matrix
Vec :: resvec ! residual vector
Vec :: xvec ! results
KSP :: ksp ! linear solver context
PC :: pc ! preconditioner
SNES :: snes ! nonlinear solver context
integer :: ierr ! error flag
logical :: adjoint_calc = .false. ! adjoint calculation
contains
!===============================================================================
! CMFD_SNES_EXECUTE
!===============================================================================
subroutine cmfd_snes_execute(adjoint)
logical, optional :: adjoint ! adjoint calculation
! check for adjoint
if (present(adjoint)) adjoint_calc = adjoint
! seed with power iteration
call cmfd_power_execute(k_tol=1.E-3_8, s_tol=1.E-3_8, adjoint=adjoint_calc)
! initialize data
call init_data()
! initialize solver
call init_solver()
! precondition
call precondition_matrix()
! solve the system
call SNESSolve(snes, PETSC_NULL_DOUBLE, xvec, ierr)
! extracts results to cmfd object
call extract_results()
! deallocate all slepc data
call finalize()
end subroutine cmfd_snes_execute
!===============================================================================
! INIT_DATA allocates matrices vectors for CMFD solution
!===============================================================================
subroutine init_data()
use constants, only: ONE
use global, only: cmfd, n_procs_cmfd, rank
integer :: k ! implied do counter
integer :: n ! problem size
integer :: row_start ! local row start
integer :: row_end ! local row end
real(8), pointer :: xptr(:) ! solution pointer
! set up operator matrices
call init_M_operator(ctx%loss)
call init_F_operator(ctx%prod)
call init_J_operator(jac_prec, ctx)
! get problem size
n = jac_prec%n
! create PETSc vectors
call VecCreateMPI(PETSC_COMM_WORLD, jac_prec%localn, PETSC_DECIDE, &
resvec, ierr)
call VecCreateMPI(PETSC_COMM_WORLD, jac_prec%localn, PETSC_DECIDE, &
xvec, ierr)
! get the local dimensions for each process
call VecGetOwnershipRange(xvec, row_start, row_end, ierr)
if (rank == n_procs_cmfd - 1) row_end = n
! set flux in guess
if (adjoint_calc) then
call VecSetValues(xvec, row_end-row_start, (/(k,k=row_start,row_end-1)/), &
cmfd%adj_phi(row_start+1:row_end), INSERT_VALUES, ierr)
else
call VecSetValues(xvec, row_end-row_start, (/(k,k=row_start,row_end-1)/), &
cmfd%phi(row_start+1:row_end), INSERT_VALUES, ierr)
end if
call VecAssemblyBegin(xvec, ierr)
call VecAssemblyEnd(xvec, ierr)
! set keff in guess
if (rank == n_procs_cmfd - 1) then
call VecGetArrayF90(xvec, xptr, ierr)
if (adjoint_calc) then
xptr(size(xptr)) = ONE/cmfd%adj_keff
else
xptr(size(xptr)) = ONE/cmfd%keff
end if
call VecRestoreArrayF90(xvec, xptr, ierr)
end if
! nullify pointers
if (associated(xptr)) nullify(xptr)
end subroutine init_data
!===============================================================================
! INIT_SOLVER
!===============================================================================
subroutine init_solver()
use global, only: cmfd_snes_monitor, n_procs_cmfd
! turn on mf_operator option
call PetscOptionsSetValue("-snes_mf_operator", "TRUE", ierr)
! create SNES context
call SNESCreate(PETSC_COMM_WORLD, snes, ierr)
! set the residual function
call SNESSetFunction(snes, resvec, compute_nonlinear_residual, &
PETSC_NULL_DOUBLE, ierr)
! set GMRES solver
call SNESGetKSP(snes, ksp, ierr)
call KSPSetType(ksp, KSPGMRES, ierr)
! create matrix free jacobian
call MatCreateSNESMF(snes, jac, ierr)
! set matrix free finite difference
call SNESSetJacobian(snes, jac, jac_prec, build_jacobian_matrix, ctx, ierr)
! set lags
call SNESSetLagJacobian(snes, -2, ierr)
call SNESSetLagPreconditioner(snes, -1, ierr)
! set convergence
! call SNESSetTolerances(snes, 1.0e-8_8, 1.0e-8_8, 1.0e-8_8, 50, 10000, ierr)
! set SNES options
call SNESSetFromOptions(snes, ierr)
! turn off line searching
! call SNESLineSearchSet(snes, SNESLineSearchNo, PETSC_NULL_OBJECT, ierr)
end subroutine init_solver
!===============================================================================
! PRECONDITION_MATRIX
!===============================================================================
subroutine precondition_matrix()
use global, only: cmfd_snes_monitor, cmfd_solver_type, &
n_procs_cmfd, cmfd_ilu_levels, master
use string, only: to_str
use, intrinsic :: ISO_FORTRAN_ENV
character(LEN=20) :: snestype, ksptype, pctype
! set up preconditioner
call KSPGetPC(ksp, pc, ierr)
if (n_procs_cmfd == 1) then
call PCSetType(pc, PCILU, ierr)
call PCFactorSetLevels(pc, cmfd_ilu_levels, ierr)
else
call PetscOptionsSetValue("-pc_type", "bjacobi", ierr)
call PetscOptionsSetValue("-sub_pc_type", "ilu", ierr)
call PetscOptionsSetValue("-sub_pc_factor_levels", &
trim(to_str(cmfd_ilu_levels)), ierr)
end if
! get options
call PCSetFromOptions(pc, ierr)
call KSPSetFromOptions(ksp, ierr)
! finalize ksp setup
! call KSPSetUp(ksp, ierr)
! get all types and print
call SNESGetType(snes, snestype, ierr)
call KSPGetType(ksp, ksptype, ierr)
call PCGetType(pc, pctype, ierr)
! display solver info to user
if (cmfd_snes_monitor .and. master) then
write(OUTPUT_UNIT,'(A)') ''
write(OUTPUT_UNIT,'(A)') '########################################################'
write(OUTPUT_UNIT,'(A)') '################ JFNK Nonlinear Solver ################'
write(OUTPUT_UNIT,'(A)') '########################################################'
write(OUTPUT_UNIT,'(A)') ''
write(OUTPUT_UNIT,101) 'NONLINEAR SOLVER: ', snestype
write(OUTPUT_UNIT,101) 'LINEAR SOLVER: ', ksptype
write(OUTPUT_UNIT,101) 'PRECONDITIONER: ', pctype
write(OUTPUT_UNIT,100) 'ILU levels: ', cmfd_ilu_levels
write(OUTPUT_UNIT,'(A)') ''
write(OUTPUT_UNIT,'(A)') '---------------------------------------------'
write(OUTPUT_UNIT,'(A)') ''
end if
100 FORMAT(A,1X,I0)
101 FORMAT(A,1X,A)
end subroutine precondition_matrix
!===============================================================================
! COMPUTE_NONLINEAR_RESIDUAL
!===============================================================================
subroutine compute_nonlinear_residual(snes, x, res, ierr)
use global, only: n_procs_cmfd, cmfd_write_matrices, &
cmfd_adjoint_type, rank
SNES :: snes ! nonlinear solver context
Vec :: x ! independent vector
Vec :: res ! residual vector
integer :: ierr ! error flag
Vec :: phi ! flux vector
Vec :: rphi ! flux part of residual
Vec :: phiM ! M part of residual flux calc
integer :: n ! problem size
real(8) :: lambda ! eigenvalue
real(8) :: reslamb ! residual for lambda
real(8), pointer :: xptr(:) ! pointer to solution vector
real(8), pointer :: rptr(:) ! pointer to residual vector
PetscViewer :: viewer
logical :: physical_adjoint = .false.
! check for physical adjoint
if (adjoint_calc .and. trim(cmfd_adjoint_type) == 'physical') &
physical_adjoint = .true.
! create operators
call build_loss_matrix(ctx%loss, adjoint = physical_adjoint)
call build_prod_matrix(ctx%prod, adjoint = physical_adjoint)
! check for adjoint
if (adjoint_calc .and. trim(cmfd_adjoint_type) == 'math') &
call compute_adjoint()
! get problem size
n = ctx%loss%n
! get pointers to vectors
call VecGetArrayF90(x, xptr, ierr)
call VecGetArrayF90(res, rptr, ierr)
! create petsc vector for flux
call VecCreateMPI(PETSC_COMM_WORLD, ctx%loss%localn, &
PETSC_DECIDE, phi, ierr)
call VecCreateMPI(PETSC_COMM_WORLD, ctx%loss%localn, &
PETSC_DECIDE, rphi, ierr)
! extract flux and place in petsc vector
call VecPlaceArray(phi, xptr, ierr)
call VecPlaceArray(rphi, rptr, ierr)
! extract eigenvalue and broadcast (going to want to make this more general in future)
if (rank == n_procs_cmfd - 1) lambda = xptr(size(xptr))
call MPI_BCAST(lambda, 1, MPI_REAL8, n_procs_cmfd-1, PETSC_COMM_WORLD, ierr)
! create new petsc vectors to perform math
call VecCreateMPI(PETSC_COMM_WORLD, ctx%loss%localn, PETSC_DECIDE, &
phiM, ierr)
! calculate flux part of residual vector
call MatMult(ctx%loss%M, phi, phiM, ierr)
call MatMult(ctx%prod%F, phi, rphi, ierr)
call VecAYPX(rphi, -lambda, phiM, ierr)
! set eigenvalue part of residual vector
call VecDot(phi, phi, reslamb, ierr)
! map to ptr
if (rank == n_procs_cmfd - 1) rptr(size(rptr)) = 0.5_8 - 0.5_8*reslamb
! reset arrays that are not used
call VecResetArray(phi, ierr)
call VecResetArray(rphi, ierr)
! restore arrays for residual and solution
call VecRestoreArrayF90(x, xptr, ierr)
call VecRestoreArrayF90(res, rptr, ierr)
! destroy all temp vectors
call VecDestroy(phi, ierr)
call VecDestroy(phiM, ierr)
call VecDestroy(rphi, ierr)
! nullify all pointers
if (associated(xptr)) nullify(xptr)
if (associated(rptr)) nullify(rptr)
! write out matrix in binary file (debugging)
if (cmfd_write_matrices) then
if (adjoint_calc) then
call PetscViewerBinaryOpen(PETSC_COMM_WORLD, 'adj_res.bin', &
FILE_MODE_WRITE, viewer, ierr)
else
call PetscViewerBinaryOpen(PETSC_COMM_WORLD, 'res.bin', &
FILE_MODE_WRITE, viewer, ierr)
end if
call VecView(res, viewer, ierr)
call PetscViewerDestroy(viewer, ierr)
if (adjoint_calc) then
call PetscViewerBinaryOpen(PETSC_COMM_WORLD, 'adj_x.bin', &
FILE_MODE_WRITE, viewer, ierr)
else
call PetscViewerBinaryOpen(PETSC_COMM_WORLD, 'x.bin', &
FILE_MODE_WRITE, viewer, ierr)
end if
call VecView(x, viewer, ierr)
call PetscViewerDestroy(viewer, ierr)
end if
end subroutine compute_nonlinear_residual
!===============================================================================
! COMPUTE_ADJOINT
!===============================================================================
subroutine compute_adjoint()
use global, only: cmfd_write_matrices
PetscViewer :: viewer
! transpose matrices
call MatTranspose(ctx%loss%M, MAT_REUSE_MATRIX, ctx%loss%M, ierr)
call MatTranspose(ctx%prod%F, MAT_REUSE_MATRIX, ctx%prod%F, ierr)
! write out matrix in binary file (debugging)
if (cmfd_write_matrices) then
call PetscViewerBinaryOpen(PETSC_COMM_WORLD, 'adj_lossmat.bin', &
FILE_MODE_WRITE, viewer, ierr)
call MatView(ctx%loss%M, viewer, ierr)
call PetscViewerDestroy(viewer, ierr)
call PetscViewerBinaryOpen(PETSC_COMM_WORLD, 'adj_prodmat.bin', &
FILE_MODE_WRITE, viewer, ierr)
call MatView(ctx%prod%F, viewer, ierr)
call PetscViewerDestroy(viewer, ierr)
end if
end subroutine compute_adjoint
!===============================================================================
! EXTRACT_RESULTS
!===============================================================================
subroutine extract_results()
use constants, only: ZERO, ONE
use global, only: cmfd, n_procs_cmfd, rank
integer :: n ! problem size
integer :: row_start ! local row start
integer :: row_end ! local row end
real(8) :: keff ! keff of problem
real(8),allocatable :: mybuf(:) ! temp buffer
PetscScalar, pointer :: xptr(:) ! pointer to eigenvector info
! get problem size
n = ctx%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
if (.not. allocated(mybuf)) allocate(mybuf(n))
! get ownership range
call VecGetOwnershipRange(xvec, row_start, row_end, ierr)
! resize the last proc
if (rank == n_procs_cmfd - 1) row_end = row_end - 1
! convert petsc phi_object to cmfd_obj
call VecGetArrayF90(xvec, xptr, ierr)
if (adjoint_calc) then
cmfd%adj_phi(row_start+1:row_end) = xptr(1:row_end - row_start)
else
cmfd%phi(row_start+1:row_end) = xptr(1:row_end - row_start)
end if
! reduce result to all
mybuf = ZERO
if (adjoint_calc) then
call MPI_ALLREDUCE(cmfd%adj_phi, mybuf, n, MPI_REAL8, MPI_SUM, &
PETSC_COMM_WORLD, ierr)
else
call MPI_ALLREDUCE(cmfd%phi, mybuf, n, MPI_REAL8, MPI_SUM, &
PETSC_COMM_WORLD, ierr)
end if
! move buffer to object and deallocate
cmfd%phi = mybuf
if(allocated(mybuf)) deallocate(mybuf)
! save eigenvalue
if(rank == n_procs_cmfd - 1) keff = ONE / xptr(size(xptr))
if (adjoint_calc) then
cmfd%adj_keff = keff
call MPI_BCAST(cmfd%adj_keff, 1, MPI_REAL8, n_procs_cmfd-1, &
PETSC_COMM_WORLD, ierr)
else
cmfd%keff = keff
call MPI_BCAST(cmfd%keff, 1, MPI_REAL8, n_procs_cmfd-1, &
PETSC_COMM_WORLD, ierr)
end if
call VecRestoreArrayF90(xvec, xptr, ierr)
! nullify pointers and deallocate local variables
if (associated(xptr)) nullify(xptr)
if (allocated(mybuf)) deallocate(mybuf)
end subroutine extract_results
!===============================================================================
! FINALIZE
!===============================================================================
subroutine finalize()
! finalize data objects
call destroy_M_operator(ctx%loss)
call destroy_F_operator(ctx%prod)
call destroy_J_operator(jac_prec)
call VecDestroy(xvec, ierr)
call VecDestroy(resvec, ierr)
call MatDestroy(jac, ierr)
! finalize solver objects
call SNESDestroy(snes, ierr)
end subroutine finalize
# endif
end module cmfd_snes_solver

View file

@ -11,7 +11,7 @@ module constants
integer, parameter :: VERSION_RELEASE = 2
! Revision numbers for binary files
integer, parameter :: REVISION_STATEPOINT = 9
integer, parameter :: REVISION_STATEPOINT = 10
integer, parameter :: REVISION_PARTICLE_RESTART = 1
! Binary file types
@ -362,9 +362,8 @@ module constants
integer, parameter :: UNIT_TALLY = 12 ! unit # for writing tally file
integer, parameter :: UNIT_PLOT = 13 ! unit # for writing plot file
integer, parameter :: UNIT_XS = 14 ! unit # for writing xs summary file
integer, parameter :: CMFD_BALANCE = 15 ! unit # for writing cmfd balance file
integer, parameter :: UNIT_PARTICLE = 16 ! unit # for writing particle restart
integer, parameter :: UNIT_OUTPUT = 17 ! unit # for writing output
integer, parameter :: UNIT_PARTICLE = 15 ! unit # for writing particle restart
integer, parameter :: UNIT_OUTPUT = 16 ! unit # for writing output
!=============================================================================
! CMFD CONSTANTS

View file

@ -44,7 +44,7 @@ contains
if (master) call header("K EIGENVALUE SIMULATION", level=1)
! Display column titles
call print_columns()
if(master) call print_columns()
! Turn on inactive timer
call time_inactive % start()

View file

@ -1,8 +1,5 @@
module finalize
# ifdef PETSC
use cmfd_output, only: finalize_cmfd
# endif
use global
use output, only: print_runtime, print_results, &
print_overlap_check, write_tallies
@ -42,11 +39,11 @@ contains
end if
#ifdef PETSC
! finalize cmfd
if (cmfd_run) call finalize_cmfd()
! Finalize PETSc
if (cmfd_run) call PetscFinalize(mpi_err)
#endif
! stop timers and show timing statistics
! Stop timers and show timing statistics
call time_finalize % stop()
call time_total % stop()
if (master .and. (run_mode /= MODE_PLOTTING .and. &
@ -56,7 +53,7 @@ contains
if (check_overlaps) call print_overlap_check()
end if
! deallocate arrays
! Deallocate arrays
call free_memory()
#ifdef HDF5

View file

@ -15,10 +15,6 @@ module global
use tally_header, only: TallyObject, TallyMap, TallyResult
use timer_header, only: Timer
#ifdef MPI
use mpi
#endif
#ifdef HDF5
use hdf5_interface, only: HID_T
#endif
@ -300,80 +296,64 @@ module global
logical :: cmfd_run = .false.
! Timing objects
type(Timer) :: time_cmfd ! timer for whole cmfd calculation
type(Timer) :: time_solver ! timer for solver
type(Timer) :: time_cmfd ! timer for whole cmfd calculation
type(Timer) :: time_cmfdbuild ! timer for matrix build
type(Timer) :: time_cmfdsolve ! timer for solver
! Flag for CMFD only
logical :: cmfd_only = .false.
! Flag for coremap accelerator
! Flag for active core map
logical :: cmfd_coremap = .false.
! number of processors for cmfd
integer :: n_procs_cmfd
! reset dhats to zero
! Flag to reset dhats to zero
logical :: dhat_reset = .false.
! activate neutronic feedback
! Flag to activate neutronic feedback via source weights
logical :: cmfd_feedback = .false.
! activate auto-balance of tallies (2grp only)
! logical :: cmfd_balance = .false.
! User-defined tally information
integer :: n_cmfd_meshes = 1 ! # of structured meshes
integer :: n_cmfd_tallies = 3 ! # of user-defined tallies
! calculate effective downscatter
! logical :: cmfd_downscatter = .false.
! user-defined tally information
integer :: n_cmfd_meshes = 1 ! # of structured meshes
integer :: n_cmfd_tallies = 3 ! # of user-defined tallies
! overwrite with 2grp xs
logical :: cmfd_run_2grp = .false.
! hold cmfd weight adjustment factors
! Flag to hold cmfd weight adjustment factors
logical :: cmfd_hold_weights = .false.
! eigenvalue solver type
! Eigenvalue solver type
character(len=10) :: cmfd_solver_type = 'power'
! adjoint method type
! Adjoint method type
character(len=10) :: cmfd_adjoint_type = 'physical'
! number of incomplete ilu factorization levels
! Number of incomplete ilu factorization levels
integer :: cmfd_ilu_levels = 1
! batch to begin cmfd
! Batch to begin cmfd
integer :: cmfd_begin = 1
! when and how long to flush cmfd tallies during inactive batches
! When and how long to flush cmfd tallies during inactive batches
integer :: cmfd_inact_flush(2) = (/9999,1/)
! batch to last flush before active batches
! Batch to last flush before active batches
integer :: cmfd_act_flush = 0
! compute effective downscatter cross section
! Compute effective downscatter cross section
logical :: cmfd_downscatter = .false.
! convergence monitoring
! Convergence monitoring
logical :: cmfd_snes_monitor = .false.
logical :: cmfd_ksp_monitor = .false.
logical :: cmfd_power_monitor = .false.
! cmfd output
logical :: cmfd_write_balance = .false.
! Cmfd output
logical :: cmfd_write_matrices = .false.
logical :: cmfd_write_hdf5 = .false.
! run an adjoint calculation (last batch only)
! Run an adjoint calculation (last batch only)
logical :: cmfd_run_adjoint = .false.
! cmfd run logicals
! CMFD run logicals
logical :: cmfd_on = .false.
logical :: cmfd_tally_on = .true.
! tolerance on keff to run cmfd
real(8) :: cmfd_keff_tol = 0.005_8
! CMFD display info
character(len=25) :: cmfd_display
! Information about state points to be written
integer :: n_state_points = 0
@ -463,7 +443,7 @@ contains
! Deallocate array of work indices
if (allocated(work_index)) deallocate(work_index)
! Deallocate cmfd
! Deallocate CMFD
call deallocate_cmfd(cmfd)
! Deallocate tally node lists

View file

@ -575,15 +575,13 @@ contains
end do
end if
! check for cmfd run
! Check for cmfd run
call lower_case(run_cmfd_)
if (run_cmfd_ == 'true' .or. run_cmfd_ == '1') then
cmfd_run = .true.
#ifndef PETSC
if (master) then
message = 'CMFD is not available, compile OpenMC with PETSc'
call fatal_error()
end if
message = 'CMFD is not available, recompile OpenMC with PETSc'
call fatal_error()
#endif
end if

358
src/matrix_header.F90 Normal file
View file

@ -0,0 +1,358 @@
module matrix_header
#ifdef PETSC
use petscmat
#endif
implicit none
private
type, public :: Matrix
integer :: n ! number of rows/cols in matrix
integer :: nnz ! number of nonzeros in matrix
integer :: n_count ! counter for length of matrix
integer :: nz_count ! counter for number of non zeros
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
procedure :: add_value => matrix_add_value
procedure :: new_row => matrix_new_row
procedure :: assemble => matrix_assemble
procedure :: get_row => matrix_get_row
procedure :: get_col => matrix_get_col
procedure :: setup_petsc => matrix_setup_petsc
procedure :: write_petsc_binary => matrix_write_petsc_binary
procedure :: transpose => matrix_transpose
procedure :: vector_multiply => matrix_vector_multiply
end type matrix
integer :: petsc_err
contains
!===============================================================================
! MATRIX_CREATE allocates CSR vectors
!===============================================================================
subroutine matrix_create(self, n, nnz)
integer, intent(in) :: n ! dimension of matrix
integer, intent(in) :: nnz ! number of nonzeros
class(Matrix), intent(inout) :: self ! matrix instance
! Preallocate vectors
if (.not.allocated(self % row)) allocate(self % row(n+1))
if (.not.allocated(self % col)) allocate(self % col(nnz))
if (.not.allocated(self % val)) allocate(self % val(nnz))
! Set counters to 1
self % n_count = 1
self % nz_count = 1
! Set n and nnz
self % n = n
self % nnz = nnz
! Set PETSc active by default to false
self % petsc_active = .false.
end subroutine matrix_create
!===============================================================================
! MATRIX_DESTROY deallocates all space associated with a matrix
!===============================================================================
subroutine matrix_destroy(self)
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)
end subroutine matrix_destroy
!===============================================================================
! MATRIX_ADD_VALUE adds a value to the matrix
!===============================================================================
subroutine matrix_add_value(self, col, val)
integer, intent(in) :: col ! col location in matrix
real(8), intent(in) :: val ! value to store in matrix
class(Matrix), intent(inout) :: self ! matrix instance
! Record the data
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
end subroutine matrix_add_value
!===============================================================================
! MATRIX_NEW_ROW adds a new row by saving the column position of the new row
!===============================================================================
subroutine matrix_new_row(self)
class(Matrix), intent(inout) :: self ! matrix instance
! 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
end subroutine matrix_new_row
!===============================================================================
! MATRIX_ASSEMBLE main rountine to sort all the columns in CSR matrix
!===============================================================================
subroutine matrix_assemble(self)
class(Matrix), intent(inout) :: self ! matrix instance
integer :: i
integer :: first
integer :: last
! Loop around row vector
do i = 1, self % n
! Get bounds
first = self % row(i)
last = self % row(i+1) - 1
! Sort a row
call sort_csr(self % col, self % val, first, last)
end do
end subroutine matrix_assemble
!===============================================================================
! SORT_CSR main routine that performs a sort on a CSR col/val subvector
!===============================================================================
recursive subroutine sort_csr(col, val, first, last)
integer :: col(:) ! column vector to sort
integer :: first ! first value in sort
integer :: last ! last value in sort
real(8) :: val(:) ! value vector to be sorted like columns
integer :: mid ! midpoint value
if (first < last) then
call split(col, val, first, last, mid)
call sort_csr(col, val, first, mid-1)
call sort_csr(col, val, mid+1, last)
end if
end subroutine sort_csr
!===============================================================================
! SPLIT bisects the search space for the sorting routine
!===============================================================================
subroutine split(col, val, low, high, mid)
integer :: col(:) ! column vector to sort
integer :: low ! low index to sort
integer :: high ! high index to sort
integer :: mid ! middle of sort
real(8) :: val(:) ! value vector to be sorted like columns
integer :: left ! contains left value in sort
integer :: right ! contains right value in sort
integer :: iswap ! temporary interger swap
integer :: pivot ! pivotting variable for columns
real(8) :: rswap ! temporary real swap
real(8) :: val0 ! pivot for value vector
left = low
right = high
pivot = col(low)
val0 = val(low)
! Repeat the following while left and right havent met
do while (left < right)
! Scan right to left to find element < pivot
do while (left < right .and. col(right) >= pivot)
right = right - 1
end do
! Scan left to right to find element > pivot
do while (left < right .and. col(left) <= pivot)
left = left + 1
end do
! If left and right havent met, exchange the items
if (left < right) then
iswap = col(left)
col(left) = col(right)
col(right) = iswap
rswap = val(left)
val(left) = val(right)
val(right) = rswap
end if
end do
! Switch the element in split position with pivot
col(low) = col(right)
col(right) = pivot
mid = right
val(low) = val(right)
val(right) = val0
end subroutine split
!===============================================================================
! MATRIX_GET_ROW is a method to get row and checks for PETSc C indexing use
!===============================================================================
function matrix_get_row(self, i) result(row)
class(Matrix), intent(in) :: self ! matrix instance
integer, intent(in) :: i ! row to get
integer :: row ! index in col where row begins
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
!===============================================================================
function matrix_get_col(self, i) result(col)
class(Matrix), intent(in) :: self ! matrix instance
integer, intent(in) :: i ! index from row vector
integer :: col ! the actual column
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
!===============================================================================
subroutine matrix_setup_petsc(self)
class(Matrix), intent(inout) :: self ! matrix instance
! change indices to c notation
self % row = self % row - 1
self % col = self % col - 1
! Link to petsc
#ifdef PETSC
call MatCreateSeqAIJWithArrays(PETSC_COMM_WORLD, self % n, self % n, &
self % row, self % col, self % val, self % petsc_mat, petsc_err)
#endif
! Petsc is now active
self % petsc_active = .true.
end subroutine matrix_setup_petsc
!===============================================================================
! MATRIX_WRITE_PETSC_BINARY writes a PETSc matrix binary file
!===============================================================================
subroutine matrix_write_petsc_binary(self, filename)
character(*), intent(in) :: filename ! file name to write to
class(Matrix), intent(in) :: self ! matrix instance
#ifdef PETSC
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)
#endif
end subroutine matrix_write_petsc_binary
!===============================================================================
! MATRIX_TRANSPOSE uses PETSc to transpose a matrix
!===============================================================================
subroutine matrix_transpose(self)
class(Matrix), intent(inout) :: self ! matrix instance
#ifdef PETSC
call MatTranspose(self % petsc_mat, MAT_REUSE_MATRIX, self % petsc_mat, &
petsc_err)
#endif
end subroutine matrix_transpose
!===============================================================================
! MATRIX_VECTOR_MULTIPLY allow a vector to multiply the matrix
!===============================================================================
subroutine matrix_vector_multiply(self, vec_in, vec_out)
use constants, only: ZERO
use vector_header, only: Vector
class(Matrix), intent(in) :: self ! matrix instance
type(Vector), intent(in) :: vec_in ! vector to multiply matrix against
type(Vector), intent(inout) :: vec_out ! resulting vector
integer :: i ! row iteration counter
integer :: j ! column iteration counter
! Begin loop around rows
ROWS: do i = 1, self % n
! Initialize target location in vector
vec_out % val(i) = ZERO
! Begin loop around columns
COLS: do j = self % get_row(i), self % get_row(i + 1) - 1
vec_out % val(i) = vec_out % val(i) + self % val(j) * &
vec_in % val(self % get_col(j))
end do COLS
end do ROWS
end subroutine matrix_vector_multiply
end module matrix_header

View file

@ -1223,37 +1223,40 @@ contains
subroutine print_columns()
if (entropy_on) then
if (cmfd_run) then
message = " Bat./Gen. k Entropy Average k CMFD k CMFD Ent"
call write_message(1)
message = " ========= ======== ======== ==================== ======== ========"
call write_message(1)
else
message = " Bat./Gen. k Entropy Average k"
call write_message(1)
message = " ========= ======== ======== ===================="
call write_message(1)
end if
else
if (cmfd_run) then
message = " Bat./Gen. k Average k CMFD k"
call write_message(1)
message = " ========= ======== ==================== ========"
call write_message(1)
else
message = " Bat./Gen. k Average k"
call write_message(1)
message = " ========= ======== ===================="
call write_message(1)
end if
write(UNIT=ou, FMT='(2X,A9,3X)', ADVANCE='NO') "Bat./Gen."
write(UNIT=ou, FMT='(A8,3X)', ADVANCE='NO') " k "
if (entropy_on) write(UNIT=ou, FMT='(A8,3X)', ADVANCE='NO') "Entropy "
write(UNIT=ou, FMT='(A20,3X)', ADVANCE='NO') " Average k "
if (cmfd_run) then
write(UNIT=ou, FMT='(A8,3X)', ADVANCE='NO') " CMFD k "
select case(trim(cmfd_display))
case('entropy')
write(UNIT=ou, FMT='(A8,3X)', ADVANCE='NO') "CMFD Ent"
case('balance')
write(UNIT=ou, FMT='(A8,3X)', ADVANCE='NO') "RMS Bal "
case('source')
write(UNIT=ou, FMT='(A8,3X)', ADVANCE='NO') "RMS Src "
case('dominance')
write(UNIT=ou, FMT='(A8,3X)', ADVANCE='NO') "Dom Rat "
end select
end if
write(UNIT=ou, FMT=*)
write(UNIT=ou, FMT='(2X,A9,3X)', ADVANCE='NO') "========="
write(UNIT=ou, FMT='(A8,3X)', ADVANCE='NO') "========"
if (entropy_on) write(UNIT=ou, FMT='(A8,3X)', ADVANCE='NO') "========"
write(UNIT=ou, FMT='(A20,3X)', ADVANCE='NO') "===================="
if (cmfd_run) then
write(UNIT=ou, FMT='(A8,3X)', ADVANCE='NO') "========"
if (cmfd_display /= '') &
write(UNIT=ou, FMT='(A8,3X)', ADVANCE='NO') "========"
end if
write(UNIT=ou, FMT=*)
end subroutine print_columns
!===============================================================================
! PRINT_GENERATION displays information for a generation of neutrons. For now,
! if the user has entropy on, it will print out the entropy
! PRINT_GENERATION displays information for a generation of neutrons.
!===============================================================================
subroutine print_generation()
@ -1303,14 +1306,25 @@ contains
write(UNIT=OUTPUT_UNIT, FMT='(23X)', ADVANCE='NO')
end if
! write out cmfd keff if it is active
if (cmfd_on) write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') &
cmfd % keff
! write out cmfd entopy
if (cmfd_on .and. entropy_on) write(UNIT=OUTPUT_UNIT, &
FMT='(3X, F8.5)', ADVANCE='NO') cmfd % entropy
! write out cmfd keff if it is active and other display info
if (cmfd_on) then
write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') &
cmfd % k_cmfd(current_batch)
select case(trim(cmfd_display))
case('entropy')
write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') &
cmfd % entropy(current_batch)
case('balance')
write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') &
cmfd % balance(current_batch)
case('source')
write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') &
cmfd % src_cmp(current_batch)
case('dominance')
write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') &
cmfd % dom(current_batch)
end select
end if
! next line
write(UNIT=OUTPUT_UNIT, FMT=*)
@ -1405,13 +1419,17 @@ contains
write(ou,100) "Total time in simulation", time_inactive % elapsed + &
time_active % elapsed
write(ou,100) " Time in transport only", time_transport % elapsed
if(cmfd_run) write(ou,100) "Total CMFD time", time_cmfd % elapsed
write(ou,100) " Time in inactive batches", time_inactive % elapsed
write(ou,100) " Time in active batches", time_active % elapsed
write(ou,100) " Time synchronizing fission bank", time_bank % elapsed
write(ou,100) " Sampling source sites", time_bank_sample % elapsed
write(ou,100) " SEND/RECV source sites", time_bank_sendrecv % elapsed
write(ou,100) " Time accumulating tallies", time_tallies % elapsed
if (cmfd_run) write(ou,100) " Time in CMFD", time_cmfd % elapsed
if (cmfd_run) write(ou,100) " Building matrices", &
time_cmfdbuild % elapsed
if (cmfd_run) write(ou,100) " Solving matrices", &
time_cmfdsolve % elapsed
write(ou,100) "Total time for finalization", time_finalize % elapsed
write(ou,100) "Total time elapsed", time_total % elapsed

294
src/solver_interface.F90 Normal file
View file

@ -0,0 +1,294 @@
module solver_interface
use error, only: fatal_error
use global, only: message
use matrix_header, only: Matrix
use vector_header, only: Vector
#ifdef PETSC
use petscksp
use petscsnes
#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
#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)
call KSPSetTolerances(self % ksp_, rtol, atol, &
PETSC_DEFAULT_DOUBLE_PRECISION, PETSC_DEFAULT_INTEGER, petsc_err)
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
call KSPSetOperators(self % ksp_, mat_in % petsc_mat, prec_mat % petsc_mat, &
SAME_NONZERO_PATTERN, petsc_err)
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
!===============================================================================
subroutine petsc_jfnk_compute_jacobian(snes_, x, jac_mf, jac_prec, flag, &
ctx, ierr)
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.
integer, intent(inout) :: flag ! unused madatory flag
type(Jfnk_ctx), intent(inout) :: ctx ! JFNK context instance
integer, intent(inout) :: ierr ! error code
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

@ -20,6 +20,10 @@ module state_point
use output_interface
use tally_header, only: TallyObject
#ifdef MPI
use mpi
#endif
implicit none
type(BinaryOutput) :: sp ! statepoint/source output file
@ -96,6 +100,28 @@ contains
call sp % write_data(k_col_tra, "k_col_tra")
call sp % write_data(k_abs_tra, "k_abs_tra")
call sp % write_data(k_combined, "k_combined", length=2)
! Write out CMFD info
if (cmfd_on) then
call sp % write_data(1, "cmfd_on")
call sp % write_data(cmfd % indices, "indicies", length=4, group="cmfd")
call sp % write_data(cmfd % k_cmfd, "k_cmfd", length=current_batch, &
group="cmfd")
call sp % write_data(cmfd % cmfd_src, "cmfd_src", &
length=(/cmfd % indices(4), cmfd % indices(1), &
cmfd % indices(2), cmfd % indices(3)/), &
group="cmfd")
call sp % write_data(cmfd % entropy, "cmfd_entropy", &
length=current_batch, group="cmfd")
call sp % write_data(cmfd % balance, "cmfd_balance", &
length=current_batch, group="cmfd")
call sp % write_data(cmfd % dom, "cmfd_dominance", &
length = current_batch, group="cmfd")
call sp % write_data(cmfd % src_cmp, "cmfd_srccmp", &
length = current_batch, group="cmfd")
else
call sp % write_data(0, "cmfd_on")
end if
end if
! Write number of meshes
@ -507,6 +533,28 @@ contains
! Take maximum of statepoint n_inactive and input n_inactive
n_inactive = max(n_inactive, int_array(1))
! Read in to see if CMFD was on
call sp % read_data(int_array(1), "cmfd_on")
! Write out CMFD info
if (int_array(1) == 1) then
call sp % read_data(cmfd % indices, "indicies", length=4, group="cmfd")
call sp % read_data(cmfd % k_cmfd, "k_cmfd", length=restart_batch, &
group="cmfd")
call sp % read_data(cmfd % cmfd_src, "cmfd_src", &
length=(/cmfd % indices(4), cmfd % indices(1), &
cmfd % indices(2), cmfd % indices(3)/), &
group="cmfd")
call sp % read_data(cmfd % entropy, "cmfd_entropy", &
length=restart_batch, group="cmfd")
call sp % read_data(cmfd % balance, "cmfd_balance", &
length=restart_batch, group="cmfd")
call sp % read_data(cmfd % dom, "cmfd_dominance", &
length = restart_batch, group="cmfd")
call sp % read_data(cmfd % src_cmp, "cmfd_srccmp", &
length = restart_batch, group="cmfd")
end if
end if
! Read number of meshes

View file

@ -13,27 +13,22 @@
<component name="energy" type="double-array" />
</typedef>
<variable name="mesh_" tag="mesh" type="mesh_xml" />
<variable name="norm_" tag="norm" type="double" default="1.0" />
<variable name="feedback_" tag="feedback" type="word" length="5" />
<variable name="n_cmfd_procs_" tag="n_procs" type="integer" default="1" />
<variable name="reset_" tag="reset" type="word" length="5" />
<variable name="balance_" tag="balance" type="word" length="5" />
<variable name="downscatter_" tag="downscatter" type="word" length="5" />
<variable name="run_2grp_" tag="run_2grp" type="word" length="5" />
<variable name="solver_" tag="solver" type="word" default="'power'" length="250" />
<variable name="snes_monitor_" tag="snes_monitor" type="word" length="5" />
<variable name="ksp_monitor_" tag="ksp_monitor" type="word" length="5" />
<variable name="power_monitor_" tag="power_monitor" type="word" length="5" />
<variable name="write_balance_" tag="write_balance" type="word" length="5" />
<variable name="write_matrices_" tag="write_matrices" type="word" length="5" />
<variable name="run_adjoint_" tag="run_adjoint" type="word" length="5" />
<variable name="write_hdf5_" tag="write_hdf5" type="word" length="5" />
<variable name="begin_" tag="begin" type="integer" default="1" />
<variable name="inactive_" tag="inactive" type="word" length="5" />
<variable name="active_flush_" tag="active_flush" type="integer" default="0" />
<variable name="keff_tol_" tag="keff_tol" type="double" default="0.005" />
<variable name="inactive_flush_" tag="inactive_flush" type="integer" default="9999" />
<variable name="num_flushes_" tag="num_flushes" type="integer" default="1" />
<variable name="mesh_" tag="mesh" type="mesh_xml" />
<variable name="norm_" tag="norm" type="double" default="1.0" />
<variable name="feedback_" tag="feedback" type="word" length="5" />
<variable name="reset_" tag="reset" type="word" length="5" />
<variable name="downscatter_" tag="downscatter" type="word" length="5" />
<variable name="solver_" tag="solver" type="word" default="'power'" length="25" />
<variable name="snes_monitor_" tag="snes_monitor" type="word" length="5" />
<variable name="ksp_monitor_" tag="ksp_monitor" type="word" length="5" />
<variable name="power_monitor_" tag="power_monitor" type="word" length="5" />
<variable name="write_matrices_" tag="write_matrices" type="word" length="5" />
<variable name="run_adjoint_" tag="run_adjoint" type="word" length="5" />
<variable name="begin_" tag="begin" type="integer" default="1" />
<variable name="inactive_" tag="inactive" type="word" length="5" />
<variable name="active_flush_" tag="active_flush" type="integer" default="0" />
<variable name="inactive_flush_" tag="inactive_flush" type="integer" default="9999" />
<variable name="num_flushes_" tag="num_flushes" type="integer" default="1" />
<variable name="display_" tag="display" type="word" default="''" length="25" />
</template>

View file

@ -11,7 +11,8 @@ for src in glob.iglob('*.F90'):
d = re.findall(r'\n\s*use\s+(\w+)',
open(src, 'r').read())
for name in d:
if name in ['mpi', 'hdf5', 'h5lt']:
if name in ['mpi', 'hdf5', 'h5lt', 'petscsys', 'petscmat', 'petscksp',
'petscsnes', 'petscvec']:
continue
if name.startswith('xml_data_'):
name = name.replace('xml_data_', 'templates/')
@ -20,7 +21,7 @@ for src in glob.iglob('*.F90'):
dependencies[module] = sorted(list(deps))
for module in dependencies.keys():
for module in sorted(dependencies.keys()):
for dep in dependencies[module]:
print("{0}.o: {1}.o".format(module, dep))
print('')

View file

@ -152,6 +152,8 @@ class StatePoint(object):
# Read statepoint revision
self.revision = self._get_int(path='revision')[0]
if self.revision != 10:
raise Exception('Statepoint Revision is not consistent.')
# Read OpenMC version
if self._hdf5:
@ -186,11 +188,29 @@ class StatePoint(object):
self.current_batch*self.gen_per_batch, path='k_generation')
self.entropy = self._get_double(
self.current_batch*self.gen_per_batch, path='entropy')
if self.revision >= 8:
self.k_col_abs = self._get_double(path='k_col_abs')[0]
self.k_col_tra = self._get_double(path='k_col_tra')[0]
self.k_abs_tra = self._get_double(path='k_abs_tra')[0]
self.k_combined = self._get_double(2, path='k_combined')
self.k_col_abs = self._get_double(path='k_col_abs')[0]
self.k_col_tra = self._get_double(path='k_col_tra')[0]
self.k_abs_tra = self._get_double(path='k_abs_tra')[0]
self.k_combined = self._get_double(2, path='k_combined')
# Read CMFD information
cmfd_present = self._get_int(path='cmfd_on')[0]
if cmfd_present == 1:
self.cmfd_indices = self._get_int(4, path='cmfd/indicies')
self.k_cmfd = self._get_double(self.current_batch,
path='cmfd/k_cmfd')
self.cmfd_src = self._get_double_array(np.product(self.cmfd_indices),
path='cmfd/cmfd_src')
self.cmfd_src = np.reshape(self.cmfd_src,
tuple(self.cmfd_indices), order='F')
self.cmfd_entropy = self._get_double(self.current_batch,
path='cmfd/cmfd_entropy')
self.cmfd_balance = self._get_double(self.current_batch,
path='cmfd/cmfd_balance')
self.cmfd_dominance = self._get_double(self.current_batch,
path='cmfd/cmfd_dominance')
self.cmfd_srccmp = self._get_double(self.current_batch,
path='cmfd/cmfd_srccmp')
# Read number of meshes
n_meshes = self._get_int(path='tallies/n_meshes')[0]
@ -581,6 +601,12 @@ class StatePoint(object):
else:
return [float(v) for v in self._get_data(n, 'd', 8)]
def _get_double_array(self, n=1, path=None):
if self._hdf5:
return self._f[path].value
else:
return self._get_data(n, 'd', 8)
def _get_string(self, n=1, path=None):
if self._hdf5:
return str(self._f[path].value)

126
src/vector_header.F90 Normal file
View file

@ -0,0 +1,126 @@
module vector_header
use constants, only: ZERO
#ifdef PETSC
use petscvec
#endif
implicit none
private
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 :: setup_petsc => vector_setup_petsc
procedure :: write_petsc_binary => vector_write_petsc_binary
end type Vector
integer :: petsc_err ! petsc error code
contains
!===============================================================================
! VECTOR_CREATE allocates and initializes a vector
!===============================================================================
subroutine vector_create(self, n)
integer, intent(in) :: n ! size of vector
class(Vector), intent(inout), target :: self ! vector instance
! Preallocate vector
if (.not.allocated(self % data)) allocate(self % data(n))
self % val => self % data(1:n)
! Set n
self % n = n
! Initialize to zero
self % val = ZERO
! Petsc is default not active
self % petsc_active = .false.
end subroutine vector_create
!===============================================================================
! VECTOR_DESTROY deallocates all space associated with a vector
!===============================================================================
subroutine vector_destroy(self)
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)
end subroutine vector_destroy
!===============================================================================
! VECTOR_ADD_VALUE adds a value to the vector
!===============================================================================
subroutine vector_add_value(self, idx, val)
integer, intent(in) :: idx ! index location in vector
real(8), intent(in) :: val ! value to add
class(Vector), intent(inout) :: self ! vector instance
self % val(idx) = val
end subroutine vector_add_value
!===============================================================================
! VECTOR_SETUP_PETSC links the data to a PETSc vector
!===============================================================================
subroutine vector_setup_petsc(self)
class(Vector), intent(inout) :: self ! vector instance
! Link to PETSc
#ifdef PETSC
call VecCreateSeqWithArray(PETSC_COMM_WORLD, 1, self % n, self % val, &
self % petsc_vec, petsc_err)
#endif
! Set that PETSc is now active
self % petsc_active = .true.
end subroutine vector_setup_petsc
!===============================================================================
! VECTOR_WRITE_PETSC_BINARY writes the PETSc vector to a binary file
!===============================================================================
subroutine vector_write_petsc_binary(self, filename)
character(*), intent(in) :: filename ! name of file to write to
class(Vector), intent(in) :: self ! vector instance
#ifdef PETSC
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)
#endif
end subroutine vector_write_petsc_binary
end module vector_header

View file

@ -22,7 +22,7 @@ def run_compile():
os.remove(exe)
# run compile test
result = nose.run(argv=['run_tests.py', 'test_compile'] + flags)
result = nose.run(argv=['nosetests', 'test_compile'] + flags)
if not result:
print('Did not pass compile tests.')
results.append(('compile', result))
@ -35,7 +35,7 @@ def run_suite(name=None, mpi=False):
# Set arguments list. Note that the first argument is a dummy argument (the
# script name). It's not actually recursively calling run_tests.py
argv = ['run_tests.py', '--exclude', 'test_compile'] + flags
argv = ['nosetests', '--exclude', 'test_compile'] + flags
# Add MPI plugin if set
if mpi:
@ -48,11 +48,13 @@ def run_suite(name=None, mpi=False):
os.chdir(pwd)
os.rename(pwd + '/../src/openmc-' + name, pwd + '/../src/openmc')
result = nose.run(argv=argv, addplugins=plugins)
finally:
os.rename(pwd + '/../src/openmc', pwd + '/../src/openmc-' + name)
if not result:
print('Did not pass ' + name + ' tests')
results.append((name, result))
except OSError:
result = False
print('No OpenMC executable found for ' + name + ' tests')
if not result:
print('Did not pass ' + name + ' tests')
results.append((name, result))
# set mpiexec path
if 'COMPILER' in os.environ:
@ -67,23 +69,77 @@ sys.path.append(pwd)
# Set list of tests, either default or from command line
flags = []
tests = ['compile', 'normal', 'debug', 'optimize', 'mpi', 'omp', 'hdf5',
'petsc', 'mpi-omp', 'omp-hdf5', 'phdf5', 'phdf5-petsc',
'phdf5-petsc-optimize', 'omp-phdf5-petsc-optimize']
tests = ['compile', 'normal', 'debug', 'optimize',
'omp', 'omp-debug', 'omp-optimize',
'hdf5', 'hdf5-debug', 'hdf5-optimize',
'omp-hdf5', 'omp-hdf5-debug', 'omp-hdf5-optimize',
'mpi', 'mpi-debug', 'mpi-optimize',
'mpi-omp', 'mpi-omp-debug', 'mpi-omp-optimize',
'phdf5', 'phdf5-debug', 'phdf5-optimize',
'phdf5-omp', 'phdf5-omp-debug', 'phdf5-omp-optimize',
'petsc', 'petsc-debug', 'petsc-optimize',
'phdf5-petsc', 'phdf5-petsc-debug', 'phdf5-petsc-optimize',
'omp-phdf5-petsc', 'omp-phdf5-petsc-debug',
'omp-phdf5-petsc-optimize']
if len(sys.argv) > 1:
flags = [i for i in sys.argv[1:] if i.startswith('-')]
tests_ = [i for i in sys.argv[1:] if not i.startswith('-')]
tests = tests_ if tests_ else tests
# Check for special subsets of tests
tests__ = []
for i in tests_:
# All tests will run all the tests except for compile unless
# it is also specified on the command line. Note that specifying
# compile all-tests is the same as not specifying any args
if i == 'all-tests':
tests__ = tests
try:
idx = tests_.index('compile') # check for compile test
except ValueError:
del tests__[0]
finally:
break # don't need to check for anything else
# This checks for any subsets of tests. The string after
# all-XXXX will be used to search through all tests.
# Specifying XXXX=normal will run tests that don't contain
# debug or optimize substring.
if i.startswith('all-'):
suffix = i.split('all-')[1]
if suffix == 'normal':
for j in tests:
if j.rfind('debug') == -1 and \
j.rfind('optimize') == -1:
tests__.append(j)
else:
for j in tests:
if j.rfind(suffix) != -1:
if suffix == 'omp' and j == 'compile':
continue
tests__.append(j)
else:
tests__.append(i) # append specific test (e.g., mpi-debug)
tests = tests__ if tests__ else tests
# Run tests
results = []
for name in tests:
if name == 'compile':
run_compile()
elif name in ['normal', 'debug', 'optimize', 'omp', 'hdf5', 'omp-hdf5']:
elif name in ['normal', 'debug', 'optimize',
'hdf5', 'hdf5-debug', 'hdf5-optimize',
'omp', 'omp-debug', 'omp-optimize',
'omp-hdf5', 'omp-hdf5-debug', 'omp-hdf5-optimize']:
run_suite(name=name)
elif name in ['mpi', 'mpi-omp' 'phdf5', 'petsc', 'phdf5-petsc',
'phdf5-petsc-optimize', 'omp-phdf5-petsc-optimize']:
elif name in ['mpi', 'mpi-debug', 'mpi-optimize',
'mpi-omp', 'mpi-omp-debug', 'mpi-omp-optimize',
'phdf5', 'phdf5-debug', 'phdf5-optimize',
'phdf5-omp', 'phdf5-omp-debug', 'phdf5-omp-optimize',
'petsc', 'petsc-debug', 'petsc-optimize',
'phdf5-petsc', 'phdf5-petsc-debug', 'phdf5-petsc-optimize',
'omp-phdf5-petsc', 'omp-phdf5-petsc-debug',
'omp-phdf5-petsc-optimize']:
run_suite(name=name, mpi=True)
# print out summary of results

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_created_statepoint():

View file

@ -34,7 +34,7 @@ results4 = np.reshape(results4, size4)
# 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])
@ -42,17 +42,41 @@ 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 += "{0:12.6E}\n".format(item)
outstr += 'tally 2:\n'
for item in results2:
outstr += "{0:12.6E}\n".format(item)
outstr += "{0:12.6E}\n".format(item)
outstr += 'tally 3:\n'
for item in results3:
outstr += "{0:12.6E}\n".format(item)
outstr += "{0:12.6E}\n".format(item)
outstr += 'tally 4:\n'
for item in results4:
outstr += "{0:12.6E}\n".format(item)
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:
with open('results_test.dat', 'w') as fh:
fh.write(outstr)

View file

@ -1,14 +1,14 @@
k-combined:
1.167124E+00 1.217343E-02
1.167124E+00 1.217344E-02
tally 1:
1.126891E+01
1.275755E+01
2.086556E+01
4.369501E+01
1.275756E+01
2.086557E+01
4.369504E+01
2.853629E+01
8.169559E+01
3.404506E+01
1.165019E+02
8.169563E+01
3.404507E+01
1.165020E+02
3.723908E+01
1.389881E+02
3.760048E+01
@ -16,86 +16,86 @@ tally 1:
3.455469E+01
1.197186E+02
2.837503E+01
8.086291E+01
2.151789E+01
4.648315E+01
8.086287E+01
2.151788E+01
4.648312E+01
1.194541E+01
1.432331E+01
1.432330E+01
tally 2:
2.292173E+01
2.650792E+01
1.598183E+01
1.288692E+01
2.194143E+00
2.510557E-01
4.217816E+01
8.938583E+01
2.292174E+01
2.650794E+01
1.598184E+01
1.288693E+01
2.194144E+00
2.510560E-01
4.217817E+01
8.938588E+01
2.974466E+01
4.447038E+01
3.939269E+00
7.882413E-01
5.742289E+01
4.447040E+01
3.939271E+00
7.882419E-01
5.742290E+01
1.654599E+02
4.067581E+01
8.308454E+01
5.566955E+00
1.570072E+00
6.810586E+01
4.067582E+01
8.308457E+01
5.566956E+00
1.570073E+00
6.810587E+01
2.331019E+02
4.831543E+01
4.831544E+01
1.174673E+02
6.244999E+00
1.967938E+00
6.245000E+00
1.967939E+00
7.258131E+01
2.646756E+02
5.185146E+01
1.351425E+02
6.622488E+00
1.351426E+02
6.622489E+00
2.231695E+00
7.246115E+01
2.641295E+02
5.159542E+01
5.159541E+01
1.339441E+02
6.702239E+00
2.277825E+00
6.716397E+01
6.716396E+01
2.266815E+02
4.765063E+01
1.141627E+02
6.384098E+00
4.765062E+01
1.141626E+02
6.384097E+00
2.063760E+00
5.549542E+01
1.545361E+02
5.549541E+01
1.545360E+02
3.940819E+01
7.794880E+01
5.183375E+00
1.369483E+00
4.155685E+01
8.679439E+01
2.926744E+01
4.307145E+01
3.938852E+00
7.877949E-01
2.335180E+01
2.758241E+01
1.622683E+01
1.331996E+01
2.248921E+00
2.594934E-01
7.794877E+01
5.183373E+00
1.369482E+00
4.155684E+01
8.679434E+01
2.926743E+01
4.307143E+01
3.938851E+00
7.877944E-01
2.335179E+01
2.758239E+01
1.622682E+01
1.331995E+01
2.248920E+00
2.594933E-01
tally 3:
1.537807E+01
1.193746E+01
1.044189E+00
5.627797E-02
2.861614E+01
4.118322E+01
4.118324E+01
1.846208E+00
1.729007E-01
3.921462E+01
7.723576E+01
3.921463E+01
7.723579E+01
2.427723E+00
2.983153E-01
4.653207E+01
4.653208E+01
1.089818E+02
3.128394E+00
4.945409E-01
@ -111,15 +111,15 @@ tally 3:
1.061019E+02
2.936424E+00
4.356307E-01
3.783939E+01
7.188489E+01
3.783938E+01
7.188486E+01
2.485680E+00
3.117618E-01
2.819358E+01
3.997588E+01
3.997585E+01
1.875908E+00
1.772000E-01
1.562421E+01
1.562420E+01
1.236020E+01
1.013140E+00
5.260420E-02
@ -160,8 +160,8 @@ tally 4:
0.000000E+00
0.000000E+00
0.000000E+00
3.126365E+00
4.930682E-01
3.126366E+00
4.930686E-01
0.000000E+00
0.000000E+00
0.000000E+00
@ -208,10 +208,10 @@ tally 4:
0.000000E+00
0.000000E+00
0.000000E+00
5.456410E+00
1.499416E+00
2.727755E+00
3.788867E-01
5.456412E+00
1.499417E+00
2.727756E+00
3.788869E-01
0.000000E+00
0.000000E+00
0.000000E+00
@ -256,10 +256,10 @@ tally 4:
0.000000E+00
0.000000E+00
0.000000E+00
7.453480E+00
2.790159E+00
5.240103E+00
1.386895E+00
7.453481E+00
2.790160E+00
5.240104E+00
1.386896E+00
0.000000E+00
0.000000E+00
0.000000E+00
@ -304,10 +304,10 @@ tally 4:
0.000000E+00
0.000000E+00
0.000000E+00
8.741563E+00
3.830225E+00
7.137968E+00
2.555042E+00
8.741565E+00
3.830226E+00
7.137969E+00
2.555043E+00
0.000000E+00
0.000000E+00
0.000000E+00
@ -354,8 +354,8 @@ tally 4:
0.000000E+00
9.221731E+00
4.266573E+00
8.474116E+00
3.604647E+00
8.474118E+00
3.604648E+00
0.000000E+00
0.000000E+00
0.000000E+00
@ -401,9 +401,9 @@ tally 4:
0.000000E+00
0.000000E+00
9.054937E+00
4.124537E+00
9.066119E+00
4.133833E+00
4.124536E+00
9.066120E+00
4.133834E+00
0.000000E+00
0.000000E+00
0.000000E+00
@ -448,8 +448,8 @@ tally 4:
0.000000E+00
0.000000E+00
0.000000E+00
8.362253E+00
3.517245E+00
8.362251E+00
3.517244E+00
9.086548E+00
4.156118E+00
0.000000E+00
@ -496,10 +496,10 @@ tally 4:
0.000000E+00
0.000000E+00
0.000000E+00
6.871573E+00
2.369190E+00
8.399584E+00
3.549784E+00
6.871571E+00
2.369189E+00
8.399583E+00
3.549783E+00
0.000000E+00
0.000000E+00
0.000000E+00
@ -544,10 +544,10 @@ tally 4:
0.000000E+00
0.000000E+00
0.000000E+00
4.998103E+00
1.258194E+00
7.145296E+00
2.568115E+00
4.998101E+00
1.258193E+00
7.145294E+00
2.568114E+00
0.000000E+00
0.000000E+00
0.000000E+00
@ -592,10 +592,10 @@ tally 4:
0.000000E+00
0.000000E+00
0.000000E+00
2.716076E+00
3.753568E-01
5.427593E+00
1.484154E+00
2.716075E+00
3.753566E-01
5.427591E+00
1.484153E+00
0.000000E+00
0.000000E+00
0.000000E+00
@ -642,8 +642,8 @@ tally 4:
0.000000E+00
0.000000E+00
0.000000E+00
3.056645E+00
4.702349E-01
3.056644E+00
4.702345E-01
0.000000E+00
0.000000E+00
0.000000E+00
@ -652,3 +652,124 @@ tally 4:
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.158665E+00
1.173348E+00
1.180514E+00
1.168646E+00
1.152485E+00
1.154989E+00
1.152285E+00
1.152275E+00
1.149476E+00
1.154912E+00
1.164393E+00
1.170145E+00
1.169161E+00
1.167966E+00
1.170153E+00
1.170989E+00
cmfd entropy
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
3.229801E+00
3.226763E+00
3.225011E+00
3.230316E+00
3.232043E+00
3.228854E+00
3.227439E+00
3.228206E+00
3.226749E+00
3.222635E+00
3.218882E+00
3.219196E+00
3.218057E+00
3.220324E+00
3.219099E+00
3.219847E+00
cmfd balance
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
6.244292E-03
5.210436E-03
4.212224E-03
3.321578E-03
3.067254E-03
3.076634E-03
2.601637E-03
2.666441E-03
2.490817E-03
2.571411E-03
3.013819E-03
2.802836E-03
2.725356E-03
2.475281E-03
2.279123E-03
2.110816E-03
cmfd dominance ratio
0.000E+00
0.000E+00
0.000E+00
0.000E+00
5.509E-01
5.473E-01
5.461E-01
5.519E-01
5.540E-01
5.514E-01
5.505E-01
5.499E-01
5.486E-01
5.450E-01
5.438E-01
5.422E-01
5.417E-01
5.437E-01
5.430E-01
5.428E-01
cmfd openmc source comparison
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
1.136193E-02
1.059214E-02
9.675100E-03
7.296572E-03
5.812484E-03
6.138879E-03
6.070297E-03
5.791570E-03
6.020032E-03
7.287183E-03
7.068423E-03
6.380521E-03
5.828010E-03
4.741191E-03
4.848114E-03
4.686772E-03
cmfd source
4.367444E-02
7.866919E-02
1.050799E-01
1.366581E-01
1.421824E-01
1.393453E-01
1.269612E-01
1.060419E-01
7.920472E-02
4.218277E-02

View file

@ -14,28 +14,31 @@ pwd = os.path.dirname(__file__)
skipAll = False
def setup():
os.putenv('PWD', pwd)
os.chdir(pwd)
global skipAll
skipAll = False
def test_run():
openmc_path = pwd + '/../../src/openmc'
if int(NoseMPI.mpi_np) > 0:
proc = Popen([NoseMPI.mpi_exec, '-np', NoseMPI.mpi_np, openmc_path],
stderr=STDOUT, stdout=PIPE)
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
output = proc.communicate()[0]
print(output)
returncode = proc.returncode
if 'CMFD is not available' in output:
global skipAll
skipAll = True
raise SkipTest
assert returncode == 0
def test_created_statepoint():
if skipAll:
raise SkipTest
@ -43,11 +46,13 @@ def test_created_statepoint():
assert len(statepoint) == 1
assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5')
def test_output_exists():
if skipAll:
raise SkipTest
assert os.path.exists(pwd + '/tallies.out')
def test_results():
if skipAll:
raise SkipTest
@ -55,9 +60,10 @@ def test_results():
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')
os.rename('results_test.dat', 'results_error.dat')
assert compare
def teardown():
output = glob.glob(pwd + '/statepoint.20.*')
output.append(pwd + '/tallies.out')
@ -65,4 +71,3 @@ def teardown():
for f in output:
if os.path.exists(f):
os.remove(f)

View file

@ -0,0 +1,16 @@
<?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

@ -0,0 +1,43 @@
<?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

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

View file

@ -0,0 +1,82 @@
#!/usr/bin/env python
import sys
import numpy as np
# import statepoint
sys.path.append('../../src/utils')
import statepoint
# read in statepoint file
if len(sys.argv) > 1:
sp = statepoint.StatePoint(sys.argv[1])
else:
sp = statepoint.StatePoint('statepoint.20.binary')
sp.read_results()
# extract tally results and convert to vector
results1 = sp.tallies[0].results
shape1 = results1.shape
size1 = (np.product(shape1))
results1 = np.reshape(results1, size1)
results2 = sp.tallies[1].results
shape2 = results2.shape
size2 = (np.product(shape2))
results2 = np.reshape(results2, size2)
results3 = sp.tallies[2].results
shape3 = results3.shape
size3 = (np.product(shape3))
results3 = np.reshape(results3, size3)
results4 = sp.tallies[3].results
shape4 = results4.shape
size4 = (np.product(shape4))
results4 = np.reshape(results4, size4)
# 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

@ -0,0 +1,775 @@
k-combined:
1.167124E+00 1.217344E-02
tally 1:
1.126891E+01
1.275756E+01
2.086557E+01
4.369504E+01
2.853629E+01
8.169563E+01
3.404507E+01
1.165020E+02
3.723908E+01
1.389881E+02
3.760048E+01
1.416691E+02
3.455469E+01
1.197186E+02
2.837503E+01
8.086287E+01
2.151788E+01
4.648312E+01
1.194541E+01
1.432330E+01
tally 2:
2.292174E+01
2.650794E+01
1.598184E+01
1.288693E+01
2.194144E+00
2.510560E-01
4.217817E+01
8.938588E+01
2.974466E+01
4.447040E+01
3.939271E+00
7.882419E-01
5.742290E+01
1.654599E+02
4.067582E+01
8.308457E+01
5.566956E+00
1.570073E+00
6.810587E+01
2.331019E+02
4.831544E+01
1.174673E+02
6.245000E+00
1.967939E+00
7.258131E+01
2.646756E+02
5.185146E+01
1.351426E+02
6.622489E+00
2.231695E+00
7.246115E+01
2.641295E+02
5.159541E+01
1.339441E+02
6.702239E+00
2.277825E+00
6.716396E+01
2.266815E+02
4.765062E+01
1.141626E+02
6.384097E+00
2.063760E+00
5.549541E+01
1.545360E+02
3.940819E+01
7.794877E+01
5.183373E+00
1.369482E+00
4.155684E+01
8.679434E+01
2.926743E+01
4.307143E+01
3.938851E+00
7.877944E-01
2.335179E+01
2.758239E+01
1.622682E+01
1.331995E+01
2.248920E+00
2.594933E-01
tally 3:
1.537807E+01
1.193746E+01
1.044189E+00
5.627797E-02
2.861614E+01
4.118324E+01
1.846208E+00
1.729007E-01
3.921463E+01
7.723579E+01
2.427723E+00
2.983153E-01
4.653208E+01
1.089818E+02
3.128394E+00
4.945409E-01
4.989896E+01
1.251503E+02
3.242191E+00
5.330743E-01
4.967378E+01
1.241657E+02
3.184751E+00
5.164524E-01
4.592989E+01
1.061019E+02
2.936424E+00
4.356307E-01
3.783938E+01
7.188486E+01
2.485680E+00
3.117618E-01
2.819358E+01
3.997585E+01
1.875908E+00
1.772000E-01
1.562420E+01
1.236020E+01
1.013140E+00
5.260420E-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.126366E+00
4.930686E-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.456412E+00
1.499417E+00
2.727756E+00
3.788869E-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.453481E+00
2.790160E+00
5.240104E+00
1.386896E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
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.741565E+00
3.830226E+00
7.137969E+00
2.555043E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
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.221731E+00
4.266573E+00
8.474118E+00
3.604648E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
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.054937E+00
4.124536E+00
9.066120E+00
4.133834E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
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.362251E+00
3.517244E+00
9.086548E+00
4.156118E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
6.871571E+00
2.369189E+00
8.399583E+00
3.549783E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
4.998101E+00
1.258193E+00
7.145294E+00
2.568114E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
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.716075E+00
3.753566E-01
5.427591E+00
1.484153E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
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.056644E+00
4.702345E-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.158665E+00
1.173348E+00
1.180514E+00
1.168646E+00
1.152485E+00
1.154989E+00
1.152285E+00
1.152275E+00
1.149476E+00
1.154912E+00
1.164393E+00
1.170145E+00
1.169161E+00
1.167966E+00
1.170153E+00
1.170989E+00
cmfd entropy
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
3.229801E+00
3.226763E+00
3.225011E+00
3.230316E+00
3.232043E+00
3.228854E+00
3.227439E+00
3.228206E+00
3.226749E+00
3.222635E+00
3.218882E+00
3.219196E+00
3.218057E+00
3.220324E+00
3.219099E+00
3.219847E+00
cmfd balance
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
6.244292E-03
5.210436E-03
4.212224E-03
3.321578E-03
3.067254E-03
3.076634E-03
2.601637E-03
2.666441E-03
2.490817E-03
2.571411E-03
3.013819E-03
2.802836E-03
2.725356E-03
2.475281E-03
2.279123E-03
2.110816E-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
1.136193E-02
1.059214E-02
9.675100E-03
7.296571E-03
5.812485E-03
6.138879E-03
6.070297E-03
5.791569E-03
6.020032E-03
7.287183E-03
7.068423E-03
6.380521E-03
5.828010E-03
4.741191E-03
4.848115E-03
4.686773E-03
cmfd source
4.367444E-02
7.866919E-02
1.050799E-01
1.366581E-01
1.421824E-01
1.393453E-01
1.269612E-01
1.060419E-01
7.920472E-02
4.218277E-02

View file

@ -0,0 +1,32 @@
<?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

@ -0,0 +1,16 @@
<?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

@ -0,0 +1,73 @@
#!/usr/bin/env python
import os
from subprocess import Popen, STDOUT, PIPE, call
import filecmp
import glob
from nose.plugins.skip import SkipTest
from nose_mpi import NoseMPI
pwd = os.path.dirname(__file__)
skipAll = False
def setup():
os.putenv('PWD', pwd)
os.chdir(pwd)
global skipAll
skipAll = False
def test_run():
openmc_path = pwd + '/../../src/openmc'
if int(NoseMPI.mpi_np) > 0:
proc = Popen([NoseMPI.mpi_exec, '-np', NoseMPI.mpi_np, openmc_path],
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
output = proc.communicate()[0]
print(output)
returncode = proc.returncode
if 'CMFD is not available' in output:
global skipAll
skipAll = True
raise SkipTest
assert returncode == 0
def test_created_statepoint():
if skipAll:
raise SkipTest
statepoint = glob.glob(pwd + '/statepoint.20.*')
assert len(statepoint) == 1
assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5')
def test_output_exists():
if skipAll:
raise SkipTest
assert os.path.exists(pwd + '/tallies.out')
def test_results():
if skipAll:
raise SkipTest
statepoint = glob.glob(pwd + '/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
def teardown():
output = glob.glob(pwd + '/statepoint.20.*')
output.append(pwd + '/tallies.out')
output.append(pwd + '/results_test.dat')
for f in output:
if os.path.exists(f):
os.remove(f)

View file

@ -34,7 +34,7 @@ results4 = np.reshape(results4, size4)
# 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])
@ -42,17 +42,41 @@ 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 += "{0:12.6E}\n".format(item)
outstr += 'tally 2:\n'
for item in results2:
outstr += "{0:12.6E}\n".format(item)
outstr += "{0:12.6E}\n".format(item)
outstr += 'tally 3:\n'
for item in results3:
outstr += "{0:12.6E}\n".format(item)
outstr += "{0:12.6E}\n".format(item)
outstr += 'tally 4:\n'
for item in results4:
outstr += "{0:12.6E}\n".format(item)
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:
with open('results_test.dat', 'w') as fh:
fh.write(outstr)

View file

@ -652,3 +652,124 @@ tally 4:
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.158665E+00
1.177323E+00
1.179341E+00
1.169660E+00
1.179827E+00
1.180220E+00
1.177699E+00
1.179321E+00
1.180728E+00
1.180996E+00
1.186680E+00
1.181455E+00
1.178393E+00
1.176470E+00
1.172593E+00
1.167042E+00
cmfd entropy
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
3.229801E+00
3.227279E+00
3.224715E+00
3.226705E+00
3.226404E+00
3.227329E+00
3.227407E+00
3.228601E+00
3.230423E+00
3.228093E+00
3.226358E+00
3.226226E+00
3.227324E+00
3.228594E+00
3.229513E+00
3.229820E+00
cmfd balance
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
6.244292E-03
5.300885E-03
5.198276E-03
3.759109E-03
3.994559E-03
3.064086E-03
2.749389E-03
2.388483E-03
2.304428E-03
2.188627E-03
2.269076E-03
2.016776E-03
1.622043E-03
1.470984E-03
1.367144E-03
1.402525E-03
cmfd dominance ratio
0.000E+00
0.000E+00
0.000E+00
0.000E+00
5.509E-01
5.474E-01
5.434E-01
5.469E-01
5.461E-01
5.474E-01
5.486E-01
5.485E-01
5.491E-01
5.474E-01
5.467E-01
5.471E-01
5.480E-01
5.486E-01
5.500E-01
5.502E-01
cmfd openmc source comparison
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
1.136193E-02
1.007349E-02
8.776995E-03
6.729971E-03
6.921055E-03
5.867589E-03
4.366276E-03
4.051899E-03
3.406846E-03
4.043983E-03
4.864903E-03
3.711687E-03
3.124413E-03
2.675579E-03
2.523549E-03
2.724884E-03
cmfd source
4.534769E-02
8.593793E-02
1.059054E-01
1.317554E-01
1.358070E-01
1.352224E-01
1.290227E-01
1.093443E-01
7.943056E-02
4.222668E-02

View file

@ -13,28 +13,31 @@ pwd = os.path.dirname(__file__)
skipAll = False
def setup():
os.putenv('PWD', pwd)
os.chdir(pwd)
global skipAll
skipAll = False
def test_run():
openmc_path = pwd + '/../../src/openmc'
if int(NoseMPI.mpi_np) > 0:
proc = Popen([NoseMPI.mpi_exec, '-np', NoseMPI.mpi_np, openmc_path],
stderr=STDOUT, stdout=PIPE)
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
output = proc.communicate()[0]
print(output)
returncode = proc.returncode
if 'CMFD is not available' in output:
global skipAll
skipAll = True
raise SkipTest
assert returncode == 0
def test_created_statepoint():
if skipAll:
raise SkipTest
@ -42,11 +45,13 @@ def test_created_statepoint():
assert len(statepoint) == 1
assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5')
def test_output_exists():
if skipAll:
raise SkipTest
assert os.path.exists(pwd + '/tallies.out')
def test_results():
if skipAll:
raise SkipTest
@ -54,9 +59,10 @@ def test_results():
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')
os.rename('results_test.dat', 'results_error.dat')
assert compare
def teardown():
output = glob.glob(pwd + '/statepoint.20.*')
output.append(pwd + '/tallies.out')
@ -64,4 +70,3 @@ def teardown():
for f in output:
if os.path.exists(f):
os.remove(f)

View file

@ -60,6 +60,20 @@ def test_mpi():
shutil.move('openmc', 'openmc-mpi')
def test_mpi_debug():
returncode = run(['make', 'distclean'])
returncode = run(['make', compiler, 'MPI=yes', 'DEBUG=yes'])
assert returncode == 0
shutil.move('openmc', 'openmc-mpi-debug')
def test_mpi_optimize():
returncode = run(['make', 'distclean'])
returncode = run(['make', compiler, 'MPI=yes', 'OPTIMIZE=yes'])
assert returncode == 0
shutil.move('openmc', 'openmc-mpi-optimize')
def test_omp():
returncode = run(['make', 'distclean'])
returncode = run(['make', compiler, 'OPENMP=yes'])
@ -67,6 +81,20 @@ def test_omp():
shutil.move('openmc', 'openmc-omp')
def test_omp_debug():
returncode = run(['make', 'distclean'])
returncode = run(['make', compiler, 'OPENMP=yes', 'DEBUG=yes'])
assert returncode == 0
shutil.move('openmc', 'openmc-omp-debug')
def test_omp_optimize():
returncode = run(['make', 'distclean'])
returncode = run(['make', compiler, 'OPENMP=yes', 'OPTIMIZE=yes'])
assert returncode == 0
shutil.move('openmc', 'openmc-omp-optimize')
def test_hdf5():
returncode = run(['make', 'distclean'])
returncode = run(['make', compiler, 'HDF5=yes'])
@ -74,6 +102,41 @@ def test_hdf5():
shutil.move('openmc', 'openmc-hdf5')
def test_hdf5_debug():
returncode = run(['make', 'distclean'])
returncode = run(['make', compiler, 'HDF5=yes', 'DEBUG=yes'])
assert returncode == 0
shutil.move('openmc', 'openmc-hdf5-debug')
def test_hdf5_optimize():
returncode = run(['make', 'distclean'])
returncode = run(['make', compiler, 'HDF5=yes', 'OPTIMIZE=yes'])
assert returncode == 0
shutil.move('openmc', 'openmc-hdf5-optimize')
def test_omp_hdf5():
returncode = run(['make', 'distclean'])
returncode = run(['make', compiler, 'OPENMP=yes', 'HDF5=yes'])
assert returncode == 0
shutil.move('openmc', 'openmc-omp-hdf5')
def test_omp_hdf5_debug():
returncode = run(['make', 'distclean'])
returncode = run(['make', compiler, 'OPENMP=yes', 'HDF5=yes', 'DEBUG=yes'])
assert returncode == 0
shutil.move('openmc', 'openmc-omp-hdf5-debug')
def test_omp_hdf5_optimize():
returncode = run(['make', 'distclean'])
returncode = run(['make', compiler, 'OPENMP=yes', 'HDF5=yes', 'OPTIMIZE=yes'])
assert returncode == 0
shutil.move('openmc', 'openmc-omp-hdf5-optimize')
def test_petsc():
returncode = run(['make', 'distclean'])
returncode = run(['make', compiler, 'MPI=yes', 'PETSC=yes'])
@ -81,6 +144,21 @@ def test_petsc():
shutil.move('openmc', 'openmc-petsc')
def test_petsc_debug():
returncode = run(['make', 'distclean'])
returncode = run(['make', compiler, 'MPI=yes', 'PETSC=yes', 'DEBUG=yes'])
assert returncode == 0
shutil.move('openmc', 'openmc-petsc-debug')
def test_petsc_optimize():
returncode = run(['make', 'distclean'])
returncode = run(['make', compiler, 'MPI=yes', 'PETSC=yes',
'OPTIMIZE=yes'])
assert returncode == 0
shutil.move('openmc', 'openmc-petsc-optimize')
def test_mpi_omp():
returncode = run(['make', 'distclean'])
returncode = run(['make', compiler, 'MPI=yes', 'OPENMP=yes'])
@ -88,11 +166,18 @@ def test_mpi_omp():
shutil.move('openmc', 'openmc-mpi-omp')
def test_omp_hdf5():
def test_mpi_omp_debug():
returncode = run(['make', 'distclean'])
returncode = run(['make', compiler, 'MPI=yes', 'OPENMP=yes', 'HDF5=yes'])
returncode = run(['make', compiler, 'MPI=yes', 'OPENMP=yes', 'DEBUG=yes'])
assert returncode == 0
shutil.move('openmc', 'openmc-omp-hdf5')
shutil.move('openmc', 'openmc-mpi-omp-debug')
def test_mpi_omp_optimize():
returncode = run(['make', 'distclean'])
returncode = run(['make', compiler, 'MPI=yes', 'OPENMP=yes', 'OPTIMIZE=yes'])
assert returncode == 0
shutil.move('openmc', 'openmc-mpi-omp-optimize')
def test_mpi_hdf5():
@ -102,6 +187,43 @@ def test_mpi_hdf5():
shutil.move('openmc', 'openmc-phdf5')
def test_mpi_hdf5_debug():
returncode = run(['make', 'distclean'])
returncode = run(['make', compiler, 'MPI=yes', 'HDF5=yes', 'DEBUG=yes'])
assert returncode == 0
shutil.move('openmc', 'openmc-phdf5-debug')
def test_mpi_hdf5_optimize():
returncode = run(['make', 'distclean'])
returncode = run(['make', compiler, 'MPI=yes', 'HDF5=yes', 'OPTIMIZE=yes'])
assert returncode == 0
shutil.move('openmc', 'openmc-phdf5-optimize')
def test_mpi_omp_hdf5():
returncode = run(['make', 'distclean'])
returncode = run(['make', compiler, 'MPI=yes', 'OPENMP=yes', 'HDF5=yes'])
assert returncode == 0
shutil.move('openmc', 'openmc-phdf5-omp')
def test_mpi_omp_hdf5_debug():
returncode = run(['make', 'distclean'])
returncode = run(['make', compiler, 'MPI=yes', 'OPENMP=yes', 'HDF5=yes',
'DEBUG=yes'])
assert returncode == 0
shutil.move('openmc', 'openmc-phdf5-omp-debug')
def test_mpi_omp_hdf5_optimize():
returncode = run(['make', 'distclean'])
returncode = run(['make', compiler, 'MPI=yes', 'OPENMP=yes', 'HDF5=yes',
'OPTIMIZE=yes'])
assert returncode == 0
shutil.move('openmc', 'openmc-phdf5-omp-optimize')
def test_mpi_hdf5_petsc():
returncode = run(['make', 'distclean'])
returncode = run(['make', compiler, 'MPI=yes', 'HDF5=yes', 'PETSC=yes'])
@ -109,6 +231,14 @@ def test_mpi_hdf5_petsc():
shutil.move('openmc', 'openmc-phdf5-petsc')
def test_mpi_hdf5_petsc_debug():
returncode = run(['make', 'distclean'])
returncode = run(['make', compiler, 'MPI=yes', 'HDF5=yes', 'PETSC=yes',
'DEBUG=yes'])
assert returncode == 0
shutil.move('openmc', 'openmc-phdf5-petsc-debug')
def test_mpi_hdf5_petsc_optimize():
returncode = run(['make', 'distclean'])
returncode = run(['make', compiler, 'MPI=yes', 'HDF5=yes', 'PETSC=yes',
@ -117,6 +247,22 @@ def test_mpi_hdf5_petsc_optimize():
shutil.move('openmc', 'openmc-phdf5-petsc-optimize')
def test_mpi_omp_hdf5_petsc():
returncode = run(['make', 'distclean'])
returncode = run(['make', compiler, 'MPI=yes', 'OMP=yes', 'HDF5=yes',
'PETSC=yes'])
assert returncode == 0
shutil.move('openmc', 'openmc-omp-phdf5-petsc')
def test_mpi_omp_hdf5_petsc_debug():
returncode = run(['make', 'distclean'])
returncode = run(['make', compiler, 'MPI=yes', 'OMP=yes', 'HDF5=yes',
'PETSC=yes', 'DEBUG=yes'])
assert returncode == 0
shutil.move('openmc', 'openmc-omp-phdf5-petsc-debug')
def test_mpi_omp_hdf5_petsc_optimize():
returncode = run(['make', 'distclean'])
returncode = run(['make', compiler, 'MPI=yes', 'OMP=yes', 'HDF5=yes',
@ -127,8 +273,8 @@ def test_mpi_omp_hdf5_petsc_optimize():
def run(commands):
proc = Popen(commands, stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
return returncode

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_created_statepoint():

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_created_statepoint():

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_created_statepoint():

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_created_statepoint():

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_created_statepoint():

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_created_statepoint():

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_created_statepoint():

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_created_statepoint():

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_created_statepoint():

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_created_statepoint():

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_created_statepoint():

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_created_statepoint():

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_created_statepoint():

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_created_statepoint():

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_created_statepoint():

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_created_statepoint():

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_created_statepoint():

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_created_statepoint():

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_created_statepoint():

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_created_statepoint():

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_created_statepoint():

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_created_statepoint():

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_summary_exists():

View file

@ -17,8 +17,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path, '-p'], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_plot_exists():

View file

@ -17,8 +17,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path, '-p'], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_plots_exists():

View file

@ -17,8 +17,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path, '-p'], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_plot_exists():

View file

@ -17,8 +17,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path, '-p'], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_plots_exists():

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_created_statepoint():

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_created_statepoint():

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_created_statepoint():

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_created_statepoint():

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_created_statepoint():

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_created_statepoint():

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_created_statepoint():

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_created_statepoint():

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_created_statepoint():

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_created_statepoint():

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_created_statepoint():

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_created_statepoint():

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_created_statepoint():

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_created_statepoint():

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_created_statepoint():

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_created_statepoint():

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_created_statepoint():

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_created_statepoint():

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_created_statepoint():

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_created_statepoint():

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_created_statepoint():

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_created_statepoint():

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_created_statepoint():

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_created_statepoint():

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_created_statepoint():

View file

@ -19,8 +19,8 @@ def test_run():
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0
def test_created_statepoint():

Some files were not shown because too many files have changed in this diff Show more