tstXS.f90

./Manager/wrap/tests/tstXS.f90

1 #include "ScaleDBCF/DBCF.h"
2 module myhelpers
3 use origen_librarywrapper_m
4 use dbcf_m
5 use origen_definitions_m
6 
7 #include "ScaleSTL/FortranTestMacros.h"
8 
9 implicit none
10 
11 contains
12 
13 !!search for the daughter [ITOT] corresponding to A(N) with N in [NON]
14 function search_daughter(this,n) result(i)
15  type(origen_librarywrapper),intent(inout) :: this
16  integer, intent(in) :: n
17  integer :: i,last_reaction
18  !keep searching until we go past
19  do i=1,this%get_itot()
20  if( this%non0(i) > n )then
21  exit
22  end if
23  end do
24  if( i>this%get_itot() )then
25  write(error_unit,*)'error finding this daughter'
26  stop 7
27  end if
28 
29 end function search_daughter
30 
31 
32 subroutine super_daughter_check(this,parent_id)
33  implicit none
34  type(origen_librarywrapper),intent(inout) :: this
35  integer :: parent_id
36 
37  integer :: idaughter,adiff,z
38  logical :: was_alpha_decay,was_fp_decay
39  real(C_DOUBLE) :: myfsum,mydsum,mycsum,factor
40  real(C_DOUBLE) :: afy,libc,libd,libf
41  integer :: i,ilite
42  integer,allocatable :: nucl(:), mt(:),loc(:),typ(:)
43  real(C_DOUBLE),allocatable :: a(:),fiss(:),tocap(:)
44  real,allocatable :: dis(:)
45 
46  ilite=this%get_ilite()
47  call this%get_mt(mt)
48  call this%get_nucl(nucl)
49  call this%get_loc(loc)
50  call this%get_typ_nuc(typ)
51  call this%get_a(a)
52  call this%get_fiss(fiss)
53  call this%get_tocap(tocap)
54  call this%get_dis(dis)
55 
56 ! write(ERROR_UNIT,'(3a12,a12,2a6,a5)')'id','type','mt','coeff','alpha','fpdec','mult'
57  myfsum=0
58  mydsum=0
59  mycsum=0
60  do i=1,this%get_non()
61  if( loc(i)==parent_id )then
62  idaughter=search_daughter(this,i)
63  !in order for the sums to be correct, we need to assign a factor
64  !to reactions which produce more than one "nuclide"
65  !for example, alpha-particles appear to be accumulated as id 40020 type=2 (actinide)
66  !thus there needs to be a 0.5 factor to account for alpha + residual heavy nuclide
67 
68  was_alpha_decay=.false.
69  was_fp_decay=.false.
70  factor=1
71 
72  if( mt(i)==18 )then
73  z=floor(nucl(idaughter)/10000.d0)
74  if( z==1 .OR. z==2 )then
75  was_fp_decay=.true.
76  else
77  factor=0.5d0
78  myfsum=myfsum+a(i)*factor
79  end if
80  else if ( mt(i)==-1 )then
81  adiff=nucl(loc(i))-nucl(idaughter)
82  !write(ERROR_UNIT,*)'adiff=',adiff
83  if( nucl(idaughter)==20040 .OR. adiff==20040 )then
84  was_alpha_decay=.true.
85  factor=0.5d0
86  else
87  factor=1
88  end if
89  mydsum=mydsum+factor*a(i)
90  else
91 
92  select case(mt(i))
93  case(16,17,28,37,51,102,103,111); factor = 1; mycsum=mycsum+a(i)*factor
94  case(22,24,25,32,33,34,104,105,106,107,112); factor = 1.d0/2.d0; mycsum=mycsum+a(i)*factor
95  case(29,108); factor = 1.d0/3.d0; mycsum=mycsum+a(i)*factor
96  case default
97  write(error_unit,*)'error in reaction handling with unknown mt=',mt(i)
98  stop 6
99  end select
100 
101  end if
102 ! write(ERROR_UNIT,'(3i12,es12.5,2g6.6,f5.2)')nucl(idaughter),typ(idaughter),mt(i),a(i),was_alpha_decay,was_fp_decay,factor
103  end if
104  end do
105 
106  libd=dis(parent_id)
107  libc=tocap(parent_id)-fiss(parent_id-ilite)
108  libf=fiss(parent_id-ilite)
109  afy=myfsum/(0.5d0*fiss(parent_id-ilite))
110 
111 ! write(ERROR_UNIT,*)'********************'
112 ! write(ERROR_UNIT,*)'LIBRARY CONSISTENCY'
113 ! write(ERROR_UNIT,*)'********************'
114 ! write(ERROR_UNIT,'(a10,3a12)')'source|','disint.','capture','fission'
115 ! write(ERROR_UNIT,'(a10,3es12.5)')'summed|',mydsum,mycsum,myfsum
116 ! write(ERROR_UNIT,'(a10,3es12.5)')'stored|',libd,libc,libf
117 ! write(ERROR_UNIT,*)
118 ! write(ERROR_UNIT,'(a,es12.5)')'average fission yield: ',afy
119 ! write(ERROR_UNIT,*)'********************'
120 
121  expect_near(libd, mydsum, 1.0d-4)
122  expect_near(libc, mycsum, 1.0d-5)
123  expect_near(libf, myfsum, 1.0d-5)
124  expect_near(2.d0, afy, 1.d-5)
125 
126 end subroutine
127 
128 end module myhelpers
129 
130 
131 program tstxs
132 
133 use origen_librarywrapper_m
134 use origen_testpaths_m
135 use dbcf_m
136 use integer2dlist_m
137 use integerlist_m
138 use origen_definitions_m, only: approx_eq
139 
140 #include "ScaleSTL/FortranTestMacros.h"
141 
142 use myhelpers !local module at the bottom
143 !
144 use nemesis_comm, only: build_types, initialize, finalize, node
145 
146 implicit none
147 
148 type(origen_librarywrapper) :: lib
149 
150 character(256) :: rbuffer(256)
151 integer :: rstatus
152 
153 integer :: itot,ilite,iact,ifp,nmd
154 integer,allocatable :: myids(:)
155 integer :: kfound
156 integer,allocatable :: mt_set(:)
157 integer :: nmts
158 
159 integer,allocatable :: nucl(:), typ(:),mt(:),loc(:),kd(:),non0(:),kd_(:),non0_(:)
160 real(C_DOUBLE),allocatable :: tocap(:),fiss(:),a(:)
161 real,allocatable :: dis(:)
162 
163 integer :: i,j,decay_start,decay_end,reaction_start,reaction_end,u235
164 real(C_DOUBLE) :: myfis,xs1,xs2,xs3
165 integer,allocatable :: parents(:),daughters(:),mts_p(:),mts_d(:),tinds_p(:),tinds_d(:)
166 
167 integer,allocatable :: mymts(:)
168 real(C_DOUBLE),allocatable :: myxs(:)
169 type(integer2dlist) :: tinds_2d
170 
171 ! CHARACTER(len=255) :: cwd
172 ! CALL getcwd(cwd)
173 ! write(ERROR_UNIT,*)TRIM(cwd)
174 
175 call initialize
176 call build_types
177 if ( node() /= 0 ) call finalize
178 
179 call lib%initialize(pwrlib_filepath)
180 
181 !check the MTS present on this library
182 call lib%find_all_mts(mt_set)
183 !write(ERROR_UNIT,*)'MT numbers present on this library:'
184 !write(ERROR_UNIT,'(a5,x,a5)')'i','mt'
185 !do i=1,size(mt_set)
186 ! write(ERROR_UNIT,'(i5,x,i5)')i,mt_set(i)
187 !end do
188 insist(mt_set(1)==16,'present: Direct (n,2n) cross section MT=16')
189 insist(mt_set(2)==17,'present: (n,3n) cross section MT=17')
190 insist(mt_set(3)==18,'present: Total fission cross section (MT=18, sum of MT = 19, 20, 21, 38)')
191 insist(mt_set(4)==22,"present: (n,n'ALPHA) cross section MT=22")
192 insist(mt_set(5)==23,"present: (n,2n'ALPHA) cross section MT=23")
193 insist(mt_set(6)==24,"present: (n,2n'ALPHA) cross section MT=24")
194 insist(mt_set(7)==25,"present: (n,3n'ALPHA) cross section MT=25")
195 insist(mt_set(8)==28,"present: (n,n'p) cross section MT=28")
196 insist(mt_set(9)==29,"present: (n,nN2'ALPHA) cross section MT=29")
197 insist(mt_set(10)==32,"present: (n,n'd) cross section MT=32")
198 insist(mt_set(11)==33,"present: (n,n't) cross section MT=33")
199 insist(mt_set(12)==34,"present: (n,n'He-3) MT=34")
200 insist(mt_set(13)==37,"present: (n,4n) cross section MT=37")
201 insist(mt_set(14)==51,"present: (n,n') to the 1st excited state MT=51")
202 insist(mt_set(15)==102,"present: (n,GAMMA) radiative capture cross section MT=102")
203 insist(mt_set(16)==103,"present: (n,p) cross section MT=103")
204 insist(mt_set(17)==104,"present: (n,d) cross section MT=104")
205 insist(mt_set(18)==105,"present: (n,t) cross section MT=105")
206 insist(mt_set(19)==106,"present: (n,He-3) cross section MT=106")
207 insist(mt_set(20)==107,"present: (n,ALPHA) cross section MT=107")
208 insist(mt_set(21)==108,"present: (n,2ALPHA) cross section MT=108")
209 insist(mt_set(22)==111,"present: (n,2p) cross section MT=111")
210 insist(mt_set(23)==112,"present: (n,pALPHA) cross section MT=112")
211 
212 !check the nuclide guess functionality
213 allocate(myids(2))
214 myids(1) = lib%find_nuclide_guess(922350)
215 myids(2) = lib%find_nuclide_guess(922380)
216 
217 !get a bunch of data
218 call lib%get_nucl(nucl)
219 call lib%get_typ_nuc(typ)
220 call lib%get_fiss(fiss)
221 call lib%get_tocap(tocap)
222 call lib%get_dis(dis)
223 call lib%get_kd(kd) !last index of decay parents
224 call lib%get_non0(non0) !last index of parents
225 call lib%get_non0_(non0_) !total number of parents
226 call lib%get_kd_(kd_) !total number of decay parents
227 call lib%get_a(a)
228 call lib%get_loc(loc)
229 call lib%get_mt(mt)
230 ilite=lib%get_ilite()
231 iact =lib%get_iact()
232 
233 !you could change this loop to any ids (just populate myids)
234 do j=1,size(myids)
235 
236  i=myids(j)
237 
238  if( i>=(1+ilite) .AND. i<=(ilite+iact) )then
239  if( typ(i)/=origen_sublib_2ac )then
240  write(error_unit,*)'index i=',i,' should have been an actinide?'
241  stop 10
242  end if
243  myfis=fiss(i-ilite)
244  else
245  myfis=0
246  end if
247 
248  decay_start=kd(i)-kd_(i)+1
249  decay_end=kd(i)
250  reaction_start=kd(i)+1
251  reaction_end=non0(i)
252 
253  expect_eq(decay_start, lib%first_transition(i))
254  expect_eq(decay_end, lib%last_transition(i,0.0d0))
255  expect_eq(reaction_end, lib%last_transition(i,1.0d0))
256 
257  expect_eq(reaction_end-decay_start+1, non0_(i))
258 
259  call super_daughter_check(lib,i)
260 
261 end do
262 
263 ! check num_massive daughters
264 nmd = lib%num_massive_daughters(51);
265 expect_eq(1,nmd)
266 
267 !! check actual getting/setting capability
268 u235=myids(1)
269 
270 expect_near(31.250144958496094_c_double, lib%get_xs(u235,18), 1d-5)
271 expect_near(7.8026990890502930_c_double, lib%get_xs(u235,102), 1.0d-5)
272 expect_near(7.8109970092773438_c_double, lib%get_xs(u235,101), 1.0d-5)
273 expect_near(39.061141967773438_c_double, lib%get_xs(u235,27), 1.0d-5)
274 
275 call lib%set_xs(u235,102,10.d0)
276 expect_eq(10.d0, lib%get_xs(u235,102))
277 
278 call lib%set_xs(u235,18,60.d0)
279 expect_eq(60.d0, lib%get_xs(u235,18))
280 expect_near(70.008297920227051_c_double, lib%get_xs(u235,27), 1d-5)
281 
282 call lib%find_parents(u235,parents,mts_p,tinds_p)
283 call lib%find_daughters(u235,daughters,mts_d,tinds_d)
284 
285 !set nuclide/mt pairs (use previous myids as nuclide indices)
286 allocate(mymts(2))
287 mymts(1)=18
288 mymts(2)=102
289 call tinds_2d%initialize()
290 call lib%setup_fast_transition_access(myids,mymts,tinds_2d)
291 
292 !deliver new xs values
293 allocate(myxs(2))
294 myxs(1)=34.3434
295 myxs(2)=6.0102
296 call lib%fast_set_xs(myids,mymts,tinds_2d,myxs)
297 
298 deallocate(parents,mts_p,tinds_p,daughters,mts_d,tinds_d)
299 
300 !destroy
301 call lib%destroy()
302 call finalize
303 
304 end program
305 
306 
Generated by   doxygen 1.8.10