Gretl submodule

.gitmodules

@@ -22,3 +22,6 @@
 [submodule "ContinuationSolvers"]
 	path = ContinuationSolvers
 	url = https://github.com/LLNL/ContinuationSolvers
+[submodule "gretl"]
+	path = gretl
+	url = https://github.com/LLNL/gretl

CMakeLists.txt

@@ -198,7 +198,8 @@ set(exported_targets
     smith_infrastructure
     smith_mesh_utils
     smith_numerics
-    smith_physics)
+    smith_physics
+    smith_differentiable_numerics)
 
 add_library(smith INTERFACE)
 target_link_libraries(smith INTERFACE ${exported_targets})
@@ -220,7 +221,7 @@ if (SMITH_ENABLE_CODEVELOP)
             )
 endif()
 
-install(EXPORT smith-targets 
+install(EXPORT smith-targets
         NAMESPACE smith::
         DESTINATION lib/cmake
         )

cmake/SmithBasics.cmake

@@ -37,6 +37,8 @@ if (NOT SMITH_BASICS_SETUP)
     cmake_dependent_option(SMITH_ENABLE_MPI "Enables Smith with MPI support" ON "ENABLE_MPI" OFF)
     cmake_dependent_option(SMITH_ENABLE_OPENMP "Enables Smith with OPENMP support" ON "ENABLE_OPENMP" OFF)
 
+    option(SMITH_ENABLE_GRETL "Enables Smith with Gretl Support" ON)
+
     #------------------------------------------------------------------------------
     # Profiling options
     #------------------------------------------------------------------------------

cmake/SmithConfigHeader.cmake

@@ -36,7 +36,7 @@ message(STATUS "Configuring Smith version ${SMITH_VERSION_FULL}")
 #------------------------------------------------------------------------------
 # Create variable for every TPL
 #------------------------------------------------------------------------------
-set(TPL_DEPS ADIAK AXOM CALIPER CAMP CONDUIT CONTINUATION CUDA ENZYME FMT HDF5 HIP LUA MFEM MPI PETSC RAJA SLEPC STRUMPACK SUNDIALS TRIBOL UMPIRE)
+set(TPL_DEPS ADIAK AXOM CALIPER CAMP CONDUIT CONTINUATION CUDA ENZYME FMT GRETL HDF5 HIP LUA MFEM MPI PETSC RAJA SLEPC STRUMPACK SUNDIALS TRIBOL UMPIRE)
 foreach(dep ${TPL_DEPS})
     if( ${dep}_FOUND OR ENABLE_${dep} )
         set(SMITH_USE_${dep} TRUE)

cmake/smith-config.cmake.in

@@ -40,6 +40,7 @@ if(NOT SMITH_FOUND)
   set(SMITH_USE_CONTINUATION   @SMITH_USE_CONTINUATION@)
   set(SMITH_USE_DFEM           @SMITH_USE_DFEM@)
   set(SMITH_USE_ENZYME         @SMITH_USE_ENZYME@)
+  set(SMITH_USE_GRETL          @SMITH_USE_GRETL@)
   set(SMITH_USE_HDF5           @SMITH_USE_HDF5@)
   set(SMITH_USE_MFEM           @SMITH_USE_MFEM@)
   set(SMITH_USE_MPI            @SMITH_USE_MPI@)
@@ -58,6 +59,7 @@ if(NOT SMITH_FOUND)
   set(SMITH_CONDUIT_DIR        "@CONDUIT_DIR@")
   set(SMITH_CONTINUATION_DIR   "@CONTINUATION_DIR@")
   set(SMITH_ENZYME_DIR         "@ENZYME_DIR@")
+  set(SMITH_GRETL_DIR          "@SMITH_INSTALL_PREFIX@")
   set(SMITH_HDF5_DIR           "@HDF5_DIR@")
   set(SMITH_MFEM_DIR           "@MFEM_DIR@")
   set(SMITH_PETSC_DIR          "@PETSC_DIR@")
@@ -170,6 +172,11 @@ if(NOT SMITH_FOUND)
     find_dependency(Enzyme REQUIRED)
   endif()
 
+  # Gretl
+  if(SMITH_USE_GRETL)
+    find_dependency(gretl REQUIRED PATHS ${SMITH_GRETL_DIR}/lib/cmake)
+  endif()
+
   #----------------------------------------------------------------------------
   # Load BLT and Smith targets
   #----------------------------------------------------------------------------

cmake/thirdparty/SetupSmithThirdParty.cmake

@@ -413,6 +413,7 @@ if (NOT SMITH_THIRD_PARTY_LIBRARIES_FOUND)
 
         set(MFEM_BUILT_WITH_CMAKE TRUE)
     endif()
+
     #------------------------------------------------------------------------------
     # ContinuationSolvers
     #------------------------------------------------------------------------------
@@ -451,6 +452,7 @@ if (NOT SMITH_THIRD_PARTY_LIBRARIES_FOUND)
         set(CONTINUATION_FOUND TRUE)
         add_subdirectory("${CONTINUATION_SOURCE_DIR}" ${CMAKE_BINARY_DIR}/ContinuationSolvers)
     endif()
+
     #------------------------------------------------------------------------------
     # Axom
     #------------------------------------------------------------------------------
@@ -553,6 +555,36 @@ if (NOT SMITH_THIRD_PARTY_LIBRARIES_FOUND)
         set(ENABLE_FORTRAN ON CACHE BOOL "" FORCE)
     endif()
 
+    #------------------------------------------------------------------------------
+    # Gretl
+    #------------------------------------------------------------------------------
+    message(STATUS "Smith Enable Gretl: ${SMITH_ENABLE_GRETL}")
+
+    if(SMITH_ENABLE_GRETL)
+        # Allow gretl as a non-submodule
+        if (DEFINED GRETL_SOURCE_DIR)
+            if(NOT EXISTS "${GRETL_SOURCE_DIR}/CMakeLists.txt")
+                message(FATAL_ERROR "Given GRETL_SOURCE_DIR does not contain CMakeLists.txt")
+            endif()
+        else()
+            set(GRETL_SOURCE_DIR "${PROJECT_SOURCE_DIR}/gretl" CACHE PATH "")
+
+            if (NOT EXISTS "${GRETL_SOURCE_DIR}/CMakeLists.txt")
+                message(FATAL_ERROR
+                    "The gretl repo is not present. "
+                    "Either run the following command in your git repository: \n"
+                    "    git submodule update --init --recursive\n"
+                    "Or add -DGRETL_SOURCE_DIR=/path/to/gretl to your CMake command." )
+            endif()
+        endif()
+
+        set(GRETL_FOUND TRUE)
+        message(STATUS "Gretl dir = ${GRETL_SOURCE_DIR}")
+        add_subdirectory("${GRETL_SOURCE_DIR}" ${CMAKE_BINARY_DIR}/gretl)
+
+    endif()
+
+
     #------------------------------------------------------------------------------
     # Tribol
     #------------------------------------------------------------------------------
@@ -678,6 +710,7 @@ if (NOT SMITH_THIRD_PARTY_LIBRARIES_FOUND)
     # Restore cleared Adiak/Caliper directories, reason at top of file.
     set(ADIAK_DIR ${_adiak_dir} CACHE PATH "" FORCE)
     set(CALIPER_DIR ${_caliper_dir} CACHE PATH "" FORCE)
+
     #------------------------------------------------------------------------------
     # Adiak
     #------------------------------------------------------------------------------

examples/buckling/CMakeLists.txt

@@ -7,7 +7,7 @@
 smith_add_executable( NAME        buckling_cylinder
                       SOURCES     cylinder.cpp
                       OUTPUT_DIR ${EXAMPLE_OUTPUT_DIRECTORY}
-                      DEPENDS_ON  smith_physics
+                      DEPENDS_ON  smith
                       )
 
 if(SMITH_ENABLE_HIP)

examples/conduction/CMakeLists.txt

@@ -7,7 +7,7 @@
 smith_add_executable( NAME        simple_conduction
                       SOURCES     simple_conduction.cpp
                       OUTPUT_DIR ${EXAMPLE_OUTPUT_DIRECTORY}
-                      DEPENDS_ON  smith_physics
+                      DEPENDS_ON  smith
                       )
 
 install(

examples/contact/CMakeLists.txt

@@ -19,7 +19,7 @@ if(TRIBOL_FOUND AND STRUMPACK_DIR)
         smith_add_executable(NAME        contact_${example_name}
                              SOURCES     ${filename}
                              OUTPUT_DIR  ${EXAMPLE_OUTPUT_DIRECTORY}
-                             DEPENDS_ON  smith_physics smith_mesh_utils)
+                             DEPENDS_ON  smith)
 
         smith_add_example_test(NAME           contact_${example_name}
                                COMMAND        contact_${example_name}

examples/explicit/CMakeLists.txt

@@ -16,7 +16,7 @@ if(SMITH_USE_DFEM)
         blt_add_executable(NAME        explicit_${example_name}
                            SOURCES     ${filename}
                            OUTPUT_DIR  ${EXAMPLE_OUTPUT_DIRECTORY}
-                           DEPENDS_ON  smith_physics smith_mesh_utils)
+                           DEPENDS_ON  smith)
 
         smith_add_example_test(NAME           explicit_${example_name}
                                COMMAND        explicit_${example_name}

examples/explicit/dfem_force_gpu.cpp

@@ -203,7 +203,7 @@ int main(int argc, char* argv[])
 
   axom::utilities::Timer timer(true);
   for (size_t step = 0; step < num_steps; ++step) {
-    auto no_mass_resid = solid_dfem_weak_form->residual(time, dt, &u_pred, pred_states);
+    auto no_mass_resid = solid_dfem_weak_form->residual(smith::TimeInfo(time, dt), &u_pred, pred_states);
     time += dt;
   }
   timer.stop();

examples/inertia_relief/CMakeLists.txt

@@ -8,7 +8,7 @@ if (SMITH_ENABLE_CONTINUATION AND STRUMPACK_DIR)
     smith_add_executable( NAME        inertia_relief_example
                           SOURCES     inertia_relief_example.cpp
                           OUTPUT_DIR ${EXAMPLE_OUTPUT_DIRECTORY}
-                          DEPENDS_ON  smith_physics continuation_solver
+                          DEPENDS_ON  smith continuation_solver
                           )
     smith_add_example_test(NAME            inertia_relief
                            COMMAND         inertia_relief_example

examples/inertia_relief/inertia_relief_example.cpp

@@ -127,6 +127,8 @@ auto createParaviewOutput(const mfem::ParMesh& mesh, const std::vector<smith::Fi
  * and thus the NLMC problem is guaranteed to be semi-monotone.
  */
 class InertialReliefProblem : public EqualityConstrainedHomotopyProblem {
+  InertialReliefProblem() : time_info_(0.0, 0.0, 0) {}
+
  protected:
   mfem::HypreParMatrix* drdu_ = nullptr;
   mfem::HypreParMatrix* dcdu_ = nullptr;
@@ -141,8 +143,7 @@ class InertialReliefProblem : public EqualityConstrainedHomotopyProblem {
   std::unique_ptr<smith::FiniteElementState> shape_disp_;
   std::shared_ptr<smith::Mesh> mesh_;
   std::vector<std::shared_ptr<smith::ScalarObjective>> constraints_;
-  double time_ = 0.0;
-  double dt_ = 0.0;
+  smith::TimeInfo time_info_;
   std::vector<double> jacobian_weights_ = {1.0, 0.0, 0.0, 0.0};
 
  public:
@@ -257,6 +258,7 @@ int main(int argc, char* argv[])
 
   double time = 0.0;
   double dt = 1.0;
+  smith::TimeInfo time_info(time, dt, 0);
   auto all_states = getConstFieldPointers(states, params);
   auto objective_states = {all_states[DISP], all_states[DENSITY]};
 
@@ -266,7 +268,7 @@ int main(int argc, char* argv[])
 
   mass_objective.addBodyIntegral(smith::DependsOn<1>{}, mesh->entireBodyName(),
                                  [](double /*t*/, auto /*X*/, auto RHO) { return get<smith::VALUE>(RHO); });
-  double mass = mass_objective.evaluate(time, dt, shape_disp.get(), objective_states);
+  double mass = mass_objective.evaluate(time_info, shape_disp.get(), objective_states);
 
   smith::tensor<double, dim> initial_cg;
 
@@ -278,7 +280,7 @@ int main(int argc, char* argv[])
                                       auto X, auto U, auto RHO) {
                                     return (get<smith::VALUE>(X)[i] + get<smith::VALUE>(U)[i]) * get<smith::VALUE>(RHO);
                                   });
-    initial_cg[i] = cg_objective->evaluate(time, dt, shape_disp.get(), objective_states) / mass;
+    initial_cg[i] = cg_objective->evaluate(time_info, shape_disp.get(), objective_states) / mass;
 
     constraints.push_back(cg_objective);
   }
@@ -343,7 +345,7 @@ InertialReliefProblem::InertialReliefProblem(std::vector<smith::FiniteElementSta
                                              std::shared_ptr<smith::Mesh> mesh,
                                              std::shared_ptr<SolidWeakFormT> weak_form,
                                              std::vector<std::shared_ptr<smith::ScalarObjective>> constraints)
-    : EqualityConstrainedHomotopyProblem()
+    : EqualityConstrainedHomotopyProblem(), time_info_(0.0, 0.0, 0)
 {
   weak_form_ = weak_form;
   mesh_ = mesh;
@@ -380,7 +382,7 @@ InertialReliefProblem::InertialReliefProblem(std::vector<smith::FiniteElementSta
 mfem::Vector InertialReliefProblem::residual(const mfem::Vector& u) const
 {
   obj_states_[DISP]->Set(1.0, u);
-  auto res_vector = weak_form_->residual(time_, dt_, shape_disp_.get(), smith::getConstFieldPointers(all_states_));
+  auto res_vector = weak_form_->residual(time_info_, shape_disp_.get(), smith::getConstFieldPointers(all_states_));
   return res_vector;
 }
 
@@ -402,7 +404,7 @@ mfem::Vector InertialReliefProblem::constraintJacobianTvp(const mfem::Vector& u,
 
   for (size_t i = 0; i < constraints_.size(); i++) {
     mfem::Vector grad_temp =
-        constraints_[i]->gradient(time_, dt_, shape_disp_.get(), smith::getConstFieldPointers(obj_states_), DISP);
+        constraints_[i]->gradient(time_info_, shape_disp_.get(), smith::getConstFieldPointers(obj_states_), DISP);
     constraint_gradient.Set(1.0, grad_temp);
     output_vec.Add(multipliers[i], constraint_gradient);
   }
@@ -414,7 +416,7 @@ mfem::HypreParMatrix* InertialReliefProblem::residualJacobian(const mfem::Vector
 {
   obj_states_[DISP]->Set(1.0, u);
   auto drdu_unique =
-      weak_form_->jacobian(time_, dt_, shape_disp_.get(), getConstFieldPointers(all_states_), jacobian_weights_);
+      weak_form_->jacobian(time_info_, shape_disp_.get(), getConstFieldPointers(all_states_), jacobian_weights_);
 
   if (drdu_) {
     delete drdu_;
@@ -437,7 +439,7 @@ mfem::Vector InertialReliefProblem::constraint(const mfem::Vector& u) const
     SLIC_ERROR_ROOT_IF(i2 != i, "Constraint index is out of range, bad cast from size_t to int");
 
     double constraint_i =
-        constraints_[i]->evaluate(time_, dt_, shape_disp_.get(), smith::getConstFieldPointers(obj_states_));
+        constraints_[i]->evaluate(time_info_, shape_disp_.get(), smith::getConstFieldPointers(obj_states_));
     if (dimc_ > 0) {
       output_vec(idx) = constraint_i;
     }
@@ -466,7 +468,7 @@ mfem::HypreParMatrix* InertialReliefProblem::constraintJacobian(const mfem::Vect
     SLIC_ERROR_ROOT_IF(i2 != i, "Constraint index is out of range, bad cast from size_t to int");
     mfem::HypreParVector gradVector(MPI_COMM_WORLD, dimuglb_, uOffsets_);
     gradVector.Set(
-        1.0, constraints_[i]->gradient(time_, dt_, shape_disp_.get(), smith::getConstFieldPointers(obj_states_), DISP));
+        1.0, constraints_[i]->gradient(time_info_, shape_disp_.get(), smith::getConstFieldPointers(obj_states_), DISP));
     mfem::Vector* globalGradVector = gradVector.GlobalVector();
     if (myid == 0) {
       dcdumat.SetRow(idx, cols, *globalGradVector);

examples/uniaxial/CMakeLists.txt

@@ -8,7 +8,7 @@ if (STRUMPACK_DIR)
     smith_add_executable( NAME        uniaxial_tension
                           SOURCES     uniaxial_tension.cpp
                           OUTPUT_DIR ${EXAMPLE_OUTPUT_DIRECTORY}
-                          DEPENDS_ON  smith_physics
+                          DEPENDS_ON  smith
                           )
 
     smith_add_example_test(NAME            uniaxial_tension

src/smith/CMakeLists.txt

@@ -8,3 +8,4 @@ add_subdirectory(infrastructure)
 add_subdirectory(numerics)
 add_subdirectory(physics)
 add_subdirectory(mesh_utils)
+add_subdirectory(differentiable_numerics)

src/smith/differentiable_numerics/CMakeLists.txt

@@ -0,0 +1,55 @@
+# Copyright (c) Lawrence Livermore National Security, LLC and
+# other Smith Project Developers. See the top-level LICENSE file for
+# details.
+#
+# SPDX-License-Identifier: (BSD-3-Clause)
+
+set(differentiable_numerics_sources
+    field_state.cpp
+    differentiable_utils.cpp
+    differentiable_physics.cpp
+    lumped_mass_explicit_newmark_state_advancer.cpp
+    )
+
+set(differentiable_numerics_headers
+    field_state.hpp
+    state_advancer.hpp
+    timestep_estimator.hpp
+    explicit_dynamic_solve.hpp
+    lumped_mass_weak_form.hpp
+    differentiable_utils.hpp
+    differentiable_physics.hpp
+    lumped_mass_explicit_newmark_state_advancer.hpp
+    )
+
+set(differentiable_numerics_depends
+    gretl
+    smith_state
+    smith_boundary_conditions
+    smith_contact
+    smith_numerics
+    smith_physics
+    )
+
+blt_add_library(
+    NAME        smith_differentiable_numerics
+    SOURCES     ${differentiable_numerics_sources}
+    HEADERS     ${differentiable_numerics_headers}
+    DEPENDS_ON  ${differentiable_numerics_depends}
+    )
+
+smith_write_unified_header(
+    NAME    differentiable_numerics
+    HEADERS ${differentiable_numerics_headers}
+    )
+
+install(FILES ${differentiable_numerics_headers} DESTINATION include/smith/differentiable_numerics )
+
+install(TARGETS              smith_differentiable_numerics
+        EXPORT               smith-targets
+        DESTINATION          lib
+        )
+
+if(SMITH_ENABLE_TESTS)
+    add_subdirectory(tests)
+endif()