Updated explicit time integration (idaholab#30335)

Author: maxnezdyur

modules/contact/test/tests/explicit_dynamics/deep_impact.i

@@ -257,6 +257,7 @@
     type = DirectCentralDifference
     mass_matrix_tag = 'mass'
     use_constant_mass = true
+    second_order_vars = 'disp_x disp_y disp_z'
   []
   skip_exception_check = true
 []

modules/contact/test/tests/explicit_dynamics/highvel.i

@@ -345,6 +345,7 @@
   [TimeIntegrator]
     type = DirectCentralDifference
     mass_matrix_tag = mass
+    second_order_vars = 'disp_x disp_y disp_z'
   []
 []
 

modules/contact/test/tests/explicit_dynamics/settlement.i

@@ -345,6 +345,7 @@
   [TimeIntegrator]
     type = DirectCentralDifference
     mass_matrix_tag = 'mass'
+    second_order_vars = 'disp_x disp_y disp_z'
   []
 []
 

modules/contact/test/tests/explicit_dynamics/test_balance.i

@@ -305,6 +305,7 @@
   [TimeIntegrator]
     type = DirectCentralDifference
     mass_matrix_tag = 'mass'
+    second_order_vars = 'disp_x disp_y disp_z'
   []
 []
 

modules/contact/test/tests/explicit_dynamics/test_balance_optimized.i

@@ -263,6 +263,7 @@
     type = DirectCentralDifference
     mass_matrix_tag = 'mass'
     use_constant_mass = true
+    second_order_vars = 'disp_x disp_y disp_z'
   []
 []
 [Outputs]

modules/solid_mechanics/include/bcs/DirectDirichletBCBase.h

@@ -9,6 +9,7 @@
 
 #pragma once
 
+#include "DirectCentralDifference.h"
 #include "NodalBC.h"
 #include "libmesh/numeric_vector.h"
 
@@ -30,7 +31,7 @@ class DirectDirichletBCBase : public NodalBC
 
 protected:
   virtual Real computeQpResidual() override;
-
+  virtual void timestepSetup() override;
   /**
    * Compute the value of the DirichletBC at the current quadrature point
    */
@@ -41,4 +42,10 @@ class DirectDirichletBCBase : public NodalBC
 
   const Real & _u_old;
   const Real & _u_dot_old;
+
+  /// explicit time integrator
+  const DirectCentralDifference * _exp_integrator;
+
+  /// variable time order need for computing BC
+  DirectCentralDifference::TimeOrder _var_time_order;
 };

modules/solid_mechanics/include/timeintegrators/DirectCentralDifference.h

@@ -10,6 +10,7 @@
 #pragma once
 
 #include "ExplicitTimeIntegrator.h"
+#include "libmesh/string_to_enum.h"
 
 // Forward declarations
 namespace libMesh
@@ -42,6 +43,19 @@ class DirectCentralDifference : public ExplicitTimeIntegrator
   {
     mooseError("NOT SUPPORTED");
   }
+  virtual void init() override;
+
+  enum TimeOrder
+  {
+    FIRST,
+    SECOND
+  };
+
+  /**
+   * Retrieve the order of the highest time derivative of a variable.
+   * @return Returns the time order enum of this variable.
+   */
+  virtual TimeOrder findVariableTimeOrder(unsigned int var_num) const;
 
 protected:
   virtual TagID massMatrixTagID() const override;
@@ -52,6 +66,19 @@ class DirectCentralDifference : public ExplicitTimeIntegrator
 
   /// The older solution
   const NumericVector<Number> & _solution_older;
+
+  // Variables that forward euler time integration will be used
+  std::unordered_set<unsigned int> & _vars_first;
+
+  // local dofs that will have forward euler time integration
+  std::vector<dof_id_type> & _local_first_order_indicies;
+
+  // Variables that central difference time integration will be used
+  std::unordered_set<unsigned int> & _vars_second;
+
+  // local dofs that will have central difference time integration
+  std::vector<dof_id_type> & _local_second_order_indicies;
+
   /**
    * Helper function that actually does the math for computing the time derivative
    */

modules/solid_mechanics/src/bcs/DirectDirichletBCBase.C

@@ -7,8 +7,12 @@
 //* Licensed under LGPL 2.1, please see LICENSE for details
 //* https://www.gnu.org/licenses/lgpl-2.1.html
 
+#include "DirectCentralDifference.h"
 #include "DirectDirichletBCBase.h"
+#include "MooseError.h"
 #include "NonlinearSystemBase.h"
+#include <iostream>
+#include <ostream>
 
 InputParameters
 DirectDirichletBCBase::validParams()
@@ -21,21 +25,40 @@ DirectDirichletBCBase::DirectDirichletBCBase(const InputParameters & parameters)
   : NodalBC(parameters),
     _mass_diag(_sys.getVector("mass_matrix_diag_inverted")),
     _u_old(_var.nodalValueOld()),
-    _u_dot_old(_var.nodalValueDotOld())
+    _u_dot_old(_var.nodalValueDotOld()),
+    _exp_integrator(
+        dynamic_cast<const DirectCentralDifference *>(&_sys.getTimeIntegrator(_var.number())))
 {
+  if (!_exp_integrator)
+    mooseError("Time integrator for the variable is not of the right type.");
 }
 
 Real
 DirectDirichletBCBase::computeQpResidual()
 {
   // Get dof for current var
   const auto dofnum = _variable->nodalDofIndex();
+  Real resid = 0;
+  // Compute residual to enforce BC based on time order
+  switch (_var_time_order)
+  {
+    case DirectCentralDifference::FIRST:
+      resid = (computeQpValue() - _u_old) / _dt;
+      resid /= -_mass_diag(dofnum);
+      break;
 
-  // Compute residual to enforce BC
-  // This is the force required to enforce the BC in a central difference scheme
-  Real avg_dt = (_dt + _dt_old) / 2;
-  Real resid = (computeQpValue() - _u_old) / (avg_dt * _dt) - (_u_dot_old) / avg_dt;
-  resid /= -_mass_diag(dofnum);
-
+    case DirectCentralDifference::SECOND:
+      Real avg_dt = (_dt + _dt_old) / 2;
+      resid = (computeQpValue() - _u_old) / (avg_dt * _dt) - (_u_dot_old) / avg_dt;
+      resid /= -_mass_diag(dofnum);
+      break;
+  }
   return resid;
 }
+
+void
+DirectDirichletBCBase::timestepSetup()
+{
+  // Now is the point that the time integrator has the variable time orders setup
+  _var_time_order = _exp_integrator->findVariableTimeOrder(_var.number());
+}