Tupek/gretl field state no error check

CMakeLists.txt

@@ -198,7 +198,9 @@ set(exported_targets
     serac_infrastructure
     serac_mesh_utils
     serac_numerics
-    serac_physics)
+    serac_physics
+    serac_gretl
+    serac_differentiable_numerics)
 
 add_library(serac INTERFACE)
 target_link_libraries(serac INTERFACE ${exported_targets})

examples/inertia_relief/inertia_relief_example.cpp

@@ -127,6 +127,8 @@ auto createParaviewOutput(const mfem::ParMesh& mesh, const std::vector<serac::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<serac::FiniteElementState> shape_disp_;
   std::shared_ptr<serac::Mesh> mesh_;
   std::vector<std::shared_ptr<serac::ScalarObjective>> constraints_;
-  double time_ = 0.0;
-  double dt_ = 0.0;
+  serac::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;
+  serac::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(serac::DependsOn<1>{}, mesh->entireBodyName(),
                                  [](double /*t*/, auto /*X*/, auto RHO) { return get<serac::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);
 
   serac::tensor<double, dim> initial_cg;
 
@@ -278,7 +280,7 @@ int main(int argc, char* argv[])
                                       auto X, auto U, auto RHO) {
                                     return (get<serac::VALUE>(X)[i] + get<serac::VALUE>(U)[i]) * get<serac::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<serac::FiniteElementSta
                                              std::shared_ptr<serac::Mesh> mesh,
                                              std::shared_ptr<SolidWeakFormT> weak_form,
                                              std::vector<std::shared_ptr<serac::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<serac::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(), serac::getConstFieldPointers(all_states_));
+  auto res_vector = weak_form_->residual(time_info_, shape_disp_.get(), serac::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(), serac::getConstFieldPointers(obj_states_), DISP);
+        constraints_[i]->gradient(time_info_, shape_disp_.get(), serac::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(), serac::getConstFieldPointers(obj_states_));
+        constraints_[i]->evaluate(time_info_, shape_disp_.get(), serac::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(), serac::getConstFieldPointers(obj_states_), DISP));
+        1.0, constraints_[i]->gradient(time_info_, shape_disp_.get(), serac::getConstFieldPointers(obj_states_), DISP));
     mfem::Vector* globalGradVector = gradVector.GlobalVector();
     if (myid == 0) {
       dcdumat.SetRow(idx, cols, *globalGradVector);

src/serac/CMakeLists.txt

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

src/serac/differentiable_numerics/CMakeLists.txt

@@ -0,0 +1,55 @@
+# Copyright (c) Lawrence Livermore National Security, LLC and
+# other Serac Project Developers. See the top-level LICENSE file for
+# details.
+#
+# SPDX-License-Identifier: (BSD-3-Clause)
+
+set(differentiable_numerics_sources
+    state_advancer.cpp
+    field_state.cpp
+    differentiable_utils.cpp
+    differentiable_physics.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
+    )
+
+set(differentiable_numerics_depends
+    serac_state
+    serac_boundary_conditions
+    serac_contact
+    serac_numerics
+    serac_physics
+    serac_gretl
+    )
+
+blt_add_library(
+    NAME        serac_differentiable_numerics
+    SOURCES     ${differentiable_numerics_sources}
+    HEADERS     ${differentiable_numerics_headers}
+    DEPENDS_ON  ${differentiable_numerics_depends}
+    )
+
+serac_write_unified_header(
+    NAME    differentiable_numerics
+    HEADERS ${differentiable_numerics_headers}
+    )
+
+install(FILES ${differentiable_numerics_headers} DESTINATION include/serac/differentiable_numerics )
+
+install(TARGETS              serac_differentiable_numerics
+        EXPORT               serac-targets
+        DESTINATION          lib
+        )
+
+if(SERAC_ENABLE_TESTS)
+    add_subdirectory(tests)
+endif()
+