base_physics.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 <cstddef>
#include <format>
#include <memory>
#include <optional>
#include <string>
#include <unordered_map>
#include <vector>
 
#include "mpi.h"
#include "mfem.hpp"
#include "axom/sidre.hpp"
#include "axom/slic.hpp"
 
#include "smith/physics/boundary_conditions/boundary_condition_manager.hpp"
#include "smith/numerics/equation_solver.hpp"
#include "smith/physics/state/finite_element_state.hpp"
#include "smith/physics/state/finite_element_dual.hpp"
#include "smith/physics/state/state_manager.hpp"
#include "smith/physics/common.hpp"
 
namespace smith {
 
class Mesh;
 
namespace detail {
std::string addPrefix(const std::string& prefix, const std::string& target);
 
std::string removePrefix(const std::string& prefix, const std::string& target);
 
}  // namespace detail
 
class BasePhysics {
 public:
  BasePhysics(std::string physics_name, std::shared_ptr<smith::Mesh> mesh, int cycle = 0, double time = 0.0,
              bool checkpoint_to_disk = false);
 
  BasePhysics(BasePhysics&& other) = default;
 
  virtual double time() const;
 
  virtual int cycle() const;
 
  virtual double maxTime() const;
 
  virtual double minTime() const;
 
  virtual int maxCycle() const;
 
  virtual int minCycle() const;
 
  bool isQuasistatic() const { return is_quasistatic_; }
 
  virtual const std::vector<double>& timesteps() const;
 
  virtual void resetStates(int cycle = 0, double time = 0.0) = 0;
 
  virtual void resetAdjointStates()
  {
    time_ = max_time_;
    cycle_ = max_cycle_;
  }
 
  virtual void completeSetup() = 0;
 
  virtual const FiniteElementState& state(const std::string& state_name) const = 0;
 
  virtual void setState(const std::string&, const FiniteElementState&) = 0;
 
  virtual std::vector<std::string> stateNames() const = 0;
 
  virtual const FiniteElementState& adjoint(const std::string& adjoint_name) const = 0;
 
  virtual std::vector<std::string> adjointNames() const { return {}; }
 
  virtual std::vector<std::string> dualNames() const { return {}; }
 
  virtual const FiniteElementDual& dual(const std::string& dual_name) const
  {
    SLIC_ERROR_ROOT(
        std::format("dual '{}' requested from physics module '{}' which does not support duals", dual_name, name_));
    return *duals_[0];
  }
 
  virtual const FiniteElementState& dualAdjoint(const std::string& dual_name) const
  {
    SLIC_ERROR_ROOT(std::format("dualAdjoint '{}' requested from physics module '{}' which does not support duals",
                                dual_name, name_));
    return *dual_adjoints_[0];
  }
 
  virtual const FiniteElementState& shapeDisplacement() const;
 
  virtual const FiniteElementState& parameter(const std::string& parameter_name) const
  {
    std::string appended_name = detail::addPrefix(name_, parameter_name);
 
    for (auto& parameter : parameters_) {
      if (appended_name == parameter.state->name()) {
        return *parameter.state;
      }
    }
 
    SLIC_ERROR_ROOT(
        std::format("Parameter {} requested from physics module {}, but it doesn't exist.", parameter_name, name_));
 
    return *states_[0];
  }
 
  virtual const FiniteElementState& parameter(std::size_t parameter_index) const
  {
    SLIC_ERROR_ROOT_IF(parameter_index >= parameters_.size(),
                       std::format("Parameter index {} requested, but only {} parameters exist in physics module {}.",
                                   parameter_index, parameters_.size(), name_));
 
    return *parameters_[parameter_index].state;
  }
 
  virtual std::vector<std::string> parameterNames() const
  {
    std::vector<std::string> parameter_names;
 
    for (auto& parameter : parameters_) {
      parameter_names.emplace_back(detail::removePrefix(name_, parameter.state->name()));
    }
 
    return parameter_names;
  }
 
  virtual void setParameter(const size_t parameter_index, const FiniteElementState& parameter_state);
 
  virtual void setShapeDisplacement(const FiniteElementState& shape_displacement);
 
  virtual FiniteElementDual computeTimestepSensitivity(size_t parameter_index);
 
  virtual const FiniteElementDual& computeTimestepShapeSensitivity();
 
  virtual const std::unordered_map<std::string, const smith::FiniteElementDual&> computeInitialConditionSensitivity()
      const
  {
    SLIC_WARNING_ROOT(std::format("Initial condition sensitivities not enabled in physics module {}", name_));
    return {};
  }
 
  virtual void advanceTimestep(double dt) = 0;
 
  virtual void setAdjointLoad(std::unordered_map<std::string, const smith::FiniteElementDual&> string_to_dual)
  {
    if (!string_to_dual.empty()) {
      SLIC_ERROR_ROOT(
          std::format("Failed to setAdjointLoad.  Adjoint analysis not defined for physics module {}", name_));
    }
  }
 
  virtual void setDualAdjointBcs(std::unordered_map<std::string, const smith::FiniteElementState&> string_to_bc)
  {
    if (!string_to_bc.empty()) {
      SLIC_ERROR_ROOT(
          std::format("Failed to setDualAdjointBCs.  Adjoint analysis not defined for physics module {}", name_));
    }
  }
 
  virtual void reverseAdjointTimestep()
  {
    SLIC_ERROR_ROOT(std::format("Adjoint analysis not defined for physics module {}", name_));
  }
 
  virtual void outputStateToDisk(std::optional<std::string> paraview_output_dir = {}) const;
 
  void loadCheckpointedStatesFromDisk(int cycle);
 
  virtual FiniteElementState loadCheckpointedState(const std::string& state_name, int cycle);
 
  virtual FiniteElementDual loadCheckpointedDual([[maybe_unused]] const std::string& state_name,
                                                 [[maybe_unused]] int cycle)
  {
    SLIC_ERROR_ROOT(std::format("loadCheckpointedDual not enabled in physics module {}", name_));
    return *duals_[0];
  }
 
  virtual double getCheckpointedTimestep(int cycle) const;
 
  virtual void initializeSummary(axom::sidre::DataStore& datastore, const double t_final, const double dt) const;
 
  virtual void saveSummary(axom::sidre::DataStore& datastore, const double t) const;
 
  virtual ~BasePhysics() = default;
 
  const smith::Mesh& mesh() const;
 
  const mfem::ParMesh& mfemParMesh() const;
 
  mfem::ParMesh& mfemParMesh();
 
  std::string name() const { return name_; }
 
 protected:
  const FiniteElementDual& shapeDisplacementSensitivity() const;
 
  void CreateParaviewDataCollection() const;
 
  void UpdateParaviewDataCollection(const std::string& paraview_output_dir) const;
 
  void initializeBasePhysicsStates(int cycle, double time);
 
  std::unordered_map<std::string, FiniteElementState> getCheckpointedStates(int cycle);
 
  std::string name_ = {};
 
  std::shared_ptr<smith::Mesh> mesh_;
 
  MPI_Comm comm_;
 
  std::vector<const smith::FiniteElementState*> states_;
 
  std::vector<const smith::FiniteElementState*> adjoints_;
 
  std::vector<const smith::FiniteElementDual*> duals_;
 
  std::vector<const smith::FiniteElementState*> dual_adjoints_;
 
  struct ParameterInfo {
    template <typename FunctionSpace>
    ParameterInfo(mfem::ParMesh& mesh, FunctionSpace space, const std::string& name = "")
    {
      state = std::make_unique<FiniteElementState>(mesh, space, name);
      previous_state = std::make_unique<FiniteElementState>(mesh, space, "previous_" + name);
      sensitivity = std::make_unique<FiniteElementDual>(mesh, space, name + "_sensitivity");
      StateManager::storeState(*state);
      StateManager::storeDual(*sensitivity);
    }
 
    std::unique_ptr<smith::FiniteElementState> state;
 
    std::unique_ptr<smith::FiniteElementState> previous_state;
 
    std::unique_ptr<smith::FiniteElementDual> sensitivity;
  };
 
  std::vector<ParameterInfo> parameters_;
 
  mutable std::unordered_map<std::string, std::vector<smith::FiniteElementState>> checkpoint_states_;
 
  mutable std::unordered_map<std::string, smith::FiniteElementState> cached_checkpoint_states_;
 
  mutable std::optional<int> cached_checkpoint_cycle_;
 
  bool is_quasistatic_ = true;
 
  static constexpr int FLOAT_PRECISION_ = 8;
 
  double time_;
 
  double dt_;
 
  double max_time_;
 
  double min_time_;
 
  std::vector<double> timesteps_;
 
  int cycle_;
 
  int max_cycle_;
 
  int min_cycle_;
 
  double ode_time_point_;
 
  int mpi_rank_;
 
  int mpi_size_;
 
  mutable std::unique_ptr<mfem::ParaViewDataCollection> paraview_dc_;
 
  mutable std::unordered_map<std::string, std::unique_ptr<mfem::ParGridFunction>> paraview_dual_grid_functions_;
 
  mutable std::unique_ptr<mfem::ParGridFunction> shape_sensitivity_grid_function_;
 
  smith::FiniteElementState shape_displacement_;
 
  smith::FiniteElementDual shape_displacement_dual_;
 
  BoundaryConditionManager bcs_;
 
  bool checkpoint_to_disk_;
};
 
}  // namespace smith