tstXSChange.f90

./Manager/wrap/tests/tstXSChange.f90

program tstxschange

use origen_casewrapper_m
use origen_librarywrapper_m

use origen_testpaths_m
use dbcf_m
!
use nemesis_comm, only: build_types, initialize, finalize, node


#include "ScaleSTL/FortranTestMacros.h"

implicit none

type(origen_casewrapper)    :: casew
type(origen_librarywrapper) :: libw


integer :: itot,nsteps
real(C_DOUBLE), allocatable :: times(:)
real(C_DOUBLE),allocatable :: conc_buffer1(:),conc_buffer2(:,:)
real(C_DOUBLE), allocatable :: fluxes(:) !will depend on number of steps
integer,parameter :: nnucl0=2
real(C_DOUBLE) :: conc0(nnucl0)
integer :: nucl0(nnucl0), typ0(nnucl0)
integer,allocatable :: myids(:)
real(C_DOUBLE),allocatable :: myxs(:)
integer :: kfound
logical :: mt_set(200)

integer,allocatable :: all_nucl(:), all_typ(:),all_mt(:),all_loc(:)
real(C_DOUBLE),allocatable :: conc(:,:),all_tocap(:),all_fiss(:)

integer :: i,j,ilite,iact,ifp
real(C_DOUBLE) :: dt,myfis
! CHARACTER(len=255) :: cwd
!             CALL getcwd(cwd)
!             write(ERROR_UNIT,*)TRIM(cwd)

call initialize
call build_types
if( node() /= 0 ) call finalize
!initialize input
nsteps=3
allocate(times(0:nsteps),fluxes(1:nsteps))
dt=100
times(0)=0.d0
do i=1,nsteps
    times(i)=times(i-1)+dt
end do
times=times*86400.d0 !units from days to seconds
fluxes=1.e14;

nucl0=[922350,922380]
conc0=[0.05d0,0.95d0]
typ0=[origen_sublib_2ac,origen_sublib_2ac] !Origen_SUBLIB_1LT,Origen_SUBLIB_3FP also available



!get data
call libw % initialize(pwrlib_filepath) !library must be initialized and loaded OUTSIDE the loop
itot=libw % get_itot()
allocate( conc(0:nsteps,1:itot) )
allocate(conc_buffer1(itot))
allocate(myids(2),myxs(2))
myids(1) = libw % find_nuclide_guess(922350)
myids(2) = libw % find_nuclide_guess(922380)
myxs(1) = libw % get_xs(myids(1),18)

call test_times_array(libw)

if( .false. )then
!  write(ERROR_UNIT,*)'my U-235 guess=',myids(1)
!  write(ERROR_UNIT,*)'my U-238 guess=',myids(2)
!  write(ERROR_UNIT,*)'my U-235 fission xs=',myxs(1),' kfound='

!  write(ERROR_UNIT,*)'print_transitions'
  call libw % print_transitions_for(922350,18);
!  write(ERROR_UNIT,*)'print info'
  call libw % print_info();
!  write(ERROR_UNIT,*)'print library'
  call libw % library_printout(.true.,.true.,.true.);
end if

call libw % get_nucl(all_nucl)
call libw % get_typ_nuc(all_typ)
call libw % get_fiss(all_fiss)
call libw % get_tocap(all_tocap)
!write(ERROR_UNIT,*)'size(all_nucl)=',size(all_nucl)
!write(ERROR_UNIT,*)'size(all_typ)=',size(all_typ)

!write(ERROR_UNIT,*)'for all isotopes:'
!write(ERROR_UNIT,'(a5,4(x,a12))')'n','nucl','type','fiss','tocap'
!1:ilite Origen_SUBLIB_1LT
!ilite+1:ilite+iact Origen_SUBLIB_2AC
!ilite+iact+1:ilite+iact+ifp Origen_SUBLIB_3FP
ilite=libw % get_ilite()
iact=libw % get_iact()
ifp=libw % get_ifp()
do i=1,itot
  if( i>=(1+ilite) .AND. i<=(ilite+iact) )then
    expect_eq(origen_sublib_2ac, all_typ(i))  ! Should be actinide
    myfis=all_fiss(i-ilite)
  else
    myfis=0
  end if
!  write(ERROR_UNIT,'(i5,2(x,i12),2(x,es12.5))')i,all_nucl(i),all_typ(i),myfis,all_tocap(i)
end do

call libw % get_mt(all_mt)
call libw % get_loc(all_loc)
!write(ERROR_UNIT,*)'size(all_mt)=',size(all_mt)
!write(ERROR_UNIT,*)'size(all_loc)=',size(all_loc)

if( .true. )then
!  write(ERROR_UNIT,*)'for all nonzeros in the matrix:'
!  write(ERROR_UNIT,'(a5,3(x,a12))')'n','mt','loc','parent'
  do i=1,libw % get_non()
!    write(ERROR_UNIT,'(i5,3(x,i12))')i,all_mt(i),all_loc(i),all_nucl(all_loc(i))
    mt_set(abs(all_mt(i)))=.true.
  end do
  mt_set=.false.
  do i=1,200
  if(mt_set(i))write(error_unit,*)mt_set(i)
  end do
end if

!stop

do i=1,nsteps
   call casew % initialize(libw)
   expect_eq(itot,casew % total_nuclides())
   call casew % set_times(times(i-1:i)) !only 1 step
   call casew % set_fluxes(fluxes(i:i))

   !we need a special condition for the first step we will pass in the initial concentrations
   !which are a much shorter list and we need to get the full nuclide ids and types lists
   !for future steps
   if( i==1 )then
      call casew % set_initial_concentrations(nnucl0,nucl0,conc0,typ0,origen_concentrationunit_gatoms)

   !now we are in the general looping regime
   else
      conc_buffer1=conc(i-1,1:itot) !pass in beginning of time step concentrations
      call casew % set_initial_concentrations_by_vector(conc_buffer1,origen_concentrationunit_gatoms)
   endif

   !the big run
   call casew % run()

   !extract beginning and end-of step
   call casew % get_concentrations(conc_buffer2)
   conc(i-1:i,1:itot)=conc_buffer2

end do

!destroy
call casew % destroy()
call libw % destroy()

call finalize

contains

subroutine test_times_array(libwt)
type(origen_librarywrapper) :: libwt

type(origen_casewrapper) :: casewt
real(C_DOUBLE),allocatable :: times(:)
real(C_DOUBLE) :: times1a(0:1)
real(C_DOUBLE) :: times2a(1:2)
real(C_DOUBLE) :: times1b(0:2)
real(C_DOUBLE) :: times2b(1:3)


call casewt%initialize(libwt)

!1a
times1a = [0d0,1d0]
call casewt%set_times(times1a)
expect_eq(1,casewt%nsteps())
call casewt%get_times(times)
expect_eq(0.0,times(0))
expect_eq(1.0,times(1))
expect_eq(2,size(times,1))

!2a (same result from shifted times)
times2a = [0d0,1d0]
call casewt%set_times(times2a)
expect_eq(1,casewt%nsteps())
call casewt%get_times(times)
expect_eq(0.0,times(0))
expect_eq(1.0,times(1))
expect_eq(2,size(times,1))

!1b
times1b = [0d0,1d0,2d0]
call casewt%set_times(times1b)
expect_eq(2,casewt%nsteps())
call casewt%get_times(times)
expect_eq(0.0,times(0))
expect_eq(1.0,times(1))
expect_eq(2.0,times(2))
expect_eq(3,size(times,1))

!2b (same result from shifted times)
times2b = [0d0,1d0,2d0]
call casewt%set_times(times2b)
expect_eq(2,casewt%nsteps())
call casewt%get_times(times)
expect_eq(0.0,times(0))
expect_eq(1.0,times(1))
expect_eq(2.0,times(2))
expect_eq(3,size(times,1))

call casewt%destroy()
end subroutine

end program