tstExternalMultistep.f90

./Manager/wrap/tests/tstExternalMultistep.f90

program tstexternalmultistep

use origen_casewrapper_m
use origen_librarywrapper_m

#include "ScaleSTL/FortranTestMacros.h"
use origen_definitions_m, only: approx_eq

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

implicit none

type(origen_casewrapper)    :: casew1,casew2
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)

!important variables
integer,allocatable :: all_nucl_in(:), all_typ_in(:)
integer,allocatable :: all_nucl_ex(:), all_typ_ex(:)
real(C_DOUBLE),allocatable :: conc_in(:,:),conc_ex(:,:)

integer :: i,j
real(C_DOUBLE) :: dt

logical :: file_exists

call initialize
call build_types
if ( node() /=0 ) call finalize

inquire(file=ce14_e15_filepath,exist=file_exists)

expect_true(file_exists)

!! initialize times and fluxes.
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;

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


!! Option 1: Let solver make multi-steps.
call libw%initialize(ce14_e15_filepath)
call libw%set_pos(1)
call casew1%initialize(libw)
itot = casew1 % total_nuclides()
call casew1%set_times(times)
call casew1%set_fluxes(fluxes)
call casew1%set_initial_concentrations(nnucl0,nucl0,conc0,typ0,origen_concentrationunit_gatoms)
call casew1%libw%get_nucl(all_nucl_in)
call casew1%libw%get_typ_nuc(all_typ_in)
call casew1%run()
call casew1%get_concentrations(conc_in)



!! Option 2: Make our own external multi-steps.
allocate( conc_ex(0:nsteps,1:itot) )
allocate(conc_buffer1(itot))
call libw%initialize(ce14_e15_filepath) !library must be initialized and loaded OUTSIDE the loop
do i=1,nsteps
   call libw%set_pos(1) !get new cross sections!
   expect_eq(1,libw%get_pos())
   call casew2%initialize(libw)
   itot = casew2 % total_nuclides()
   call casew2%set_times(times(i-1:i)) !only 1 step
   call casew2%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 casew2%set_initial_concentrations(nnucl0,nucl0,conc0,typ0,origen_concentrationunit_gatoms)
      call casew2%libw%get_nucl(all_nucl_ex)
      call casew2%libw%get_typ_nuc(all_typ_ex)

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

   !the big run
   call casew2%run()

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

end do

!destroy
call casew1%destroy()
call casew2%destroy()
call libw%destroy()

!! Check consistency (but not exact because we change cross sections each time)
do j=1,itot
   expect_eq(all_nucl_in(j), all_nucl_ex(j))
   expect_eq(all_typ_in(j), all_typ_ex(j))
   do i = 0, nsteps
      expect_eq(conc_in(i,j), conc_ex(i,j))
   end do
end do

call finalize


1001 format('*** Error: nuclide ids do not match for nuclide number j=',i5.5,&
'internal (',i7,') vs. external (',i7,')')
1002 format('*** Error: nuclide types do not match for nuclide number j=',i5.5,&
'internal (',i7,') vs. external (',i7,')')
1003 format('*** Error: nuclide concentrations do not match for time step i=',i5.5,' nuclide number j=',i5.5,&
'internal (',i7,') vs. external (',i7,')')
end program