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];
  }
};
 
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
  {
    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 (2.0 / t.dt()) * (field_new - field_old) - velo_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
  {
    auto dt = t.dt();
    return (4.0 / (dt * dt)) * (field_new - field_old) - (4.0 / dt) * velo_old - accel_old;
  }
 
  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