odes.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 <functional>
#include <memory>
 
#include "mfem.hpp"
 
#include "smith/physics/boundary_conditions/boundary_condition_manager.hpp"
#include "smith/numerics/equation_solver.hpp"
#include "smith/numerics/solver_config.hpp"
 
namespace smith::mfem_ext {
 
class SecondOrderODE : public mfem::SecondOrderTimeDependentOperator {
 public:
  static constexpr double epsilon = 0.0001;
 
  struct State {
    double& time;
 
    double& c0;
 
    double& c1;
 
    mfem::Vector& u;
 
    mfem::Vector& du_dt;
 
    mfem::Vector& d2u_dt2;
  };
 
  SecondOrderODE(int n, State&& state, const EquationSolver& solver, const BoundaryConditionManager& bcs);
 
  void Mult(const mfem::Vector& u, const mfem::Vector& du_dt, mfem::Vector& d2u_dt2) const override
  {
    Solve(t, 0.0, 0.0, u, du_dt, d2u_dt2);
  }
 
  void ImplicitSolve(const double c0, const double c1, const mfem::Vector& u, const mfem::Vector& du_dt,
                     mfem::Vector& d2u_dt2) override
  {
    Solve(t, c0, c1, u, du_dt, d2u_dt2);
  }
 
  void ImplicitSolve(const double dt, const mfem::Vector& u, mfem::Vector& du_dt) override;
 
  void SetEnforcementMethod(const DirichletEnforcementMethod method) { enforcement_method_ = method; }
 
  void SetTimestepper(const smith::TimestepMethod timestepper);
 
  void Step(mfem::Vector& x, mfem::Vector& dxdt, double& time, double& dt);
 
  const State& GetState() { return state_; }
 
  TimestepMethod GetTimestepper() { return timestepper_; }
 
 private:
  void Solve(const double time, const double c0, const double c1, const mfem::Vector& u, const mfem::Vector& du_dt,
             mfem::Vector& d2u_dt2) const;
 
  State state_;
  DirichletEnforcementMethod enforcement_method_ = smith::DirichletEnforcementMethod::RateControl;
  const EquationSolver& solver_;
  std::unique_ptr<mfem::SecondOrderODESolver> second_order_ode_solver_;
 
  std::unique_ptr<mfem::ODESolver> first_order_system_ode_solver_;
 
  const BoundaryConditionManager& bcs_;
  mfem::Vector zero_;
 
  mutable mfem::Vector U_minus_;
  mutable mfem::Vector U_;
  mutable mfem::Vector U_plus_;
  mutable mfem::Vector dU_dt_;
  mutable mfem::Vector d2U_dt2_;
 
  smith::TimestepMethod timestepper_;
};
 
class FirstOrderODE : public mfem::TimeDependentOperator {
 public:
  static constexpr double epsilon = 0.000001;
 
  struct State {
    double& time;
 
    mfem::Vector& u;
 
    double& dt;
 
    mfem::Vector& du_dt;
 
    double& previous_dt;
  };
 
  FirstOrderODE(int n, FirstOrderODE::State&& state, const EquationSolver& solver, const BoundaryConditionManager& bcs);
 
  void Mult(const mfem::Vector& u, mfem::Vector& du_dt) const { Solve(t, 0.0, u, du_dt); }
 
  void ImplicitSolve(const double dt, const mfem::Vector& u, mfem::Vector& du_dt) { Solve(t, dt, u, du_dt); }
 
  void SetEnforcementMethod(const DirichletEnforcementMethod method) { enforcement_method_ = method; }
 
  void SetTimestepper(const smith::TimestepMethod timestepper);
 
  void Step(mfem::Vector& x, double& time, double& dt)
  {
    if (ode_solver_) {
      ode_solver_->Step(x, time, dt);
    } else {
      SLIC_ERROR("ode_solver_ unspecified");
    }
  }
 
  TimestepMethod GetTimestepper() { return timestepper_; }
 
 private:
  virtual void Solve(const double time, const double dt, const mfem::Vector& u, mfem::Vector& du_dt) const;
 
  FirstOrderODE::State state_;
 
  DirichletEnforcementMethod enforcement_method_ = smith::DirichletEnforcementMethod::RateControl;
  const EquationSolver& solver_;
  std::unique_ptr<mfem::ODESolver> ode_solver_;
  const BoundaryConditionManager& bcs_;
  mfem::Vector zero_;
 
  mutable mfem::Vector U_minus_;
  mutable mfem::Vector U_;
  mutable mfem::Vector U_plus_;
  mutable mfem::Vector dU_dt_;
 
  TimestepMethod timestepper_;
};
 
}  // namespace smith::mfem_ext