Write explicit diffusion simulation

_code/pset_05/cpp/CMakeLists.txt

@@ -13,6 +13,11 @@ target_include_directories(sim PUBLIC
 )
 
 add_executable(waves
-    main.cpp
+    waves.cpp
 )
 target_link_libraries(waves sim clara)
+
+add_executable(diffusion
+    diffusion.cpp
+)
+target_link_libraries(diffusion sim)

_code/pset_05/cpp/diffusion.cpp

+#include <iostream>
+#include "glut_grapher.h"
+#include "scalar.h"
+#include "vector.h"
+
+//--------------------------------------------------------------------------------------------------
+GlutGrapher grapher;
+void display() { grapher.display(); }
+void idle() { grapher.idle(); }
+void mouse(int, int, int, int) { exit(0); }
+
+//--------------------------------------------------------------------------------------------------
+class DiffusionSim : public SimulationInterface {
+public:
+    DiffusionSim();
+    void initialize();
+    void step() final;
+    VectorX<Scalar> const& state() const final { return state_n; }
+
+    Scalar dt() const { return delta_t; }
+    Scalar dx() const { return delta_x; }
+
+private:
+    static constexpr uint32_t n = 1000;
+    Scalar const D = 0.0001;
+    Scalar const delta_t = 0.001;
+    Scalar const delta_x = 0.001;
+    Scalar const a1 = D * delta_t / (delta_x * delta_x);
+    Scalar const a2 = 1.0 - 2 * a1;
+
+    Scalar const boundary_0 = 0;
+    Scalar const boundary_1 = 0;
+
+    VectorX<Scalar> state_next_n;
+    VectorX<Scalar> state_n;
+    VectorX<Scalar> state_prev_n;
+};
+
+//..................................................................................................
+DiffusionSim::DiffusionSim() {
+    state_next_n.resize(n);
+    state_n.resize(n);
+    state_prev_n.resize(n);
+    state_next_n.setZero();
+    state_n.setZero();
+    state_prev_n.setZero();
+}
+
+//..................................................................................................
+void DiffusionSim::initialize() {
+    state_n.setZero();
+
+    // delta function
+    state_n[n/2] += 1.0;
+
+    // ensure boundary conditions are met
+    state_n[0] = boundary_0;
+    state_n[n - 1] = boundary_1;
+    state_prev_n = state_n;
+}
+
+//..................................................................................................
+void DiffusionSim::step() {
+    for (uint32_t i = 1; i < n - 1; ++i) {
+        state_next_n[i] = a2 * state_n[i] + a1 * (state_n[i - 1] + state_n[i + 1]);
+    }
+
+    state_prev_n.swap(state_n);
+    state_n.swap(state_next_n);
+}
+
+//--------------------------------------------------------------------------------------------------
+int main() {
+    DiffusionSim sim;
+    sim.initialize();
+    grapher.initialize(
+        &sim,
+        BoundingBox<Vector2<Scalar>>(
+            {0, -0.1},
+            {sim.dx() * sim.state().size(), 1.1}
+        )
+    );
+    grapher.init_glut();
+    glutDisplayFunc(display);
+    glutMouseFunc(mouse);
+    glutIdleFunc(idle);
+    grapher.run();
+
+    return 0;
+}