tstCRAMInterface.f90

./Manager/wrap/tests/tstCRAMInterface.f90

!
! Unit test for using the CRAM solver in Origen.
! Author: Aarno Isotalo
!
program tstcraminterface

#include "ScaleSTL/FortranTestMacros.h"

  use origen_casewrapper_m
  use origen_librarywrapper_m
  !2. Add dynamic path info for ARP library.
  use origen_testpaths_m
  use origen_definitions_m, only: approx_eq
  !
  use nemesis_comm, only: build_types, initialize, finalize, node

  implicit none

  ! timing
  integer :: count1, count2, count_rate

  !instance of origen calc class
  !instance of origen library class
  type(origen_casewrapper)    :: casew1, casew2
  type(origen_librarywrapper) :: libw

  !some local variables
  integer :: itot,nsteps
  real(C_DOUBLE), allocatable :: times(:)
  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 ones
  real(C_DOUBLE),allocatable :: conc_fw(:,:), conc_ad(:,:)

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


  !solver options
  integer :: solver, nterm, nshrt, cram_order, internal_steps
  real(C_DOUBLE) :: abstol, reltol

  real(C_DOUBLE) :: dot1, dot2

  call initialize
  call build_types 
  if(node() /= 0 ) call finalize
  
  nsteps=1
  allocate(times(0:nsteps),fluxes(1:nsteps))
  times(0)=0.d0
  times(1)=100.0
  times=times*86400.d0 !units from days to seconds
  fluxes=1.e14;

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

  call libw%initialize(pwrlib_filepath)
  call libw%set_pos(1)

  call casew1%initialize(libw)
  call casew1%set_times(times)
  call casew1%set_fluxes(fluxes)
  call casew1%set_initial_concentrations(nnucl0,nucl0,conc0,typ0,origen_concentrationunit_gatoms)

  itot = casew1 % libw % get_itot()

  ! Set the solver options (order 14 cram, normally use order 16 (default))
  ! Also set the use of internal substeps (4)
  call casew1%b_case%set_solver_options(2, -1, -2, -3.0d0, -4.0d0, 14, 4)

  ! Check that the options were correctly set
  call casew1%b_case%get_solver_options(solver, nterm, nshrt, abstol, reltol, cram_order, internal_steps)

  expect_eq(2, solver)
  expect_eq(-1, nterm)
  expect_eq(-2, nshrt)
  expect_eq(-3.0d0, abstol)
  expect_eq(-4.0d0, reltol)
  expect_eq(14, cram_order)
  expect_eq(4, internal_steps)

  ! run and time origen
  call system_clock(count1, count_rate)

  call casew1%run()
  call casew1%get_concentrations(conc_fw)

  call system_clock(count2, count_rate)

  Write(*,*) "ORIGEN RUN TOOK (ms):", dble(count2-count1)/dble(count_rate)*1000.0

  ! check the results for a few nuclides to see that the calculation actually
  ! used CRAM

  expect_near(1.1397904580875401d-5, conc_fw(1,libw%find_nuclide(origen_sublib_2ac, 942400)), 1d-10)
  expect_near(3.1582984913113634d-7, conc_fw(1,libw%find_nuclide(origen_sublib_3fp,  541350)), 1d-10)

  ! write(*,*) conc_fw(1,libw%find_nuclide(Origen_SUBLIB_2AC, 942400))
  ! write(*,*) conc_fw(1,libw%find_nuclide(Origen_SUBLIB_3FP,  541350))

  ! Test adjoint calculation

  nucl0=[942390,942400]
  conc0=[2.00d0,1.5d0]
  typ0=[origen_sublib_2ac,origen_sublib_2ac]

  call casew2%initialize(libw,.true.) !use adjoint!
  call casew2%set_times(times)
  call casew2%set_fluxes(fluxes)
  call casew2%set_initial_concentrations(nnucl0,nucl0,conc0,typ0,origen_concentrationunit_gatoms)
  call casew2%b_case%set_solver_options(2, -1, -2, -3.0d0, -4.0d0, 14, 2)
  call casew2%run()
  call casew2%get_concentrations(conc_ad)

  dot1 = 0.0d0
  dot2 = 0.0d0
  do i=1,itot
    dot1 = dot1 + conc_fw(0,i)*conc_ad(1,i);
    dot2 = dot2 + conc_fw(1,i)*conc_ad(0,i);
  end do

  expect_near(1d0, dot1/dot2, 1d-10)

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

  deallocate(times,fluxes)

  call finalize
  
end program tstcraminterface