base_physics.hpp

Go to the documentation of this file.

// Copyright (c) 2019-2023, Lawrence Livermore National Security, LLC and
// other Serac 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 "axom/sidre.hpp"
 
#include "serac/physics/boundary_conditions/boundary_condition_manager.hpp"
#include "serac/numerics/equation_solver.hpp"
#include "serac/physics/state/finite_element_state.hpp"
#include "serac/physics/state/finite_element_dual.hpp"
#include "serac/physics/state/state_manager.hpp"
#include "serac/physics/common.hpp"
 
namespace serac {
 
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::string mesh_tag, 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 std::vector<double> timesteps() const;
 
  virtual void resetStates(int cycle = 0, double time = 0.0) = 0;
 
  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 {}; }
 
  const FiniteElementState& shapeDisplacement() const { return shape_displacement_; }
 
  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(axom::fmt::format("Parameter {} requested from physics module {}, but it doesn't exist.",
                                      parameter_name, name_));
 
    return *states_[0];
  }
 
  const FiniteElementState& parameter(std::size_t parameter_index) const
  {
    SLIC_ERROR_ROOT_IF(
        parameter_index >= parameters_.size(),
        axom::fmt::format("Parameter index {} requested, but only {} parameters exist in physics module {}.",
                          parameter_index, parameters_.size(), name_));
 
    return *parameters_[parameter_index].state;
  }
 
  std::vector<std::string> parameterNames()
  {
    std::vector<std::string> parameter_names;
 
    for (auto& parameter : parameters_) {
      parameter_names.emplace_back(detail::removePrefix(name_, parameter.state->name()));
    }
 
    return parameter_names;
  }
 
  void setParameter(const size_t parameter_index, const FiniteElementState& parameter_state);
 
  void setShapeDisplacement(const FiniteElementState& shape_displacement);
 
  virtual const FiniteElementDual& computeTimestepSensitivity(size_t /* parameter_index */)
  {
    SLIC_ERROR_ROOT(axom::fmt::format("Parameter sensitivities not enabled in physics module {}", name_));
    return *parameters_[0].sensitivity;
  }
 
  virtual const FiniteElementDual& computeTimestepShapeSensitivity()
  {
    SLIC_ERROR_ROOT(axom::fmt::format("Shape sensitivities not enabled in physics module {}", name_));
    return *shape_displacement_sensitivity_;
  }
 
  virtual const std::unordered_map<std::string, const serac::FiniteElementDual&> computeInitialConditionSensitivity()
  {
    SLIC_ERROR_ROOT(axom::fmt::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 serac::FiniteElementDual&>)
  {
    SLIC_ERROR_ROOT(axom::fmt::format("Adjoint analysis not defined for physics module {}", name_));
  }
 
  virtual void reverseAdjointTimestep()
  {
    SLIC_ERROR_ROOT(axom::fmt::format("Adjoint analysis not defined for physics module {}", name_));
  }
 
  virtual void outputStateToDisk(std::optional<std::string> paraview_output_dir = {}) const;
 
  FiniteElementState loadCheckpointedState(const std::string& state_name, int cycle) const;
 
  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 mfem::ParMesh& mesh() const { return mesh_; }
  mfem::ParMesh& mesh() { return mesh_; }
 
protected:
  void CreateParaviewDataCollection() const;
 
  void UpdateParaviewDataCollection(const std::string& paraview_output_dir) const;
 
  void initializeBasePhysicsStates(int cycle, double time);
 
  virtual std::unordered_map<std::string, FiniteElementState> getCheckpointedStates(int cycle) const;
 
  std::string name_ = {};
 
  std::string mesh_tag_ = {};
 
  mfem::ParMesh& mesh_;
 
  MPI_Comm comm_;
 
  std::vector<const serac::FiniteElementState*> states_;
 
  std::vector<const serac::FiniteElementState*> adjoints_;
 
  std::vector<const serac::FiniteElementDual*> duals_;
 
  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<serac::FiniteElementState> state;
 
    std::unique_ptr<serac::FiniteElementState> previous_state;
 
    std::unique_ptr<serac::FiniteElementDual> sensitivity;
  };
 
  std::vector<ParameterInfo> parameters_;
 
  FiniteElementState& shape_displacement_;
 
  std::unique_ptr<FiniteElementDual> shape_displacement_sensitivity_;
 
  mutable std::unordered_map<std::string, std::vector<serac::FiniteElementState>> checkpoint_states_;
 
  mutable std::unordered_map<std::string, serac::FiniteElementState> cached_checkpoint_states_;
 
  mutable std::optional<int> cached_checkpoint_cycle_;
 
  bool is_quasistatic_ = true;
 
  static constexpr int FLOAT_PRECISION_ = 8;
 
  double time_;
 
  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_;
 
  BoundaryConditionManager bcs_;
 
  bool checkpoint_to_disk_;
};
 
}  // namespace serac