diff --git a/src/clib/Make.config.objects b/src/clib/Make.config.objects index 8229772a..26dc4649 100644 --- a/src/clib/Make.config.objects +++ b/src/clib/Make.config.objects @@ -20,11 +20,9 @@ OBJS_CONFIG_LIB = \ calculate_gamma.lo \ calculate_pressure.lo \ calculate_temperature.lo \ - calc_temp1d_cloudy_g.lo \ calc_temp_cloudy_g.lo \ calc_tdust_1d_g.lo \ calc_tdust_3d_g.lo \ - cool1d_cloudy_g.lo \ cool1d_cloudy_old_tables_g.lo \ cool1d_multi_g.lo \ cool_multi_time_g.lo \ @@ -34,11 +32,14 @@ OBJS_CONFIG_LIB = \ initialize_chemistry_data.lo \ initialize_cloudy_data.lo \ initialize_UVbackground_data.lo \ - interpolators_g.lo \ set_default_chemistry_parameters.lo \ solve_chemistry.lo \ solve_rate_cool_g.lo \ update_UVbackground_rates.lo \ rate_functions.lo \ initialize_rates.lo \ - utils.lo \ No newline at end of file + utils.lo \ + interop/shims_g.lo \ + interop/calc_temp1d_cloudy_g.lo \ + interop/cool1d_cloudy_g.lo \ + interop/interpolators_g.lo \ No newline at end of file diff --git a/src/clib/Makefile b/src/clib/Makefile index 74b46054..1f97e6a8 100644 --- a/src/clib/Makefile +++ b/src/clib/Makefile @@ -138,7 +138,7 @@ verbose: VERBOSE = 1 @rm -f $@ @(if [ $(VERBOSE) -eq 0 ]; then \ echo "Compiling $<" ; \ - $(LIBTOOL) --mode=compile --tag=FC $(FC) -c $(FFLAGS) $(DEFINES) $*.F >& $(OUTPUT) ; \ + $(LIBTOOL) --mode=compile --tag=FC $(FC) -c $(FFLAGS) $(DEFINES) $*.F -o $@ >& $(OUTPUT) ; \ if [ ! -e $@ ]; then \ echo; \ echo "Compiling $< failed!"; \ @@ -147,7 +147,7 @@ verbose: VERBOSE = 1 exit 1; \ fi ; \ else \ - $(LIBTOOL) --mode=compile --tag=FC $(FC) -c $(FFLAGS) $(DEFINES) $*.F; \ + $(LIBTOOL) --mode=compile --tag=FC $(FC) -c $(FFLAGS) $(DEFINES) $*.F -o $@; \ if [ ! -e $@ ]; then \ exit 1; \ fi ; \ @@ -157,7 +157,7 @@ verbose: VERBOSE = 1 @rm -f $@ @(if [ $(VERBOSE) -eq 0 ]; then \ echo "Compiling $<" ; \ - $(LIBTOOL) --mode=compile --tag=CXX $(CXX) -c $(DEFINES) $(CXXFLAGS) $(INCLUDES) $*.C \ + $(LIBTOOL) --mode=compile --tag=CXX $(CXX) -c $(DEFINES) $(CXXFLAGS) $(INCLUDES) $*.C -o $@\ >& $(OUTPUT) ; \ if [ ! -e $@ ]; then \ echo; \ @@ -167,7 +167,7 @@ verbose: VERBOSE = 1 exit 1; \ fi ; \ else \ - $(LIBTOOL) --mode=compile --tag=CXX $(CXX) -c $(DEFINES) $(CXXFLAGS) $(INCLUDES) $*.C; \ + $(LIBTOOL) --mode=compile --tag=CXX $(CXX) -c $(DEFINES) $(CXXFLAGS) $(INCLUDES) $*.C -o $@; \ if [ ! -e $@ ]; then \ exit 1; \ fi ; \ @@ -177,7 +177,7 @@ verbose: VERBOSE = 1 @rm -f $@ @(if [ $(VERBOSE) -eq 0 ]; then \ echo "Compiling $<" ; \ - $(LIBTOOL) --mode=compile --tag=CC $(CC) -c $(DEFINES) $(CFLAGS) $(INCLUDES) $*.c \ + $(LIBTOOL) --mode=compile --tag=CC $(CC) -c $(DEFINES) $(CFLAGS) $(INCLUDES) $*.c -o $@ \ >& $(OUTPUT) ; \ if [ ! -e $@ ]; then \ echo; \ @@ -187,7 +187,7 @@ verbose: VERBOSE = 1 exit 1; \ fi ; \ else \ - $(LIBTOOL) --mode=compile --tag=CC $(CC) -c $(DEFINES) $(CFLAGS) $(INCLUDES) $*.c; \ + $(LIBTOOL) --mode=compile --tag=CC $(CC) -c $(DEFINES) $(CFLAGS) $(INCLUDES) $*.c -o $@; \ if [ ! -e $@ ]; then \ exit 1; \ fi ; \ @@ -307,7 +307,10 @@ install: #----------------------------------------------------------------------- clean: - -@rm -f *.la .libs/* *.o *.lo DEPEND.bak *~ $(OUTPUT) grackle_float.h *.exe auto_show*.c DEPEND out.make.DEPEND + -@rm -f *.la .libs/* *.o *.lo \ + interop/.libs/* interop/*.o interop/*.lo \ + DEPEND.bak *~ $(OUTPUT) grackle_float.h *.exe auto_show*.c \ + DEPEND out.make.DEPEND -@touch DEPEND #----------------------------------------------------------------------- diff --git a/src/clib/calc_temp_cloudy_g.F b/src/clib/calc_temp_cloudy_g.F index ac80f6f2..8e935452 100644 --- a/src/clib/calc_temp_cloudy_g.F +++ b/src/clib/calc_temp_cloudy_g.F @@ -59,8 +59,10 @@ subroutine calc_temp_cloudy_g(d, e, metal, temperature, ! !----------------------------------------------------------------------- + USE ISO_C_BINDING implicit NONE #include "grackle_fortran_types.def" +#include "interop/interop_funcs.def" #ifdef _OPENMP #include "omp_lib.h" #endif @@ -102,7 +104,7 @@ subroutine calc_temp_cloudy_g(d, e, metal, temperature, ! Iteration mask - logical itmask(in) + integer*4 itmask(in) ! !\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\///////////////////////////////// !======================================================================= @@ -145,7 +147,7 @@ subroutine calc_temp_cloudy_g(d, e, metal, temperature, ! Initialize iteration mask to true for all cells. do i = is+1, ie+1 - itmask(i) = .true. + itmask(i) = 1 rhoH(i) = fh * d(i,j,k) enddo diff --git a/src/clib/cool1d_cloudy_old_tables_g.F b/src/clib/cool1d_cloudy_old_tables_g.F index 2e95caea..26261889 100644 --- a/src/clib/cool1d_cloudy_old_tables_g.F +++ b/src/clib/cool1d_cloudy_old_tables_g.F @@ -56,8 +56,10 @@ subroutine cool1D_cloudy_old_tables_g(d, de, rhoH, metallicity, ! !----------------------------------------------------------------------- + USE ISO_C_BINDING implicit NONE #include "grackle_fortran_types.def" +#include "interop/interop_funcs.def" ! General Arguments diff --git a/src/clib/cool1d_multi_g.F b/src/clib/cool1d_multi_g.F index d8762b3f..ecf9d564 100644 --- a/src/clib/cool1d_multi_g.F +++ b/src/clib/cool1d_multi_g.F @@ -237,7 +237,7 @@ subroutine cool1d_multi_g( enddo endif - call calc_temp1d_cloudy_g(d, metal, e, rhoH, + call calc_temp1d_cloudy_shim_g(d, metal, e, rhoH, & in, jn, kn, is, ie, j, k, & tgas, mmw, dom, zr, & temstart, temend, @@ -941,7 +941,7 @@ subroutine cool1d_multi_g( iZscale = 0 mycmbTfloor = 0 - call cool1d_cloudy_g(d, rhoH, metallicity, + call cool1d_cloudy_shim_g(d, rhoH, metallicity, & in, jn, kn, is, ie, j, k, & logtem, edot, comp2, dom, zr, & mycmbTfloor, iClHeat, iZscale, @@ -1094,7 +1094,7 @@ subroutine cool1d_multi_g( if (clnew .eq. 1) then iZscale = 1 - call cool1d_cloudy_g(d, rhoH, metallicity, + call cool1d_cloudy_shim_g(d, rhoH, metallicity, & in, jn, kn, is, ie, j, k, & logtem, edot, comp2, dom, zr, & icmbTfloor, iClHeat, iZscale, diff --git a/src/clib/grackle_macros.h b/src/clib/grackle_macros.h index 4bce918f..44a093e7 100644 --- a/src/clib/grackle_macros.h +++ b/src/clib/grackle_macros.h @@ -19,6 +19,22 @@ / ************************************************************************/ +// define the GR_RESTRICT macro be the restrict keyword introduced in C99 +// -> restrict isn't technically part of the C++ standard, but it is commonly +// provided by c++ compilers (so we try to make use of those alternatives) +#ifdef CONFIG_NO_RESTRICT + #define GR_RESTRICT /* ... */ +#elif !defined(__cplusplus) /* simple case (we are compiling C) */ + #define GR_RESTRICT restrict +#elif __GNUC__ + // C++ compilers other than g++ define this macro. To my knowledge, all of + // them (e.g. clang++, the new & old intel c++ compilers) define the same + // the restrict-extension in the same way + #define GR_RESTRICT __restrict +#else + #define GR_RESTRICT /* ... */ +#endif /* GR_RESTRICT */ + #define GRACKLE_FREE(p) \ { \ if (p != NULL) { \ diff --git a/src/clib/interop/README.md b/src/clib/interop/README.md new file mode 100644 index 00000000..caa5b119 --- /dev/null +++ b/src/clib/interop/README.md @@ -0,0 +1,67 @@ +## Overview + +Files in this directory are used to define "interoperable functions." These are C functions that may need to be called from Fortran subroutines (usually these C functions define functionality that was previously transcribed from Fortran subroutines). + +Currently, the C function-prototypes for these "interoperable functions" are defined in [interop_funcs.h](./interop_funcs.h) while the corresponding Fortran interfaces are declared within [interop_funcs.def](./interop_funcs.def). + +## Calling These Functions from Fortran + +To be able to call the C functions from Fortran, we make use of Fortran 2003's ```ISO_C_BINDING`` module (which is backported on most Fortran compilers). + +When you're writing a Fortran subroutine (or function) that needs to call one or more of these "interoperable functions", it's critical that your function/routine properly include the interface. To accomplish this, you **MUST** make sure your subroutine includes lines with the following directives: + +1. It must contain a line that reads ``USE ISO_C_BINDING`` (this usually comes just before the line that reads ``implicit NONE``) + +2. After the line that reads ``implicit NONE``, you should insert your include directives. There are 2 relevant include directives (order matters!): + + 1. First, you should have an include-directive for ``grackle_fortran_types.def`` file[^1] + + 2. Next, you should have an include-directive for the ``interop_funcs.def`` file + +> [!CAUTION] +> If you try to call an interoperable function without following the above instructions, the problems will arise. Specifically, the Fortran compiler may not make proper guesses about name-mangling **AND** it has no way of knowing how to pass certain arguments into the function (whether it should be passed by value or by reference...) + +For the sake of concreteness, the following codeblock illustrates the general structure that a Fortran subroutine, defined in `src/clib/`, would need to have in order for it to call one of these interoperable functions: + +```fortran + subroutine my_subroutine( arg1, arg2, arg3 ) + +! < docstring sumarizing purpose and describing args ... > + + USE ISO_C_BINDING + implicit NONE +#include "grackle_fortran_types.def" +#include "interop/interop_funcs.def" + +! < declare arg types and then declare local variables (with their types) ... > + +! < do actual work work... > + + return + end +``` + + +## Fortran Shims + +There are some issues with the datatype of the iteration mask and interoperability. At it's core, this issue boils down to the fact that the Fortran doesn't guarantee any compatibility between its default ``LOGICAL`` type and any C type. + +As a short-term solution, we introduce shims to convert the logical mask into a temporary integer mask. This is somewhat undesirable given that memory is allocated each time we call a shim. + +There are essentially 2 long-term solutions: + +1. the simplest solution is to replace all occurences of the ``LOGICAL`` type in the Fortran code with ``logical(C_BOOL)``. + + - In more detail, ``C_BOOL`` is a "KIND type parameters" provided by the ``ISO_C_BINDING`` module. Using ``logical(C_BOOL)``, ensures that the datatype is compatible with the ``_Bool`` datatype introduced in C99.[^2] + + - However, this may produce problems if we want to compile the code with a C++ compiler (this is essentially required to support GPUs). In more detail, C++'s ``bool`` type is [NOT guaranteed to be compatible](https://stackoverflow.com/q/40020423) with C's ``_Bool`` type (in practice, this probably isn't much of an issue). + +2. The alternative, is to define some consistent integer datatype that is used in to hold the mask values in the Fortran and C layers. This approach might look like the following: + + - inserting something like ``#define gr_mask_int int32_t`` into ``grackle_macros.h`` and ``#define GR_MASK_INT integer*4`` into ``grackle_fortran_types.def`` + + - replacing all usage of ``LOGICAL`` with ``GR_MASK_INT`` (all boolean operations would need to be replaced with comparisons against 0) + +[^1]: (this may not be strictly necessary right now, but it will probably be necessary in the future). + +[^2]: This is footnote is provided for the sake of clarity since people often get confused about C's boolean type. The 1999 C standard introduced ``_Bool`` as a standard type for holding booleans (i.e. any compliant C compiler always defines this type, even when no headers are included). The type is only capable of holding ``1`` or ``0``. The 1999 standandard also introduced the [](https://en.cppreference.com/w/c/types/boolean) component to the "standard library" to define useful macros. These macros include ``bool``, which expands to ``_Bool``, as well as the ``true`` and ``false`` macros, which expand to the integer constants ``1`` and ``0``. (As an aside, things are a little different when a compiler the 2023 standard, but that's not relevant to us) \ No newline at end of file diff --git a/src/clib/interop/calc_temp1d_cloudy_g.c b/src/clib/interop/calc_temp1d_cloudy_g.c new file mode 100644 index 00000000..57bcbcf5 --- /dev/null +++ b/src/clib/interop/calc_temp1d_cloudy_g.c @@ -0,0 +1,128 @@ +/*********************************************************************** +/ +/ Calculate temperature field of a 1D slice from a cloudy table +/ +/ +/ Copyright (c) 2013, Enzo/Grackle Development Team. +/ +/ Distributed under the terms of the Enzo Public Licence. +/ +/ The full license is in the file LICENSE, distributed with this +/ software. +************************************************************************/ + +#include +#include // int32_t +#include +#include // log, log10 +#include "../grackle_chemistry_data.h" +#include "../phys_constants.h" +#include "./interop_funcs.h" + +void calc_temp1d_cloudy_g( + const gr_float* d, const gr_float* metal, // 3D arrays + const gr_float* e, // 3D array + const double* rhoH, // 1D array + int in, int jn, int kn, int is, int ie, int j, int k, + double * GR_RESTRICT tgas, double * GR_RESTRICT mmw, // 1D array + double dom, double zr, double temstart, double temend, double gamma, + double utem, int imetal, + long long clGridRank, + const long long* clGridDim, + const double* clPar1, + const double* clPar2, + const double* clPar3, + long long clDataSize, + const double* clMMW, + const int32_t* itmask) +{ + const double mu_metal = 16.0; + const int ti_max = 20; + + const double inv_log10 = 1.0 / log(10.0); + + // Calculate parameter value slopes + const double dclPar[3] = { + ((clPar1[clGridDim[0]-1] - clPar1[0]) / (double)(clGridDim[0] - 1)), + + (clGridRank > 1) + ? ((clPar2[clGridDim[1]-1] - clPar2[0]) / (double)(clGridDim[1] - 1)) : 0, + + (clGridRank > 2) + ? ((clPar3[clGridDim[2]-1] - clPar3[0]) / (double)(clGridDim[2] - 1)) : 0 + }; + + // Calculate index for redshift dimension - intentionally kept 1-indexed + const long long zindex = find_zindex(zr, clGridRank, clGridDim, clPar2); + const gr_int64 end_int = ((clGridRank > 2) && (zindex == clGridDim[1])); + + for (int i = is + 1; i <= (ie + 1); i++) { + if ( !itmask[i-1] ) { continue; } + + const double log_n_h = log10(rhoH[i-1] * dom); // Calculate proper log(n_H) + + const long long ind_3D = (i-1) + in *( (j-1) + jn * (k-1)); + + double munew = 1.0; + double muold; + for (int ti = 1; ti <= ti_max; ti++) { + muold = munew; + + tgas[i-1] = max((gamma - 1.0) * e[ind_3D] * munew * utem, temstart); + // the original version doesn't use log10*(tgas[i-1]) either + const double log10tem = log(tgas[i-1]) * inv_log10; + + // Call interpolation functions to get mmw + if (clGridRank == 1) { // Interpolate over temperature. + interpolate_1d_g(log10tem, clGridDim, + clPar1, dclPar[0], + clDataSize, clMMW, &munew); + } else if ( clGridRank == 2) { // Interpolate over density & temperature. + interpolate_2d_g(log_n_h, log10tem, + clGridDim, + clPar1, dclPar[0], + clPar2, dclPar[1], + clDataSize, clMMW, &munew); + } else if (clGridRank == 3) { // Interpolate over density, redshift, + // & temperature. + interpolate_3dz_g(log_n_h, zr, log10tem, + clGridDim, + clPar1, dclPar[0], + clPar2, zindex, + clPar3, dclPar[2], + clDataSize, clMMW, + end_int, &munew); + } else { + // no need to be within an openmp critical section + printf("Maximum mmw data grid rank is 3!\n"); + } + + munew = 0.5 * (munew + muold); + tgas[i-1] = tgas[i-1] * munew / muold; + + if (fabs((munew/muold) - 1.0) < 1.e-2) { + muold = munew; + + if (imetal == 1){ + munew = d[ind_3D] / ((d[ind_3D] - metal[ind_3D])/ munew + + metal[ind_3D] / mu_metal); + tgas[i-1] = tgas[i-1] * munew / muold; + } + + mmw[i-1] = munew; + goto completed_T_calc; // TODO: refactor to remove this goto statement! + } + } + + mmw[i-1] = munew; + + // no need to put this in an omp critical statement + fprintf(stderr, + "Mean molecular weight not converged! %#.16g, %#.16g, %#.16g\n", + munew, muold, fabs((munew/muold) - 1.0)); + + completed_T_calc: + continue; + + } +} diff --git a/src/clib/interop/cool1d_cloudy_g.c b/src/clib/interop/cool1d_cloudy_g.c new file mode 100644 index 00000000..eb663fc6 --- /dev/null +++ b/src/clib/interop/cool1d_cloudy_g.c @@ -0,0 +1,162 @@ +/*********************************************************************** +/ +/ Calculate temperature field of a 1D slice from a cloudy table +/ +/ +/ Copyright (c) 2013, Enzo/Grackle Development Team. +/ +/ Distributed under the terms of the Enzo Public Licence. +/ +/ The full license is in the file LICENSE, distributed with this +/ software. +************************************************************************/ + +#include +#include // int32_t +#include +#include // log, log10 +#include "../grackle_chemistry_data.h" +#include "../phys_constants.h" +#include "./interop_funcs.h" + +void cool1d_cloudy_g( + const gr_float* d, // 3D arrays + const double* rhoH, const gr_float* metallicity, // 1D array + int in, int jn, int kn, int is, int ie, int j, int k, + const double* logtem, double * GR_RESTRICT edot, // 1D array + double comp2, double dom, double zr, + int icmbTfloor, int iClHeat, int iZscale, + long long clGridRank, + const long long* clGridDim, + const double* clPar1, const double* clPar2, const double* clPar3, + long long clDataSize, const double* clCooling, const double* clHeating, + const int32_t* itmask) +{ + const double inv_log10 = 1.0 / log(10.0); + const double log10_tCMB = log10(comp2); + + // Calculate parameter value slopes + const double dclPar[3] = { + ((clPar1[clGridDim[0]-1] - clPar1[0]) / (double)(clGridDim[0] - 1)), + + (clGridRank > 1) + ? ((clPar2[clGridDim[1]-1] - clPar2[0]) / (double)(clGridDim[1] - 1)) : 0, + + (clGridRank > 2) + ? ((clPar3[clGridDim[2]-1] - clPar3[0]) / (double)(clGridDim[2] - 1)) : 0 + }; + + // Calculate index for redshift dimension - intentionally kept 1-indexed + const long long zindex = find_zindex(zr, clGridRank, clGridDim, clPar2); + const gr_int64 end_int = ((clGridRank > 2) && (zindex == clGridDim[1])); + const int get_heat = ((clGridRank > 2) && (zindex == clGridDim[1])) + ? 0 : iClHeat; + + // do work + for (int i = is + 1; i <= (ie + 1); i++) { + if ( !itmask[i-1] ) { continue; } + + const double log10tem = logtem[i-1] * inv_log10; + const double log_n_h = log10(rhoH[i-1] * dom); // Calculate proper log(n_H) + + double edot_met = 0.0; + double log_cool = 0.0; + double log_cool_cmb = 0.0; + double log_heat = 0.0; + + // Call interpolation functions to get heating and cooling + if (clGridRank == 1) { // Interpolate over temperature. + interpolate_1d_g(log10tem, clGridDim, + clPar1, dclPar[0], + clDataSize, clCooling, &log_cool); + edot_met = -1.0 * pow(10.0,log_cool); + + // Ignore CMB term if T >> T_CMB + if ((icmbTfloor == 1) && ((log10tem - log10_tCMB) < 2.0)) { + interpolate_1d_g(log10_tCMB, clGridDim, + clPar1, dclPar[0], + clDataSize, clCooling, &log_cool_cmb); + edot_met += pow(10.0,log_cool_cmb); + } + + if (get_heat == 1){ + interpolate_1d_g(log10tem, clGridDim, + clPar1, dclPar[0], + clDataSize, clHeating, &log_heat); + edot_met += pow(10.0,log_heat); + } + + } else if ( clGridRank == 2) { // Interpolate over density & temperature. + interpolate_2d_g(log_n_h, log10tem, + clGridDim, + clPar1, dclPar[0], + clPar2, dclPar[1], + clDataSize, clCooling, &log_cool); + edot_met = -1.0 * pow(10.0,log_cool); + + // Ignore CMB term if T >> T_CMB + if ((icmbTfloor == 1) && ((log10tem - log10_tCMB) < 2.0)) { + interpolate_2d_g(log_n_h, log10_tCMB, + clGridDim, + clPar1, dclPar[0], + clPar2, dclPar[1], + clDataSize, clCooling, &log_cool_cmb); + edot_met += pow(10.0,log_cool_cmb); + } + + if (get_heat == 1){ + interpolate_2d_g(log_n_h, log10tem, + clGridDim, + clPar1, dclPar[0], + clPar2, dclPar[1], + clDataSize, clHeating, &log_heat); + edot_met += pow(10.0,log_heat); + } + + } else if ( clGridRank == 3) { // Interpolate over density, redshift, + // & temperature. + interpolate_3dz_g(log_n_h, zr, log10tem, + clGridDim, + clPar1, dclPar[0], + clPar2, zindex, + clPar3, dclPar[2], + clDataSize, clCooling, + end_int, &log_cool); + edot_met = -1.0 * pow(10.0,log_cool); + + // Ignore CMB term if T >> T_CMB + if ((icmbTfloor == 1) && ((log10tem - log10_tCMB) < 2.0)) { + interpolate_3dz_g(log_n_h, zr, log10_tCMB, + clGridDim, + clPar1, dclPar[0], + clPar2, zindex, + clPar3, dclPar[2], + clDataSize, clCooling, + end_int, &log_cool_cmb); + edot_met += pow(10.0,log_cool_cmb); + } + + if (get_heat == 1){ + interpolate_3dz_g(log_n_h, zr, log10tem, + clGridDim, + clPar1, dclPar[0], + clPar2, zindex, + clPar3, dclPar[2], + clDataSize, clHeating, + end_int, &log_heat); + edot_met += pow(10.0,log_heat); + } + + } else { + // no need to be within an openmp critical section + printf("Maximum cooling data grid rank is 3!\n"); + } + + // scale cooling by metallicity + if (iZscale == 1) { + edot_met = edot_met * metallicity[i-1]; + } + + edot[i-1] += (edot_met * rhoH[i-1] * rhoH[i-1]); + } +} diff --git a/src/clib/interop/interop_funcs.def b/src/clib/interop/interop_funcs.def new file mode 100644 index 00000000..03e38628 --- /dev/null +++ b/src/clib/interop/interop_funcs.def @@ -0,0 +1,250 @@ +!======================================================================= +! +! +! fortran interfaces of internal interoperable functions +! +! +! Copyright (c) 2013, Enzo/Grackle Development Team. +! +! Distributed under the terms of the Enzo Public Licence. +! +! The full license is in the file LICENSE, distributed with this +! software. +!======================================================================= + +c Here we define the Fortran 90 interfaces to the C routines using the 2003 +c Fortran C interoperability features. Note that some care must be taken +c defining matching types. + +c before including this file, be sure that: +c 1. you have invoked ``USE ISO_C_BINDING`` +c 2. grackle_fortran_types.def has already been included + + interface + subroutine calc_temp1d_cloudy_g(d, metal, e, rhoH, + & in, jn, kn, is, ie, j, k, + & tgas, mmw, dom, zr, + & temstart, temend, + & gamma, utem, imetal, + & clGridRank, clGridDim, + & clPar1, clPar2, clPar3, + & clDataSize, clMMW, + & itmask) bind(C) + IMPORT + R_PREC, INTENT(IN) :: d(*) + R_PREC, INTENT(IN) :: metal(*) + R_PREC, INTENT(IN) :: e(*) + REAL(C_DOUBLE), INTENT(IN) :: rhoH(*) + INTEGER(C_INT), VALUE, INTENT(IN) :: in + INTEGER(C_INT), VALUE, INTENT(IN) :: jn + INTEGER(C_INT), VALUE, INTENT(IN) :: kn + INTEGER(C_INT), VALUE, INTENT(IN) :: is + INTEGER(C_INT), VALUE, INTENT(IN) :: ie + INTEGER(C_INT), VALUE, INTENT(IN) :: j + INTEGER(C_INT), VALUE, INTENT(IN) :: k + REAL(C_DOUBLE), INTENT(OUT) :: tgas(*) + REAL(C_DOUBLE), INTENT(OUT) :: mmw(*) + REAL(C_DOUBLE), VALUE, INTENT(IN) :: dom + REAL(C_DOUBLE), VALUE, INTENT(IN) :: zr + REAL(C_DOUBLE), VALUE, INTENT(IN) :: temstart + REAL(C_DOUBLE), VALUE, INTENT(IN) :: temend + REAL(C_DOUBLE), VALUE, INTENT(IN) :: gamma + REAL(C_DOUBLE), VALUE, INTENT(IN) :: utem + INTEGER(C_INT), VALUE, INTENT(IN) :: imetal + INTEGER(C_LONG_LONG), VALUE, INTENT(IN) :: clGridRank + INTEGER(C_LONG_LONG), INTENT(IN) :: clGridDim(*) + REAL(C_DOUBLE), INTENT(IN) :: clPar1(*) + REAL(C_DOUBLE), INTENT(IN) :: clPar2(*) + REAL(C_DOUBLE), INTENT(IN) :: clPar3(*) + INTEGER(C_LONG_LONG), VALUE, INTENT(IN) :: clDataSize + REAL(C_DOUBLE), INTENT(IN) :: clMMW(*) + INTEGER(C_INT32_T), INTENT(IN) :: itmask(*) + end subroutine calc_temp1d_cloudy_g + end interface + + interface + subroutine cool1d_cloudy_g(d, rhoH, metallicity, + & in, jn, kn, is, ie, j, k, + & logtem, edot, comp2, dom, zr, + & icmbTfloor, iClHeat, iZscale, + & clGridRank, clGridDim, + & clPar1, clPar2, clPar3, + & clDataSize, clCooling, clHeating, + & itmask) bind(C) + IMPORT + R_PREC, INTENT(IN) :: d(*) + REAL(C_DOUBLE), INTENT(IN) :: rhoH(*) + REAL(C_DOUBLE), INTENT(IN) :: metallicity(*) + INTEGER(C_INT), VALUE, INTENT(IN) :: in + INTEGER(C_INT), VALUE, INTENT(IN) :: jn + INTEGER(C_INT), VALUE, INTENT(IN) :: kn + INTEGER(C_INT), VALUE, INTENT(IN) :: is + INTEGER(C_INT), VALUE, INTENT(IN) :: ie + INTEGER(C_INT), VALUE, INTENT(IN) :: j + INTEGER(C_INT), VALUE, INTENT(IN) :: k + REAL(C_DOUBLE), INTENT(IN) :: logtem(*) + REAL(C_DOUBLE), INTENT(OUT) :: edot(*) + REAL(C_DOUBLE), VALUE, INTENT(IN) :: comp2 + REAL(C_DOUBLE), VALUE, INTENT(IN) :: dom + REAL(C_DOUBLE), VALUE, INTENT(IN) :: zr + INTEGER(C_INT), VALUE, INTENT(IN) :: icmbTfloor + INTEGER(C_INT), VALUE, INTENT(IN) :: iClHeat + INTEGER(C_INT), VALUE, INTENT(IN) :: iZscale + INTEGER(C_LONG_LONG), VALUE, INTENT(IN) :: clGridRank + INTEGER(C_LONG_LONG), INTENT(IN) :: clGridDim(*) + REAL(C_DOUBLE), INTENT(IN) :: clPar1(*) + REAL(C_DOUBLE), INTENT(IN) :: clPar2(*) + REAL(C_DOUBLE), INTENT(IN) :: clPar3(*) + INTEGER(C_LONG_LONG), VALUE, INTENT(IN) :: clDataSize + REAL(C_DOUBLE), INTENT(IN) :: clCooling(*) + REAL(C_DOUBLE), INTENT(IN) :: clHeating(*) + INTEGER(C_INT32_T), INTENT(IN) :: itmask(*) + end subroutine cool1d_cloudy_g + end interface + + interface + subroutine interpolate_1D_g( + & input1, gridDim, gridPar1, dgridPar1, + & dataSize, dataField, value) bind(C) + IMPORT + REAL(C_DOUBLE), VALUE, INTENT(IN) :: input1 + INTEGER(C_LONG_LONG), INTENT(IN) :: gridDim(*) + REAL(C_DOUBLE), INTENT(IN) :: gridPar1(*) + REAL(C_DOUBLE), VALUE, INTENT(IN) :: dgridPar1 + INTEGER(C_LONG_LONG), VALUE, INTENT(IN) :: dataSize + REAL(C_DOUBLE), INTENT(IN) :: dataField(*) + REAL(C_DOUBLE), INTENT(OUT) :: value + end subroutine interpolate_1D_g + end interface + + interface + subroutine interpolate_2D_g( + & input1, input2, gridDim, + & gridPar1, dgridPar1, + & gridPar2, dgridPar2, + & dataSize, dataField, value) bind(C) + IMPORT + REAL(C_DOUBLE), VALUE, INTENT(IN) :: input1 + REAL(C_DOUBLE), VALUE, INTENT(IN) :: input2 + INTEGER(C_LONG_LONG), INTENT(IN) :: gridDim(*) + REAL(C_DOUBLE), INTENT(IN) :: gridPar1(*) + REAL(C_DOUBLE), VALUE, INTENT(IN) :: dgridPar1 + REAL(C_DOUBLE), INTENT(IN) :: gridPar2(*) + REAL(C_DOUBLE), VALUE, INTENT(IN) :: dgridPar2 + INTEGER(C_LONG_LONG), VALUE, INTENT(IN) :: dataSize + REAL(C_DOUBLE), INTENT(IN) :: dataField(*) + REAL(C_DOUBLE), INTENT(OUT) :: value + end subroutine interpolate_2D_g + end interface + + interface + subroutine interpolate_3D_g( + & input1, input2, input3, gridDim, + & gridPar1, dgridPar1, + & gridPar2, dgridPar2, + & gridPar3, dgridPar3, + & dataSize, dataField, value) bind(C) + IMPORT + REAL(C_DOUBLE), VALUE, INTENT(IN) :: input1 + REAL(C_DOUBLE), VALUE, INTENT(IN) :: input2 + REAL(C_DOUBLE), VALUE, INTENT(IN) :: input3 + INTEGER(C_LONG_LONG), INTENT(IN) :: gridDim(*) + REAL(C_DOUBLE), INTENT(IN) :: gridPar1(*) + REAL(C_DOUBLE), VALUE, INTENT(IN) :: dgridPar1 + REAL(C_DOUBLE), INTENT(IN) :: gridPar2(*) + REAL(C_DOUBLE), VALUE, INTENT(IN) :: dgridPar2 + REAL(C_DOUBLE), INTENT(IN) :: gridPar3(*) + REAL(C_DOUBLE), VALUE, INTENT(IN) :: dgridPar3 + INTEGER(C_LONG_LONG), VALUE, INTENT(IN) :: dataSize + REAL(C_DOUBLE), INTENT(IN) :: dataField(*) + REAL(C_DOUBLE), INTENT(OUT) :: value + end subroutine interpolate_3D_g + end interface + + interface + subroutine interpolate_3Dz_g( + & input1, input2, input3, gridDim, + & gridPar1, dgridPar1, + & gridPar2, index2, + & gridPar3, dgridPar3, + & dataSize, dataField, + & end_int, value) bind(C) + IMPORT + REAL(C_DOUBLE), VALUE, INTENT(IN) :: input1 + REAL(C_DOUBLE), VALUE, INTENT(IN) :: input2 + REAL(C_DOUBLE), VALUE, INTENT(IN) :: input3 + INTEGER(C_LONG_LONG), INTENT(IN) :: gridDim(*) + REAL(C_DOUBLE), INTENT(IN) :: gridPar1(*) + REAL(C_DOUBLE), VALUE, INTENT(IN) :: dgridPar1 + REAL(C_DOUBLE), INTENT(IN) :: gridPar2(*) + INTEGER(C_LONG_LONG), VALUE, INTENT(IN) :: index2 + REAL(C_DOUBLE), INTENT(IN) :: gridPar3(*) + REAL(C_DOUBLE), VALUE, INTENT(IN) :: dgridPar3 + INTEGER(C_LONG_LONG), VALUE, INTENT(IN) :: dataSize + REAL(C_DOUBLE), INTENT(IN) :: dataField(*) + INTEGER(C_LONG_LONG), VALUE, INTENT(IN) :: end_int + REAL(C_DOUBLE), INTENT(OUT) :: value + end subroutine interpolate_3Dz_g + end interface + + interface + subroutine interpolate_4D_g( + & input1, input2, input3, input4, + & gridDim, + & gridPar1, dgridPar1, + & gridPar2, dgridPar2, + & gridPar3, dgridPar3, + & gridPar4, dgridPar4, + & dataSize, dataField, value) bind(C) + IMPORT + REAL(C_DOUBLE), VALUE, INTENT(IN) :: input1 + REAL(C_DOUBLE), VALUE, INTENT(IN) :: input2 + REAL(C_DOUBLE), VALUE, INTENT(IN) :: input3 + REAL(C_DOUBLE), VALUE, INTENT(IN) :: input4 + INTEGER(C_LONG_LONG), INTENT(IN) :: gridDim(*) + REAL(C_DOUBLE), INTENT(IN) :: gridPar1(*) + REAL(C_DOUBLE), VALUE, INTENT(IN) :: dgridPar1 + REAL(C_DOUBLE), INTENT(IN) :: gridPar2(*) + REAL(C_DOUBLE), VALUE, INTENT(IN) :: dgridPar2 + REAL(C_DOUBLE), INTENT(IN) :: gridPar3(*) + REAL(C_DOUBLE), VALUE, INTENT(IN) :: dgridPar3 + REAL(C_DOUBLE), INTENT(IN) :: gridPar4(*) + REAL(C_DOUBLE), VALUE, INTENT(IN) :: dgridPar4 + INTEGER(C_LONG_LONG), VALUE, INTENT(IN) :: dataSize + REAL(C_DOUBLE), INTENT(IN) :: dataField(*) + REAL(C_DOUBLE), INTENT(OUT) :: value + end subroutine interpolate_4D_g + end interface + + interface + subroutine interpolate_5D_g( + & input1, input2, input3, input4, input5, + & gridDim, + & gridPar1, dgridPar1, + & gridPar2, dgridPar2, + & gridPar3, dgridPar3, + & gridPar4, dgridPar4, + & gridPar5, dgridPar5, + & dataSize, dataField, value) bind(C) + IMPORT + REAL(C_DOUBLE), VALUE, INTENT(IN) :: input1 + REAL(C_DOUBLE), VALUE, INTENT(IN) :: input2 + REAL(C_DOUBLE), VALUE, INTENT(IN) :: input3 + REAL(C_DOUBLE), VALUE, INTENT(IN) :: input4 + REAL(C_DOUBLE), VALUE, INTENT(IN) :: input5 + INTEGER(C_LONG_LONG), INTENT(IN) :: gridDim(*) + REAL(C_DOUBLE), INTENT(IN) :: gridPar1(*) + REAL(C_DOUBLE), VALUE, INTENT(IN) :: dgridPar1 + REAL(C_DOUBLE), INTENT(IN) :: gridPar2(*) + REAL(C_DOUBLE), VALUE, INTENT(IN) :: dgridPar2 + REAL(C_DOUBLE), INTENT(IN) :: gridPar3(*) + REAL(C_DOUBLE), VALUE, INTENT(IN) :: dgridPar3 + REAL(C_DOUBLE), INTENT(IN) :: gridPar4(*) + REAL(C_DOUBLE), VALUE, INTENT(IN) :: dgridPar4 + REAL(C_DOUBLE), INTENT(IN) :: gridPar5(*) + REAL(C_DOUBLE), VALUE, INTENT(IN) :: dgridPar5 + INTEGER(C_LONG_LONG), VALUE, INTENT(IN) :: dataSize + REAL(C_DOUBLE), INTENT(IN) :: dataField(*) + REAL(C_DOUBLE), INTENT(OUT) :: value + end subroutine interpolate_5D_g + end interface \ No newline at end of file diff --git a/src/clib/interop/interop_funcs.h b/src/clib/interop/interop_funcs.h new file mode 100644 index 00000000..04ef0ffb --- /dev/null +++ b/src/clib/interop/interop_funcs.h @@ -0,0 +1,218 @@ +/*********************************************************************** +/ +/ C interfaces of internal interop functions (and helper functions) +/ +/ +/ Copyright (c) 2013, Enzo/Grackle Development Team. +/ +/ Distributed under the terms of the Enzo Public Licence. +/ +/ The full license is in the file LICENSE, distributed with this +/ software. +************************************************************************/ +#ifndef __INTEROP_FUNCS_H +#define __INTEROP_FUNCS_H + +#include // int32_t + +#include "grackle_macros.h" +#include "grackle_types.h" + +typedef long long gr_int64; + +void interpolate_1d_g(double input1, + const gr_int64 * GR_RESTRICT gridDim, // 1 elements + const double * GR_RESTRICT gridPar1, double dgridPar1, + gr_int64 dataSize, const double * GR_RESTRICT dataField, + double * GR_RESTRICT value); + +void interpolate_2d_g(double input1, double input2, + const gr_int64 * GR_RESTRICT gridDim, // 2 elements + const double * GR_RESTRICT gridPar1, double dgridPar1, + const double * GR_RESTRICT gridPar2, double dgridPar2, + gr_int64 dataSize, const double* dataField, + double * GR_RESTRICT value); + +void interpolate_3dz_g(double input1, double input2, double input3, + const gr_int64 * GR_RESTRICT gridDim, // 3 elements + const double * GR_RESTRICT gridPar1, double dgridPar1, + const double * GR_RESTRICT gridPar2, gr_int64 index2, + const double * GR_RESTRICT gridPar3, double dgridPar3, + gr_int64 dataSize, const double * GR_RESTRICT dataField, + gr_int64 end_int, double * GR_RESTRICT value); + +void interpolate_3d_g(double input1, double input2, double input3, + const gr_int64 * GR_RESTRICT gridDim, // 3 elements + const double * GR_RESTRICT gridPar1, double dgridPar1, + const double * GR_RESTRICT gridPar2, double dgridPar2, + const double * GR_RESTRICT gridPar3, double dgridPar3, + gr_int64 dataSize, const double * GR_RESTRICT dataField, + double* value); + +void interpolate_4d_g(double input1, double input2, double input3, + double input4, + const gr_int64 * GR_RESTRICT gridDim, // 4 elements + const double * GR_RESTRICT gridPar1, double dgridPar1, + const double * GR_RESTRICT gridPar2, double dgridPar2, + const double * GR_RESTRICT gridPar3, double dgridPar3, + const double * GR_RESTRICT gridPar4, double dgridPar4, + gr_int64 dataSize, const double * GR_RESTRICT dataField, + double * GR_RESTRICT value); + +void interpolate_5d_g(double input1, double input2, double input3, + double input4, double input5, + const gr_int64 * GR_RESTRICT gridDim, // 5 elements + const double * GR_RESTRICT gridPar1, double dgridPar1, + const double * GR_RESTRICT gridPar2, double dgridPar2, + const double * GR_RESTRICT gridPar3, double dgridPar3, + const double * GR_RESTRICT gridPar4, double dgridPar4, + const double * GR_RESTRICT gridPar5, double dgridPar5, + gr_int64 dataSize, const double * GR_RESTRICT dataField, + double * GR_RESTRICT value); + + +// helper function that retrieves index for redshift dimension (of cloudy +// tables) via bisection +// - the index is one-indexed +// - the names of variables have not been changed for backwards compatibility +// (it may seem counter-intuitive that clGridDim[1] gives the length of +// clPar2, but that's because in Fortran you would access clGridDim(2) ) +// - NOTE: since we define this function in a header, we must declare it as +// static inline (in C++ we could just declare it as inline) +static inline long long find_zindex(double zr, long long clGridRank, + const long long* clGridDim, + const double* clPar2){ + if (clGridRank > 2){ + long long zindex; + if (zr <= clPar2[0]) { + zindex = 1; + } else if (zr >= clPar2[clGridDim[1]-2]) { + zindex = clGridDim[1]; + } else if (zr >= clPar2[clGridDim[1]-3]) { + zindex = clGridDim[1] - 2; + } else { + zindex = 1; + long long zhighpt = clGridDim[1] - 2; + while ((zhighpt - zindex) > 1) { + long long zmidpt = (long long)((zhighpt + zindex) / 2); + if (zr >= clPar2[zmidpt-1]){ + zindex = zmidpt; + } else { + zhighpt = zmidpt; + } + } + } + return zindex; + } else { + return 1; + } +} + + +/// Calculate temperature and mean molecular weight for tabulated cooling. +/// @details This performs a calculation for a 1D slice from a 3D field. +/// +/// @author Britton Smith +/// @author Matthew Abruzzo +/// @date May, 2015 +/// @remark This was originally written by Britton Smith in May 2015 in +/// Fortran. It was transcribed to C by Matthew Abruzzo in April 2023. +/// +/// @param[in] d 3D density field +/// @param[in] metal 3D metal density field +/// @param[in] e 3D specific internal energy field +/// @param[in] rhoH (precomputed) total H mass density. This only holds +/// values for the 1D slice +/// @param[in] in,jn,kn dimensions of 3D fields (1D array hold ``in`` items) +/// @param[in] is,ie start and (inclusive) end indices of active region +/// (zero-based) +/// @param[in] j,k indices along other dimensions (one-based) +/// @param[out] tgas 1D array to store output temperature values +/// @param[out] mmw 1D array to store output mean molecular weight values +/// @param[in] dom unit conversion to proper number density in code units +/// @param[in] zr current redshift +/// @param[in] temstart, temend start and end of temperature range for rate +/// table +/// @param[in] gamma adiabatic index +/// @param[in] utem temperature units +/// @param[in] imetal flag if metal field is active (0 = no, 1 = yes) +/// @param[in] clGridRank rank of cloudy cooling data grid +/// @param[in] clGridDim array containing dimensions of cloudy data +/// @param[in] clPar1, clPar2, clPar3 arrays containing cloudy grid parameter +/// values. +/// @param[in] clDataSize total size of flattened mmw data array +/// @param[in] clMMW cloudy mmw data +/// @param[in] itmask iteration mask +/// +/// @note +/// None of the pointers are allowed to alias. Other than ``clPar2`` (when +/// ``clGridRank < 2``) or ``clPar3`` (when ``clGridRank < 3``), no parameters +/// should be passed ``NULL`` +void calc_temp1d_cloudy_g( + const gr_float* d, const gr_float* metal, const gr_float* e, + const double* rhoH, + int in, int jn, int kn, int is, int ie, int j, int k, + double * GR_RESTRICT tgas, double * GR_RESTRICT mmw, + double dom, double zr, double temstart, double temend, double gamma, + double utem, int imetal, + long long clGridRank, + const long long* clGridDim, + const double* clPar1, + const double* clPar2, + const double* clPar3, + long long clDataSize, + const double* clMMW, + const int32_t* itmask); + +/// Solve cloudy cooling by interpolating from data tables. +/// @details This performs a calculation for a 1D slice from a 3D field. +/// +/// @author Britton Smith +/// @author Matthew Abruzzo +/// @date September, 2009 +/// @remark This was originally written by Britton Smith in September 2009 in +/// Fortran. It was transcribed to C by Matthew Abruzzo in May 2023. +/// +/// @param[in] d 3D density field +/// @param[in] rhoH (precomputed) total H mass density. This only holds +/// values for the 1D slice +/// @param[in] metallicity (precomputed) metallicity. This only holds +/// values for the 1D slice +/// @param[in] in,jn,kn dimensions of 3D fields (1D array hold ``in`` +/// items) +/// @param[in] is,ie start and (inclusive) end indices of active region +/// (zero-based) +/// @param[in] j,k indices along other dimensions (one-based) +/// @param[in] logtem (precomputed) natural log of temperature values. +/// This only holds values for the 1D slice +/// @param[in,out] edot the heating/cooling contributions computed in this +/// function are used to update this array (the net contributions are added +/// to prexisting values. This only holds values for the 1D slice +/// @param[in] comp2 Temperature of the CMB at the redshift given by ``zr`` +/// @param[in] dom unit conversion to proper number density in code units +/// @param[in] zr current redshift +/// @param[in] icmbTfloor flag to include temperature floor from cmb +/// @param[in] iClHeat flag to include cloudy heating +/// @param[in] iZscale flag to scale cooling by metallicity +/// @param[in] clGridRank rank of cloudy cooling data grid +/// @param[in] clGridDim array containing dimensions of cloudy data +/// @param[in] clPar1, clPar2, clPar3 arrays containing cloudy grid +/// parameter values. +/// @param[in] clDataSize total size of flattened cloudy data arrays +/// @param[in] clCooling cloudy cooling data +/// @param[in] clHeating cloudy heating data +/// @param[in] itmask iteration mask +void cool1d_cloudy_g( + const gr_float* d, // 3D arrays + const double* rhoH, const gr_float* metallicity, // 1D array + int in, int jn, int kn, int is, int ie, int j, int k, + const double* logtem, double * GR_RESTRICT edot, // 1D array + double comp2, double dom, double zr, + int icmbTfloor, int iClHeat, int iZscale, + long long clGridRank, + const long long* clGridDim, + const double* clPar1, const double* clPar2, const double* clPar3, + long long clDataSize, const double* clCooling, const double* clHeating, + const int32_t* itmask); + +#endif /* __INTEROP_FUNCS_H */ diff --git a/src/clib/interop/interpolators_g.c b/src/clib/interop/interpolators_g.c new file mode 100644 index 00000000..2778cbca --- /dev/null +++ b/src/clib/interop/interpolators_g.c @@ -0,0 +1,342 @@ +/*********************************************************************** +/ +/ Implement interpolation routines +/ +/ +/ Copyright (c) 2013, Enzo/Grackle Development Team. +/ +/ Distributed under the terms of the Enzo Public Licence. +/ +/ The full license is in the file LICENSE, distributed with this +/ software. +************************************************************************/ + +#include // log + +#include "interop_funcs.h" + +// There is a huge potential for optimization in these functions. For example: +// - we could take advantage of the regular spacing between parameter values +// and precompute 1/dgridPar1, 1/dgridPar2, 1/dgridPar3, ... ahead of time. +// This would allow us to avoid a lot of expensive division operations +// - we could make these functions directly return the result (rather than +// returning through a pointer argument). This change generally helps with +// compiler optimizations. +// - we could restructure the ordering of the tables (maybe when we read them +// in?) so that we can further reduce the number of calls to the log +// function, which is generally very slow. + + +/// helper function that determines the 1-indexed interpolation index +/// +/// This assumes that parameter is evenly spaced on the grid +static inline gr_int64 get_index_(double input, gr_int64 parLen, + const double *gridPar, double dgridPar) +{ + gr_int64 index = (gr_int64)((input-gridPar[0])/dgridPar)+1; + gr_int64 index_with_floor = (index > 1) ? index : 1; + return ((parLen - 1) < index_with_floor) ? (parLen - 1) : index_with_floor; +}; + +/// helper function to interpolate along a single dimension +/// +/// @note +/// we could probably make this a 3 argument function where the second & third +/// arguments each expect a pointer to a pair of values. The main reason I have +/// avoided this right now is that it may interfere with vectorization +/// +/// @note +/// We may want to use some compiler specific extensions to force inlining of +/// this function (OR define it as a macro) +static inline double interp_(double x, + double xref0, double xref1, + double yref0, double yref1) +{ + double slope = (yref1 - yref0) / (xref1 - xref0); + return (x - xref0) * slope + yref0; +} + +void interpolate_1d_g(double input1, + const gr_int64 * GR_RESTRICT gridDim, // 1 elements + const double * GR_RESTRICT gridPar1, double dgridPar1, + gr_int64 dataSize, const double * GR_RESTRICT dataField, + double * GR_RESTRICT value) +{ + const gr_int64 index1 = get_index_(input1, gridDim[0], gridPar1, dgridPar1); + + // interpolate over parameter 1 + *value = interp_(input1, gridPar1[index1-1], gridPar1[index1], + dataField[index1-1], dataField[index1]); +} + +void interpolate_2d_g(double input1, double input2, + const gr_int64 * GR_RESTRICT gridDim, // 2 elements + const double * GR_RESTRICT gridPar1, double dgridPar1, + const double * GR_RESTRICT gridPar2, double dgridPar2, + gr_int64 dataSize, const double * GR_RESTRICT dataField, + double * GR_RESTRICT value) +{ + double value2[2]; + + const gr_int64 index1 = get_index_(input1, gridDim[0], gridPar1, dgridPar1); + const gr_int64 index2 = get_index_(input2, gridDim[1], gridPar2, dgridPar2); + + for (gr_int64 q=0; q < 2; q++) { + + // interpolate over parameter 2 + gr_int64 int_index = (q+index1-1) * gridDim[1] + index2; + + value2[q] = interp_(input2, gridPar2[index2-1], gridPar2[index2], + dataField[int_index-1], dataField[int_index]); + } + + *value = interp_(input1, gridPar1[index1-1], gridPar1[index1], + value2[0], value2[1]); +} + +void interpolate_3d_g(double input1, double input2, double input3, + const gr_int64 * GR_RESTRICT gridDim, // 3 elements + const double * GR_RESTRICT gridPar1, double dgridPar1, + const double * GR_RESTRICT gridPar2, double dgridPar2, + const double * GR_RESTRICT gridPar3, double dgridPar3, + gr_int64 dataSize, const double * GR_RESTRICT dataField, + double * GR_RESTRICT value) +{ + double value3[2], value2[2]; + + const gr_int64 index1 = get_index_(input1, gridDim[0], gridPar1, dgridPar1); + const gr_int64 index2 = get_index_(input2, gridDim[1], gridPar2, dgridPar2); + const gr_int64 index3 = get_index_(input3, gridDim[2], gridPar3, dgridPar3); + + for (gr_int64 q = 0; q < 2; q++) { + for (gr_int64 w = 0; w < 2; w++) { + + // interpolate over parameter 3 + gr_int64 int_index = ((q+index1-1) * gridDim[1] + + (w+index2-1)) * gridDim[2] + index3; + + value3[w] = interp_(input3, gridPar3[index3-1], gridPar3[index3], + dataField[int_index-1], dataField[int_index]); + } + + // interpolate over parameter 2 + value2[q] = interp_(input2, gridPar2[index2-1], gridPar2[index2], + value3[0], value3[1]); + } + + // interpolate over parameter 1 + *value = interp_(input1, gridPar1[index1-1], gridPar1[index1], + value2[0], value2[1]); +} + +void interpolate_4d_g(double input1, double input2, double input3, + double input4, + const gr_int64 * GR_RESTRICT gridDim, // 4 elements + const double * GR_RESTRICT gridPar1, double dgridPar1, + const double * GR_RESTRICT gridPar2, double dgridPar2, + const double * GR_RESTRICT gridPar3, double dgridPar3, + const double * GR_RESTRICT gridPar4, double dgridPar4, + gr_int64 dataSize, const double * GR_RESTRICT dataField, + double * GR_RESTRICT value) +{ + + double value4[2], value3[2], value2[2]; + + const gr_int64 index1 = get_index_(input1, gridDim[0], gridPar1, dgridPar1); + const gr_int64 index2 = get_index_(input2, gridDim[1], gridPar2, dgridPar2); + const gr_int64 index3 = get_index_(input3, gridDim[2], gridPar3, dgridPar3); + const gr_int64 index4 = get_index_(input4, gridDim[3], gridPar4, dgridPar4); + + for (gr_int64 q = 0; q < 2; q++) { + for (gr_int64 w = 0; w < 2; w++) { + for (gr_int64 e = 0; e < 2; e++) { + + // interpolate over parameter 4 + gr_int64 int_index = (((q+index1-1) * gridDim[1] + + (w+index2-1)) * gridDim[2] + + (e+index3-1)) * gridDim[3] + index4; + + value4[e] = interp_(input4, gridPar4[index4-1], gridPar4[index4], + dataField[int_index-1], dataField[int_index]); + } + + // interpolate over parameter 3 + value3[w] = interp_(input3, gridPar3[index3-1], gridPar3[index3], + value4[0], value4[1]); + } + + // interpolate over parameter 2 + value2[q] = interp_(input2, gridPar2[index2-1], gridPar2[index2], + value3[0], value3[1]); + } + + // interpolate over parameter 1 + *value = interp_(input1, gridPar1[index1-1], gridPar1[index1], + value2[0], value2[1]); +} + +void interpolate_5d_g(double input1, double input2, double input3, + double input4, double input5, + const gr_int64 * GR_RESTRICT gridDim, // 5 elements + const double * GR_RESTRICT gridPar1, double dgridPar1, + const double * GR_RESTRICT gridPar2, double dgridPar2, + const double * GR_RESTRICT gridPar3, double dgridPar3, + const double * GR_RESTRICT gridPar4, double dgridPar4, + const double * GR_RESTRICT gridPar5, double dgridPar5, + gr_int64 dataSize, const double * GR_RESTRICT dataField, + double * GR_RESTRICT value) +{ + double value5[2], value4[2], value3[2], value2[2]; + + const gr_int64 index1 = get_index_(input1, gridDim[0], gridPar1, dgridPar1); + const gr_int64 index2 = get_index_(input2, gridDim[1], gridPar2, dgridPar2); + const gr_int64 index3 = get_index_(input3, gridDim[2], gridPar3, dgridPar3); +#define INDEX_4_BISECTION +#ifdef INDEX_4_BISECTION + // get index 4 with bisection, since not evenly spaced + gr_int64 index4; + if (input4 <= gridPar4[0]) { + index4 = 1; + } else if (input4 >= gridPar4[gridDim[3]-2]) { // -2 isn't a typo + index4 = gridDim[3] - 1; + } else { + index4 = 1; + gr_int64 highPt = gridDim[3]; + while ((highPt - index4) > 1) { + gr_int64 midPt = (gr_int64)((highPt + index4) / 2); + if (input4 >= gridPar4[midPt-1]) { + index4 = midPt; + } else { + highPt = midPt; + } + } + } +#else + gr_int64 index4 = get_index_(input4, gridDim[3], gridPar4, dgridPar4); +#endif /* INDEX_4_BISECTION */ + const gr_int64 index5 = get_index_(input5, gridDim[4], gridPar5, dgridPar5); + + for (gr_int64 q = 0; q < 2; q++) { + for (gr_int64 w = 0; w < 2; w++) { + for (gr_int64 e = 0; e < 2; e++) { + for (gr_int64 r = 0; r < 2; r++) { + + // interpolate over parameter 5 + gr_int64 int_index = ((((q+index1-1) * gridDim[1] + + (w+index2-1)) * gridDim[2] + (e+index3-1)) * + gridDim[3] + (r+index4-1)) * gridDim[4] + + index5; + + value5[r] = interp_(input5, gridPar5[index5-1], gridPar5[index5], + dataField[int_index-1], dataField[int_index]); + } + + // interpolate over parameter 4 + value4[e] = interp_(input4, gridPar4[index4-1], gridPar4[index4], + value5[0], value5[1]); + } + + // interpolate over parameter 3 + value3[w] = interp_(input3, gridPar3[index3-1], gridPar3[index3], + value4[0], value4[1]); + } + + // interpolate over parameter 2 + value2[q] = interp_(input2, gridPar2[index2-1], gridPar2[index2], + value3[0], value3[1]); + } + + *value = interp_(input1, gridPar1[index1-1], gridPar1[index1], + value2[0], value2[1]); +} + + +// OTHER FUNCTIONS +// Interpolation in 2 dimensions but with a 3D grid. +// This is used for interpolating from just the last +// slice in the datacube before the redshift where +// the UV background turns on. +static inline double interpolate_2Df3D_g(double input1, double input3, + const gr_int64* GR_RESTRICT gridDim, // 3 elements + const double* GR_RESTRICT gridPar1, double dgridPar1, + gr_int64 index2, + const double* GR_RESTRICT gridPar3, double dgridPar3, + gr_int64 dataSize, + const double* dataField) +{ + double value3[2]; + + // Calculate interpolation indices + const gr_int64 index1 = get_index_(input1, gridDim[0], gridPar1, dgridPar1); + const gr_int64 index3 = get_index_(input3, gridDim[2], gridPar3, dgridPar3); + + for (gr_int64 q = 0; q < 2; q++) { // interpolate over parameter 3 + + gr_int64 int_index = ((q+index1-1) * gridDim[1] + (index2-1)) * + gridDim[2] + index3; + + value3[q] = interp_(input3, gridPar3[index3-1], gridPar3[index1], + dataField[int_index-1], dataField[int_index]); + } + + // interpolate over parameter 1 + return interp_(input1, gridPar1[index1-1], gridPar1[index1], + value3[0], value3[1]); +} + + +void interpolate_3dz_g(double input1, double input2, double input3, + const gr_int64 * GR_RESTRICT gridDim, // 3 elements + const double * GR_RESTRICT gridPar1, double dgridPar1, + const double * GR_RESTRICT gridPar2, gr_int64 index2, + const double * GR_RESTRICT gridPar3, double dgridPar3, + gr_int64 dataSize, const double * GR_RESTRICT dataField, + gr_int64 end_int, double * GR_RESTRICT value) +{ + if (end_int == 1) { + *value = interpolate_2Df3D_g(input1, input3, + gridDim, + gridPar1, dgridPar1, + index2, + gridPar3, dgridPar3, + dataSize, dataField); + return; + } + + double value3[2], value2[2]; + + // Calculate interpolation indices + const gr_int64 index1 = get_index_(input1, gridDim[0], gridPar1, dgridPar1); + const gr_int64 index3 = get_index_(input3, gridDim[2], gridPar3, dgridPar3); + + + // it turns out that precomputing the following 2 variables reduces runtime + // appreciably (because the C compiler can't automatically hoist these + // calculations out of the loop) + const double par2_slope_denom = log((1+gridPar2[index2]) / + (1+gridPar2[index2-1])); + const double par2_offset_from_grid = log((1+input2)/(1+gridPar2[index2-1])); + + + // preliminary testing on gcc 9.4 suggests that unrolling the outer loop + // could speed this function up by ~10% + for (gr_int64 q = 0; q < 2; q++) { + for (gr_int64 w = 0; w < 2; w++) { + + // interpolate over parameter 3 + gr_int64 int_index = ((q+index1-1) * gridDim[1] + + (w+index2-1)) * gridDim[2] + index3; + + value3[w] = interp_(input3, gridPar3[index3-1], gridPar3[index3], + dataField[int_index-1], dataField[int_index]); + } + + // interpolate over parameter 2 + double slope = (value3[1] - value3[0]) / par2_slope_denom; + value2[q] = par2_offset_from_grid * slope + value3[0]; + } + + // interpolate over parameter 1 + *value = interp_(input1, gridPar1[index1-1], gridPar1[index1], + value2[0], value2[1]); +} diff --git a/src/clib/interop/shims_g.F b/src/clib/interop/shims_g.F new file mode 100644 index 00000000..19d27915 --- /dev/null +++ b/src/clib/interop/shims_g.F @@ -0,0 +1,220 @@ +!======================================================================= +!////////////////// SUBROUTINE CALC_TEMP1D_CLOUDY_SHIM_G \\\\\\\\\\\\\\\\\\ + + subroutine calc_temp1d_cloudy_shim_g(d, metal, e, rhoH, + & in, jn, kn, is, ie, j, k, + & tgas, mmw, dom, zr, + & temstart, temend, + & gamma, utem, imetal, + & clGridRank, clGridDim, + & clPar1, clPar2, clPar3, + & clDataSize, clMMW, + & itmask) + +! +! SHIM TO CALCULATE TEMPERATURE AND MEAN MOLECULAR WEIGHT FROM CLOUDY TABLES +! +! written by: Matthew Abruzzo +! date: April, 2023 +! +! PURPOSE: +! Calculate temperature and mean molecular weight for tabulated cooling. +! This converts itmask from a logical dtype to integer*4 since logical +! dtypes can't be trivially passed to a C function. +! +! INPUTS: +! in,jn,kn - dimensions of 3D fields +! +! d - density field +! metal - metal density field +! e - internal energy field +! +! rhoH - total H mass density +! +! is,ie - start and end indices of active region (zero based) +! tgas - temperature values +! mmw - mean molecular weight values +! +! dom - unit conversion to proper number density in code units +! zr - current redshift +! temstart, temend - start and end of temperature range for rate table +! gamma - adiabatic index +! utem - temperature units +! imetal - flag if metal field is active (0 = no, 1 = yes) +! +! clGridRank - rank of cloudy cooling data grid +! clGridDim - array containing dimensions of cloudy data +! clPar1, clPar2, clPar3 - arrays containing cloudy grid parameter values +! clDataSize - total size of flattened 1D cooling data array +! clMMW - cloudy mmw data +! +! itmask - iteration mask +! +! OUTPUTS: +! fills temperature and mmw arrays +! +! PARAMETERS: +! +!----------------------------------------------------------------------- + + USE ISO_C_BINDING + implicit NONE +#include "../grackle_fortran_types.def" +#include "interop_funcs.def" + +! General Arguments + + integer in, jn, kn, is, ie, j, k, + & imetal + + real*8 dom, zr, gamma, temstart, temend, utem + R_PREC d(in,jn,kn), metal(in,jn,kn), e(in,jn,kn) + real*8 rhoH(in), tgas(in), + & mmw(in), logtem(in) + +! Cloudy parameters and data + + integer*8 clGridRank, clDataSize, + & clGridDim(clGridRank) + real*8 clPar1(clGridDim(1)), clPar2(clGridDim(2)), + & clPar3(clGridDim(3)), + & clMMW(clDataSize) + +! Iteration mask + + logical itmask(in) + +! Local variable + integer*4 temp_itmask(in) + integer i + + + do i=is+1, ie+1 + if ( itmask(i) ) then + temp_itmask(i) = 1 + else + temp_itmask(i) = 0 + endif + enddo + + call calc_temp1d_cloudy_g(d, metal, e, rhoH, + & in, jn, kn, is, ie, j, k, + & tgas, mmw, dom, zr, + & temstart, temend, + & gamma, utem, imetal, + & clGridRank, clGridDim, + & clPar1, clPar2, clPar3, + & clDataSize, clMMW, + & temp_itmask) + + return + end + +!======================================================================= +!////////////////// SUBROUTINE COOL1D_CLOUDY_SHIM_G \\\\\\\\\\\\\\\\\\ + + subroutine cool1d_cloudy_shim_g(d, rhoH, metallicity, + & in, jn, kn, is, ie, j, k, + & logtem, edot, comp2, dom, zr, + & icmbTfloor, iClHeat, iZscale, + & clGridRank, clGridDim, + & clPar1, clPar2, clPar3, + & clDataSize, clCooling, clHeating, + & itmask) + +! +! SHIM TO SOLVE CLOUDY METAL COOLING +! +! written by: Matthew Abruzzo +! date: May, 2023 +! +! PURPOSE: +! Solve cloudy cooling by interpolating from the data. This converts itmask +! itmask from a logical dtype to integer*4 since logical +! dtypes can't be trivially passed to a C function. +! +! INPUTS: +! in,jn,kn - dimensions of 3D fields +! +! d - total density field +! +! rhoH - total H mass density +! metallicity - metallicity +! +! is,ie - start and end indices of active region (zero based) +! logtem - natural log of temperature values +! +! dom - unit conversion to proper number density in code units +! zr - current redshift +! +! icmbTfloor - flag to include temperature floor from cmb +! iClHeat - flag to include cloudy heating +! iZscale - flag to scale cooling by metallicity +! clGridRank - rank of cloudy cooling data grid +! clGridDim - array containing dimensions of cloudy data +! clPar1, clPar2, clPar3 - arrays containing cloudy grid parameter values +! clDataSize - total size of flattened 1D cooling data array +! clCooling - cloudy cooling data +! clHeating - cloudy heating data +! +! itmask - iteration mask +! +! OUTPUTS: +! update edot with heating/cooling contributions from metals +! +! PARAMETERS: +! +!----------------------------------------------------------------------- + + USE ISO_C_BINDING + implicit NONE +#include "../grackle_fortran_types.def" +#include "interop_funcs.def" + +! General Arguments + + integer in, jn, kn, is, ie, j, k + + real*8 comp2, dom, zr + R_PREC d(in,jn,kn) + real*8 rhoH(in), metallicity(in), logtem(in), + & edot(in) + +! Cloudy parameters and data + + integer icmbTfloor, iClHeat, iZscale + integer*8 clGridRank, clDataSize, + & clGridDim(clGridRank) + real*8 clPar1(clGridDim(1)), clPar2(clGridDim(2)), + & clPar3(clGridDim(3)), + & clCooling(clDataSize), clHeating(clDataSize) + +! Iteration mask + + logical itmask(in) + + +! Local variable + integer*4 temp_itmask(in) + integer i + + + do i=is+1, ie+1 + if ( itmask(i) ) then + temp_itmask(i) = 1 + else + temp_itmask(i) = 0 + endif + enddo + + call cool1d_cloudy_g(d, rhoH, metallicity, + & in, jn, kn, is, ie, j, k, + & logtem, edot, comp2, dom, zr, + & icmbTfloor, iClHeat, iZscale, + & clGridRank, clGridDim, + & clPar1, clPar2, clPar3, + & clDataSize, clCooling, clHeating, + & temp_itmask) + + return + end diff --git a/src/example/Makefile b/src/example/Makefile index 81382f75..01ddfd42 100644 --- a/src/example/Makefile +++ b/src/example/Makefile @@ -186,6 +186,72 @@ fortran_example: $(MODULES) fortran_example.o echo "Failed! See $(OUTPUT) for error messages"; \ fi) +#----------------------------------------------------------------------- +# Unit test for calc_temp1d_cloudy +#----------------------------------------------------------------------- + +TEST_CALC_TEMP1D_CLOUDY_OBJ = \ + utest_calc_temp1d_cloudy/test_calc_temp1d.o \ + utest_calc_temp1d_cloudy/calc_temp1d_cloudy_g.o + +test_calc_temp1d_cloudy: $(TEST_CALC_TEMP1D_CLOUDY_OBJ) utest_helpers.hpp + @rm -f $@ + @echo "Linking" + @($(CXX) $(LDFLAGS) -o test_calc_temp1d_cloudy $(TEST_CALC_TEMP1D_CLOUDY_OBJ) $(GRACKLE_LIB) $(LIBS) -lm) >& $(OUTPUT) + @(if [ -e $@ ]; then \ + echo "Success!"; \ + else \ + echo "$(CXX) $(LDFLAGS) -o test_calc_temp1d_cloudy $(TEST_CALC_TEMP1D_CLOUDY_OBJ) $(GRACKLE_LIB) $(LIBS) -lm" >> temp1; \ + cat temp1 $(OUTPUT) > temp2; \ + rm -f temp1; \ + mv -f temp2 $(OUTPUT); \ + echo "Failed! See $(OUTPUT) for error messages"; \ + fi) + +#----------------------------------------------------------------------- +# Unit test for cool1d_cloudy +#----------------------------------------------------------------------- + +TEST_COOL1D_CLOUDY_OBJ = \ + utest_cool1d_cloudy/test_cool1d_cloudy.o \ + utest_cool1d_cloudy/cool1d_cloudy_g.o + +test_cool1d_cloudy: $(TEST_COOL1D_CLOUDY_OBJ) utest_helpers.hpp + @rm -f $@ + @echo "Linking" + @($(CXX) $(LDFLAGS) -o test_cool1d_cloudy $(TEST_COOL1D_CLOUDY_OBJ) $(GRACKLE_LIB) $(LIBS) -lm) >& $(OUTPUT) + @(if [ -e $@ ]; then \ + echo "Success!"; \ + else \ + echo "$(CXX) $(LDFLAGS) -o test_cool1d_cloudy $(TEST_COOL1D_CLOUDY_OBJ) $(GRACKLE_LIB) $(LIBS) -lm" >> temp1; \ + cat temp1 $(OUTPUT) > temp2; \ + rm -f temp1; \ + mv -f temp2 $(OUTPUT); \ + echo "Failed! See $(OUTPUT) for error messages"; \ + fi) + +#----------------------------------------------------------------------- +# Unit test for interpolators +#----------------------------------------------------------------------- + +TEST_INTERPOLATORS_OBJ = \ + utest_interpolators/test_interpolators.o \ + utest_interpolators/interpolators_g.o + +test_interpolators: $(TEST_INTERPOLATORS_OBJ) utest_helpers.hpp + @rm -f $@ + @echo "Linking" + @($(CXX) $(LDFLAGS) -o test_interpolators $(TEST_INTERPOLATORS_OBJ) $(GRACKLE_LIB) $(LIBS) -lm) >& $(OUTPUT) + @(if [ -e $@ ]; then \ + echo "Success!"; \ + else \ + echo "$(CXX) $(LDFLAGS) -o test_cool1d_cloudy $(TEST_COOL1D_CLOUDY_OBJ) $(GRACKLE_LIB) $(LIBS) -lm" >> temp1; \ + cat temp1 $(OUTPUT) > temp2; \ + rm -f temp1; \ + mv -f temp2 $(OUTPUT); \ + echo "Failed! See $(OUTPUT) for error messages"; \ + fi) + #----------------------------------------------------------------------- # Implicit rules #----------------------------------------------------------------------- @@ -257,7 +323,7 @@ help: #----------------------------------------------------------------------- clean: - -@rm -f *.o *.mod *.f *.f90 *~ *.exe $(OUTPUT) cxx_example cxx_omp_example cxx_grid_example c_example c_local_example fortran_example + -@rm -f *.o *.mod *.f *.f90 *~ *.exe utest_calc_temp1d_cloudy/*.o utest_cool1d_cloudy/*.o utest_interpolators/*.o $(OUTPUT) cxx_example cxx_omp_example cxx_grid_example c_example c_local_example fortran_example test_calc_temp1d_cloudy test_cool1d_cloudy test_interpolators #----------------------------------------------------------------------- # Include configuration targets diff --git a/src/clib/calc_temp1d_cloudy_g.F b/src/example/utest_calc_temp1d_cloudy/calc_temp1d_cloudy_g.F similarity index 93% rename from src/clib/calc_temp1d_cloudy_g.F rename to src/example/utest_calc_temp1d_cloudy/calc_temp1d_cloudy_g.F index 0495d955..c561a755 100644 --- a/src/clib/calc_temp1d_cloudy_g.F +++ b/src/example/utest_calc_temp1d_cloudy/calc_temp1d_cloudy_g.F @@ -1,7 +1,10 @@ !======================================================================= !////////////////// SUBROUTINE CALC_TEMP1D_CLOUDY_G \\\\\\\\\\\\\\\\\\ - subroutine calc_temp1d_cloudy_g(d, metal, e, rhoH, +! NOTE: "_fortran" has been appended to the name of this subroutine to avoid +! name collisions when we include interop_funcs.def file + + subroutine calc_temp1d_cloudy_g_fortran(d, metal, e, rhoH, & in, jn, kn, is, ie, j, k, & tgas, mmw, dom, zr, & temstart, temend, @@ -55,8 +58,10 @@ subroutine calc_temp1d_cloudy_g(d, metal, e, rhoH, ! !----------------------------------------------------------------------- + USE ISO_C_BINDING implicit NONE -#include "grackle_fortran_types.def" +#include "../../clib/grackle_fortran_types.def" +#include "../../clib/interop/interop_funcs.def" ! General Arguments @@ -78,7 +83,7 @@ subroutine calc_temp1d_cloudy_g(d, metal, e, rhoH, ! Iteration mask - logical itmask(in) + integer*4 itmask(in) ! Parameters @@ -93,7 +98,7 @@ subroutine calc_temp1d_cloudy_g(d, metal, e, rhoH, integer*8 zindex, zmidpt, zhighpt real*8 dclPar(clGridRank), inv_log10, & muold, munew - logical end_int + integer*8 end_int ! Slice locals @@ -102,7 +107,7 @@ subroutine calc_temp1d_cloudy_g(d, metal, e, rhoH, !\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\///////////////////////////////// !======================================================================= - end_int = .false. + end_int = 0 inv_log10 = 1._DKIND / log(10._DKIND) @@ -129,7 +134,7 @@ subroutine calc_temp1d_cloudy_g(d, metal, e, rhoH, zindex = 1 else if (zr .ge. clPar2(clGridDim(2)-1)) then zindex = clGridDim(2) - end_int = .true. + end_int = 1 else if (zr .ge. clPar2(clGridDim(2)-2)) then zindex = clGridDim(2) - 2 else @@ -148,14 +153,14 @@ subroutine calc_temp1d_cloudy_g(d, metal, e, rhoH, endif do i=is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. 0 ) then ! Calculate proper log(n_H) log_n_h(i) = log10(rhoH(i) * dom) endif enddo do i=is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. 0 ) then munew = 1._DKIND do ti = 1, ti_max muold = munew diff --git a/src/example/utest_calc_temp1d_cloudy/test_calc_temp1d.C b/src/example/utest_calc_temp1d_cloudy/test_calc_temp1d.C new file mode 100644 index 00000000..1b7a808f --- /dev/null +++ b/src/example/utest_calc_temp1d_cloudy/test_calc_temp1d.C @@ -0,0 +1,164 @@ +#include +#include + +#include "../utest_helpers.hpp" + +extern "C" { + #include "../../clib/interop/interop_funcs.h" + + // legacy version of the function + extern void FORTRAN_NAME(calc_temp1d_cloudy_g_fortran)( + const gr_float* d, const gr_float* metal, // 3D array + const gr_float* e, // 3D array + const double* rhoH, // 1D array + const int* in, const int* jn, const int* kn, + const int* is, const int* ie, const int* j, const int* k, + double* tgas, double*mmw, + const double* dom, const double* zr, + const double* temstart, const double* temend, const double* gamma, + const double* utem, const int* imetal, + const long long* clGridRank, + const long long* clGridDim, + const double* clPar1, + const double* clPar2, + const double* clPar3, + const long long* clDataSize, + const double* clMMW, + const int32_t* itmask); +} + + +struct calc_temp_outputs{ + std::vector tgas; + std::vector mmw; +}; + +calc_temp_outputs run_test(DummyGrackleConfig& config, + bool use_fortran, bool slc_from_3D_arr = true){ + + /* Calculate temperature units. */ + const double temperature_units = get_temperature_units(&(config.units)); + + // initialize physical quantities! + const gr_float Temp0 = 1.0e3; // background temperature (in K) + + const std::size_t length = 60; + std::vector density(length); + std::vector metal_density(length); + std::vector eint(length); + std::vector rhoH(length); // double is NOT a typo + for (std::size_t i = 0; i < length; i++) { + density[i] = gr_float(i+1); + metal_density[i] = gr_float(config.chem_data.SolarMetalFractionByMass * + density[i]); + rhoH[i] = config.chem_data.HydrogenFractionByMass * density[i]; + eint[i] = gr_float(Temp0 / temperature_units); + } + + // prepare some other arguments + // -> when running in 3D, we aren't + const int in = (!slc_from_3D_arr) ? int(length) : 5; + const int jn = (!slc_from_3D_arr) ? 1 : 4; + const int kn = (!slc_from_3D_arr) ? 1 : 3; + + if (std::size_t(in * jn * kn) > length) { error("something is wrong\n"); } + + const int j = jn; // remember, these are currently 1-indexed + const int k = kn; + + const int is = 0; // not a typo! (even though other indices are 1-indexed) + const int ie = in - 1; // not a typo! + + const double aye = config.units.a_value; // expansion factor (in code units) + const double dom = config.units.density_units * pow(aye,3) / mh; + const double zr = 1.0/(aye * config.units.a_units) - 1.; + const int imetal = 1; + + // prepare the output arguments + std::vector tgas(in); + std::vector mmw(in); + + const std::vector itmask(std::size_t(in), int32_t(1)); + + if (use_fortran) { + FORTRAN_NAME(calc_temp1d_cloudy_g_fortran) + (density.data(), metal_density.data(), eint.data(), + rhoH.data(), + &in, &jn, &kn, &is, &ie, &j, &k, + tgas.data(), mmw.data(), + &dom, &zr, + &config.chem_data.TemperatureStart, &config.chem_data.TemperatureEnd, + &config.chem_data.Gamma, &temperature_units, &imetal, + &config.chem_rates.cloudy_primordial.grid_rank, // clGridRank + config.chem_rates.cloudy_primordial.grid_dimension, // clGridDim + config.chem_rates.cloudy_primordial.grid_parameters[0], // clPar1 + config.chem_rates.cloudy_primordial.grid_parameters[1], // clPar2 + config.chem_rates.cloudy_primordial.grid_parameters[2], // clPar3 + &config.chem_rates.cloudy_primordial.data_size, // clDataSize + config.chem_rates.cloudy_primordial.mmw_data, // clMMW + itmask.data()); + + } else { + calc_temp1d_cloudy_g + (density.data(), metal_density.data(), eint.data(), + rhoH.data(), + in, jn, kn, is, ie, j, k, + tgas.data(), mmw.data(), + dom, zr, + config.chem_data.TemperatureStart, config.chem_data.TemperatureEnd, + config.chem_data.Gamma, temperature_units, imetal, + config.chem_rates.cloudy_primordial.grid_rank, // clGridRank + config.chem_rates.cloudy_primordial.grid_dimension, // clGridDim + config.chem_rates.cloudy_primordial.grid_parameters[0], // clPar1 + config.chem_rates.cloudy_primordial.grid_parameters[1], // clPar2 + config.chem_rates.cloudy_primordial.grid_parameters[2], // clPar3 + config.chem_rates.cloudy_primordial.data_size, // clDataSize + config.chem_rates.cloudy_primordial.mmw_data, // clMMW + itmask.data()); + } + + return {tgas, mmw}; +} + +int main(void){ + for (int n_tab_dims = 1; n_tab_dims < 4; n_tab_dims++) { + printf("\nConsidering a %dD table of MMW values\n", n_tab_dims); + + std::vector z_vals {0.0}; + if (n_tab_dims == 3) { + z_vals.push_back(0.13242); + z_vals.push_back(0.5); + z_vals.push_back(10); + + DummyGrackleConfig conf(n_tab_dims,0.0); + double* clPar2 = conf.chem_rates.cloudy_primordial.grid_parameters[1]; + long long* clGridDim = conf.chem_rates.cloudy_primordial.grid_dimension; + z_vals.push_back(clPar2[clGridDim[1]-3]); + z_vals.push_back(clPar2[clGridDim[1]-2]); + z_vals.push_back(clPar2[clGridDim[1]-1]); + z_vals.push_back(clPar2[clGridDim[1]-1] * 1.01); + } + + for (const auto& z_val : z_vals) { + DummyGrackleConfig config(n_tab_dims,z_val); + + for (bool slc_from_3D_arr : {false, true}) { + const char* descr = (slc_from_3D_arr) ? "slice from 3D arr" : "1D arr"; + printf("-> z = %g, comparing %s\n", z_val, descr); + + //printf("using the c version:\n"); fflush(stdout); + calc_temp_outputs actual = run_test(config, false, slc_from_3D_arr); + + //printf("using the fortran version:\n"); fflush(stdout); + calc_temp_outputs reference = run_test(config, true, slc_from_3D_arr); + + compare_values(actual.mmw, reference.mmw, 0.0, 0.0, + "**Error during comparison of mmw**"); + compare_values(actual.tgas, reference.tgas, 0.0, 0.0, + "**Error during comparison of tgas**"); + } + } + } + + return 0; +} diff --git a/src/clib/cool1d_cloudy_g.F b/src/example/utest_cool1d_cloudy/cool1d_cloudy_g.F similarity index 94% rename from src/clib/cool1d_cloudy_g.F rename to src/example/utest_cool1d_cloudy/cool1d_cloudy_g.F index 6d77b318..de7f8f86 100644 --- a/src/clib/cool1d_cloudy_g.F +++ b/src/example/utest_cool1d_cloudy/cool1d_cloudy_g.F @@ -1,7 +1,10 @@ !======================================================================= !//////////////////// SUBROUTINE COOL1D_CLOUDY_G \\\\\\\\\\\\\\\\\\\\\ - subroutine cool1d_cloudy_g(d, rhoH, metallicity, +! NOTE: "_fortran" has been appended to the name of this subroutine to avoid +! name collisions when we include interop_funcs.def file + + subroutine cool1d_cloudy_g_fortran(d, rhoH, metallicity, & in, jn, kn, is, ie, j, k, & logtem, edot, comp2, dom, zr, & icmbTfloor, iClHeat, iZscale, @@ -52,8 +55,10 @@ subroutine cool1d_cloudy_g(d, rhoH, metallicity, ! !----------------------------------------------------------------------- + USE ISO_C_BINDING implicit NONE -#include "grackle_fortran_types.def" +#include "../../clib/grackle_fortran_types.def" +#include "../../clib/interop/interop_funcs.def" ! General Arguments @@ -75,7 +80,7 @@ subroutine cool1d_cloudy_g(d, rhoH, metallicity, ! Iteration mask - logical itmask(in) + integer*4 itmask(in) ! Parameters @@ -84,7 +89,7 @@ subroutine cool1d_cloudy_g(d, rhoH, metallicity, integer i, get_heat integer*8 zindex, zmidpt, zhighpt real*8 dclPar(clGridRank), inv_log10, log10_tCMB - logical end_int + integer*8 end_int ! Slice locals @@ -95,7 +100,7 @@ subroutine cool1d_cloudy_g(d, rhoH, metallicity, !\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\///////////////////////////////// !======================================================================= - end_int = .false. + end_int = 0 get_heat = iClHeat inv_log10 = 1._DKIND / log(10._DKIND) @@ -115,7 +120,7 @@ subroutine cool1d_cloudy_g(d, rhoH, metallicity, endif do i=is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. 0 ) then log10tem(i) = logtem(i) * inv_log10 @@ -133,7 +138,7 @@ subroutine cool1d_cloudy_g(d, rhoH, metallicity, zindex = 1 else if (zr .ge. clPar2(clGridDim(2)-1)) then zindex = clGridDim(2) - end_int = .true. + end_int = 1 get_heat = 0 else if (zr .ge. clPar2(clGridDim(2)-2)) then zindex = clGridDim(2) - 2 diff --git a/src/example/utest_cool1d_cloudy/test_cool1d_cloudy.C b/src/example/utest_cool1d_cloudy/test_cool1d_cloudy.C new file mode 100644 index 00000000..6ced4a91 --- /dev/null +++ b/src/example/utest_cool1d_cloudy/test_cool1d_cloudy.C @@ -0,0 +1,244 @@ +#include +#include + +#include "../utest_helpers.hpp" + +extern "C" { + #include "../../clib/interop/interop_funcs.h" + + extern void FORTRAN_NAME(cool1d_cloudy_g_fortran)( + const gr_float* d, // 3D arrays + const double* rhoH, const gr_float* metallicity, // 1D array + const int* in, const int* jn, const int* kn, + const int* is, const int* ie, const int* j, const int* k, + const double* logtem, double* edot, // 1D array + const double* comp2, const double* dom, const double* zr, + const int* icmbTfloor, const int* iClHeat, const int* iZscale, + const long long* clGridRank, + const long long* clGridDim, + const double* clPar1, const double* clPar2, const double* clPar3, + const long long* clDataSize, + const double* clCooling, const double* clHeating, + const int32_t* itmask); +} + +std::vector run_test(DummyGrackleConfig& config, + bool use_metal_table, + bool use_fortran, + bool slc_from_3D_arr = true){ + + /* Calculate temperature units. */ + const double temperature_units = get_temperature_units(&(config.units)); + + // initialize physical quantities! + const gr_float Temp0 = 1.0e3; // background temperature (in K) + + const double z_solar = config.chem_data.SolarMetalFractionByMass; + + const std::size_t length = 60; + std::vector density(length); + std::vector metal_density(length); + std::vector eint(length); + std::vector rhoH(length); // double is NOT a typo + for (std::size_t i = 0; i < length; i++) { + density[i] = gr_float(i+1); + metal_density[i] = gr_float(z_solar * density[i]); + rhoH[i] = config.chem_data.HydrogenFractionByMass * density[i]; + eint[i] = gr_float(Temp0 / temperature_units); + } + + // prepare some other arguments + // -> when running in 3D, we aren't + const int in = (!slc_from_3D_arr) ? int(length) : 5; + const int jn = (!slc_from_3D_arr) ? 1 : 4; + const int kn = (!slc_from_3D_arr) ? 1 : 3; + + if (std::size_t(in * jn * kn) > length) { error("something is wrong\n"); } + + const int j = jn; // remember, these are currently 1-indexed + const int k = kn; + + const int is = 0; // not a typo! (even though other indices are 1-indexed) + const int ie = in - 1; // not a typo! + + const double aye = config.units.a_value; // expansion factor (in code units) + const double dom = config.units.density_units * pow(aye,3) / mh; + const double zr = 1.0/(aye * config.units.a_units) - 1.; + const int imetal = 1; + + // prepare the output arguments + std::vector tgas(in); + std::vector mmw(in); + + const std::vector itmask(std::size_t(in), int32_t(1)); + + // calculate temperature (as part of the setup) + calc_temp1d_cloudy_g + (density.data(), metal_density.data(), eint.data(), + rhoH.data(), + in, jn, kn, is, ie, j, k, + tgas.data(), mmw.data(), + dom, zr, + config.chem_data.TemperatureStart, config.chem_data.TemperatureEnd, + config.chem_data.Gamma, temperature_units, imetal, + config.chem_rates.cloudy_primordial.grid_rank, // clGridRank + config.chem_rates.cloudy_primordial.grid_dimension, // clGridDim + config.chem_rates.cloudy_primordial.grid_parameters[0], // clPar1 + config.chem_rates.cloudy_primordial.grid_parameters[1], // clPar2 + config.chem_rates.cloudy_primordial.grid_parameters[2], // clPar3 + config.chem_rates.cloudy_primordial.data_size, // clDataSize + config.chem_rates.cloudy_primordial.mmw_data, // clMMW + itmask.data()); + + // prepare a couple other arguments! + std::vector metallicity(length); // double is NOT a typo + for (std::size_t i = 0; i < length; i++) { + metallicity[i] = metal_density[i] / density[i] / z_solar; + } + + // these use natural logarithms (intentional!) + const double logtem0 = log(config.chem_data.TemperatureStart); + const double logtem9 = log(config.chem_data.TemperatureEnd); + + std::vector logtem(length); // double is NOT a typo + for (std::size_t i = 0; i < length; i++) { + if (tgas[i] <= 0) { logtem[i] = logtem0; } // tgas wasn't computed here + logtem[i] = fmin(fmax(log(tgas[i]), logtem0), logtem9); + } + + const int iClHeat = config.chem_data.UVbackground; + + const cloudy_data& cloudy_data = (use_metal_table) ? + config.chem_rates.cloudy_metal : config.chem_rates.cloudy_primordial; + + const int iZscale = (use_metal_table) ? 1 : 0; + const int icmbTfloor = (use_metal_table) ? + config.chem_data.cmb_temperature_floor : 0; + + // things to vary (in the future): + // metal cooling: icmbTfloor (when using metal cooling) + // in both cases, vary iClHeat + + // compton cooling coefficient + const double comp2 = 2.73 * (1.0 + zr); + + std::vector edot(length, 0.0); // all outputs are set to 0 + // now we actually compute the cooling time + if (use_fortran) { + FORTRAN_NAME(cool1d_cloudy_g_fortran)( + density.data(), // 3D arrays + rhoH.data(), metallicity.data(), // 1D array + &in, &jn, &kn, &is, &ie, &j, &k, + logtem.data(), edot.data(), // 1D array + &comp2, &dom, &zr, + &icmbTfloor, &iClHeat, &iZscale, + &cloudy_data.grid_rank, // clGridRank + cloudy_data.grid_dimension, // clGridDim + cloudy_data.grid_parameters[0], // clPar1 + cloudy_data.grid_parameters[1], // clPar2 + cloudy_data.grid_parameters[2], // clPar3 + &cloudy_data.data_size, // clDataSize + cloudy_data.cooling_data, // clCooling + cloudy_data.heating_data, // clHeating + itmask.data()); + } else { + cool1d_cloudy_g( + density.data(), // 3D arrays + rhoH.data(), metallicity.data(), // 1D array + in, jn, kn, is, ie, j, k, + logtem.data(), edot.data(), // 1D array + comp2, dom, zr, + icmbTfloor, iClHeat, iZscale, + cloudy_data.grid_rank, // clGridRank + cloudy_data.grid_dimension, // clGridDim + cloudy_data.grid_parameters[0], // clPar1 + cloudy_data.grid_parameters[1], // clPar2 + cloudy_data.grid_parameters[2], // clPar3 + cloudy_data.data_size, // clDataSize + cloudy_data.cooling_data, // clCooling + cloudy_data.heating_data, // clHeating + itmask.data()); + } + + return edot; +} + + +int main(void){ + for (int n_tab_dims = 1; n_tab_dims < 4; n_tab_dims++) { + printf("\nConsidering a %dD table of values\n", n_tab_dims); + + std::vector z_vals {0.0}; + if (n_tab_dims == 3) { + z_vals.push_back(0.13242); + z_vals.push_back(0.5); + z_vals.push_back(10); + + DummyGrackleConfig conf(n_tab_dims,0.0); + double* clPar2 = conf.chem_rates.cloudy_primordial.grid_parameters[1]; + long long* clGridDim = conf.chem_rates.cloudy_primordial.grid_dimension; + z_vals.push_back(clPar2[clGridDim[1]-3]); + z_vals.push_back(clPar2[clGridDim[1]-2]); + z_vals.push_back(clPar2[clGridDim[1]-1]); + z_vals.push_back(clPar2[clGridDim[1]-1] * 1.01); + } + + for (const auto& z_val : z_vals) { + for (int j = 0; j < 3; j++){ + std::string descr; + bool use_metal_table; + bool cmb_temperature_floor; + + // there's no reason to apply the cmb_temperature_floor when + // primordial chemistry is used + + if (j == 0) { + use_metal_table = true; + descr = "metals, cmb floor"; + cmb_temperature_floor = true; + } else if (j == 1) { + use_metal_table = true; + descr = "metals, no cmb floor"; + cmb_temperature_floor = false; + } else if (j == 2) { + use_metal_table = false; + descr = "primordial, no cmb floor"; + cmb_temperature_floor = true; + } else { + error("SOMETHING IS WRONG"); + } + + for (bool use_UVbackground : {false, true}) { + std::string UV_descr; + if (!use_UVbackground) { + UV_descr = "cool-only"; + } else if (n_tab_dims == 3) { + UV_descr = "heat & cool"; + } else { + continue; // can't model UV background + } + + DummyGrackleConfig config(n_tab_dims,z_val, use_UVbackground, + cmb_temperature_floor); + + for (bool slc_from_3D_arr : {false, true}) { + const char* descr2 = (slc_from_3D_arr) ? + "slice from 3D arr" : "1D arr"; + printf("-> z = %g, %s, %s, comparing %s\n", + z_val, descr.c_str(), UV_descr.c_str(), descr2); + + std::vector actual = run_test(config, use_metal_table, + false, slc_from_3D_arr); + std::vector reference = run_test(config, use_metal_table, + true, slc_from_3D_arr); + + //printf("%.15e, %.15e\n", actual[0], reference[0]); + + compare_values(actual, reference, 0.0, 0.0, + "**Error during comparison of edot**"); + } + } + } + } + } +} diff --git a/src/example/utest_helpers.hpp b/src/example/utest_helpers.hpp new file mode 100644 index 00000000..491e4814 --- /dev/null +++ b/src/example/utest_helpers.hpp @@ -0,0 +1,307 @@ +#include // abort +#include // stderr, vfprintf +#include // va_list, va_start, va_end + +#include + +#include +#include + +#define mh 1.67262171e-24 + +#ifndef OMIT_LEGACY_INTERNAL_GRACKLE_FUNC +#define OMIT_LEGACY_INTERNAL_GRACKLE_FUNC +#endif + +extern "C" { + #include + #include "../clib/grackle_macros.h" +} + +#undef max +#undef min + +[[noreturn]] void error(const char *fmt, ...){ + + // parse argument list + va_list arg_list; + va_start(arg_list, fmt); // access variable arguments after fmt + vfprintf(stderr, fmt, arg_list); // print to stderr + va_end(arg_list); // cleanup variable arguments + abort(); // exit program with nonzero exit code +} + +inline std::string vec_to_string(const std::vector& vec) { + std::string out = "{"; + + std::size_t len = vec.size(); + + std::size_t pause_start; + std::size_t pause_stop; + + if (len > 30){ + pause_start = 3; + pause_stop = len - 3; + } else { + pause_start = len *2; + pause_stop = pause_start; + } + + for (std::size_t i = 0; i < len; i++) { + if ((i > pause_start) && (i < pause_stop)) { continue; } + + if (i == pause_stop) { + out += ", ... "; + } else if (i != 0) { + out += ", "; + } + + char buf[30]; + sprintf(buf, "%g", vec[i]); + out += buf; + } + return out + "}"; +} + +inline void compare_values(const std::vector& actual, + const std::vector& desired, + double rtol = 0.0, double atol = 0.0, + std::string err_msg = "") +{ + if (actual.size() != desired.size()){ + error("the compared arrays have different lengths\n"); + } + + std::size_t num_mismatches = 0; + double max_absDiff = 0.0; + std::size_t max_absDiff_ind = 0; + double max_relDiff = 0.0; + std::size_t max_relDiff_ind = 0; + bool has_nan_mismatch = false; + + for (std::size_t i = 0; i < actual.size(); i++) { + double cur_absDiff = fabs(actual[i]-desired[i]); + + bool isnan_actual = isnan(actual[i]); + bool isnan_desired = isnan(desired[i]); + + if ( (cur_absDiff > (atol + rtol * fabs(desired[i]))) || + (isnan_actual != isnan_desired) ) { + + num_mismatches++; + if (isnan_actual != isnan_desired){ + has_nan_mismatch = true; + max_absDiff = NAN; + max_absDiff_ind = i; + max_relDiff = NAN; + max_relDiff_ind = i; + } else if (!has_nan_mismatch) { + if (cur_absDiff > max_absDiff){ + max_absDiff = cur_absDiff; + max_absDiff_ind = i; + } + + if ( cur_absDiff > (max_relDiff * fabs(desired[i])) ) { + max_relDiff = cur_absDiff / fabs(desired[i]); + max_relDiff_ind = i; + } + } + } + } + + if (num_mismatches == 0) { return; } + + std::string actual_vec_str = vec_to_string(actual); + std::string ref_vec_str = vec_to_string(desired); + + error + (("arrays are unequal for the tolerance: rtol = %g, atol = %g\n" + "%s\n" // custom error message + "Mismatched elements: %d / %d\n" + "Max absolute difference: %g, ind = %d, actual = %g, reference = %g\n" + "Max relative difference: %g, ind = %d, actual = %g, reference = %g\n" + "actual: %s\n" + "desired: %s\n"), + rtol, atol, + err_msg.c_str(), + (int)num_mismatches, (int)actual.size(), + max_absDiff, (int)max_absDiff_ind, actual[max_absDiff_ind], + desired[max_absDiff_ind], + max_relDiff, (int)max_relDiff_ind, actual[max_relDiff_ind], + desired[max_relDiff_ind], + actual_vec_str.c_str(), ref_vec_str.c_str()); +} + +/// This function takes an input cloudy-cooling data table and cuts it down so +/// that it is a 1D table (that only varies in terms of temperature) +inline void cut_down_to_1D_table(cloudy_data* ptr){ + if (ptr == nullptr) { + error("something is wrong"); + } else if (ptr->grid_rank != 2) { + error("only tested while cutting a 2D table down to 1D"); + } + + // - for a rank N table, ptr->grid_dimension[N-1] holds values corresponding + // to temperature + // - for a 2D table, ptr->grid_dimension[0] & ptr->grid_dimension[1] holds + // values that respectively corresponding to mass density & temperature + + // table over rho and Temperature + // -> I confirmed from the mmw table that temperature axis is contiguous + long long num_T_vals = ptr->grid_dimension[ptr->grid_rank - 1]; + + cloudy_data newObj; + for (long long i = 0LL; i < GRACKLE_CLOUDY_TABLE_MAX_DIMENSION; i++){ + newObj.grid_dimension[i] = 0LL; + newObj.grid_parameters[i] = NULL; + } + + // set rank, dimensions, & parameter vals + newObj.grid_rank = 1; + newObj.grid_dimension[0] = num_T_vals; + newObj.grid_parameters[0] = ptr->grid_parameters[ptr->grid_rank-1]; + + // since Temperature is the fast-index, we are going to just reuse the + // original pointers (even though they will have a bunch of unused data) + newObj.heating_data = ptr->heating_data; + newObj.cooling_data = ptr->cooling_data; + newObj.mmw_data = ptr->mmw_data; + newObj.data_size = num_T_vals; + + // clean up unused data originally stored within ptr + for (long long i = 0LL; i < ptr->grid_rank; i++){ + if ((i+1) == ptr->grid_rank) { + continue; // don't free the temperature values (we reuse them later) + } else { + GRACKLE_FREE(ptr->grid_parameters[0]); + } + } + + (*ptr) = newObj; +} + +struct DummyGrackleConfig{ + // the central purpose here is to hold grackle configuration options for + // tabulated solver + // + // do NOT try to use this as a general purpose C++ interface + + chemistry_data chem_data; + chemistry_data_storage chem_rates; + code_units units; + + DummyGrackleConfig(int n_tab_dims, double radiation_redshit, + bool UVbackground = false, + bool cmb_temperature_floor = true) { + // radiation_redshift is meaningless when n_tab_dims isn't 3 + + // setup units! + code_units my_units; + my_units.comoving_coordinates = 0; // 1 if cosmological sim, 0 if not + my_units.density_units = mh; // = mass_H/cm^3 (in g/cm^3) + my_units.length_units = 1.0 / 3.24077929e-25; // = 1 Mpc (in cm) + my_units.time_units = 3.15576e13; // = 1 Myr (in seconds) + my_units.a_units = 1.0; // units for the expansion factor + // Set expansion factor to 1 for non-cosmological simulation. + my_units.a_value = 1. / (1. + radiation_redshit) / my_units.a_units; + set_velocity_units(&my_units); + + // setup runtime parameters + chemistry_data my_chem; + local_initialize_chemistry_parameters(&my_chem); + my_chem.use_grackle = 1; // chemistry on + my_chem.with_radiative_cooling = 1; // cooling on + my_chem.primordial_chemistry = 0; + my_chem.metal_cooling = 1; // metal cooling on + my_chem.UVbackground = (UVbackground) ? 1 : 0; // UV background on + my_chem.cmb_temperature_floor = (cmb_temperature_floor) ? 1 : 0; + + if (UVbackground && (n_tab_dims != 3)) { + error("can't enable UVbackground without the redshift dependence"); + } + + if ((n_tab_dims <= 0) || (n_tab_dims > 3)) { + error("n_tab_dims must be 1, 2, or 3\n"); + } else if (n_tab_dims <= 2) { + my_chem.grackle_data_file + = "../../grackle_data_files/input/CloudyData_noUVB.h5"; + } else { + my_chem.grackle_data_file = + "../../grackle_data_files/input/CloudyData_UVB=HM2012.h5"; + } + + chemistry_data_storage my_rates; + local_initialize_chemistry_data(&my_chem, &my_rates, &my_units); + + if (n_tab_dims == 1) { + cut_down_to_1D_table(&(my_rates.cloudy_metal)); + cut_down_to_1D_table(&(my_rates.cloudy_primordial)); + //error("Can't currently support n_tab_dims == 1\n"); + } + + this->chem_data = my_chem; + this->chem_rates = my_rates; + this->units = my_units; + } + + ~DummyGrackleConfig() { + local_free_chemistry_data(&(this->chem_data), &(this->chem_rates)); + // don't deallocate (this->chem_data).grackle_data_file, that member is + // assigned a string-literal (it's not heap-allocated) + } + + DummyGrackleConfig(const DummyGrackleConfig&) = delete; + DummyGrackleConfig(DummyGrackleConfig&&) = delete; + DummyGrackleConfig& operator=(const DummyGrackleConfig&) = delete; + DummyGrackleConfig& operator=(DummyGrackleConfig&&) = delete; + +}; + + +/// This provides an interface for a very crude timer +/// -> the main reason this is a struct is so that we can easily swap out the +/// internal functionality +struct Timer { + Timer() { + active_ = false; + running_total_ = {0,0}; + } + + void start() { + if (active_) error("TIMER ALREADY ACTIVE!"); + active_ = true; + clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &last_start_); + } + + void stop() { + timespec stop; + clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &stop); + if(!active_) error("TIMER ALREADY STOPPED!"); + active_ = false; + + timespec elapsed; + if (last_start_.tv_nsec > stop.tv_nsec) { + stop.tv_sec -= time_t(1); + stop.tv_nsec += long(1.0e9); + } + elapsed.tv_sec = stop.tv_sec - last_start_.tv_sec; + elapsed.tv_nsec = stop.tv_nsec - last_start_.tv_nsec; + + running_total_.tv_sec += elapsed.tv_sec; + running_total_.tv_nsec += elapsed.tv_nsec; + if (running_total_.tv_nsec > long(1.0e9)) { + running_total_.tv_sec++; + running_total_.tv_nsec -= long(1.0e9); + } + } + + double get() const { + return running_total_.tv_sec + (running_total_.tv_nsec / 1.0e9); + } + +private: + bool active_; + timespec running_total_; + timespec last_start_; + +}; diff --git a/src/clib/interpolators_g.F b/src/example/utest_interpolators/interpolators_g.F similarity index 97% rename from src/clib/interpolators_g.F rename to src/example/utest_interpolators/interpolators_g.F index 3bf04ef4..998ce843 100644 --- a/src/clib/interpolators_g.F +++ b/src/example/utest_interpolators/interpolators_g.F @@ -6,7 +6,7 @@ subroutine interpolate_1D_g( & dataSize, dataField, value) implicit NONE -#include "grackle_fortran_types.def" +#include "../../clib/grackle_fortran_types.def" ! General Arguments @@ -49,7 +49,7 @@ subroutine interpolate_2D_g( & dataSize, dataField, value) implicit NONE -#include "grackle_fortran_types.def" +#include "../../clib/grackle_fortran_types.def" ! General Arguments @@ -111,7 +111,7 @@ subroutine interpolate_3D_g( & dataSize, dataField, value) implicit NONE -#include "grackle_fortran_types.def" +#include "../../clib/grackle_fortran_types.def" ! General Arguments @@ -192,7 +192,7 @@ subroutine interpolate_3Dz_g( & end_int, value) implicit NONE -#include "grackle_fortran_types.def" +#include "../../clib/grackle_fortran_types.def" ! General Arguments integer*8 dataSize, index2 @@ -202,10 +202,10 @@ subroutine interpolate_3Dz_g( & gridPar2(gridDim(2)), & gridPar3(gridDim(3)), dgridPar3 real*8 dataField(dataSize) + integer*8 end_int ! Locals - logical end_int integer*8 index1, index3, int_index integer q, w real*8 slope, value3(2), value2(2) @@ -213,7 +213,7 @@ subroutine interpolate_3Dz_g( !\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\///////////////////////////////// !======================================================================= - if (end_int) then + if (end_int .eq. 1) then call interpolate_2Df3D_g(input1, & input3, gridDim, & gridPar1, dgridPar1, @@ -285,7 +285,7 @@ subroutine interpolate_2Df3D_g( & value) implicit NONE -#include "grackle_fortran_types.def" +#include "../../clib/grackle_fortran_types.def" ! General Arguments @@ -350,7 +350,7 @@ subroutine interpolate_4D_g( & dataSize, dataField, value) implicit NONE -#include "grackle_fortran_types.def" +#include "../../clib/grackle_fortran_types.def" ! General Arguments @@ -444,7 +444,7 @@ subroutine interpolate_5D_g( & dataSize, dataField, value) implicit NONE -#include "grackle_fortran_types.def" +#include "../../clib/grackle_fortran_types.def" ! General Arguments diff --git a/src/example/utest_interpolators/test_interpolators.C b/src/example/utest_interpolators/test_interpolators.C new file mode 100644 index 00000000..ce1c03e2 --- /dev/null +++ b/src/example/utest_interpolators/test_interpolators.C @@ -0,0 +1,775 @@ +#include // std::minstd_rand +#include // std::int32_t +#include // std::printf + +#include + + +#include "../utest_helpers.hpp" + + +typedef long long gr_int64; + +extern "C" { + +#include "../../clib/interop/interop_funcs.h" + +void FORTRAN_NAME(interpolate_1d_g)( + const double* input1, const gr_int64* gridDim, + const double* gridPar1, const double* dgridPar1, + const gr_int64* dataSize, const double* dataField, + double* value); + +void FORTRAN_NAME(interpolate_2d_g)( + const double* input1, const double* input2, + const gr_int64* gridDim, + const double* gridPar1, const double* dgridPar1, + const double* gridPar2, const double* dgridPar2, + const gr_int64* dataSize, const double* dataField, + double* value); + +void FORTRAN_NAME(interpolate_3dz_g)( + const double* input1, const double* input2, const double* input3, + const gr_int64* gridDim, + const double* gridPar1, const double* dgridPar1, + const double* gridPar2, const gr_int64* index2, + const double* gridPar3, const double* dgridPar3, + const gr_int64* dataSize, const double* dataField, + const gr_int64* end_int, + double* value); + +void FORTRAN_NAME(interpolate_3d_g)( + const double* input1, const double* input2, const double* input3, + const gr_int64* gridDim, + const double* gridPar1, const double* dgridPar1, + const double* gridPar2, const double* dgridPar2, + const double* gridPar3, const double* dgridPar3, + const gr_int64* dataSize, const double* dataField, + double* value); + +void FORTRAN_NAME(interpolate_4d_g)( + const double* input1, const double* input2, const double* input3, + const double* input4, + const gr_int64* gridDim, + const double* gridPar1, const double* dgridPar1, + const double* gridPar2, const double* dgridPar2, + const double* gridPar3, const double* dgridPar3, + const double* gridPar4, const double* dgridPar4, + const gr_int64* dataSize, const double* dataField, + double* value); + +void FORTRAN_NAME(interpolate_5d_g)( + const double* input1, const double* input2, const double* input3, + const double* input4, const double* input5, + const gr_int64* gridDim, + const double* gridPar1, const double* dgridPar1, + const double* gridPar2, const double* dgridPar2, + const double* gridPar3, const double* dgridPar3, + const double* gridPar4, const double* dgridPar4, + const double* gridPar5, const double* dgridPar5, + const gr_int64* dataSize, const double* dataField, + double* value); + +} + + +/// The idea here is to encapsulate an DataTable that can be used to execute +/// one of the interpolation functions +class InterpTable { + +public: + // delete default constructor + InterpTable() = delete; + + InterpTable(std::vector> param_vals, + std::vector dataField) + { + if ((param_vals.size() == 0) || (param_vals.size() > 5)){ + error("InterpTable::InterpTable()", + "param_vals must have 1 to 5 entries"); + } + + std::vector gridDim; + for (std::size_t i = 0; i < param_vals.size(); i++) { + std::size_t cur_dim_size = param_vals[i].size(); + if (cur_dim_size < 2) { + error("InterpTable::InterpTable()", + "each parameter must have 2+ entries"); + } + gridDim.push_back(cur_dim_size); + } + + this->gridDim_ = gridDim; + this->param_vals_ = param_vals; + this->dataField_ = dataField; + + // sanity check: + if (this->dataSize() != dataField_.size()) { + error("InterpTable::InterpTable()", + "dataField has the wrong number of entries."); + } + } + + const gr_int64* gridDim() const + { return gridDim_.data(); } + + const double* gridField() const + { return dataField_.data(); } + + gr_int64 dataSize() const { + gr_int64 out = 1; + for (auto elem: this->gridDim_) { out *= elem; } + return out; + } + + int ndim() const + { return static_cast(gridDim_.size()); } + + const double* gridPar_0Indexed(int dim) const { + if (dim < 0) { + error("InterpTable::gridPar_0Indexed", "can't have negative dimension"); + } else if (dim >= ndim()) { + error("InterpTable::gridPar_0Indexed", "dim is too large"); + } + return param_vals_[dim].data(); + } + + double dgridPar_0Indexed(int dim) const { + const double* gridPar = gridPar_0Indexed(dim); + gr_int64 gridPar_len = this->gridDim_[dim]; + + return ((gridPar[gridPar_len - 1] - gridPar[0]) / + static_cast(gridPar_len - 1)); + } + +private: + // these shouldn't be mutated after construction + std::vector gridDim_; + // in principle param_vals_ should monotonically increase (and in almost all + // cases the spacing between values should be constant) + std::vector> param_vals_; + std::vector dataField_; +}; + + +Timer run_interp_helper(const double* val_vec, + const InterpTable& table, + std::size_t length, + bool use_fortran, double * out){ + + const int rank = table.ndim(); + + const double* gridPar1 = table.gridPar_0Indexed(0); + const double* gridPar2 = (rank > 1) ? table.gridPar_0Indexed(1) : nullptr; + const double* gridPar3 = (rank > 2) ? table.gridPar_0Indexed(2) : nullptr; + const double* gridPar4 = (rank > 3) ? table.gridPar_0Indexed(3) : nullptr; + const double* gridPar5 = (rank > 4) ? table.gridPar_0Indexed(4) : nullptr; + + const double dgridPar1 = table.dgridPar_0Indexed(0); + const double dgridPar2 = (rank > 1) ? table.dgridPar_0Indexed(1) : 0.0; + const double dgridPar3 = (rank > 2) ? table.dgridPar_0Indexed(2) : 0.0; + const double dgridPar4 = (rank > 3) ? table.dgridPar_0Indexed(3) : 0.0; + const double dgridPar5 = (rank > 4) ? table.dgridPar_0Indexed(4) : 0.0; + + const gr_int64* gridDim = table.gridDim(); + const gr_int64 dataSize = table.dataSize(); + const double* dataField = table.gridField(); + + Timer t; + + if ((rank == 1) && (use_fortran)) { + + t.start(); + for (std::size_t i = 0; i < length; i++) { + FORTRAN_NAME(interpolate_1d_g)(&val_vec[i], gridDim, + gridPar1, &dgridPar1, + &dataSize, dataField, + out+i); + } + t.stop(); + + } else if (rank == 1){ + + t.start(); + for (std::size_t i = 0; i < length; i++) { + interpolate_1d_g(val_vec[i], gridDim, + gridPar1, dgridPar1, + dataSize, dataField, + out+i); + } + t.stop(); + + } else if ((rank == 2) && use_fortran){ + + t.start(); + for (std::size_t i = 0; i < length; i++) { + FORTRAN_NAME(interpolate_2d_g)(&val_vec[i*2], &val_vec[i*2+1], + gridDim, + gridPar1, &dgridPar1, + gridPar2, &dgridPar2, + &dataSize, dataField, + out+i); + } + t.stop(); + + } else if ((rank == 2)){ + + t.start(); + for (std::size_t i = 0; i < length; i++) { + interpolate_2d_g(val_vec[i*2], val_vec[i*2+1], + gridDim, + gridPar1, dgridPar1, + gridPar2, dgridPar2, + dataSize, dataField, + out+i); + } + t.stop(); + + } else if ((rank == 3) && use_fortran){ + + t.start(); + for (std::size_t i = 0; i < length; i++) { + FORTRAN_NAME(interpolate_3d_g)(&val_vec[i*3], &val_vec[i*3+1], + &val_vec[i*3+2], + gridDim, + gridPar1, &dgridPar1, + gridPar2, &dgridPar2, + gridPar3, &dgridPar3, + &dataSize, dataField, + out+i); + } + t.stop(); + + } else if (rank == 3) { + + t.start(); + for (std::size_t i = 0; i < length; i++) { + interpolate_3d_g(val_vec[i*3], val_vec[i*3+1], val_vec[i*3+2], + gridDim, + gridPar1, dgridPar1, + gridPar2, dgridPar2, + gridPar3, dgridPar3, + dataSize, dataField, + out+i); + } + t.stop(); + + } else if ((rank == 4) && use_fortran){ + + t.start(); + for (std::size_t i = 0; i < length; i++) { + FORTRAN_NAME(interpolate_4d_g)(&val_vec[i*4], &val_vec[i*4+1], + &val_vec[i*4+2], &val_vec[i*4+3], + gridDim, + gridPar1, &dgridPar1, + gridPar2, &dgridPar2, + gridPar3, &dgridPar3, + gridPar4, &dgridPar4, + &dataSize, dataField, + out+i); + } + t.stop(); + + } else if (rank == 4) { + + t.start(); + for (std::size_t i = 0; i < length; i++) { + interpolate_4d_g(val_vec[i*4], val_vec[i*4+1], val_vec[i*4+2], + val_vec[i*4+3], + gridDim, + gridPar1, dgridPar1, + gridPar2, dgridPar2, + gridPar3, dgridPar3, + gridPar4, dgridPar4, + dataSize, dataField, + out+i); + } + t.stop(); + + } else if (use_fortran) { + + t.start(); + for (std::size_t i = 0; i < length; i++) { + FORTRAN_NAME(interpolate_5d_g)(&val_vec[i*5], &val_vec[i*5+1], + &val_vec[i*5+2], &val_vec[i*5+3], + &val_vec[i*5+4], + gridDim, + gridPar1, &dgridPar1, + gridPar2, &dgridPar2, + gridPar3, &dgridPar3, + gridPar4, &dgridPar4, + gridPar5, &dgridPar5, + &dataSize, dataField, + out+i); + } + t.stop(); + + } else { + + t.start(); + for (std::size_t i = 0; i < length; i++) { + interpolate_5d_g(val_vec[i*5], val_vec[i*5+1], val_vec[i*5+2], + val_vec[i*5+3], val_vec[i*5+4], + gridDim, + gridPar1, dgridPar1, + gridPar2, dgridPar2, + gridPar3, dgridPar3, + gridPar4, dgridPar4, + gridPar5, dgridPar5, + dataSize, dataField, + out+i); + } + t.stop(); + + } + return t; +} + + +struct interp_3dz_arg_{ + double par1; + double zr; + double par3; + gr_int64 zindex; + gr_int64 end_int; +}; + +Timer run_interp_3dz_(const double* vals, + const InterpTable& table, + std::size_t length, + bool use_fortran, + double * out) +{ + const int rank = table.ndim(); + if (rank != 3) error("run_interp_3dz", "requires a rank 3 table"); + + const double* gridPar1 = table.gridPar_0Indexed(0); + const double* gridPar2 = table.gridPar_0Indexed(1); + const double* gridPar3 = table.gridPar_0Indexed(2); + + const double dgridPar1 = table.dgridPar_0Indexed(0); + const double dgridPar3 = table.dgridPar_0Indexed(2); + + const gr_int64* gridDim = table.gridDim(); + const gr_int64 dataSize = table.dataSize(); + const double* dataField = table.gridField(); + + // need to do some conversions on the arguments (we try to do this ahead of + // time to minimize impact on the timings). Ideally, we'd do this before + // calling this function so that we could reuse the same array for both + // versions of the function + std::vector inputs_vec(length); + for (std::size_t i = 0; i < length; i++) { + const double zr = vals[i*3+1]; // z redshift + + // Calculate index for redshift dimension - intentionally kept 1-indexed + const long long zindex = find_zindex(zr, rank, gridDim, gridPar2); + const gr_int64 end_int = ((rank > 2) && (zindex == gridDim[1])); + + inputs_vec[i] = {vals[i*3], zr, vals[i*3+2], zindex, end_int}; + } + + const interp_3dz_arg_ * inputs = inputs_vec.data(); + + Timer timer; + + if (use_fortran) { + + timer.start(); + for (std::size_t i = 0; i < length; i++) { + FORTRAN_NAME(interpolate_3dz_g)(&(inputs[i].par1), &(inputs[i].zr), + &(inputs[i].par3), + gridDim, + gridPar1, &dgridPar1, + gridPar2, &(inputs[i].zindex), + gridPar3, &dgridPar3, + &dataSize, dataField, + &(inputs[i].end_int), out+i); + } + timer.stop(); + + } else { + + timer.start(); + for (std::size_t i = 0; i < length; i++) { + interpolate_3dz_g(inputs[i].par1, inputs[i].zr, inputs[i].par3, + gridDim, + gridPar1, dgridPar1, + gridPar2, inputs[i].zindex, + gridPar3, dgridPar3, + dataSize, dataField, + inputs[i].end_int, out+i); + } + timer.stop(); + + } + + return timer; +} + +std::vector run_interp(const std::vector>& vals, + const InterpTable& table, bool use_fortran, + bool use_3dz, double * elapsed = nullptr) +{ + const std::size_t ndim = table.ndim(); + std::vector flat_vec(vals.size() * ndim); + for (std::size_t i = 0; i < vals.size(); i++) { + for (std::size_t j = 0; j < ndim; j++) { + flat_vec[i*ndim+j] = vals[i][j]; + } + } + + std::vector out(vals.size()); + Timer t; + if (use_3dz) { + t = run_interp_3dz_(flat_vec.data(), table, vals.size(), + use_fortran, out.data()); + } else { + t = run_interp_helper(flat_vec.data(), table, vals.size(), + use_fortran, out.data()); + } + + if (elapsed != nullptr) *elapsed = t.get(); + return out; +} + +void compare_interp_impls_(const std::vector> &inputs, + const InterpTable& table, bool use_3dz, + bool show_timing = false) +{ + std::vector ref, actual; + double fortran_duration_sec, c_duration_sec; + + ref = run_interp(inputs, table, true, use_3dz); // run an extra-time for + // more-fair timing + ref = run_interp(inputs, table, true, use_3dz, &fortran_duration_sec); + actual = run_interp(inputs, table, false, use_3dz, &c_duration_sec); + if (show_timing) { + printf(" ...Timing Compare: Fortran = %g sec C = %g sec\n", + fortran_duration_sec, c_duration_sec); + } + + if (false) { // debugging statement: + std::string ref_s = vec_to_string(ref); + std::string actual_s = vec_to_string(actual); + std::printf(" reference: %s\n", ref_s.c_str()); + std::printf(" actual: %s\n", actual_s.c_str()); + } + compare_values(actual, ref, 0.0, 0.0, ""); +}; + +// Miscelaneous Functions +// ====================== + +// uniform distribution on the interval (0, 1] +double uniform_dist_transform_(std::minstd_rand &prng){ + + // sanity check to confirm that the largest value returned by generator is + // representable (without any loss of precision) by a double (its <= 2^53) + static_assert( (prng.max() <= 9007199254740992) & (prng.min() == 1), + "sanity check failed"); + + return ( static_cast(prng()) / static_cast(prng.max()) ); +} + +// Functions for constructing InterpTable +// ====================================== + +struct ParamProps { double min_val; double max_val; int num_vals; }; + +InterpTable build_table(std::uint32_t seed, + const std::vector& paramprop_vec) +{ + std::minstd_rand prng = std::minstd_rand(seed); + std::vector> param_vals; + std::size_t field_size = 1; + + for (const ParamProps& cur_param_prop : paramprop_vec) { + const int dim_size = cur_param_prop.num_vals; + + // initialize cur_param_vals + std::vector cur_param_vals(dim_size); + double dt = (cur_param_prop.max_val - cur_param_prop.min_val) / dim_size; + for (int i = 0; i < dim_size; i++){ + cur_param_vals[i] = cur_param_prop.min_val + i * dt; + } + cur_param_vals[dim_size - 1] = cur_param_prop.max_val; + param_vals.push_back(cur_param_vals); + + // update field_size + field_size *= dim_size; + } + + // now initialize field_data + std::vector field_data(field_size); + for (int i = 0; i < field_size; i++) { + field_data[i] = 100.0 * uniform_dist_transform_(prng) - 50.0; + } + + return InterpTable(param_vals, field_data); +} + +InterpTable get_3dz_table() +{ + DummyGrackleConfig config(3, 0.0, true, true); + + cloudy_data& cloud_data_obj = config.chem_rates.cloudy_primordial; + + const int rank = (int)cloud_data_obj.grid_rank; + + // setup param_vals + std::vector> param_vals; + for (int i = 0; i < rank; i++) { + const int len = (int)(cloud_data_obj.grid_dimension[i]); + std::vector cur(len); + for (int j = 0; j < len; j++) { + cur[j] = cloud_data_obj.grid_parameters[i][j]; + } + param_vals.push_back(cur); + } + + // setup dataField + const long long data_size = (long long)cloud_data_obj.data_size; + std::vector dataField(data_size); + for (long long i = 0; i < data_size; i++) { + dataField[i] = cloud_data_obj.heating_data[i]; + } + + return {param_vals, dataField}; +} + +// Functions to help generate locations to perform interpolations at +// ================================================================= + +// this returns a vector of vectors appropriate interpolate function for cases where 1 or more +// inputs come from special_vals (the other inputs are taken from +// ordinary_vals) +// +// EXAMPLES +// -------- +// if special_vals = {1.0,2.0} & ordinary_vals = {-1.0, -2.0}, this passes the +// following to run_interp in separate calls: +// { 1.0, 2.0}, {-1.0, 2.0}, { 1.0,-2.0} +// if special_vals = {1.0,2.0,3.0} & ordinary_vals = {-1.0, -2.0, -3.0}, this +// passes the following to run_interp in separate calls: +// { 1.0, 2.0, 3.0}, {-1.0, 2.0, 3.0}, { 1.0,-2.0, 3.0}, {-1.0,-2.0, 3.0} +// { 1.0, 2.0,-3.0}, {-1.0, 2.0,-3.0}, { 1.0,-2.0,-3.0} +// +// The equivalent in python would be the result of the following command: +// >>> list(itertools.product(*zip(special_vals, ordinary_vals)))[:-1] +std::vector> get_combinations_ +(const std::vector& special_vals, + const std::vector& ordinary_vals) +{ + std::size_t rank = special_vals.size(); + + int choice[5] = {0, 0, 0, 0, 0}; + + int out_count = pow(2, rank) - 1; + std::vector> out; + out.reserve(out_count - 1); + + while (out.size() < out_count) { + std::vector cur(rank); + for (int i = 0; i < rank; i++){ + cur[i] = (choice[i] == 0) ? special_vals[i] : ordinary_vals[i]; + } + out.push_back(cur); + + // now increment choice + for (int i = 0; i < rank; i++){ + if (choice[i] == 0) { choice[i] = 1; break; } + choice[i] = 0; + } + } + + return out; +} + +struct TableSummaryProps{ + std::vector min_param_val; + std::vector max_param_val; + std::vector unaligned_inrange; + std::vector unaligned_inrange_alt; + std::vector half_min_param_val; + std::vector double_max_param_val; +}; + +TableSummaryProps get_interp_table_summary(const InterpTable& table) { + const int rank = table.ndim(); + + // get useful values for test: + std::vector min_param_val(rank); + std::vector max_param_val(rank); + std::vector unaligned_inrange(rank); + std::vector unaligned_inrange_alt(rank); + std::vector half_min_param_val(rank); + std::vector double_max_param_val(rank); + for (int i = 0; i < rank; i++){ + const double* cur_param_vals = table.gridPar_0Indexed(i); + const int cur_param_len = table.gridDim()[i]; + + min_param_val[i] = cur_param_vals[0]; + max_param_val[i] = cur_param_vals[cur_param_len-1]; + + double delta = cur_param_vals[1] - cur_param_vals[0]; + + unaligned_inrange[i] = cur_param_vals[0] + delta * 0.3; + unaligned_inrange_alt[i] = cur_param_vals[0] + delta * 0.5; + + half_min_param_val[i] = min_param_val[i]/2; + double_max_param_val[i] = max_param_val[i]*2; + } + + TableSummaryProps out; + out.min_param_val = min_param_val; + out.max_param_val = max_param_val; + out.unaligned_inrange = unaligned_inrange; + out.unaligned_inrange_alt = unaligned_inrange_alt; + out.half_min_param_val = half_min_param_val; + out.double_max_param_val = double_max_param_val; + return out; +} + +// Define the actual tests +// ======================= + +void run_test(const InterpTable& table, const bool use_3dz = false) +{ + TableSummaryProps props = get_interp_table_summary(table); + + // define the actual tests in the test suite + + std::printf(" -> comparing some in-range values\n"); + { + std::vector> inputs = get_combinations_ + (props.unaligned_inrange, props.unaligned_inrange_alt); + compare_interp_impls_(inputs, table, use_3dz); + } + + std::printf(" -> comparing cases with 1+ inputs on left grid edge\n"); + { + std::vector> inputs = get_combinations_ + (props.min_param_val, props.unaligned_inrange); + compare_interp_impls_(inputs, table, use_3dz); + } + + std::printf(" -> comparing cases with 1+ inputs on right grid edge\n"); + { + std::vector> inputs = get_combinations_ + (props.max_param_val, props.unaligned_inrange); + compare_interp_impls_(inputs, table, use_3dz); + } + + if (!use_3dz) { + std::printf(" -> comparing cases with 1+ inputs less than min grid val\n"); + std::vector> inputs = get_combinations_ + (props.half_min_param_val, props.unaligned_inrange); + compare_interp_impls_(inputs, table, use_3dz); + } + + std::printf(" -> comparing cases with 1+ inputs greater than max grid val\n"); + { + std::vector> inputs = get_combinations_ + (props.double_max_param_val, props.unaligned_inrange); + compare_interp_impls_(inputs, table, use_3dz); + } + + if (use_3dz) { + std::printf(" -> comparing cases with some extra redshifts\n"); + + gr_int64 n_gridded_z_vals = table.gridDim()[1]; + + double temp[5] = {table.gridPar_0Indexed(1)[n_gridded_z_vals - 4], + table.gridPar_0Indexed(1)[n_gridded_z_vals - 3], + table.gridPar_0Indexed(1)[n_gridded_z_vals - 2], + table.gridPar_0Indexed(1)[n_gridded_z_vals - 1], + table.gridPar_0Indexed(1)[n_gridded_z_vals - 1] * 1.01}; + + std::vector> inputs; + for (int i = 0; i < 5; i++) { + inputs.push_back({props.unaligned_inrange[0], temp[i], + props.unaligned_inrange[2]}); + } + compare_interp_impls_(inputs, table, use_3dz); + } + + std::size_t n_random = 1000000; + std::printf(" -> comparing %d random points within grid (with timing)\n", + int(n_random)); + { + const int rank = table.ndim(); + // the fact that we use vectors of vectors as inputs makes this inefficient + std::minstd_rand prng = std::minstd_rand(353545); + std::vector> inputs(n_random); + for (std::size_t i = 0; i < n_random; i++) { + std::vector tmp(rank); + for (int j = 0; j < rank; j++) { + double max_val = props.max_param_val[j]; + double min_val = props.min_param_val[j]; + tmp[j] = uniform_dist_transform_(prng) * (max_val - min_val) + min_val; + } + inputs[i] = tmp; + } + + compare_interp_impls_(inputs, table, use_3dz, true); + } + +} + + +int main(){ + + const std::uint32_t seed = 342; + + + std::printf("comparing interpolate_1d_g:\n"); + { + InterpTable table = build_table(seed, {{-6.0, 6.0, 25}}); + run_test(table); + } + + + std::printf("\ncomparing interpolate_2d_g:\n"); + { + InterpTable table = build_table(seed, {{-6.0, 6.0, 25}, {0.0, 10.0, 11}}); + run_test(table); + } + + + std::printf("\ncomparing interpolate_3d_g:\n"); + { + InterpTable table = build_table(seed, + {{-6.0, 6.0, 25}, {0.0, 10.0, 11}, + {-1.0, 0.0, 5}}); + run_test(table); + } + + + std::printf("\ncomparing interpolate_4d_g:\n"); + { + InterpTable table = build_table(seed, + {{-6.0, 6.0, 25}, {0.0, 10.0, 11}, + {-1.0, 0.0, 5}, {0.5, 5.5, 6}}); + run_test(table); + } + + + std::printf("\ncomparing interpolate_5d_g:\n"); + { + InterpTable table = build_table(seed, + {{-6.0, 6.0, 25}, {0.0, 10.0, 11}, + {-1.0, 0.0, 5}, {0.5, 5.5, 6}, + {-10.0, 0.0, 11}}); + run_test(table); + } + + + std::printf("\ncomparing interpolate_3dz_g:\n"); + { + InterpTable table = get_3dz_table(); + run_test(table, true); + } + + return 0; +} diff --git a/src/python/tests/test_code_examples.py b/src/python/tests/test_code_examples.py index ad6c3ce4..ff2f7c99 100644 --- a/src/python/tests/test_code_examples.py +++ b/src/python/tests/test_code_examples.py @@ -25,7 +25,10 @@ "cxx_example", "cxx_omp_example", "cxx_grid_example", - "fortran_example"] + "fortran_example", + "test_calc_temp1d_cloudy", + "test_cool1d_cloudy", + "test_interpolators"] def run_command(command, cwd, env):