Commit 0993c117 authored by Messer II, Bronson's avatar Messer II, Bronson
Browse files

fixed interpolation test

parent 15823654

microphysics/InterpolationTest/Makefile

+4 −2
Original line number Diff line number Diff line
@@ -19,6 +19,7 @@ WL_OBJECTS = \ 
  $(WEAKLIB_DIR)/Distributions/EOSSource/wlEquationOfStateTableModule.o \

  $(WEAKLIB_DIR)/Distributions/EOSSource/wlEOSIOModuleHDF.o \

  $(WEAKLIB_DIR)/Distributions/OpacitySource/wlOpacityFieldsModule.o \

  $(WEAKLIB_DIR)/Distributions/OpacitySource/wlOpacityInterpolationModule.o \

  $(WEAKLIB_DIR)/Distributions/OpacitySource/wlOpacityTableModule.o \

  $(WEAKLIB_DIR)/Distributions/OpacitySource/wlOpacityTableIOModuleHDF.o
@@ -36,8 +37,9 @@ wlInterpolationTest: \ 
	wlDependentVariablesModule.o wlInterpolationModule.o \

	wlIOModuleHDF.o wlEquationOfStateTableModule.o \

	wlEOSIOModuleHDF.o wlOpacityFieldsModule.o \

	wlOpacityTableModule.o wlOpacityTableIOModuleHDF.o \

	wlInterpolationTest.o $(LIBRARIES)

	wlOpacityInterpolationModule.o wlOpacityTableModule.o \

	wlOpacityTableIOModuleHDF.o wlInterpolationTest.o \

	$(LIBRARIES)



clean:

	rm -f *.o *.mod *.ld

microphysics/InterpolationTest/wlInterpolationTest.f90

+157 −21
Original line number Diff line number Diff line
@@ -20,7 +20,11 @@ PROGRAM wlInterpolationTest 
  USE wlOpacityTableIOModuleHDF, ONLY: &

    ReadOpacityTableHDF

  USE wlInterpolationModule, ONLY: &

    LogInterpolateSingleVariable, &

    LogInterpolateSingleVariable_1D3D

  USE wlOpacityInterpolationModule, ONLY: &

    wlInterpolateOpacity_NES, &

    wlInterpolateOpacity_Pair



  USE HDF5
@@ -31,11 +35,18 @@ PROGRAM wlInterpolationTest 
  INTEGER,     PARAMETER :: iD = 1

  INTEGER,     PARAMETER :: iT = 2

  INTEGER,     PARAMETER :: iY = 3

  INTEGER,     PARAMETER :: iMu_e = 4

  INTEGER,     PARAMETER :: nSpecies = 2

  INTEGER                :: iS

  REAL(dp),    PARAMETER :: kmev = 8.61733d-11 ! [ MeV K^{-1} ]

  INTEGER                :: iX, iE, iS

  REAL(dp)               :: xR(nx), xD(nx), xT(nx), xY(nx)

  REAL(dp)               :: xMu_e(nx), xEta(nX)

  REAL(dp)               :: xOpacAbs (ne,nx,nSpecies)

  REAL(dp)               :: xOpacSct (ne,nx,nSpecies)

  REAL(dp)               :: xOpacNes (ne,nx,nSpecies)

  REAL(dp)               :: xOpacPair(ne,nx,nSpecies)

  REAL(dp)               :: xPhiNES (ne,ne,2,nSpecies,nx)

  REAL(dp)               :: xPhiPair(ne,ne,2,nSpecies,nx)

  TYPE(GridType)         :: xEnergyGrid

  TYPE(OpacityTableType) :: OPACITIES
@@ -49,8 +60,8 @@ PROGRAM wlInterpolationTest 


  xEnergyGrid % Unit      = 'MeV'

  xEnergyGrid % Name      = 'Energy Grid'

  xEnergyGrid % MinValue  = 1.0d0

  xEnergyGrid % MaxValue  = 3.0d2

  xEnergyGrid % MinValue  = 1.05d0

  xEnergyGrid % MaxValue  = 2.95d2

  xEnergyGrid % LogInterp = 1



  CALL MakeLogGrid &
@@ -66,22 +77,44 @@ PROGRAM wlInterpolationTest 
  CALL InitializeHDF( )



  CALL ReadOpacityTableHDF &

         ( OPACITIES, FileName = "wl-Op-SFHo-15-25-50-AbIs.h5", &

         ( OPACITIES, &

           FileName_EmAb_Option = 'wl-Op-SFHo-15-25-50-E40-B85-AbEm.h5', &

           FileName_Iso_Option  = 'wl-Op-SFHo-15-25-50-E40-B85-Iso.h5',  &

           FileName_NES_Option  = 'wl-Op-SFHo-15-25-50-E40-B85-NES.h5',  &

           FileName_Pair_Option = 'wl-Op-SFHo-15-25-50-E40-B85-Pair.h5', &

           Verbose_Option = .TRUE. )



  CALL FinalizeHDF( )



  ! --- Interpolate Electron Chemical Potential ---



  ASSOCIATE &

    ( EosMu_e => OPACITIES % EOSTable % DV % Variables(iMu_e) % Values, &

      OSMu_e  => OPACITIES % EOSTable % DV % Offsets  (iMu_e),          &

      EosD    => OPACITIES % EOSTable % TS % States   (iD   ) % Values, &

      EosT    => OPACITIES % EOSTable % TS % States   (iT   ) % Values, &

      EosY    => OPACITIES % EOSTable % TS % States   (iY   ) % Values )



  CALL LogInterpolateSingleVariable &

         ( xD, xT, xY, EosD, EosT, EosY, [1,1,0], OSMu_e, EosMu_e, xMu_e )



  xEta = xMu_e / ( kmev * xT )



  END ASSOCIATE



  ! --- Interpolate Opacities ---



  ! --- EmAb and Iso Scattering ---



  DO iS = 1, nSpecies



    ASSOCIATE &

      ( OpacAbs   => OPACITIES % thermEmAb % Absorptivity(iS) &

      ( OpacAbs   => OPACITIES % EmAb % Opacity(iS) &

                               % Values(:,:,:,:), &

        OpacAbsOS => OPACITIES % thermEmAb % Offsets(iS), &

        OpacSct   => OPACITIES % scatt_Iso % Kernel(iS) &

                               % Values(:,:,:,:,1), &

        OpacSctOS => OPACITIES % scatt_Iso % Offsets(iS,1), &

        OpacAbsOS => OPACITIES % EmAb % Offsets(iS), &

        OpacSct   => OPACITIES % scat_Iso % Kernel(iS) &

                               % Values(:,1,:,:,:), &

        OpacSctOS => OPACITIES % scat_Iso % Offsets(iS,1), &

        OpacE     => OPACITIES % EnergyGrid % Values(:), &

        EosD      => OPACITIES % EOSTable % TS % States(iD) % Values(:), &

        EosT      => OPACITIES % EOSTable % TS % States(iT) % Values(:), &
@@ -99,12 +132,75 @@ PROGRAM wlInterpolationTest 


  END DO



  ! --- NES Kernels ---



  ASSOCIATE &

    ( OpacNES => OPACITIES % Scat_NES  % Kernel(1) % Values(:,:,:,:,:), &

      OpacOS  => OPACITIES % Scat_NES  % Offsets(1,:), &

      OpacE   => OPACITIES % EnergyGrid % Values(:), &

      EosEta  => OPACITIES % EtaGrid % Values(:), &

      EosT    => OPACITIES % EOSTable % TS % States(iT) % Values(:) )



  CALL wlInterpolateOpacity_NES &

         ( xEnergyGrid % Values, xT * kmev, xEta, OpacE, &

           EosT * kmev, EosEta, OpacOS, OpacNES, xPhiNES )



  END ASSOCIATE



  ASSOCIATE( xE => xEnergyGrid % Values )



  DO iX = 1, nx

  DO iS = 1, nSpecies

  DO iE = 1, ne



    xOpacNES(ie,ix,iS) &

      = TRAPEZ( ne, xE, xPhiNES(:,ie,1,iS,ix) * xE**2 )



  END DO

  END DO

  END DO



  END ASSOCIATE



  ! --- Pair Kernels ---



  ASSOCIATE &

    ( OpacPair => OPACITIES % Scat_Pair % Kernel(1) % Values(:,:,:,:,:), &

      OpacOS   => OPACITIES % Scat_Pair % Offsets(1,:), &

      OpacE    => OPACITIES % EnergyGrid % Values(:), &

      EosEta   => OPACITIES % EtaGrid % Values(:), &

      EosT     => OPACITIES % EOSTable % TS % States(iT) % Values(:) )



  CALL wlInterpolateOpacity_Pair &

         ( xEnergyGrid % Values, xT * kmev, xEta, OpacE, &

           EosT * kmev, EosEta, OpacOS, OpacPair, xPhiPair )



  END ASSOCIATE



  ASSOCIATE( xE => xEnergyGrid % Values )



  DO iX = 1, nx

  DO iS = 1, nSpecies

  DO iE = 1, ne



    xOpacPair(ie,ix,iS) &

      = TRAPEZ( ne, xE, xPhiPair(:,ie,1,iS,ix) * xE**2 &

                          * EXP( - (xE(iE)+xE) / (kmev*xT(iX)) ) )



  END DO

  END DO

  END DO



  END ASSOCIATE





  CALL DeallocateOpacityTable( OPACITIES, Verbose_Option = .TRUE. )



  ASSOCIATE( xE => xEnergyGrid % Values )



  CALL WriteOpacityProfilesHDF &

         ( ne, nx, nSpecies, xE, xR, xD, xT, xY, xOpacAbs, xOpacSct )

         ( ne, nx, nSpecies, xE, xR, xD, xT, xY, xEta, &

           xOpacAbs, xOpacSct, xOpacNes, xOpacPair )



  END ASSOCIATE
@@ -140,12 +236,15 @@ CONTAINS 




  SUBROUTINE WriteOpacityProfilesHDF &

    ( nE, nX, nS, xE, xR, xD, xT, xY, xOpacAbs, xOpacSct )

    ( nE, nX, nS, xE, xR, xD, xT, xY, xEta, xOpacAbs, xOpacSct, &

      xOpacNes, xOpacPair )



    INTEGER,  INTENT(in) :: nE, nX, nS

    REAL(DP), INTENT(in) :: xE(nE), xR(nX), xD(nX), xT(nX), xY(nX)

    REAL(DP), INTENT(in) :: xE(nE), xR(nX), xD(nX), xT(nX), xY(nX), xEta(nx)

    REAL(DP), INTENT(in) :: xOpacAbs (nE,nX,nS)

    REAL(DP), INTENT(in) :: xOpacSct (nE,nX,nS)

    REAL(DP), INTENT(in) :: xOpacNes (nE,nX,nS)

    REAL(DP), INTENT(in) :: xOpacPair(nE,nX,nS)



    INTEGER(HID_T)   :: file_id, group_id

    INTEGER(HSIZE_T) :: dsize1d(1), dsize2d(2)
@@ -165,6 +264,7 @@ CONTAINS 
    CALL WriteHDF( "Density",           xD,   group_id, dsize1d )

    CALL WriteHDF( "Temperature",       xT,   group_id, dsize1d )

    CALL WriteHDF( "Electron Fraction", xY,   group_id, dsize1d )

    CALL WriteHDF( "Eta",               xEta, group_id, dsize1d )



    CALL CloseGroupHDF( group_id )
@@ -173,19 +273,39 @@ CONTAINS 
    CALL OpenGroupHDF( 'Opacities', .true., file_id, group_id )



    dsize2D = [ nE, nX ]



    CALL WriteHDF &

           ( "Absorption Opacity (Electron Neutrino)", &

             xOpacAbs(:,:,1), group_id, dsize2d )



    CALL WriteHDF &

           ( "Absorption Opacity (Electron Antineutrino)", &

             xOpacAbs(:,:,2), group_id, dsize2d )



    CALL WriteHDF &

           ( "Scattering Opacity (Electron Neutrino)", &

             xOpacSct(:,:,1), group_id, dsize2d )



    CALL WriteHDF &

           ( "Scattering Opacity (Electron Antineutrino)", &

             xOpacSct(:,:,2), group_id, dsize2d )



    CALL WriteHDF &

           ( "NES Opacity (Electron Neutrino)", &

             xOpacNes(:,:,1), group_id, dsize2d )



    CALL WriteHDF &

           ( "NES Opacity (Electron Antineutrino)", &

             xOpacNes(:,:,2), group_id, dsize2d )



    CALL WriteHDF &

           ( "Pair Opacity (Electron Neutrino)", &

             xOpacPair(:,:,1), group_id, dsize2d )



    CALL WriteHDF &

           ( "Pair Opacity (Electron Antineutrino)", &

             xOpacPair(:,:,2), group_id, dsize2d )



    CALL CloseGroupHDF( group_id )



    CALL CloseFileHDF( file_id )
@@ -194,4 +314,20 @@ CONTAINS 
  END SUBROUTINE WriteOpacityProfilesHDF





  PURE REAL(dp) FUNCTION TRAPEZ( n, x, y )



    INTEGER,  INTENT(in) :: n

    REAL(dp), INTENT(in) :: x(n), y(n)



    INTEGER :: i



    TRAPEZ = 0.0_dp

    DO i = 1, n - 1

      TRAPEZ = TRAPEZ + 0.5_dp * ( x(i+1) - x(i) ) * ( y(i) + y(i+1) )

    END DO



    RETURN

  END FUNCTION TRAPEZ





END PROGRAM wlInterpolationTest