time_integration_rule.hpp

Go to the documentation of this file.

// 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)
 
#pragma once
 
#include <tuple>
#include "smith/physics/common.hpp"
#include "smith/differentiable_numerics/field_state.hpp"
 
namespace smith {
 
class TimeIntegrationRule {
 public:
  virtual ~TimeIntegrationRule() {}
 
  virtual bool requiresInitialAccelerationSolve() const { return false; }
 
  virtual FieldState corrected_value(const TimeInfo& t, const std::vector<FieldState>& states) const = 0;
 
  virtual int num_args() const = 0;
 
  virtual FieldState corrected_dot(const TimeInfo& t, const std::vector<FieldState>& states) const = 0;
 
  virtual FieldState corrected_ddot(const TimeInfo& t, const std::vector<FieldState>& states) const = 0;
};
 
class BackwardEulerFirstOrderTimeIntegrationRule : public TimeIntegrationRule {
 public:
  static constexpr int num_states = 2;  
 
  BackwardEulerFirstOrderTimeIntegrationRule() {}
 
  int num_args() const override { return num_states; }
 
  template <typename T1, typename T2>
  SMITH_HOST_DEVICE auto value(const TimeInfo& /*t*/, const T1& field_new, const T2& /*field_old*/) const
  {
    return field_new;
  }
 
  template <typename T1, typename T2>
  SMITH_HOST_DEVICE auto dot(const TimeInfo& t, const T1& field_new, const T2& field_old) const
  {
    return (1.0 / t.dt()) * (field_new - field_old);
  }
 
  template <typename T1, typename T2>
  SMITH_HOST_DEVICE auto interpolate(const TimeInfo& t, const T1& field_new, const T2& field_old) const
  {
    return std::make_tuple(value(t, field_new, field_old), dot(t, field_new, field_old));
  }
 
  FieldState corrected_value(const TimeInfo& t, const std::vector<FieldState>& states) const override
  {
    return value(t, states[0], states[1]);
  }
 
  FieldState corrected_dot(const TimeInfo& t, const std::vector<FieldState>& states) const override
  {
    return dot(t, states[0], states[1]);
  }
 
  FieldState corrected_ddot(const TimeInfo& /*t*/, const std::vector<FieldState>& states) const override
  {
    SLIC_ERROR("BackwardEulerFirstOrderTimeIntegrationRule does not support second derivatives.");
    return states[0];
  }
};
 
class QuasiStaticRule : public TimeIntegrationRule {
 public:
  static constexpr int num_states = 1;  
 
  QuasiStaticRule() {}
 
  int num_args() const override { return num_states; }
 
  template <typename T1>
  SMITH_HOST_DEVICE auto value(const TimeInfo& /*t*/, const T1& field_new) const
  {
    return field_new;
  }
 
  template <typename T1>
  SMITH_HOST_DEVICE auto dot(const TimeInfo& /*t*/, const T1& /*field_new*/) const
  {
    return zero{};
  }
 
  template <typename T1>
  SMITH_HOST_DEVICE auto interpolate(const TimeInfo& t, const T1& field_new) const
  {
    return std::make_tuple(value(t, field_new), dot(t, field_new));
  }
 
  FieldState corrected_value(const TimeInfo& t, const std::vector<FieldState>& states) const override
  {
    return value(t, states[0]);
  }
 
  FieldState corrected_dot(const TimeInfo& /*t*/, const std::vector<FieldState>& states) const override
  {
    return zeroCopy(states[0]);
  }
 
  FieldState corrected_ddot(const TimeInfo& /*t*/, const std::vector<FieldState>& states) const override
  {
    SLIC_ERROR("QuasiStaticRule does not support second derivatives.");
    return states[0];
  }
};
 
class StaticTimeIntegrationRule : public TimeIntegrationRule {
 public:
  static constexpr int num_states = 1;  
 
  int num_args() const override { return num_states; }
 
  template <typename T1>
  SMITH_HOST_DEVICE auto value(const TimeInfo& /*t*/, const T1& field_new) const
  {
    return field_new;
  }
 
  template <typename T1>
  SMITH_HOST_DEVICE auto dot(const TimeInfo& /*t*/, const T1& /*field_new*/) const
  {
    return zero{};
  }
 
  template <typename T1>
  SMITH_HOST_DEVICE auto ddot(const TimeInfo& /*t*/, const T1& /*field_new*/) const
  {
    return zero{};
  }
 
  template <typename T1>
  SMITH_HOST_DEVICE auto interpolate(const TimeInfo& t, const T1& field_new) const
  {
    return std::make_tuple(value(t, field_new), dot(t, field_new), ddot(t, field_new));
  }
 
  FieldState corrected_value(const TimeInfo& t, const std::vector<FieldState>& states) const override
  {
    return value(t, states[0]);
  }
 
  FieldState corrected_dot(const TimeInfo& /*t*/, const std::vector<FieldState>& states) const override
  {
    return zeroCopy(states[0]);
  }
 
  FieldState corrected_ddot(const TimeInfo& /*t*/, const std::vector<FieldState>& states) const override
  {
    return zeroCopy(states[0]);
  }
};
 
using QuasiStaticFirstOrderTimeIntegrationRule = BackwardEulerFirstOrderTimeIntegrationRule;
 
struct ImplicitNewmarkSecondOrderTimeIntegrationRule : public TimeIntegrationRule {
 public:
  static constexpr int num_states = 4;  
 
  ImplicitNewmarkSecondOrderTimeIntegrationRule() {}
 
  int num_args() const override { return num_states; }
 
  bool requiresInitialAccelerationSolve() const override { return true; }
 
  template <typename T1, typename T2, typename T3, typename T4>
  SMITH_HOST_DEVICE auto value([[maybe_unused]] const TimeInfo& t, [[maybe_unused]] const T1& field_new,
                               [[maybe_unused]] const T2& field_old, [[maybe_unused]] const T3& velo_old,
                               [[maybe_unused]] const T4& accel_old) const
  {
    auto regular_value = field_new + (field_old - field_old);
    auto cycle_zero_correction = field_old - field_new;
    if (t.isCycleZeroEvaluation()) {
      return regular_value + cycle_zero_correction;
    }
    return regular_value + (cycle_zero_correction - cycle_zero_correction);
  }
 
  template <typename T1, typename T2, typename T3, typename T4>
  SMITH_HOST_DEVICE auto dot([[maybe_unused]] const TimeInfo& t, [[maybe_unused]] const T1& field_new,
                             [[maybe_unused]] const T2& field_old, [[maybe_unused]] const T3& velo_old,
                             [[maybe_unused]] const T4& accel_old) const
  {
    auto regular_dot = (2.0 / t.dt()) * (field_new - field_old) - velo_old;
    auto cycle_zero_correction = velo_old - regular_dot;
    if (t.isCycleZeroEvaluation()) {
      return regular_dot + cycle_zero_correction;
    }
    return regular_dot + (cycle_zero_correction - cycle_zero_correction);
  }
 
  template <typename T1, typename T2, typename T3, typename T4>
  SMITH_HOST_DEVICE auto ddot([[maybe_unused]] const TimeInfo& t, [[maybe_unused]] const T1& field_new,
                              [[maybe_unused]] const T2& field_old, [[maybe_unused]] const T3& velo_old,
                              [[maybe_unused]] const T4& accel_old) const
  {
    auto dt = t.dt();
    auto regular_ddot = (4.0 / (dt * dt)) * (field_new - field_old) - (4.0 / dt) * velo_old - accel_old;
    auto cycle_zero_correction = accel_old - regular_ddot;
    if (t.isCycleZeroEvaluation()) {
      return regular_ddot + cycle_zero_correction;
    }
    return regular_ddot + (cycle_zero_correction - cycle_zero_correction);
  }
 
  template <typename T1, typename T2, typename T3, typename T4>
  SMITH_HOST_DEVICE auto interpolate(const TimeInfo& t, const T1& field_new, const T2& field_old, const T3& velo_old,
                                     const T4& accel_old) const
  {
    return std::make_tuple(value(t, field_new, field_old, velo_old, accel_old),
                           dot(t, field_new, field_old, velo_old, accel_old),
                           ddot(t, field_new, field_old, velo_old, accel_old));
  }
 
  FieldState corrected_value(const TimeInfo& t, const std::vector<FieldState>& states) const override
  {
    return value(t, states[0], states[1], states[2], states[3]);
  }
 
  FieldState corrected_dot(const TimeInfo& t, const std::vector<FieldState>& states) const override
  {
    return dot(t, states[0], states[1], states[2], states[3]);
  }
 
  FieldState corrected_ddot(const TimeInfo& t, const std::vector<FieldState>& states) const override
  {
    return ddot(t, states[0], states[1], states[2], states[3]);
  }
};
 
struct QuasiStaticSecondOrderTimeIntegrationRule : public TimeIntegrationRule {
 public:
  static constexpr int num_states = 4;  
 
  QuasiStaticSecondOrderTimeIntegrationRule() {}
 
  int num_args() const override { return num_states; }
 
  bool requiresInitialAccelerationSolve() const override { return false; }
 
  template <typename T1, typename T2, typename T3, typename T4>
  SMITH_HOST_DEVICE auto value([[maybe_unused]] const TimeInfo& t, [[maybe_unused]] const T1& field_new,
                               [[maybe_unused]] const T2& field_old, [[maybe_unused]] const T3& velo_old,
                               [[maybe_unused]] const T4& accel_old) const
  {
    return field_new;
  }
 
  template <typename T1, typename T2, typename T3, typename T4>
  SMITH_HOST_DEVICE auto dot([[maybe_unused]] const TimeInfo& t, [[maybe_unused]] const T1& field_new,
                             [[maybe_unused]] const T2& field_old, [[maybe_unused]] const T3& velo_old,
                             [[maybe_unused]] const T4& accel_old) const
  {
    return (1.0 / t.dt()) * (field_new - field_old);
  }
 
  template <typename T1, typename T2, typename T3, typename T4>
  SMITH_HOST_DEVICE auto ddot([[maybe_unused]] const TimeInfo& t, [[maybe_unused]] const T1& field_new,
                              [[maybe_unused]] const T2& field_old, [[maybe_unused]] const T3& velo_old,
                              [[maybe_unused]] const T4& accel_old) const
  {
    return zero{};
  }
 
  template <typename T1, typename T2, typename T3, typename T4>
  SMITH_HOST_DEVICE auto interpolate(const TimeInfo& t, const T1& field_new, const T2& field_old, const T3& velo_old,
                                     const T4& accel_old) const
  {
    return std::make_tuple(value(t, field_new, field_old, velo_old, accel_old),
                           dot(t, field_new, field_old, velo_old, accel_old),
                           ddot(t, field_new, field_old, velo_old, accel_old));
  }
 
  FieldState corrected_value(const TimeInfo& t, const std::vector<FieldState>& states) const override
  {
    return value(t, states[0], states[1], states[2], states[3]);
  }
 
  FieldState corrected_dot(const TimeInfo& t, const std::vector<FieldState>& states) const override
  {
    return dot(t, states[0], states[1], states[2], states[3]);
  }
 
  FieldState corrected_ddot(const TimeInfo& /*t*/, const std::vector<FieldState>& states) const override
  {
    return zeroCopy(states[0]);
  }
};
 
}  // namespace smith