Access boundary conditions from model

[ahojukka5]

We need to create a method to retrieve boundary condition data from the model. The challenge is that boundary conditions are defined in the simulator, and the model does not have direct access to them, so we need to employ a workaround.

One way to do this might be

#include <iostream>

template <class Model> class Simulator {
public:
  void say_hi() { std::cout << "Hi" << std::endl; }
  void step() {
    Model model;
    post_hook(*this, model);
  }
};

class Model {};
class Tungsten : public Model {};
class Aluminum : public Model {};

// default post hook is no-op and is called for all models which does not have a
// post hook defined
template <class Model> void post_hook(Simulator<Model> &simulator, Model &model) {
  std::cout << "Default post hook" << std::endl;
}

void post_hook(Simulator<Tungsten> &simulator, Tungsten &model) {
  std::cout << "Tungsten post hook" << std::endl;
  // Tungsten ask simulator to say hi
  simulator.say_hi();
}

int main() {
  Simulator<Model> simulator1;
  simulator1.step();

  Simulator<Tungsten> simulator2;
  simulator2.step();

  Simulator<Aluminum> simulator3;
  simulator3.step();
}

[ahojukka5]

This can also be done without using templates by introducing a free function called void step(Simulator &s, Model &m):

#include <iostream>

// models.h

class Model {
public:
  virtual ~Model() = default;
  virtual void step(double dt) = 0;
};

// simulator.h

class Simulator {
private:
  double time = 0.0;
  double dt = 0.1;

public:
  Model &m_model;
  Simulator(Model &model) : m_model(model) {}
  double get_time() { return time; }
  double get_dt() { return dt; }
  Model &get_model() { return m_model; }
};

void step(Simulator &s, Model &m) {
  m.step(s.get_dt());
}

// tungsten.h

class Tungsten : public Model {
public:
  void step(double dt) override { std::cout << "Tungsten step" << std::endl; }
};

// aluminum.h

class Aluminum : public Model {
public:
  void step(double dt) override { std::cout << "Aluminum step" << std::endl; }
};

void step(Simulator &s, Aluminum &m) {
  std::cout << "Aluminum pre-step" << std::endl;
  m.step(s.get_dt());
  std::cout << "Aluminum post-step" << std::endl;
}

// main.cpp

int main() {
  Tungsten tungsten;
  Aluminum aluminum;
  Simulator simulator1(tungsten);
  Simulator simulator2(aluminum);
  step(simulator1, aluminum);
  step(simulator2, tungsten);
  return 0;
}

As I work on this implementation, I'm starting to think that perhaps boundary conditions and initial conditions should be embedded within the model itself. We may find ourselves in a situation where we want to use one simulator with multiple models that interact and possibly have some interdependencies. Therefore, it would make sense for each model to define its boundary conditions and initial conditions.

cc: @suviraj1

[ahojukka5]

#include <iostream>
#include <memory>
#include <string>
#include <vector>

// fieldmodifier.h

class FieldModifier {
public:
  virtual ~FieldModifier() = default;
  virtual void apply() = 0;
  virtual const std::string get_name() = 0;
};

class FixedBC : public FieldModifier {
public:
  void apply() override { std::cout << "Applying FixedBC" << std::endl; }
  const std::string get_name() override { return "FixedBC"; }
};

class MovingBC : public FieldModifier {
private:
  double m_xpos = 1.0;

public:
  void apply() override { std::cout << "Applying MovingBC" << std::endl; }
  const std::string get_name() override { return "MovingBC"; }
  double get_xpos() { return m_xpos; }
};

// models.h

class Model {
private:
  std::vector<std::unique_ptr<FieldModifier>> m_initial_conditions;
  std::vector<std::unique_ptr<FieldModifier>> m_boundary_conditions;

public:
  virtual ~Model() = default;
  virtual void step(double t, double dt) = 0;
  const std::vector<std::unique_ptr<FieldModifier>> &get_initial_conditions() { return m_initial_conditions; }
  const std::vector<std::unique_ptr<FieldModifier>> &get_boundary_conditions() { return m_boundary_conditions; }
  void add_initial_condition(std::unique_ptr<FieldModifier> ic) { m_initial_conditions.push_back(std::move(ic)); }
  void add_boundary_condition(std::unique_ptr<FieldModifier> bc) { m_boundary_conditions.push_back(std::move(bc)); }
  void apply_boundary_conditions(double t, double dt) {
    for (auto &bc : m_boundary_conditions) bc->apply();
  }
  FieldModifier *get_boundary_condition(const std::string &name) {
    for (auto &bc : m_boundary_conditions) {
      if (bc->get_name() == name) return bc.get();
    }
    std::cout << "Warning: boundary condition " << name << " not found" << std::endl;
    return nullptr;
  }
};

// simulator.h

class Simulator {
private:
  double m_time = 0.0;
  double m_dt = 0.1;

public:
  double get_dt() { return m_dt; }
  void set_dt(double dt) { m_dt = dt; }
  double get_time() { return m_time; }
  void set_time(double time) { m_time = time; }
};

void step(Simulator &s, Model &m) {
  double t = s.get_time();
  double dt = s.get_dt();
  m.apply_boundary_conditions(t, dt);
  m.step(t, dt);
}

// tungsten.h

class Tungsten : public Model {
public:
  void step(double, double) override { std::cout << "Perform Tungsten step" << std::endl; }
};

// aluminum.h

class Aluminum : public Model {
private:
  double x_pos = 0.0;

public:
  void step(double, double) override { std::cout << "Perform Aluminum step" << std::endl; }
  double get_x_pos() { return x_pos; }
  void set_x_pos(double x) {
    std::cout << "Aluminum set x_pos to " << x << std::endl;
    x_pos = x;
  }
  void post_step(double, double) {
    std::cout << "Aluminum post-step" << std::endl;
    auto bc = dynamic_cast<MovingBC *>(get_boundary_condition("MovingBC"));
    if (bc) set_x_pos(bc->get_xpos());
    std::cout << "Aluminum post-step done" << std::endl;
  }
};

void step(Simulator &s, Aluminum &m) {
  std::cout << "Aluminum pre-step" << std::endl;
  double t = s.get_time();
  double dt = s.get_dt();
  m.step(t, dt);
  m.apply_boundary_conditions(t, dt);
  m.post_step(t, dt);
}

// main.cpp

int main() {
  Tungsten tungsten;
  Aluminum aluminum;
  Simulator simulator;
  // add FixedBC to tungsten and MovingBC to aluminum
  tungsten.add_boundary_condition(std::make_unique<FixedBC>());
  aluminum.add_boundary_condition(std::make_unique<MovingBC>());
  step(simulator, tungsten);
  step(simulator, aluminum);
  return 0;
}

/*
Applying FixedBC
Perform Tungsten step
Aluminum pre-step
Perform Aluminum step
Applying MovingBC
Aluminum post-step
Aluminum set x_pos to 1
Aluminum post-step done
*/

[ahojukka5]

Actually, in the latter case, there's no need to override void step(Simulator &s, Aluminum &m) as post_step can be called after the regular step, thus this is not changing API. So in this particular case we could just move boundary conditions and initial conditions inside the model and make potential API changes in the 0.2 version.

[ahojukka5]

#include <iostream>
#include <memory>
#include <string>
#include <vector>

// fieldmodifier.h

class IFieldModifier {
public:
  virtual ~IFieldModifier() = default;
  virtual void apply(IModel &model) = 0;
  virtual const std::string get_name() = 0;
};

class FixedBC : public IFieldModifier {
public:
  void apply(IModel &model) override { std::cout << "Applying FixedBC to " << model.get_name() << std::endl; }
  const std::string get_name() override { return "FixedBC"; }
};

class MovingBC : public IFieldModifier {
private:
  double m_xpos = 1.0;

public:
  void apply(IModel &model) override { std::cout << "Applying MovingBC" << model.get_name() << std::endl; }
  const std::string get_name() override { return "MovingBC"; }
  double get_xpos() { return m_xpos; }
};

// models.h

class IModel {
private:
  std::vector<std::unique_ptr<IFieldModifier>> m_initial_conditions;
  std::vector<std::unique_ptr<IFieldModifier>> m_boundary_conditions;

public:
  virtual ~IModel() = default;
  virtual void step(double t, double dt) = 0;
  virtual std::string get_name() = 0;
  const std::vector<std::unique_ptr<IFieldModifier>> &get_initial_conditions() { return m_initial_conditions; }
  const std::vector<std::unique_ptr<IFieldModifier>> &get_boundary_conditions() { return m_boundary_conditions; }
  void add_initial_condition(std::unique_ptr<IFieldModifier> ic) { m_initial_conditions.push_back(std::move(ic)); }
  void add_boundary_condition(std::unique_ptr<IFieldModifier> bc) { m_boundary_conditions.push_back(std::move(bc)); }
  void apply_boundary_conditions(double t, double dt) {
    for (auto &bc : m_boundary_conditions) bc->apply(*this);
  }
  IFieldModifier *get_boundary_condition(const std::string &name) {
    for (auto &bc : m_boundary_conditions) {
      if (bc->get_name() == name) return bc.get();
    }
    std::cout << "Warning: boundary condition " << name << " not found" << std::endl;
    return nullptr;
  }
};

// simulator.h

class Simulator {
private:
  double m_time = 0.0;
  double m_dt = 0.1;

public:
  double get_dt() { return m_dt; }
  void set_dt(double dt) { m_dt = dt; }
  double get_time() { return m_time; }
  void set_time(double time) { m_time = time; }
};

void step(Simulator &s, IModel &m) {
  double t = s.get_time();
  double dt = s.get_dt();
  m.apply_boundary_conditions(t, dt);
  m.step(t, dt);
}

// tungsten.h

class Tungsten : public IModel {
public:
  std::string get_name() override { return "Tungsten"; }
  void step(double, double) override { std::cout << "Perform Tungsten step" << std::endl; }
};

// aluminum.h

class Aluminum : public IModel {
private:
  double x_pos = 0.0;

public:
  std::string get_name() override { return "Aluminum"; }
  void step(double, double) override { std::cout << "Perform Aluminum step" << std::endl; }
  double get_x_pos() { return x_pos; }
  void set_x_pos(double x) {
    std::cout << "Aluminum set x_pos to " << x << std::endl;
    x_pos = x;
  }
  void post_step(double, double) {
    std::cout << "Aluminum post-step" << std::endl;
    auto bc = static_cast<MovingBC *>(get_boundary_condition("MovingBC"));
    if (bc) set_x_pos(bc->get_xpos());
    std::cout << "Aluminum post-step done" << std::endl;
  }
};

void step(Simulator &s, Aluminum &m) {
  std::cout << "Aluminum pre-step" << std::endl;
  double t = s.get_time();
  double dt = s.get_dt();
  m.step(t, dt);
  m.apply_boundary_conditions(t, dt);
  m.post_step(t, dt);
}

// main.cpp

int main() {
  Tungsten tungsten;
  Aluminum aluminum;
  Simulator simulator;
  // add FixedBC to tungsten and MovingBC to aluminum
  tungsten.add_boundary_condition(std::make_unique<FixedBC>());
  aluminum.add_boundary_condition(std::make_unique<MovingBC>());
  step(simulator, tungsten);
  step(simulator, aluminum);
  return 0;
}

The problem here is the circular dependency between IModel and IFieldModifier. Where it indeed would be possible to solve this problem using forward declaration, I don't like that idea much...