Port calc_temp1d_cloudy_g & cool1d_cloudy_g from Fortran to C

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
+        utils.lo \
+        interop/shims_g.lo \
+        interop/calc_temp1d_cloudy_g.lo \
+        interop/cool1d_cloudy_g.lo \
+        interop/interpolators_g.lo

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
 
 #-----------------------------------------------------------------------

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

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
 

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,

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) {				\

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 [<stdbool.h>](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)

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 <stdlib.h>
+#include <stdint.h> // int32_t
+#include <stdio.h>
+#include <math.h> // 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;
+
+  }
+}