anderson_accelerator.hpp

Go to the documentation of this file.

 
#ifndef PECLET_FLOW_ANDERSON_ACCELERATOR_HPP
#define PECLET_FLOW_ANDERSON_ACCELERATOR_HPP
 
#include <array>
#include <chrono>
#include <cstddef>
#include <stdexcept>
#include <string>
#include <utility>
 
#include "flow_ibm.hpp"
#include "peclet/core/solver/anderson.hpp"
#ifdef PECLET_FLOW_MPI
#include "peclet/core/solver/anderson_mpi.hpp"
#endif
 
namespace peclet::flow {
 
template <class Grid>
class AndersonAccelerator {
 public:
  using Core = peclet::core::solver::AndersonCore;
 
  explicit AndersonAccelerator(Solver<Grid>& solver, int window = 5, double mixing = 1.0)
      : AndersonAccelerator(solver, solver.marchState(), window, mixing) {}
 
  AndersonAccelerator(const AndersonAccelerator&) = delete;
  AndersonAccelerator& operator=(const AndersonAccelerator&) = delete;
 
  void step(bool accelerate);
 
  // ---- status (forwarded from the core) ----
  const char* statusName() const { return core_.statusName(); }
  const std::string& reason() const { return core_.reason(); }
  double residual() const { return core_.residual(); }
  int numRestarts() const { return core_.numRestarts(); }
  int numResets() const { return core_.numResets(); }
  int numColumns() const { return core_.numColumns(); }
  int window() const { return core_.window(); }
  double mixing() const { return core_.mixing(); }
  std::size_t memoryBytes() const { return core_.memoryBytes(); }
  double seconds() const { return seconds_; }
  void reset() { core_.reset(); }
  void disable() { core_.disable(); }
  const Core& core() const { return core_; }
 
 private:
  AndersonAccelerator(Solver<Grid>& solver, typename Solver<Grid>::MarchState ms, int window,
                      double mixing);
  static Core makeCore(typename Solver<Grid>::MarchState& ms, int window, double mixing);
  static double now();
 
  Solver<Grid>& s_;
  std::array<double, 9> sig_;
  Core core_;
  double seconds_ = 0.0;
};
 
template <class Grid>
AndersonAccelerator<Grid>::AndersonAccelerator(Solver<Grid>& solver,
                                               typename Solver<Grid>::MarchState ms, int window,
                                               double mixing)
    : s_(solver), sig_(ms.signature), core_(makeCore(ms, window, mixing)) {}
 
template <class Grid>
typename AndersonAccelerator<Grid>::Core AndersonAccelerator<Grid>::makeCore(
    typename Solver<Grid>::MarchState& ms, int window, double mixing) {
  peclet::core::solver::AndersonState st;
  st.fields = ms.fields;
  st.roles = ms.roles;
  st.extent = {ms.e.x, ms.e.y, ms.e.z};
  st.ghost = ms.G;
#ifdef PECLET_FLOW_MPI
  if (ms.distributed)
    st.comm = peclet::core::solver::andersonComm(ms.comm);
#endif
  try {
    return Core(std::move(st), window, mixing);
  } catch (const std::runtime_error& e) {  // the history does not fit (§6.3), on every rank
    throw std::runtime_error(std::string("AndersonAccelerator: ") + e.what() +
                             "; pass accelerate=False");
  }
}
 
template <class Grid>
double AndersonAccelerator<Grid>::now() {
  return std::chrono::duration<double>(std::chrono::steady_clock::now().time_since_epoch()).count();
}
 
template <class Grid>
void AndersonAccelerator<Grid>::step(bool accelerate) {
  double t0 = now();
  // §5.3: the refusals are re-evaluated at every call (host flags only), so a feature enabled
  // after construction is refused at the next call.
  const auto ms = s_.marchState();
  // A different number of state fields is a configuration change (review R5): on the collocated
  // grid set_advection / the projected-face advection add or remove the face field uf/vf/wf. The
  // history belongs to the old map, so it is invalidated, and the core cannot mix a state it was
  // not built on, so the call is refused.
  const auto& held = core_.state().fields;
  if (held.size() != ms.fields.size()) {
    core_.invalidate();
    throw std::logic_error(
        "AndersonAccelerator: the march state changed from " + std::to_string(held.size()) +
        " to " + std::to_string(ms.fields.size()) +
        " fields since construction (a configuration change: set_advection on the collocated "
        "solver adds or removes the face field uf/vf/wf); construct a new accelerator");
  }
  // The core mixes the buffers it was built on; a redistribute reallocates them, after which the
  // accelerator would act on buffers the solver no longer reads.
  bool same = true;
  for (std::size_t f = 0; same && f < held.size(); ++f)
    same = held[f].data() == ms.fields[f].data() && held[f].extent(0) == ms.fields[f].extent(0);
  if (!same)
    throw std::logic_error(
        "AndersonAccelerator: the solver's state buffers were reallocated (a redistribute or a "
        "rebalance since construction); construct a new accelerator");
  // §4.3 step 0, the signature half: a parameter change is a new map (dt, rho, mu, F, and since
  // review R5 the advection settings, so staggered set_advection between calls resets too).
  if (core_.status() == Core::Status::Active && ms.signature != sig_) {
    core_.invalidate();
    sig_ = ms.signature;
  }
  core_.prepare(accelerate);
  Kokkos::fence();
  seconds_ += now() - t0;
  try {
    s_.step();
  } catch (const std::exception& e) {
    t0 = now();
    const bool absorbed = core_.stepFailed(e.what());
    seconds_ += now() - t0;
    if (!absorbed)
      throw;
    return;
  }
  t0 = now();
  core_.complete(s_.pressureSolveFailed());
  Kokkos::fence();
  seconds_ += now() - t0;
}
 
}  // namespace peclet::flow
 
// Compiled once per grid in src/flow_solver_{staggered,colocated}.cpp (QUALITY_PLAN G.8).
#ifndef PECLET_FLOW_INSTANTIATING
namespace peclet::flow {
extern template class AndersonAccelerator<Staggered>;
extern template class AndersonAccelerator<Colocated>;
}  // namespace peclet::flow
#endif
 
#endif  // PECLET_FLOW_ANDERSON_ACCELERATOR_HPP