diff --git a/ChangeLog b/ChangeLog index 90b9029d36..bc4be65a06 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,15 @@ +2011-03-22 Paul Romano + + * benchmarks: Added PU-MET-FAST-002 and PU-MET-FAST-005 benchmark + problems. + * mpi_routines.f90: Added check for no fission bank sites. + * output.f90: Changed error subroutine to work when run in + parallel. + * physics.f90: Fixed sampling of seconary angle for elastic + scattering. Changed sampling of reaction. Changed n_gamma to + n_absorption and now includes (n,p), etc. Reverted sample_angle + routine to arguments u,v,w,mu. + 2011-03-22 Paul Romano * cross_section.f90: Calculation of material total cross section diff --git a/benchmarks/pu-met-fast-001 b/benchmarks/pu-met-fast-001 index e14789b57c..18de4871f1 100644 --- a/benchmarks/pu-met-fast-001 +++ b/benchmarks/pu-met-fast-001 @@ -10,7 +10,7 @@ cell 1 0 1 -1 surface 1 sph 0. 0. 0. 6.3849 -source_box -4 -4 -4 4 4 4 +source_box -1 -1 -1 1 1 1 material 1 0.04029014 & 94239.03c 0.037047 & @@ -20,5 +20,5 @@ material 1 0.04029014 & xs_library endfb7 xs_data /opt/serpent/xsdata/endfb7 -criticality 150 25 1000 +criticality 3000 20 10000 verbosity 7 diff --git a/benchmarks/pu-met-fast-002 b/benchmarks/pu-met-fast-002 new file mode 100644 index 0000000000..5b6812cdb6 --- /dev/null +++ b/benchmarks/pu-met-fast-002 @@ -0,0 +1,25 @@ +# ==================================================================== +# Description: Bare Pu-240 Jezebel Benchmark +# Case: PU-MET-FAST-002 +# Written By: Paul Romano +# Date: 3/22/2011 +# ==================================================================== + +# cell no. universe material surfaces +cell 1 0 1 -1 + +surface 1 sph 0. 0. 0. 6.6595 + +source_box -1 -1 -1 1 1 1 + +material 1 0.04055292 & + 94239.03c 0.029934 & + 94240.03c 0.0078754 & + 94241.03c 0.0012146 & + 94242.03c 0.00015672 & + 31000.03c 0.0013722 + +xs_library endfb7 +xs_data /opt/serpent/xsdata/endfb7 +criticality 3000 20 10000 +verbosity 7 diff --git a/benchmarks/pu-met-fast-005 b/benchmarks/pu-met-fast-005 new file mode 100644 index 0000000000..ef094ad311 --- /dev/null +++ b/benchmarks/pu-met-fast-005 @@ -0,0 +1,32 @@ +# ==================================================================== +# Description: Unmoderated Plutonium Metal Button Array Benchmark +# Case: PU-MET-FAST-005 +# Written By: Paul Romano +# Date: 3/22/2011 +# ==================================================================== + +# cell no. universe material surfaces +cell 1 0 1 -1 +cell 2 0 2 1 -2 + +surface 1 sph 0. 0. 0. 5.0419 +surface 2 sph 0. 0. 0. 9.7409 + +source_box -1 -1 -1 1 1 1 + +material 1 0.04070346 & + 94239.03c 0.037291 & + 94240.03c 0.0019277 & + 94241.03c 0.00012196 & + 31000.03c 0.0013628 + +material 2 0.06605308 & + 74000.03c 0.051468 & + 28000.03c 0.0097124 & + 29000.03c 0.0040774 & + 40000.03c 0.00079528 + +xs_library endfb7 +xs_data /opt/serpent/xsdata/endfb7 +criticality 3000 20 10000 +verbosity 7 diff --git a/src/mpi_routines.f90 b/src/mpi_routines.f90 index bcf3c77695..001a711126 100644 --- a/src/mpi_routines.f90 +++ b/src/mpi_routines.f90 @@ -119,6 +119,12 @@ contains call MPI_BCAST(total, 1, MPI_INTEGER, n_procs - 1, & & MPI_COMM_WORLD, ierr) + ! Check if there are no fission sites + if (total == 0) then + msg = "No fission sites banked!" + call error(msg) + end if + allocate(temp_sites(2*work)) count = 0 ! Index for source_bank index = 0 ! Index for uid -- must account for all nodes diff --git a/src/output.f90 b/src/output.f90 index 4628c18976..0aa2a4ee26 100644 --- a/src/output.f90 +++ b/src/output.f90 @@ -159,20 +159,21 @@ module output integer :: i ! Only allow master to print to screen - if (.not. master) return + if (master) then + write(eu, fmt='(1X,A7)', advance='no') 'ERROR: ' - write(eu, fmt='(1X,A7)', advance='no') 'ERROR: ' - - n_lines = (len_trim(msg)-1)/72 + 1 - do i = 1, n_lines - if ( i == 1 ) then - write(eu, fmt='(A72)') msg(72*(i-1)+1:72*i) - else - write(eu, fmt='(7X,A72)') msg(72*(i-1)+1:72*i) - end if - end do - write(eu,*) + n_lines = (len_trim(msg)-1)/72 + 1 + do i = 1, n_lines + if ( i == 1 ) then + write(eu, fmt='(A72)') msg(72*(i-1)+1:72*i) + else + write(eu, fmt='(7X,A72)') msg(72*(i-1)+1:72*i) + end if + end do + write(eu,*) + end if + ! All processors abort call free_memory() end subroutine error diff --git a/src/physics.f90 b/src/physics.f90 index b7efa242e6..eaf413e72a 100644 --- a/src/physics.f90 +++ b/src/physics.f90 @@ -22,7 +22,8 @@ contains type(Particle), pointer :: p character(250) :: msg - integer :: i, surf + integer :: i + integer :: surf logical :: alive logical :: found_cell @@ -184,11 +185,11 @@ contains end do ! Get table, total xs, interpolation factor + density_i = cMaterial%atom_percent(i)*density table => xs_continuous(cMaterial%table(i)) - Sigma = Sigma_t(i) + Sigma = Sigma_t(i) / density_i IE = table%grid_index(p % IE) f = (p%E - table%energy(IE))/(table%energy(IE+1) - table%energy(IE)) - density = cMaterial%atom_percent(i)*density ! free memory deallocate(Sigma_t) @@ -198,9 +199,9 @@ contains prob = ZERO do i = 1, table%n_reaction rxn => table%reactions(i) - if (rxn%MT >= 200) cycle + if (rxn%MT >= 200 .or. rxn%MT == 4) cycle if (IE < rxn%IE) cycle - prob = prob + density * ((ONE-f)*rxn%sigma(IE-rxn%IE+1) & + prob = prob + ((ONE-f)*rxn%sigma(IE-rxn%IE+1) & & + f*(rxn%sigma(IE-rxn%IE+2))) if (r1 < prob) exit end do @@ -213,17 +214,17 @@ contains ! call appropriate subroutine select case (rxn%MT) case (2) - call elastic_scatter(p, table%awr) + call elastic_scatter(p, table, rxn) case (11, 16, 17, 24, 25, 30, 37, 41, 42) call n_xn(p, table, rxn) case (22, 23, 28, 29, 32:36, 43, 44, 51:91) call inelastic_scatter(p, table, rxn) case (18:21, 38) call n_fission(p, table, rxn) - case (102) - call n_gamma(p) + case (102:117) + call n_absorption(p) case default - call elastic_scatter(p, table%awr) + print *, 'warning: cant simulate reaction: ', rxn % MT end select ! find energy index, interpolation factor @@ -235,11 +236,13 @@ contains ! ELASTIC_SCATTER !===================================================================== - subroutine elastic_scatter(p, awr) + subroutine elastic_scatter(p, table, rxn) - type(Particle), pointer :: p - real(8), intent(in) :: awr + type(Particle), pointer :: p + type(AceContinuous), pointer :: table + type(AceReaction), pointer :: rxn + real(8) :: awr real(8) :: phi ! azimuthal angle real(8) :: mu ! cosine of polar angle real(8) :: vx, vy, vz @@ -250,12 +253,12 @@ contains integer :: IE vel = sqrt(p % E) + awr = table % awr vx = vel*p % uvw(1) vy = vel*p % uvw(2) vz = vel*p % uvw(3) - vcx = vx/(awr + ONE) vcy = vy/(awr + ONE) vcz = vz/(awr + ONE) @@ -267,13 +270,16 @@ contains vel = sqrt(vx*vx + vy*vy + vz*vz) - ! Select isotropic direcion -- this is only valid for s-wave - ! scattering - phi = TWO*PI*rang() - mu = TWO*rang() - ONE - u = mu - v = sqrt(ONE - mu**2) * cos(phi) - w = sqrt(ONE - mu**2) * sin(phi) + ! Sample scattering angle + mu = sample_angle(rxn, p % E) + + ! Determine direction cosines + u = vx/vel + v = vy/vel + w = vz/vel + + ! Change direction cosines according to mu + call rotate_angle(u, v, w, mu) vx = u*vel vy = v*vel @@ -523,6 +529,7 @@ contains real(8) :: mu ! cosine of scattering angle real(8) :: E ! outgoing energy in laboratory real(8) :: E_cm ! outgoing energy in center-of-mass + real(8) :: u,v,w ! direction cosines character(250) :: msg ! error message ! copy energy of neutron @@ -557,8 +564,14 @@ contains mu = mu * sqrt(E_cm/E) + ONE/(A+ONE) * sqrt(E_in/E) end if + ! copy directional cosines + u = p % uvw(1) + v = p % uvw(2) + w = p % uvw(3) + ! change direction of particle - call rotate_angle(p, mu) + call rotate_angle(u, v, w, mu) + p % uvw = (/ u, v, w /) ! change energy of particle p % E = E @@ -583,6 +596,7 @@ contains real(8) :: mu ! cosine of scattering angle real(8) :: E ! outgoing energy in laboratory real(8) :: E_cm ! outgoing energy in center-of-mass + real(8) :: u,v,w ! direction cosines character(250) :: msg ! error message ! copy energy of neutron @@ -617,8 +631,14 @@ contains mu = mu * sqrt(E_cm/E) + ONE/(A+ONE) * sqrt(E_in/E) end if + ! copy directional cosines + u = p % uvw(1) + v = p % uvw(2) + w = p % uvw(3) + ! change direction of particle - call rotate_angle(p, mu) + call rotate_angle(u, v, w, mu) + p % uvw = (/ u, v, w /) ! change energy of particle p % E = E @@ -632,10 +652,10 @@ contains end subroutine n_xn !===================================================================== -! N_GAMMA +! N_ABSORPTION !===================================================================== - subroutine n_gamma(p) + subroutine n_absorption(p) type(Particle), pointer :: p @@ -649,7 +669,7 @@ contains call message(msg, 10) end if - end subroutine n_gamma + end subroutine n_absorption !===================================================================== ! SAMPLE_ANGLE samples the cosine of the angle between incident and @@ -775,19 +795,22 @@ contains ! done in MCNP and SERPENT. !===================================================================== - subroutine rotate_angle(p, mu) + subroutine rotate_angle(u, v, w, mu) - type(Particle), pointer :: p - real(8), intent(in) :: mu ! cosine of angle in lab +! type(Particle), pointer :: p + real(8), intent(inout) :: u + real(8), intent(inout) :: v + real(8), intent(inout) :: w + real(8), intent(in) :: mu ! cosine of angle in lab real(8) :: phi, sinphi, cosphi real(8) :: a,b real(8) :: u0, v0, w0 ! Copy original directional cosines - u0 = p % uvw(1) - v0 = p % uvw(2) - w0 = p % uvw(3) + u0 = u + v0 = v + w0 = w ! Sample azimuthal angle in [0,2pi) phi = TWO * PI * rang() @@ -801,14 +824,14 @@ contains ! Need to treat special case where sqrt(1 - w**2) is close to zero ! by expanding about the v component rather than the w component if (b > 1e-10) then - p % uvw(1) = mu*u0 + a*(u0*w0*cosphi - v0*sinphi)/b - p % uvw(2) = mu*v0 + a*(v0*w0*cosphi + u0*sinphi)/b - p % uvw(3) = mu*w0 - a*b*cosphi + u = mu*u0 + a*(u0*w0*cosphi - v0*sinphi)/b + v = mu*v0 + a*(v0*w0*cosphi + u0*sinphi)/b + w = mu*w0 - a*b*cosphi else b = sqrt(ONE - v0*v0) - p % uvw(1) = mu*u0 + a*(u0*v0*cosphi + w0*sinphi)/b - p % uvw(2) = mu*v0 - a*b*cosphi - p % uvw(3) = mu*w0 + a*(v0*w0*cosphi - u0*sinphi)/b + u = mu*u0 + a*(u0*v0*cosphi + w0*sinphi)/b + v = mu*v0 - a*b*cosphi + w = mu*w0 + a*(v0*w0*cosphi - u0*sinphi)/b end if end subroutine rotate_angle