open-atmos · mattldawson · Sep 3, 2024 · Aug 31, 2024 · Sep 3, 2024
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -27,7 +27,6 @@ endif ()
 ######################################################################
 # options
 
-option(ENABLE_GSL "Enable GSL library for Jacobian evaluation" OFF)
 option(ENABLE_MPI "Enable MPI parallel support" OFF)
 option(ENABLE_DEBUG "Compile debugging functions" OFF)
 option(FAILURE_DETAIL "Output conditions before and after solver failures" OFF)
@@ -45,27 +44,6 @@ set(CPACK_SOURCE_IGNORE_FILES "${CPACK_SOURCE_IGNORE_FILES};/.*~$;/[.].*/;/build
 
 include(CPack)
 
-######################################################################
-# GSL
-
-if(ENABLE_GSL)
-  find_path(GSL_INCLUDE_DIR gsl/gsl_math.h
-    DOC "GSL include directory (must have gsl/ subdir)"
-    PATHS $ENV{GSL_HOME}/include /opt/local/include)
-  find_library(GSL_LIB gsl
-    DOC "GSL library"
-    PATHS $ENV{GSL_HOME}/lib /opt/local/lib)
-  find_library(GSL_CBLAS_LIB gslcblas
-    DOC "GSL CBLAS library"
-    PATHS $ENV{GSL_HOME}/lib /opt/local/lib)
-  find_library(M_LIB m
-    DOC "standard C math library")
-  set(GSL_SRC src/rand_gsl.c)
-  set(GSL_LIBS ${GSL_LIB} ${GSL_CBLAS_LIB} ${M_LIB})
-  include_directories(${GSL_INCLUDE_DIR})
-  add_definitions(-DCAMP_USE_GSL)
-endif()
-
 ######################################################################
 # MPI
 
@@ -381,13 +359,13 @@ set(CAMP_LIB_SRC
   src/solver_stats.F90
   src/debug_diff_check.F90
   ${CAMP_C_SRC} ${AEROSOL_REPS_SRC} ${SUB_MODELS_SRC} ${REACTIONS_SRC}
-  ${CAMP_CUDA_SRC} ${GSL_SRC} ${CAMP_CXX_SRC} )
+  ${CAMP_CUDA_SRC} ${CAMP_CXX_SRC} )
 
 add_library(camplib SHARED ${CAMP_LIB_SRC})
 add_library(camplib-static STATIC ${CAMP_LIB_SRC})
 
-target_link_libraries(camplib ${SUNDIALS_LIBS} ${GSL_LIBS} ${JSON_LIB})
-target_link_libraries(camplib-static ${SUNDIALS_LIBS} ${GSL_LIBS} ${JSON_LIB})
+target_link_libraries(camplib ${SUNDIALS_LIBS} ${JSON_LIB})
+target_link_libraries(camplib-static ${SUNDIALS_LIBS} ${JSON_LIB})
 
 set(MODULE_DIR "${CMAKE_BINARY_DIR}/include")
 

diff --git a/Dockerfile b/Dockerfile
@@ -4,7 +4,7 @@ RUN dnf -y update \
     && dnf -y install \
         gcc-gfortran \
         gcc-c++ \
-        gsl-devel \
+        make \
         metis-devel \
         lapack-devel \
         openblas-devel \
@@ -63,8 +63,7 @@ ENV PATH="${PATH}:/usr/local/jsonfortran-gnu-6.1.0/lib"
              -D CMAKE_C_FLAGS_DEBUG="-pg" \
              -D CMAKE_Fortran_FLAGS_DEBUG="-pg" \
              -D CMAKE_MODULE_LINKER_FLAGS="-pg" \
-             -D ENABLE_GSL:BOOL=TRUE \
              /camp \
     && make install
 
-WORKDIR /build
+WORKDIR /build
diff --git a/Dockerfile.debug b/Dockerfile.debug
@@ -4,7 +4,7 @@ RUN dnf -y update \
     && dnf -y install \
         gcc-gfortran \
         gcc-c++ \
-        gsl-devel \
+        make \
         metis-devel \
         lapack-devel \
         openblas-devel \
@@ -62,7 +62,6 @@ ENV PATH="${PATH}:/usr/local/jsonfortran-gnu-6.1.0/lib"
              -D CMAKE_C_FLAGS="-pg -Og" \
              -D CMAKE_Fortran_FLAGS="-pg -Og -fbacktrace" \
              -D CMAKE_MODULE_LINKER_FLAGS="-pg" \
-             -D ENABLE_GSL:BOOL=TRUE \
              -D ENABLE_DEBUG:BOOL=TRUE \
              /camp \
     && make install
diff --git a/Dockerfile.mpi b/Dockerfile.mpi
@@ -13,7 +13,7 @@ RUN sudo dnf -y install \
         openmpi-devel \
         gcc-gfortran \
         gcc-c++ \
-        gsl-devel \
+        make \
         metis-devel \
         lapack-devel \
         openblas-devel \
@@ -79,4 +79,4 @@ COPY . /home/test_user/camp/
     && make \
     && sudo make install
 
-WORKDIR /home/test_user/build
+WORKDIR /home/test_user/build
diff --git a/README.md b/README.md
@@ -51,12 +51,6 @@ Required dependencies:
    * SuiteSparse - <http://faculty.cse.tamu.edu/davis/SuiteSparse/SuiteSparse-5.1.0.tar.gz>
    * A modified version of CVODE 3.1.2 - provided here in `cvode-3.4-alpha.tar.gz`
 
-Optional dependencies:
-
-   * GSL for Jacobian evaluation and random number generators -
-     <http://www.gnu.org/software/gsl/>
-
-
 Installation
 ============
 
@@ -84,7 +78,6 @@ Installation
         export JSON_FORTRAN_HOME=${HOME}/opt
         export SUITE_SPARSE_HOME=${HOME}/opt
         export SUNDIALS_HOME=${HOME}/opt
-        export GSL_HOME=${HOME}/opt
 
    Of course the exact directories will depend on where the libraries
    are installed. You only need to set variables for libraries

diff --git a/src/camp_core.F90 b/src/camp_core.F90
@@ -41,7 +41,6 @@
 !! |--------------|---------|-----------------------------------------------|
 !! | SUNDIALS     | custom  | camp/cvode-3.4-alpha.tar.gz                   |
 !! | SuiteSparse  | 5.1.0   | http://faculty.cse.tamu.edu/davis/SuiteSparse/SuiteSparse-5.1.0.tar.gz |
-!! | GSL          |         | https://www.gnu.org/software/gsl/             |
 !! | json-fortran | 6.1.0   | https://github.com/jacobwilliams/json-fortran/archive/6.1.0.tar.gz |
 !!
 !! The SUNDIALS library must be built with the `ENABLE_KLU` flag set to `ON`

diff --git a/src/camp_solver.c b/src/camp_solver.c
@@ -23,23 +23,13 @@
 #ifdef CAMP_USE_GPU
 #include "cuda/camp_gpu_solver.h"
 #endif
-#ifdef CAMP_USE_GSL
-#include <gsl/gsl_deriv.h>
-#include <gsl/gsl_math.h>
-#include <gsl/gsl_roots.h>
-#endif
 #include "camp_debug.h"
 
 // Default solver initial time step relative to total integration time
 #define DEFAULT_TIME_STEP 1.0
 // State advancement factor for Jacobian element evaluation
 #define JAC_CHECK_ADV_MAX 1.0E-00
 #define JAC_CHECK_ADV_MIN 1.0E-12
-// Relative tolerance for Jacobian element evaluation against GSL absolute
-// errors
-#define JAC_CHECK_GSL_REL_TOL 1.0e-4
-// Absolute Jacobian error tolerance
-#define JAC_CHECK_GSL_ABS_TOL 1.0e-9
 // Set MAX_TIMESTEP_WARNINGS to a negative number to prevent output
 #define MAX_TIMESTEP_WARNINGS -1
 // Maximum number of steps in discreet addition guess helper
@@ -1149,144 +1139,15 @@ bool check_Jac(realtype t, N_Vector y, SUNMatrix J, N_Vector deriv,
   realtype *d_deriv = NV_DATA_S(deriv);
   bool retval = true;
 
-#ifdef CAMP_USE_GSL
-  GSLParam gsl_param;
-  gsl_function gsl_func;
-
-  // Set up gsl parameters needed during numerical differentiation
-  gsl_param.t = t;
-  gsl_param.y = tmp;
-  gsl_param.deriv = tmp1;
-  gsl_param.solver_data = (SolverData *)solver_data;
-
-  // Set up the gsl function
-  gsl_func.function = &gsl_f;
-  gsl_func.params = &gsl_param;
-#endif
-
   // Calculate the the derivative for the current state y
   if (f(t, y, deriv, solver_data) != 0) {
     printf("\n Derivative calculation failed.\n");
     return false;
   }
 
-  // Loop through the independent variables, numerically calculating
-  // the partial derivatives d_fy/d_x
-  for (int i_ind = 0; i_ind < NV_LENGTH_S(y); ++i_ind) {
-    // If GSL is available, use their numerical differentiation to
-    // calculate the partial derivatives. Otherwise, estimate them by
-    // advancing the state.
-#ifdef CAMP_USE_GSL
-
-    // Reset tmp to the initial state
-    N_VScale(ONE, y, tmp);
-
-    // Save the independent species concentration and index
-    double x = d_state[i_ind];
-    gsl_param.ind_var = i_ind;
-
-    // Skip small concentrations
-    if (x < SMALL) continue;
-
-    // Do the numerical differentiation for each potentially non-zero
-    // Jacobian element
-    for (int i_elem = SM_INDEXPTRS_S(J)[i_ind];
-         i_elem < SM_INDEXPTRS_S(J)[i_ind + 1]; ++i_elem) {
-      int i_dep = SM_INDEXVALS_S(J)[i_elem];
-
-      double abs_err;
-      double partial_deriv;
-
-      gsl_param.dep_var = i_dep;
-
-      bool test_pass = false;
-      double h, abs_tol, rel_diff, scaling;
-
-      // Evaluate the Jacobian element over a range of initial step sizes
-      for (scaling = JAC_CHECK_ADV_MIN;
-           scaling <= JAC_CHECK_ADV_MAX && test_pass == false;
-           scaling *= 10.0) {
-        // Get the current initial step size
-        h = x * scaling;
-
-        // Get the partial derivative d_fy/dx
-        if (gsl_deriv_forward(&gsl_func, x, h, &partial_deriv, &abs_err) == 1) {
-          printf("\nERROR in numerical differentiation for J[%d][%d]", i_ind,
-                 i_dep);
-        }
-
-        // Evaluate the results
-        abs_tol = 1.2 * fabs(abs_err);
-        abs_tol =
-            abs_tol > JAC_CHECK_GSL_ABS_TOL ? abs_tol : JAC_CHECK_GSL_ABS_TOL;
-        rel_diff = 1.0;
-        if (partial_deriv != 0.0)
-          rel_diff =
-              fabs((SM_DATA_S(J)[i_elem] - partial_deriv) / partial_deriv);
-        if (fabs(SM_DATA_S(J)[i_elem] - partial_deriv) < abs_tol ||
-            rel_diff < JAC_CHECK_GSL_REL_TOL)
-          test_pass = true;
-      }
-
-      // If the test does not pass with any initial step size, print out the
-      // failure, output the local derivative state and return false
-      if (test_pass == false) {
-        printf(
-            "\nError in Jacobian[%d][%d]: Got %le; expected %le"
-            "\n  difference %le is greater than error %le",
-            i_ind, i_dep, SM_DATA_S(J)[i_elem], partial_deriv,
-            fabs(SM_DATA_S(J)[i_elem] - partial_deriv), abs_tol);
-        printf("\n  relative error %le intial step size %le", rel_diff, h);
-        printf("\n  initial rate %le initial state %le", d_deriv[i_dep],
-               d_state[i_ind]);
-        printf(" scaling %le", scaling);
-        ModelData *md = &(((SolverData *)solver_data)->model_data);
-        for (int i_cell = 0; i_cell < md->n_cells; ++i_cell)
-          for (int i_spec = 0; i_spec < md->n_per_cell_state_var; ++i_spec)
-            printf("\n cell: %d species %d state_id %d conc: %le", i_cell,
-                   i_spec, i_cell * md->n_per_cell_state_var + i_spec,
-                   md->total_state[i_cell * md->n_per_cell_state_var + i_spec]);
-        retval = false;
-        output_deriv_local_state(t, y, deriv, solver_data, &f, i_dep, i_ind,
-                                 SM_DATA_S(J)[i_elem], h / 10.0);
-      }
-    }
-#endif
-  }
   return retval;
 }
 
-#ifdef CAMP_USE_GSL
-/** \brief Wrapper function for derivative calculations for numerical solving
- *
- * Wraps the f(t,y) function for use by the GSL numerical differentiation
- * functions.
- *
- * \param x Independent variable \f$x\f$ for calculations of \f$df_y/dx\f$
- * \param param Differentiation parameters
- * \return Partial derivative \f$df_y/dx\f$
- */
-double gsl_f(double x, void *param) {
-  GSLParam *gsl_param = (GSLParam *)param;
-  N_Vector y = gsl_param->y;
-  N_Vector deriv = gsl_param->deriv;
-
-  // Set the independent variable
-  NV_DATA_S(y)[gsl_param->ind_var] = x;
-
-  // Calculate the derivative
-  if (f(gsl_param->t, y, deriv, (void *)gsl_param->solver_data) != 0) {
-    printf("\nDerivative calculation failed!");
-    for (int i_spec = 0; i_spec < NV_LENGTH_S(y); ++i_spec)
-      printf("\n species %d conc: %le", i_spec, NV_DATA_S(y)[i_spec]);
-    return 0.0 / 0.0;
-  }
-
-  // Return the calculated derivative for the dependent variable
-  return NV_DATA_S(deriv)[gsl_param->dep_var];
-}
-#endif
-
 /** \brief Try to improve guesses of y sent to the linear solver
  *
  * This function checks if there are any negative guessed concentrations,

diff --git a/src/camp_solver.h b/src/camp_solver.h
@@ -68,17 +68,6 @@ static void print_jacobian(SUNMatrix M);
 static void print_derivative(N_Vector deriv);
 bool is_anything_going_on_here(SolverData *sd, realtype t_initial,
                                realtype t_final);
-#ifdef CAMP_USE_GSL
-double gsl_f(double x, void *param);
-typedef struct {
-  int ind_var;              // independent variable index
-  int dep_var;              // dependent variable index
-  realtype t;               // current model time (s)
-  N_Vector y;               // dependent variable array
-  N_Vector deriv;           // time derivative vector f(t,y)
-  SolverData *solver_data;  // solver data
-} GSLParam;
-#endif
 #e
4951
ndif
 
 #endif