Commit 4f3c8436 authored by Eirik Endeve's avatar Eirik Endeve
Browse files

nu/nuBar tables and interpolation program. Also reorganization of directory structures.

parent ee2bcf58

microphysics/InterpolationTest/Makefile

0 → 100644
+49 −0
Original line number Diff line number Diff line 
VPATH = .

WEAKLIB_DIR = /Users/ece/Documents/weaklib

FORTRAN     = mpif90

FLINKER     = mpif90

FLAGS       = -g -framework Accelerate -fbacktrace -ffpe-trap=invalid,zero -ffpe-summary=invalid,zero -fbounds-check

LIB_HDF5    = -L/usr/local/Cellar/hdf5/1.8.16_1/lib -lhdf5_fortran -lhdf5

INC_HDF5    = -I/usr/local/Cellar/hdf5/1.8.16_1/include



LIBRARIES   = $(LIB_HDF5)

INCLUDE     = $(INC_HDF5)



WL_OBJECTS = \

  $(WEAKLIB_DIR)/Distributions/Library/wlKindModule.o \

  $(WEAKLIB_DIR)/Distributions/Library/wlGridModule.o \

  $(WEAKLIB_DIR)/Distributions/Library/wlThermoStateModule.o \

  $(WEAKLIB_DIR)/Distributions/Library/wlDependentVariablesModule.o \

  $(WEAKLIB_DIR)/Distributions/Library/wlInterpolationModule.o \

  $(WEAKLIB_DIR)/Distributions/Library/wlIOModuleHDF.o \

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

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

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

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

  $(WEAKLIB_DIR)/Distributions/OpacitySource/wlOpacityTableIOModuleHDF.o



.SUFFIXES: .f95 .f90 .F90 .f



.f90.o:

	$(FORTRAN) -c $(FLAGS) $(INCLUDE) $(SUFFIX_f90) $<



all: wlInterpolationTest



wlInterpolationTest: \

	$(WL_OBJECTS) wlInterpolationTest.o

	$(FLINKER) $(FLAGS) -o wlInterpolationTest_exe \

	wlKindModule.o wlGridModule.o wlThermoStateModule.o \

	wlDependentVariablesModule.o wlInterpolationModule.o \

	wlIOModuleHDF.o wlEquationOfStateTableModule.o \

	wlEOSIOModuleHDF.o wlOpacityFieldsModule.o \

	wlOpacityTableModule.o wlOpacityTableIOModuleHDF.o \

	wlInterpolationTest.o $(LIBRARIES)



clean:

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



clobber: clean

	rm -f wlInterpolationTest_exe



wlInterpolationTest.o: \

  wlInterpolationTest.f90

microphysics/InterpolationTest/MatterProfiles.d

0 → 100644
+213 −0
Original line number Diff line number Diff line 
  2.275E+07   1.005E+08   6.909E+09   4.980E-01  

  2.248E+07   1.024E+08   6.925E+09   4.980E-01  

  2.222E+07   1.039E+08   6.941E+09   4.980E-01  

  2.196E+07   1.060E+08   6.958E+09   4.979E-01  

  2.171E+07   1.073E+08   6.976E+09   4.979E-01  

  2.145E+07   1.094E+08   6.994E+09   4.979E-01  

  2.120E+07   1.107E+08   7.010E+09   4.979E-01  

  2.095E+07   1.129E+08   7.026E+09   4.978E-01  

  2.071E+07   1.142E+08   7.042E+09   4.978E-01  

  2.046E+07   1.165E+08   7.057E+09   4.978E-01  

  2.022E+07   1.178E+08   7.073E+09   4.978E-01  

  1.998E+07   1.203E+08   7.091E+09   4.977E-01  

  1.975E+07   1.218E+08   7.111E+09   4.977E-01  

  1.951E+07   1.244E+08   7.130E+09   4.977E-01  

  1.928E+07   1.260E+08   7.152E+09   4.977E-01  

  1.905E+07   1.288E+08   7.171E+09   4.977E-01  

  1.883E+07   1.305E+08   7.195E+09   4.977E-01  

  1.860E+07   1.334E+08   7.212E+09   4.977E-01  

  1.838E+07   1.354E+08   7.275E+09   4.976E-01  

  1.816E+07   1.372E+08   7.212E+09   4.976E-01  

  1.795E+07   1.496E+08   7.412E+09   4.976E-01  

  1.773E+07   2.592E+08   8.598E+09   4.976E-01  

  1.752E+07   6.216E+08   1.165E+10   4.976E-01  

  1.731E+07   1.033E+09   1.463E+10   4.975E-01  

  1.710E+07   1.344E+09   1.521E+10   4.968E-01  

  1.689E+07   1.543E+09   1.543E+10   4.961E-01  

  1.669E+07   1.620E+09   1.564E+10   4.951E-01  

  1.649E+07   1.702E+09   1.600E+10   4.937E-01  

  1.629E+07   1.785E+09   1.625E+10   4.918E-01  

  1.609E+07   1.924E+09   1.640E+10   4.902E-01  

  1.589E+07   2.065E+09   1.654E+10   4.882E-01  

  1.570E+07   2.172E+09   1.681E+10   4.852E-01  

  1.551E+07   2.275E+09   1.715E+10   4.819E-01  

  1.532E+07   2.422E+09   1.736E+10   4.794E-01  

  1.513E+07   2.607E+09   1.750E+10   4.757E-01  

  1.495E+07   2.754E+09   1.777E+10   4.709E-01  

  1.476E+07   2.852E+09   1.822E+10   4.659E-01  

  1.458E+07   2.986E+09   1.856E+10   4.612E-01  

  1.440E+07   3.196E+09   1.873E+10   4.557E-01  

  1.422E+07   3.410E+09   1.898E+10   4.492E-01  

  1.405E+07   3.507E+09   1.955E+10   4.422E-01  

  1.387E+07   3.649E+09   2.002E+10   4.347E-01  

  1.370E+07   3.847E+09   2.040E+10   4.274E-01  

  1.353E+07   4.044E+09   2.080E+10   4.193E-01  

  1.336E+07   4.193E+09   2.144E+10   4.104E-01  

  1.319E+07   4.291E+09   2.231E+10   4.012E-01  

  1.303E+07   4.418E+09   2.312E+10   3.927E-01  

  1.286E+07   4.617E+09   2.370E+10   3.837E-01  

  1.270E+07   4.772E+09   2.447E+10   3.748E-01  

  1.254E+07   4.887E+09   2.540E+10   3.661E-01  

  1.238E+07   5.018E+09   2.627E+10   3.581E-01  

  1.222E+07   5.200E+09   2.695E+10   3.502E-01  

  1.207E+07   5.405E+09   2.757E+10   3.427E-01  

  1.191E+07   5.600E+09   2.826E+10   3.358E-01  

  1.176E+07   5.806E+09   2.889E+10   3.293E-01  

  1.161E+07   6.020E+09   2.955E+10   3.233E-01  

  1.146E+07   6.233E+09   3.021E+10   3.178E-01  

  1.131E+07   6.438E+09   3.090E+10   3.135E-01  

  1.117E+07   6.662E+09   3.155E+10   3.095E-01  

  1.102E+07   6.942E+09   3.207E+10   3.055E-01  

  1.088E+07   7.249E+09   3.256E+10   3.018E-01  

  1.074E+07   7.558E+09   3.309E+10   2.986E-01  

  1.060E+07   7.862E+09   3.366E+10   2.962E-01  

  1.046E+07   8.200E+09   3.417E+10   2.939E-01  

  1.032E+07   8.574E+09   3.466E+10   2.915E-01  

  1.019E+07   8.994E+09   3.512E+10   2.888E-01  

  1.005E+07   9.432E+09   3.558E+10   2.865E-01  

  9.920E+06   9.887E+09   3.606E+10   2.845E-01  

  9.788E+06   1.035E+10   3.656E+10   2.829E-01  

  9.658E+06   1.085E+10   3.707E+10   2.811E-01  

  9.529E+06   1.140E+10   3.757E+10   2.790E-01  

  9.402E+06   1.202E+10   3.798E+10   2.766E-01  

  9.275E+06   1.270E+10   3.841E+10   2.739E-01  

  9.150E+06   1.345E+10   3.881E+10   2.709E-01  

  9.027E+06   1.427E+10   3.921E+10   2.677E-01  

  8.904E+06   1.516E+10   3.960E+10   2.643E-01  

  8.783E+06   1.614E+10   3.997E+10   2.606E-01  

  8.664E+06   1.721E+10   4.032E+10   2.568E-01  

  8.545E+06   1.836E+10   4.068E+10   2.530E-01  

  8.428E+06   1.963E+10   4.103E+10   2.490E-01  

  8.311E+06   2.099E+10   4.141E+10   2.451E-01  

  8.196E+06   2.249E+10   4.178E+10   2.409E-01  

  8.083E+06   2.413E+10   4.215E+10   2.367E-01  

  7.970E+06   2.589E+10   4.252E+10   2.327E-01  

  7.859E+06   2.782E+10   4.289E+10   2.286E-01  

  7.749E+06   2.991E+10   4.328E+10   2.247E-01  

  7.639E+06   3.218E+10   4.369E+10   2.206E-01  

  7.531E+06   3.465E+10   4.411E+10   2.167E-01  

  7.425E+06   3.734E+10   4.454E+10   2.127E-01  

  7.319E+06   4.026E+10   4.499E+10   2.089E-01  

  7.214E+06   4.344E+10   4.544E+10   2.053E-01  

  7.111E+06   4.688E+10   4.593E+10   2.017E-01  

  7.008E+06   5.063E+10   4.645E+10   1.982E-01  

  6.907E+06   5.473E+10   4.698E+10   1.947E-01  

  6.806E+06   5.919E+10   4.752E+10   1.912E-01  

  6.707E+06   6.403E+10   4.809E+10   1.881E-01  

  6.609E+06   6.927E+10   4.869E+10   1.851E-01  

  6.511E+06   7.498E+10   4.932E+10   1.820E-01  

  6.415E+06   8.121E+10   4.998E+10   1.789E-01  

  6.320E+06   8.800E+10   5.063E+10   1.762E-01  

  6.226E+06   9.526E+10   5.135E+10   1.737E-01  

  6.132E+06   1.037E+11   5.211E+10   1.718E-01  

  6.040E+06   1.125E+11   5.285E+10   1.692E-01  

  5.949E+06   1.221E+11   5.363E+10   1.667E-01  

  5.858E+06   1.327E+11   5.446E+10   1.641E-01  

  5.769E+06   1.442E+11   5.533E+10   1.615E-01  

  5.680E+06   1.569E+11   5.621E+10   1.589E-01  

  5.593E+06   1.704E+11   5.712E+10   1.573E-01  

  5.506E+06   1.851E+11   5.809E+10   1.557E-01  

  5.420E+06   2.010E+11   5.912E+10   1.540E-01  

  5.335E+06   2.185E+11   6.017E+10   1.524E-01  

  5.251E+06   2.375E+11   6.128E+10   1.507E-01  

  5.168E+06   2.583E+11   6.242E+10   1.489E-01  

  5.085E+06   2.808E+11   6.360E+10   1.473E-01  

  5.004E+06   3.049E+11   6.484E+10   1.464E-01  

  4.923E+06   3.311E+11   6.614E+10   1.453E-01  

  4.844E+06   3.597E+11   6.748E+10   1.441E-01  

  4.765E+06   3.909E+11   6.887E+10   1.429E-01  

  4.686E+06   4.250E+11   7.031E+10   1.416E-01  

  4.609E+06   4.615E+11   7.180E+10   1.408E-01  

  4.533E+06   5.009E+11   7.335E+10   1.402E-01  

  4.457E+06   5.438E+11   7.495E+10   1.397E-01  

  4.382E+06   5.904E+11   7.662E+10   1.390E-01  

  4.308E+06   6.414E+11   7.833E+10   1.382E-01  

  4.234E+06   6.966E+11   8.008E+10   1.375E-01  

  4.162E+06   7.557E+11   8.189E+10   1.374E-01  

  4.090E+06   8.202E+11   8.377E+10   1.370E-01  

  4.019E+06   8.908E+11   8.569E+10   1.365E-01  

  3.948E+06   9.679E+11   8.766E+10   1.358E-01  

  3.879E+06   1.052E+12   8.967E+10   1.352E-01  

  3.810E+06   1.142E+12   9.173E+10   1.350E-01  

  3.741E+06   1.242E+12   9.384E+10   1.347E-01  

  3.674E+06   1.350E+12   9.600E+10   1.341E-01  

  3.607E+06   1.470E+12   9.821E+10   1.333E-01  

  3.541E+06   1.601E+12   1.004E+11   1.325E-01  

  3.475E+06   1.744E+12   1.027E+11   1.322E-01  

  3.411E+06   1.900E+12   1.051E+11   1.317E-01  

  3.347E+06   2.072E+12   1.075E+11   1.310E-01  

  3.283E+06   2.262E+12   1.100E+11   1.300E-01  

  3.220E+06   2.472E+12   1.125E+11   1.290E-01  

  3.158E+06   2.700E+12   1.151E+11   1.286E-01  

  3.097E+06   2.950E+12   1.178E+11   1.281E-01  

  3.036E+06   3.226E+12   1.206E+11   1.274E-01  

  2.976E+06   3.529E+12   1.235E+11   1.267E-01  

  2.916E+06   3.861E+12   1.265E+11   1.260E-01  

  2.857E+06   4.219E+12   1.296E+11   1.260E-01  

  2.799E+06   4.608E+12   1.330E+11   1.261E-01  

  2.741E+06   5.031E+12   1.365E+11   1.261E-01  

  2.684E+06   5.491E+12   1.403E+11   1.262E-01  

  2.627E+06   5.981E+12   1.442E+11   1.269E-01  

  2.571E+06   6.504E+12   1.484E+11   1.280E-01  

  2.516E+06   7.065E+12   1.528E+11   1.291E-01  

  2.461E+06   7.670E+12   1.574E+11   1.302E-01  

  2.407E+06   8.306E+12   1.622E+11   1.321E-01  

  2.353E+06   8.981E+12   1.673E+11   1.341E-01  

  2.300E+06   9.704E+12   1.725E+11   1.361E-01  

  2.247E+06   1.048E+13   1.779E+11   1.379E-01  

  2.195E+06   1.130E+13   1.835E+11   1.407E-01  

  2.144E+06   1.216E+13   1.892E+11   1.438E-01  

  2.093E+06   1.308E+13   1.951E+11   1.468E-01  

  2.043E+06   1.407E+13   2.011E+11   1.497E-01  

  1.993E+06   1.511E+13   2.072E+11   1.535E-01  

  1.943E+06   1.622E+13   2.134E+11   1.571E-01  

  1.894E+06   1.742E+13   2.198E+11   1.606E-01  

  1.846E+06   1.870E+13   2.261E+11   1.648E-01  

  1.798E+06   2.003E+13   2.324E+11   1.696E-01  

  1.751E+06   2.147E+13   2.389E+11   1.746E-01  

  1.704E+06   2.302E+13   2.454E+11   1.794E-01  

  1.658E+06   2.469E+13   2.517E+11   1.849E-01  

  1.612E+06   2.646E+13   2.580E+11   1.914E-01  

  1.566E+06   2.837E+13   2.641E+11   1.979E-01  

  1.521E+06   3.040E+13   2.702E+11   2.054E-01  

  1.477E+06   3.259E+13   2.760E+11   2.133E-01  

  1.433E+06   3.498E+13   2.815E+11   2.218E-01  

  1.389E+06   3.755E+13   2.864E+11   2.317E-01  

  1.346E+06   4.036E+13   2.907E+11   2.422E-01  

  1.303E+06   4.345E+13   2.943E+11   2.537E-01  

  1.261E+06   4.692E+13   2.968E+11   2.656E-01  

  1.219E+06   5.094E+13   2.980E+11   2.775E-01  

  1.178E+06   5.574E+13   2.975E+11   2.888E-01  

  1.137E+06   6.186E+13   2.948E+11   2.976E-01  

  1.096E+06   7.009E+13   2.895E+11   3.023E-01  

  1.056E+06   8.189E+13   2.808E+11   3.002E-01  

  1.017E+06   9.930E+13   2.678E+11   2.902E-01  

  9.774E+05   1.232E+14   2.508E+11   2.747E-01  

  9.386E+05   1.516E+14   2.304E+11   2.599E-01  

  9.001E+05   1.807E+14   2.088E+11   2.490E-01  

  8.621E+05   2.076E+14   1.883E+11   2.422E-01  

  8.244E+05   2.315E+14   1.702E+11   2.387E-01  

  7.871E+05   2.523E+14   1.559E+11   2.376E-01  

  7.503E+05   2.701E+14   1.450E+11   2.382E-01  

  7.138E+05   2.859E+14   1.377E+11   2.397E-01  

  6.776E+05   2.994E+14   1.336E+11   2.421E-01  

  6.419E+05   3.115E+14   1.324E+11   2.446E-01  

  6.065E+05   3.225E+14   1.323E+11   2.471E-01  

  5.715E+05   3.325E+14   1.332E+11   2.493E-01  

  5.369E+05   3.418E+14   1.342E+11   2.511E-01  

  5.026E+05   3.507E+14   1.346E+11   2.522E-01  

  4.686E+05   3.591E+14   1.337E+11   2.529E-01  

  4.351E+05   3.670E+14   1.321E+11   2.533E-01  

  4.018E+05   3.743E+14   1.305E+11   2.534E-01  

  3.690E+05   3.809E+14   1.293E+11   2.533E-01  

  3.364E+05   3.868E+14   1.286E+11   2.532E-01  

  3.042E+05   3.921E+14   1.282E+11   2.530E-01  

  2.723E+05   3.967E+14   1.280E+11   2.528E-01  

  2.408E+05   4.008E+14   1.275E+11   2.526E-01  

  2.096E+05   4.043E+14   1.270E+11   2.524E-01  

  1.787E+05   4.073E+14   1.266E+11   2.522E-01  

  1.481E+05   4.097E+14   1.263E+11   2.521E-01  

  1.179E+05   4.117E+14   1.261E+11   2.520E-01  

  8.794E+04   4.130E+14   1.260E+11   2.519E-01  

  5.832E+04   4.141E+14   1.259E+11   2.518E-01  

  2.901E+04   4.146E+14   1.260E+11   2.518E-01

microphysics/InterpolationTest/PlotInterpolatedOpacities.m

0 → 100644
+87 −0
Original line number Diff line number Diff line 
close all

clear all



fileName = './OpacityProfiles.h5';



indexRadius = 180;

indexEnergy = 20;



E = h5read( fileName, '/Profile Data/Energy' );

R = h5read( fileName, '/Profile Data/Radius' );

D = h5read( fileName, '/Profile Data/Density' );

T = h5read( fileName, '/Profile Data/Temperature' );

Y = h5read( fileName, '/Profile Data/Electron Fraction' );



Abs_NuE    =  h5read( fileName, '/Opacities/Absorption Opacity (Electron Neutrino)' );

Abs_NuEbar =  h5read( fileName, '/Opacities/Absorption Opacity (Electron Antineutrino)' );

Sct_NuE    =  h5read( fileName, '/Opacities/Scattering Opacity (Electron Neutrino)' );

Sct_NuEbar =  h5read( fileName, '/Opacities/Scattering Opacity (Electron Antineutrino)' );



fig = figure(1);

subplot( 2, 2, 1 )



loglog( R, D, '-k', 'linewidth', 2 ); hold on

loglog( R, T, '-b', 'linewidth', 2 )

plot( [R(indexRadius) R(indexRadius)], [D(indexRadius) D(indexRadius)], 'ok' )

plot( [R(indexRadius) R(indexRadius)], [T(indexRadius) T(indexRadius)], 'ob' )

axis( [1.d4 1.d8 1.d8 1.d15] )

set( gca, 'FontSize', 10 );

xlabel( 'Radius [ cm ]', 'fontsize', 10 );

ylabel( 'Mass Density [ g cm^{-3} ], Temperature [ K ]', 'fontsize', 8 );

legend_h...

  = legend...

      ( 'Mass Density',...

        'Temperature',...

        'location', 'best' );

set( legend_h, 'fontsize', 10 );



subplot( 2, 2, 3 )



semilogx( R, Y, '-k', 'linewidth', 2 ); hold on

plot( [R(indexRadius) R(indexRadius)], [Y(indexRadius) Y(indexRadius)], 'ok' )

axis( [1.d4 1.d8 0.1 0.6] )

set( gca, 'FontSize', 10 );

xlabel( 'Radius [ cm ]', 'fontsize', 10 );

ylabel( 'Electron Fraction', 'fontsize', 8 );



subplot( 2, 2, 2 )

loglog( E, Abs_NuE   (:,indexRadius),  '-k', 'linewidth', 2 ); hold on

loglog( E, Abs_NuEbar(:,indexRadius), '--m', 'linewidth', 2 )

loglog( E, Sct_NuE   (:,indexRadius),  '-b', 'linewidth', 2 )

loglog( E, Sct_NuEbar(:,indexRadius), '--c', 'linewidth', 2 )

axis( [1.d0 3.d2 1.d-6 1.d-0 ] )

xlabel( 'Energy [ MeV ]', 'fontsize', 10 );

ylabel( '1/\lambda [ cm^{-1} ]', 'fontsize', 8 );

title(['\rho = ',num2str(D(indexRadius),'%.1e'),' [ g cm^{-3} ]',...

       ', T = ', num2str(T(indexRadius),'%.1e'), ' [ K ]',...

       ', Y_{e} = ',num2str(Y(indexRadius),'%.1e')],'fontsize',7);

legend_h...

  = legend...

      ( {'Absorption $\nu_e$',...

         'Absorption $\bar{\nu}_e$',...

         'Scattering  $\nu_e$',...

         'Scattering  $\bar{\nu}_e$'},...

         'Interpreter', 'latex',...

         'location', 'northwest' );

set( legend_h, 'fontsize', 10 );



subplot( 2, 2, 4 )

loglog( R, Abs_NuE   (indexEnergy,:),  '-k', 'linewidth', 2 ); hold on

loglog( R, Abs_NuEbar(indexEnergy,:), '--m', 'linewidth', 2 )

loglog( R, Sct_NuE   (indexEnergy,:),  '-b', 'linewidth', 2 )

loglog( R, Sct_NuEbar(indexEnergy,:), '--c', 'linewidth', 2 )

axis( [1.d4 1.d8 1.d-13 1.d-2] )

xlabel( 'Radius [ cm ]', 'fontsize', 10 );

ylabel( '1/\lambda [ cm^{-1} ]', 'fontsize', 8 );

title(['Energy = ',num2str(E(indexEnergy),'%.1e'),' MeV']);

legend_h...

  = legend...

      ( {'Absorption $\nu_e$',...

         'Absorption $\bar{\nu}_e$',...

         'Scattering  $\nu_e$',...

         'Scattering  $\bar{\nu}_e$'},...

         'Interpreter', 'latex',...

         'location', 'southwest' );

set( legend_h, 'fontsize', 10 );



print( fig, '-dpng', './InterpolatedOpacities.png' )
 
 No newline at end of file

microphysics/InterpolationTest/wlInterpolationTest.f90

0 → 100644
+197 −0
Original line number Diff line number Diff line 
PROGRAM wlInterpolationTest



  USE wlKindModule, ONLY: dp

  USE wlGridModule, ONLY: &

    GridType, &

    AllocateGrid, &

    DescribeGrid, &

    MakeLogGrid

  USE wlOpacityTableModule, ONLY: &

    OpacityTableType, &

    DeallocateOpacityTable

  USE wlIOModuleHDF, ONLY: &

    InitializeHDF, &

    FinalizeHDF, &

    OpenFileHDF, &

    CloseFileHDF, &

    OpenGroupHDF, &

    CloseGroupHDF, &

    WriteHDF

  USE wlOpacityTableIOModuleHDF, ONLY: &

    ReadOpacityTableHDF

  USE wlInterpolationModule, ONLY: &

    LogInterpolateSingleVariable_1D3D



  USE HDF5



  IMPLICIT NONE



  INTEGER,     PARAMETER :: nx = 213

  INTEGER,     PARAMETER :: ne = 40

  INTEGER,     PARAMETER :: iD = 1

  INTEGER,     PARAMETER :: iT = 2

  INTEGER,     PARAMETER :: iY = 3

  INTEGER,     PARAMETER :: nSpecies = 2

  INTEGER                :: iS

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

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

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

  TYPE(GridType)         :: xEnergyGrid

  TYPE(OpacityTableType) :: OPACITIES



  ! --- Create Matter Profiles for Opacity Interpolation ---



  CALL ReadMatterProfiles( nx, xR, xD, xT, xY )



  ! --- Create Energy Grid for Opacity Interpolation ---



  CALL AllocateGrid( xEnergyGrid, ne )



  xEnergyGrid % Unit      = 'MeV'

  xEnergyGrid % Name      = 'Energy Grid'

  xEnergyGrid % MinValue  = 1.0d0

  xEnergyGrid % MaxValue  = 3.0d2

  xEnergyGrid % LogInterp = 1



  CALL MakeLogGrid &

         ( xEnergyGrid % MinValue, &

           xEnergyGrid % MaxValue, &

           xEnergyGrid % nPoints,  &

           xEnergyGrid % Values )



  CALL DescribeGrid( xEnergyGrid )



  ! --- Read Opacity Table ---



  CALL InitializeHDF( )



  CALL ReadOpacityTableHDF &

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

           Verbose_Option = .TRUE. )



  CALL FinalizeHDF( )



  ! --- Interpolate Opacities ---



  DO iS = 1, nSpecies



    ASSOCIATE &

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

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

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

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

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

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

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

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

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

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



    CALL LogInterpolateSingleVariable_1D3D &

           ( xEnergyGrid % Values, xD, xT, xY, OpacE, EosD, EosT, EosY, &

             [ 1, 1, 1, 0 ], OpacAbsOS, OpacAbs, xOpacAbs(:,:,iS) )



    CALL LogInterpolateSingleVariable_1D3D &

           ( xEnergyGrid % Values, xD, xT, xY, OpacE, EosD, EosT, EosY, &

             [ 1, 1, 1, 0 ], OpacSctOS, OpacSct, xOpacSct(:,:,iS) )



    END ASSOCIATE ! OpacAbs, etc.



  END DO



  CALL DeallocateOpacityTable( OPACITIES, Verbose_Option = .TRUE. )



  ASSOCIATE( xE => xEnergyGrid % Values )



  CALL WriteOpacityProfilesHDF &

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



  END ASSOCIATE



CONTAINS





  SUBROUTINE ReadMatterProfiles( nx, xR, xD, xT, xY )



    INTEGER,  INTENT(in)  :: nx

    REAL(dp), INTENT(out) :: xR(:), xD(:), xT(:), xY(:)



    CHARACTER(04) :: Format1 = "(I5)"

    CHARACTER(09) :: Format2 = "(4ES12.3)"

    INTEGER       :: ix

    REAL(dp)      :: buffer(4*nx)



    OPEN( 1, FILE = "MatterProfiles.d", FORM = "formatted", ACTION = 'read' )



    READ( 1, Format2 ) buffer



    DO ix = nx, 1, -1



      xR(ix) = buffer((ix-1)*4+1)

      xD(ix) = buffer((ix-1)*4+2)

      xT(ix) = buffer((ix-1)*4+3)

      xY(ix) = buffer((ix-1)*4+4)



    END DO



    CLOSE( 1, STATUS = 'keep' )



  END SUBROUTINE ReadMatterProfiles





  SUBROUTINE WriteOpacityProfilesHDF &

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



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

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

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

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



    INTEGER(HID_T)   :: file_id, group_id

    INTEGER(HSIZE_T) :: dsize1d(1), dsize2d(2)



    CALL InitializeHDF( )

    CALL OpenFileHDF( 'OpacityProfiles.h5', .TRUE., file_id )



    ! --- Write E, R, D, T, Ye ---



    CALL OpenGroupHDF( 'Profile Data', .true., file_id, group_id )



    dsize1d(1) = nE

    CALL WriteHDF( "Energy", xE, group_id, dsize1d ) 



    dsize1d(1) = nX

    CALL WriteHDF( "Radius",            xR, group_id, dsize1d )

    CALL WriteHDF( "Density",           xD, group_id, dsize1d )

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

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



    CALL CloseGroupHDF( group_id )



    ! --- Write Opacities ---



    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 CloseGroupHDF( group_id )



    CALL CloseFileHDF( file_id )

    CALL FinalizeHDF( )



  END SUBROUTINE WriteOpacityProfilesHDF





END PROGRAM wlInterpolationTest

microphysics/Tables/wl-Op-SFHo-15-25-50-AbIs.h5

0 → 100755LFS
+134 B

File added.

No diff preview for this file type.

View file