tstMassConservation.f90

./Manager/wrap/tests/tstMassConservation.f90

1 subroutine test_append_file(file,libw,ns)
2 use origen_casewrapper_m
3 use origen_librarywrapper_m
4 use origen_iofunctions_m
5 
6 #include "ScaleSTL/FortranTestMacros.h"
7 !
8 character(*),intent(in) :: file
9 integer,intent(in) :: ns
10 type(origen_librarywrapper),intent(in) :: libw
11 type(origen_casewrapper) :: casew2
12 integer :: i
13 type(origen_stateset) :: state_set
14 logical :: loaded,saved
15 real(C_DOUBLE), allocatable :: conc_end(:)
16 real(C_DOUBLE), allocatable :: conc_begin(:)
17 
18 !read the f71, load last position
19 call casew2 % initialize(libw)
20 loaded = casew2 % set_initial_concentrations_from_file(file,ns)
21 expect_eq(0,casew2%nsteps())
22 call casew2 % set_times([0.d0])
23 expect_true(loaded)
24 
25 !append to a.f71
26 saved = casew2 % save_states("a.f71",667)
27 expect_true(saved)
28 
29 write(*,*)'ns=',ns
30 !load the current a.f71 and check it
31 !{
32  call state_set%initialize()
33  loaded = origen_loadstateset(state_set,"a.f71")
34  expect_true(loaded)
35  expect_eq(ns+1,state_set%states_size())
36  call state_set%destroy()
37 !}
38 
39 !destroy
40 call casew2 % destroy()
41 
42 end subroutine
43 
44 
45 !-------------------------------------------------------------------------------
46 !-------------------------------------------------------------------------------
47 !-------------------------------------------------------------------------------
48 program tstmassconservation
49 
50 use origen_casewrapper_m
51 use origen_librarywrapper_m
52 use origen_definitions_m,only: avogadro
53 use origen_iofunctions_m
54 
55 #include "ScaleSTL/FortranTestMacros.h"
56 
57 !Add dynamic path info for ARP library
58 use origen_testpaths_m
59 use origen_definitions_m, only: approx_eq
60 
61 !
62 use nemesis_comm, only: build_types, initialize, finalize, node
63 
64 implicit none
65 
66 type(origen_casewrapper) :: casew,casew2
67 type(origen_librarywrapper) :: libw
68 real(C_DOUBLE),parameter :: rel_bal_limit=0.001d0
69 
70 !some local variables
71 integer :: nsteps
72 real(C_DOUBLE), allocatable :: times(:),conc(:,:) !will depend on number of times
73 real(C_DOUBLE), allocatable :: fluxes(:) !will depend on number of steps
74 real(C_DOUBLE), allocatable :: sum_typ(:,:) !will depend on number of steps
75 
76 type(origen_stateset) :: state_set
77 
78 real(C_DOUBLE) :: conc0(2)
79 integer :: ids(2), lib_type(2)
80 integer :: i, j, errors, itot, typ, stat
81 real(C_DOUBLE) :: dt,tdepl,tin,tout,runtime,mass0,rel_bal,mass
82 logical :: saved,loaded
83 real(C_DOUBLE), allocatable :: conc_end(:)
84 real(C_DOUBLE), allocatable :: conc_begin(:)
85 
86 logical :: file_exists
87 
88 call initialize
89 call build_types
90 if ( node() /= 0 ) call finalize
91 
92 inquire(file=ce14_e15_filepath,exist=file_exists)
93 expect_true(file_exists)
94 
95 !!! Initialize library, load the first burnup, and print some info.
96 call libw%initialize(ce14_e15_filepath)
97 !call libw% print_info()
98 
99 !!! initialize origen
100 call casew%initialize(libw)
101 
102 !!! Create time steps, set them and fluxes, and deallocate temporaries.
103 tdepl=1200
104 nsteps=5
105 allocate(times(0:nsteps),fluxes(1:nsteps))
106 dt=tdepl/nsteps
107 times(0)=0.d0
108 do i=1,nsteps
109  times(i)=times(i-1)+dt
110 end do
111 times=times*86400.d0 !units from days to seconds
112 call casew%set_times(times)
113 deallocate(times)
114 
115 fluxes=1.e14;
116 call casew%set_fluxes(fluxes)
117 deallocate(fluxes)
118 expect_false(casew%get_decay_only())
119 
120 !!! Declare nuclide ids, initial concentration (in grams), and library type.
121 ids=[922350,922380]
122 conc0=[5.d0,95.d0]
123 lib_type=[origen_sublib_2ac,origen_sublib_2ac] !Origen_SUBLIB_1LT,Origen_SUBLIB_3FP also available
124 call casew%set_initial_concentrations(2,ids,conc0,lib_type,origen_concentrationunit_grams)
125 
126 !!! Run solver.
127 call cpu_time(tin)
128 call casew%run()
129 call cpu_time(tout)
130 
131 !!! Get number of steps (according to casew).
132 expect_eq(nsteps,casew % nsteps())
133 nsteps = casew % nsteps()
134 itot = casew % total_nuclides()
135 
136 !!! Retrieve times and fluxes.
137 call casew % get_times(times)
138 call casew % get_fluxes(fluxes)
139 
140 !!! Get concentrations (in native mole units).
141 allocate( conc(0:nsteps,1:itot) )
142 call casew%get_concentrations(conc)
143 
144 !! convert concentrations (in place) back to masses
145 call convert_concentrations(casew,conc,origen_concentrationunit_grams)
146 
147 !! accumulate each type
148 allocate( sum_typ(0:nsteps,1:3) )
149 sum_typ=0.d0
150 do j = 1, itot
151  typ=casew%libw%typ_nuc(j)
152  do i = 0, nsteps
153  sum_typ(i,typ)=sum_typ(i,typ)+conc(i,j)
154  end do
155 end do
156 
157 
158 !!!
159 !write(ERROR_UNIT,*)'total MASS (Grams) by TYPE at each time step'
160 !write(ERROR_UNIT,'(6a12)',ADVANCE='yes')'time(d)','LITE','ACT','FP','TOT','REL_BAL'
161 mass0=sum(conc0)
162 do i = 0, nsteps
163  mass=sum(sum_typ(i,:))
164  rel_bal=mass/mass0-1.d0
165  expect_near(mass, mass0, rel_bal_limit)
166 end do
167 
168 !delete existing a.f71
169 open(unit=71, iostat=stat, file="a.f71", status='old')
170 if (stat == 0) close(71, status='delete')
171 
172 !write a.f71
173 saved = casew % save_states("a.f71",666)
174 expect_true(saved)
175 
176 !destroy after saving vector at end
177 call casew % get_concentrations_at_time(conc_end,nsteps)
178 
179 !load the current a.f71 and check it
180 !{
181  call state_set%initialize()
182  loaded = origen_loadstateset(state_set,"a.f71")
183  expect_true(loaded)
184  expect_eq(nsteps+1,state_set%states_size())
185  call state_set%destroy()
186 !}
187 
188 !read the f71, last position (nsteps+1)
189 call casew2 % initialize(libw)
190 loaded = casew2 % set_initial_concentrations_from_file("a.f71",nsteps+1)
191 expect_true(loaded)
192 
193 !get vector at beginning
194 call casew2 % get_concentrations0(conc_begin)
195 call casew2 % destroy()
196 
197 !check that the first vector in case2 with last vector in case1
198 do i=1,size(conc_end,1)
199  expect_eq(conc_end(i),conc_begin(i))
200 end do
201 
202 !!test file
203 call test_append_file("a.f71",libw,nsteps+1)
204 
205 !destroy
206 call casew % destroy()
207 call casew2 % destroy()
208 call libw%destroy()
209 
210 call finalize
211 contains
212 
213 subroutine convert_concentrations(this,concentrations,units)
214  use auxiliary
215 
216  !> the receiver
217  type(origen_casewrapper), intent(in) :: this
218  !> existing concentrations [0:nsteps, itot]
219  real(C_DOUBLE) :: concentrations(:,:)
220  integer, intent(in) :: units
221 
222  real(C_DOUBLE),allocatable :: conc(:,:)
223  integer:: i, n
224  real(C_DOUBLE) :: wt
225  integer,allocatable :: nucl(:)
226  real,allocatable :: dis(:),m(:)
227  !should verify bounds of concentrations are correct
228  !0:this%nsteps(),1:this%libw%get_itot()
229 
230  call this % get_concentrations(conc)
231  call this % libw % get_nucl(nucl)
232  call this % libw % get_dis(dis)
233  call this % libw % get_m(m)
234 
235  do i=1,this%libw%get_itot()
236 
237  select case(units)
238  case(origen_concentrationunit_gatoms)
239  !gram-atoms (old unit name for moles) are native units
240  wt = 1.0
241  case(origen_concentrationunit_atoms_b_cm)
242  !Atoms/barn-cm
243  wt = 1.0e24 / avogadro
244  case(origen_concentrationunit_grams)
245  !Grams
246  wt = 1.0 / m(i)
247  case(origen_concentrationunit_curies)
248  !Curies
249  wt = 1.6283e+13 * dis(i)
250  !Insist( (wt > 0), "Curie input (units=4) not allow for stable nuclides in ORIGEN");
251  case default
252  write(error_unit,*)"*** Error: Invalid concentration units!"
253  stop 1
254  end select
255 
256  !convert
257  do n=1,this%nsteps()+1
258  concentrations(n,i)=conc(n-1,i) / wt !NOW WE DIVIDE TO UNDO
259  end do
260 
261  end do
262 
263 end subroutine convert_concentrations
264 
265 end program
Generated by   doxygen 1.8.10