diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index 23832c6f64..52be20bc4c 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -13,16 +13,25 @@ ace.o: string.o ace_header.o: constants.o ace_header.o: endf_header.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_execute.o: constants.o cmfd_execute.o: cmfd_data.o cmfd_execute.o: cmfd_output.o cmfd_execute.o: cmfd_message_passing.o cmfd_execute.o: cmfd_power_solver.o cmfd_execute.o: cmfd_snes_solver.o +cmfd_execute.o: error.o +cmfd_execute.o: global.o +cmfd_execute.o: mesh_header.o +cmfd_execute.o: mesh.o +cmfd_execute.o: search.o cmfd_execute.o: tally.o cmfd_execute.o: timing.o @@ -37,6 +46,7 @@ cmfd_input.o: string.o cmfd_input.o: tally_header.o cmfd_input.o: templates/cmfd_t.o +cmfd_jacobian_operator.o: constants.o cmfd_jacobian_operator.o: global.o cmfd_jacobian_operator.o: cmfd_loss_operator.o cmfd_jacobian_operator.o: cmfd_prod_operator.o @@ -54,15 +64,26 @@ cmfd_output.o: mesh_header.o cmfd_output.o: string.o cmfd_output.o: timing.o +cmfd_power_solver.o: constants.o cmfd_power_solver.o: cmfd_loss_operator.o cmfd_power_solver.o: cmfd_prod_operator.o +cmfd_power_solver.o: global.o +cmfd_power_solver.o: string.o cmfd_prod_operator.o: constants.o cmfd_prod_operator.o: global.o +cmfd_slepc_solver.o: constants.o +cmfd_slepc_solver.o: cmfd_loss_operator.o +cmfd_slepc_solver.o: cmfd_prod_operator.o +cmfd_slepc_solver.o: global.o +cmfd_slepc_solver.o: timing.o + cmfd_snes_solver.o: cmfd_jacobian_operator.o cmfd_snes_solver.o: cmfd_loss_operator.o cmfd_snes_solver.o: cmfd_prod_operator.o +cmfd_snes_solver.o: constants.o +cmfd_snes_solver.o: global.o cross_section.o: ace_header.o cross_section.o: constants.o diff --git a/src/cmfd_data.F90 b/src/cmfd_data.F90 index 061597a8af..c2141648fc 100644 --- a/src/cmfd_data.F90 +++ b/src/cmfd_data.F90 @@ -53,16 +53,16 @@ contains subroutine compute_xs() - use constants, only: FILTER_MESH, FILTER_ENERGYIN, FILTER_ENERGYOUT& - , SURF_FILTER_ENERGYIN, SURF_FILTER_SURFACE, & - IN_RIGHT, OUT_RIGHT, IN_FRONT, OUT_FRONT, & - IN_TOP, OUT_TOP, N_FILTER_TYPES, CMFD_NOACCEL - use global, only: cmfd, message, n_user_tallies, n_tallies, & - tallies, meshes - use error, only: fatal_error - use mesh, only: mesh_indices_to_bin - use mesh_header, only: StructuredMesh - use tally_header, only: TallyObject, TallyScore + use constants, only: FILTER_MESH, FILTER_ENERGYIN, FILTER_ENERGYOUT, & + SURF_FILTER_ENERGYIN, SURF_FILTER_SURFACE, IN_RIGHT,& + OUT_RIGHT, IN_FRONT, OUT_FRONT, IN_TOP, OUT_TOP, & + N_FILTER_TYPES, CMFD_NOACCEL, ZERO, ONE + use error, only: fatal_error + use global, only: cmfd, message, n_user_tallies, n_tallies, tallies, & + meshes + use mesh, only: mesh_indices_to_bin + use mesh_header, only: StructuredMesh + use tally_header, only: TallyObject, TallyScore integer :: nx ! number of mesh cells in x direction integer :: ny ! number of mesh cells in y direction @@ -90,8 +90,8 @@ contains ng = cmfd % indices(4) ! set flux object and source distribution to all zeros - cmfd % flux = 0.0_8 - cmfd % openmc_src = 0.0_8 + cmfd % flux = ZERO + cmfd % openmc_src = ZERO ! associate tallies and mesh t => tallies(n_user_tallies + 1) @@ -165,8 +165,8 @@ contains cmfd % diffusion(h,i,j,k) = t % scores(4,score_index) % sum / flux ! calculate diffusion coefficient -! cmfd % diffcof(h,i,j,k) = 1.0_8/(3.0_8*cmfd%totalxs(h,i,j,k)) - cmfd % diffcof(h,i,j,k) = 1.0_8/(3.0_8*(cmfd % totalxs(h,i,j,k) - & +! cmfd % diffcof(h,i,j,k) = ONE/(3.0_8*cmfd%totalxs(h,i,j,k)) + cmfd % diffcof(h,i,j,k) = ONE/(3.0_8*(cmfd % totalxs(h,i,j,k) - & cmfd % p1scattxs(h,i,j,k))) ! cmfd % diffcof(h,i,j,k) = cmfd % diffusion(h,i,j,k) @@ -355,7 +355,7 @@ contains subroutine neutron_balance(uid) - use constants, only: ONE, CMFD_NOACCEL + use constants, only: ONE, ZERO, CMFD_NOACCEL use global, only: cmfd, keff integer :: nx ! number of mesh cells in x direction @@ -376,7 +376,7 @@ contains real(8) :: res ! residual of neutron balance ! check if keff is close to 0 (happens on first active batch) - if (keff - 0.0_8 < 1e-8_8) return + if (keff < 1e-8_8) return ! extract spatial and energy indices from object nx = cmfd % indices(1) @@ -405,7 +405,7 @@ contains end if ! get leakage - leakage = 0.0_8 + leakage = ZERO LEAK: do l = 1, 3 leakage = leakage + ((cmfd % current(4*l,g,i,j,k) - & @@ -419,8 +419,8 @@ contains interactions = cmfd % totalxs(g,i,j,k) * cmfd % flux(g,i,j,k) ! get scattering and fission - scattering = 0.0_8 - fission = 0.0_8 + scattering = ZERO + fission = ZERO GROUPH: do h = 1, ng scattering = scattering + cmfd % scattxs(h,g,i,j,k) * & @@ -611,7 +611,7 @@ contains subroutine compute_dhat() - use constants, only: CMFD_NOACCEL + use constants, only: CMFD_NOACCEL, ZERO use global, only: cmfd, cmfd_coremap integer :: nx ! maximum number of cells in x direction @@ -732,14 +732,14 @@ contains ! if ((abs(current(2*l)-current(2*l-1)) < 1e-8_8).and.xyz_idx/=3) then ! print *,'Zero net current interface',g,i,j,k ! print *,current(2*l),current(2*l-1),net_current -! dhat = 0.0_8 +! dhat = ZERO ! end if ! record dhat in cmfd object cmfd%dhat(l,g,i,j,k) = dhat ! check for dhat reset - if (dhat_reset) cmfd%dhat(l,g,i,j,k) = 0.0_8 + if (dhat_reset) cmfd%dhat(l,g,i,j,k) = ZERO end do LEAK @@ -780,8 +780,8 @@ contains shift_idx = -2*mod(l,2) + 1 ! shift neig by -1 or +1 ! calculate albedo - if ((shift_idx == 1 .and. (current(2*l )-0.0_8) < 1.0e-10_8) .or. & - (shift_idx == -1 .and. (current(2*l-1)-0.0_8) < 1.0e-10_8)) then + if ((shift_idx == 1 .and. current(2*l ) < 1.0e-10_8) .or. & + (shift_idx == -1 .and. current(2*l-1) < 1.0e-10_8)) then albedo = ALBEDO_REJECT cmfd_hold_weights = .true. else diff --git a/src/cmfd_execute.F90 b/src/cmfd_execute.F90 index d1b79f26b1..d101757614 100644 --- a/src/cmfd_execute.F90 +++ b/src/cmfd_execute.F90 @@ -201,7 +201,7 @@ contains subroutine calc_fission_source() - use constants, only: CMFD_NOACCEL + use constants, only: CMFD_NOACCEL, ZERO, TWO use global, only: cmfd, cmfd_coremap, master, mpi_err, entropy_on integer :: nx ! maximum number of cells in x direction @@ -229,7 +229,7 @@ contains if (.not. allocated(cmfd%cmfd_src)) allocate(cmfd%cmfd_src(ng,nx,ny,nz)) ! reset cmfd source to 0 - cmfd%cmfd_src = 0.0_8 + cmfd%cmfd_src = ZERO ! only perform for master if (master) then @@ -281,15 +281,15 @@ contains if (.not.allocated(source)) allocate(source(ng,nx,ny,nz)) ! initialize the source - source = 0.0_8 + source = ZERO ! compute log - where (cmfd%cmfd_src > 0.0_8) - source = cmfd%cmfd_src*log(cmfd%cmfd_src)/log(2.0_8) + where (cmfd%cmfd_src > ZERO) + source = cmfd%cmfd_src*log(cmfd%cmfd_src)/log(TWO) end where ! sum that source - cmfd%entropy = -1.0_8*sum(source) + cmfd%entropy = -sum(source) ! deallocate tmp array if (allocated(source)) deallocate(source) @@ -312,6 +312,7 @@ contains subroutine cmfd_reweight(new_weights) + use constants, only: ZERO, ONE use error, only: warning, fatal_error use global, only: n_particles, meshes, source_bank, work, & n_user_meshes, message, cmfd, master, mpi_err @@ -350,7 +351,7 @@ contains end if if (.not.allocated(cmfd%weightfactors)) then allocate(cmfd%weightfactors(ng,nx,ny,nz)) - cmfd % weightfactors = 1.0_8 + cmfd % weightfactors = ONE end if ! allocate energy grid and reverse cmfd energy grid @@ -361,7 +362,7 @@ contains if (new_weights) then ! zero out weights - cmfd%weightfactors = 0.0_8 + cmfd%weightfactors = ZERO ! count bank sites in mesh call count_bank_sites(m, source_bank, cmfd%sourcecounts, egrid, & @@ -369,7 +370,7 @@ contains ! have master compute weight factors if (master) then - where(cmfd%cmfd_src > 0.0_8 .and. cmfd%sourcecounts > 0.0_8) + 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 @@ -427,7 +428,7 @@ contains function get_matrix_idx(g, i, j, k, ng, nx, ny) result (matidx) - use global, only: cmfd,cmfd_coremap + use global, only: cmfd, cmfd_coremap integer :: matidx ! the index location in matrix integer :: i ! current x index diff --git a/src/cmfd_header.F90 b/src/cmfd_header.F90 index ed22530770..63c5ea44df 100644 --- a/src/cmfd_header.F90 +++ b/src/cmfd_header.F90 @@ -1,6 +1,6 @@ module cmfd_header - use constants, only: CMFD_NOACCEL + use constants, only: CMFD_NOACCEL, ZERO, ONE implicit none private @@ -57,17 +57,17 @@ module cmfd_header ! eigenvector/eigenvalue from cmfd run real(8), allocatable :: phi(:) - real(8) :: keff = 0.0_8 + real(8) :: keff = ZERO ! eigenvector/eigenvalue from adjoint run real(8), allocatable :: adj_phi(:) - real(8) :: adj_keff = 0.0_8 + real(8) :: adj_keff = ZERO ! residual for neutron balance real(8), allocatable :: resnb(:,:,:,:) ! openmc source normalization factor - real(8) :: norm = 1.0_8 + real(8) :: norm = ONE ! Shannon entropy from cmfd fission source real(8) :: entropy @@ -123,21 +123,21 @@ contains if (.not. allocated(this % weightfactors)) allocate(this % weightfactors(ng,nx,ny,nz)) ! set everthing to 0 except weight multiply factors if feedback isnt on - this % flux = 0.0_8 - this % totalxs = 0.0_8 - this % p1scattxs = 0.0_8 - this % scattxs = 0.0_8 - this % nfissxs = 0.0_8 - this % diffcof = 0.0_8 - this % diffusion = 0.0_8 - this % dtilde = 0.0_8 - this % dhat = 0.0_8 - this % hxyz = 0.0_8 - this % current = 0.0_8 - this % cmfd_src = 0.0_8 - this % openmc_src = 0.0_8 - this % sourcecounts = 0.0_8 - this % weightfactors = 1.0_8 + this % flux = ZERO + this % totalxs = ZERO + this % p1scattxs = ZERO + this % scattxs = ZERO + this % nfissxs = ZERO + this % diffcof = ZERO + this % diffusion = ZERO + this % dtilde = ZERO + this % dhat = ZERO + this % hxyz = ZERO + this % current = ZERO + this % cmfd_src = ZERO + this % openmc_src = ZERO + this % sourcecounts = ZERO + this % weightfactors = ONE end subroutine allocate_cmfd diff --git a/src/cmfd_jacobian_operator.F90 b/src/cmfd_jacobian_operator.F90 index 286b13727a..81c709a658 100644 --- a/src/cmfd_jacobian_operator.F90 +++ b/src/cmfd_jacobian_operator.F90 @@ -169,6 +169,7 @@ contains 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 @@ -213,7 +214,7 @@ contains if (.not. allocated(vals)) allocate(vals(& maxval(ctx%loss%d_nnz + ctx%loss%o_nnz))) cols = 0 - vals = 0.0_8 + vals = ZERO ! get pointers to residual vector call VecGetArrayF90(x, xptr, ierr) @@ -227,7 +228,7 @@ contains 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, -1.0_8*lambda, ctx%prod%F, & + call MatAXPY(ctx%loss%M, -lambda, ctx%prod%F, & DIFFERENT_NONZERO_PATTERN, ierr) ! create tmp petsc vector for source @@ -236,7 +237,7 @@ contains ! perform math (-F*phi --> source) call MatMult(ctx%prod%F, phi, source, ierr) - call VecScale(source, -1.0_8, ierr) + call VecScale(source, -ONE, ierr) ! get pointer to source call VecGetArrayF90(source, sptr, ierr) @@ -283,10 +284,10 @@ contains ! set values in last row of matrix if (rank == n_procs_cmfd - 1) then - phi_tmp = -1.0_8*phi_tmp ! negate the transpose + 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, 1.0_8, INSERT_VALUES, ierr) + call MatSetValue(jac_prec, n, n, ONE, INSERT_VALUES, ierr) end if ! assemble matrix diff --git a/src/cmfd_loss_operator.F90 b/src/cmfd_loss_operator.F90 index 74f392c8fe..4b5fb1196d 100644 --- a/src/cmfd_loss_operator.F90 +++ b/src/cmfd_loss_operator.F90 @@ -248,7 +248,7 @@ contains subroutine build_loss_matrix(this, adjoint) - use constants, only: CMFD_NOACCEL + use constants, only: CMFD_NOACCEL, ZERO use global, only: cmfd, cmfd_coremap, cmfd_write_matrices type(loss_operator) :: this @@ -313,7 +313,7 @@ contains if (allocated(cmfd%dhat)) then dhat = cmfd%dhat(:,g,i,j,k) else - dhat = 0.0_8 + dhat = ZERO end if ! create boundary vector diff --git a/src/cmfd_power_solver.F90 b/src/cmfd_power_solver.F90 index 7dd2903c94..dbb551e3cb 100644 --- a/src/cmfd_power_solver.F90 +++ b/src/cmfd_power_solver.F90 @@ -103,6 +103,8 @@ contains subroutine init_data() + use constants, only: ONE + integer :: n ! problem size real(8) :: guess ! initial guess @@ -122,7 +124,7 @@ contains call VecCreateMPI(PETSC_COMM_WORLD, n, PETSC_DECIDE, S_o, ierr) ! set initial guess - guess = 1.0_8 + guess = ONE call VecSet(phi_n, guess, ierr) call VecSet(phi_o, guess, ierr) k_n = guess @@ -248,9 +250,10 @@ contains subroutine execute_power_iter() + use constants, only: ONE + real(8) :: num ! numerator for eigenvalue update real(8) :: den ! denominator for eigenvalue update - real(8) :: one=1.0_8 ! one integer :: i ! iteration counter ! reset convergence flag @@ -263,7 +266,7 @@ contains call MatMult(prod%F, phi_o, S_o, ierr) ! normalize source vector - call VecScale(S_o, one/k_o, ierr) + call VecScale(S_o, ONE/k_o, ierr) ! compute new flux vector call KSPSolve(krylov, S_o, phi_n, ierr) @@ -339,6 +342,7 @@ contains subroutine extract_results() + use constants, only: ZERO use global, only: cmfd, n_procs_cmfd, cmfd_write_matrices integer :: n ! problem size @@ -379,7 +383,7 @@ contains end if ! reduce result to all - mybuf = 0.0_8 + mybuf = ZERO if (adjoint_calc) then call MPI_ALLREDUCE(cmfd%adj_phi, mybuf, n, MPI_REAL8, MPI_SUM, & PETSC_COMM_WORLD, ierr) diff --git a/src/cmfd_slepc_solver.F90 b/src/cmfd_slepc_solver.F90 index dfb25434c9..0dbe17434b 100644 --- a/src/cmfd_slepc_solver.F90 +++ b/src/cmfd_slepc_solver.F90 @@ -31,8 +31,8 @@ contains subroutine cmfd_slepc_execute() - use timing use global, only:time_cmfd,master + use timing call timer_start(time_cmfd) ! initialize data @@ -141,7 +141,8 @@ call timer_start(time_cmfd) subroutine extract_results() - use global, only: cmfd, path_input, master + use constants, only: ZERO + use global, only: cmfd, path_input, master integer :: n ! problem size integer :: i_eig = 0 ! eigenvalue to extract @@ -159,7 +160,7 @@ call timer_start(time_cmfd) if (.not. allocated(mybuf)) allocate(mybuf(n)) ! zero out cmfd object - cmfd%phi = 0.0_8 + cmfd%phi = ZERO ! extract run information call EPSGetEigenpair(eps, i_eig, keff, PETSC_NULL, phi, PETSC_NULL_OBJECT, ierr) @@ -176,7 +177,7 @@ call timer_start(time_cmfd) cmfd%keff = keff ! reduce result to master - mybuf = 0.0_8 + mybuf = ZERO call MPI_ALLREDUCE(cmfd%phi, mybuf, n, MPI_REAL8, MPI_SUM, & PETSC_COMM_WORLD, ierr) diff --git a/src/cmfd_snes_solver.F90 b/src/cmfd_snes_solver.F90 index 298b7fc9ed..728ce35485 100644 --- a/src/cmfd_snes_solver.F90 +++ b/src/cmfd_snes_solver.F90 @@ -72,7 +72,8 @@ contains subroutine init_data() - use global, only: cmfd, n_procs_cmfd, rank + use constants, only: ONE + use global, only: cmfd, n_procs_cmfd, rank integer :: k ! implied do counter integer :: n ! problem size @@ -114,9 +115,9 @@ contains if (rank == n_procs_cmfd - 1) then call VecGetArrayF90(xvec, xptr, ierr) if (adjoint_calc) then - xptr(size(xptr)) = 1.0_8/cmfd%adj_keff + xptr(size(xptr)) = ONE/cmfd%adj_keff else - xptr(size(xptr)) = 1.0_8/cmfd%keff + xptr(size(xptr)) = ONE/cmfd%keff end if call VecRestoreArrayF90(xvec, xptr, ierr) end if @@ -296,7 +297,7 @@ contains ! 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, -1.0_8*lambda, phiM, ierr) + call VecAYPX(rphi, -lambda, phiM, ierr) ! set eigenvalue part of residual vector call VecDot(phi, phi, reslamb, ierr) @@ -382,7 +383,8 @@ contains subroutine extract_results() - use global, only: cmfd, n_procs_cmfd, rank + use constants, only: ZERO, ONE + use global, only: cmfd, n_procs_cmfd, rank integer :: n ! problem size integer :: row_start ! local row start @@ -417,7 +419,7 @@ contains end if ! reduce result to all - mybuf = 0.0_8 + mybuf = ZERO if (adjoint_calc) then call MPI_ALLREDUCE(cmfd%adj_phi, mybuf, n, MPI_REAL8, MPI_SUM, & PETSC_COMM_WORLD, ierr) @@ -431,7 +433,7 @@ contains if(allocated(mybuf)) deallocate(mybuf) ! save eigenvalue - if(rank == n_procs_cmfd - 1) keff = 1.0_8 / xptr(size(xptr)) + 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, &