Implement a discretized ADR system

simucene/CMakeLists.txt

 add_library(simucene_lib
-    my_class.cpp
-    my_class.h
+    matrix_builder.cpp
+    matrix_builder.h
     vector.h
 )
 target_include_directories(simucene_lib PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})

simucene/main.cpp

-#include "my_class.h"
-#include "vector.h"
+#include "scalar.h"
+#include "matrix_builder.h"
 #include <iostream>
+#include <iomanip>
 
 using namespace simucene;
 
+//--------------------------------------------------------------------------------------------------
+template <typename T>
+void print_vector(VectorX<T> const& vec) {
+    std::cout << std::fixed << std::setprecision(3);
+    for (uint32_t i = 0; i < vec.size() - 1; ++i) {
+        std::cout << std::setw(5) << vec[i] << ", ";
+    }
+    std::cout << vec[vec.size() - 1] << '\n';
+}
+
+//--------------------------------------------------------------------------------------------------
+SparseMatrix<Scalar> d_x_first_order(int32_t const dim) {
+    return MatrixBuilder(dim).stripe(0, -1).stripe(1, 1);
+}
+
+//--------------------------------------------------------------------------------------------------
+SparseMatrix<Scalar> d_x2_first_order(int32_t const dim) {
+    return MatrixBuilder(dim).stripe(-1, 1).stripe(0, -2).stripe(1, 1);
+}
+
+//--------------------------------------------------------------------------------------------------
+SparseMatrix<Scalar> diagonal_ones(int32_t const dim) {
+    return MatrixBuilder(dim).stripe(0, 1);
+}
+
 //--------------------------------------------------------------------------------------------------
 int main(int const, char const**) {
-    MyClass my_class("world");
-    my_class.hello();
+    // primary parameters of the simulation
+    Scalar const max_x = 1;
+    int32_t const dim = 20;
+    Scalar const h = 0.001;
+    Scalar const advection_coefficient = 0;
+    Scalar const diffusion_coefficient = 1;
+    VectorX<Scalar> reaction_coefficients = VectorX<Scalar>::Zero(dim);
+
+    // derived quantities
+    Scalar const d = max_x / dim;
+    Scalar const alpha = advection_coefficient * h / d;
+    Scalar const delta = diffusion_coefficient * h / (d * d);
+
+    // transport matrices (assumes dirichlet boundary conditions identically zero)
+    SparseMatrix<Scalar> const advection_matrix = alpha * d_x_first_order(dim);
+    SparseMatrix<Scalar> const diffusion_matrix = delta * d_x2_first_order(dim);
+    SparseMatrix<Scalar> const reaction_matrix = h * diagonal_ones(dim) * reaction_coefficients.asDiagonal();
+    SparseMatrix<Scalar> const transport_matrix = advection_matrix + diffusion_matrix + reaction_matrix;
+
+    // debug: print useful things
+    std::cout << "alpha: " << alpha << '\n';
+    std::cout << "delta: " << delta << '\n';
+    std::cout << "Advection\n";
+    std::cout << advection_matrix << '\n';
+    std::cout << "Diffusion\n";
+    std::cout << diffusion_matrix << '\n';
+    std::cout << "Reaction\n";
+    std::cout << reaction_matrix << '\n';
+    std::cout << "Transport\n";
+    std::cout << transport_matrix << '\n';
+
+    // define the initial state (not boundary conditions)
+    VectorX<Scalar> state = VectorX<Scalar>::Zero(dim);
+    for (int32_t i = dim / 2 - 2; i < dim / 2 + 2; ++i) {
+        state[i] = 1;
+    }
 
-    Vector3<double> x = {0, 1, 2};
-    x = exp(x.array());
-    std::cout << x << '\n';
+    // iterate
+    for (uint32_t i = 0; i < 100; ++i) {
+        std::cout << i << ": ";
+        print_vector(state);
+        state += transport_matrix * state;
+    }
+    std::cout << 100 << ": ";
+    print_vector(state);
 
     return 0;
 }

simucene/matrix_builder.cpp

+#include "matrix_builder.h"
+#include <iostream>
+
+namespace simucene {
+
+//..................................................................................................
+MatrixBuilder& MatrixBuilder::stripe(int32_t k, Scalar value) {
+    // If k is positive, we want to iterate from 0 to n_ - d (inclusive, exclusive).
+    // If k is negative, we want to iterate from -d to n_.
+    int32_t const begin = std::max(0, -k);
+    int32_t const end = std::min(n_ - k, n_);
+    for (int32_t i = begin; i < end; ++i) {
+        triplets_.emplace_back(i, i + k, value);
+    }
+    return *this;
+}
+
+//..................................................................................................
+MatrixBuilder::operator SparseMatrix<Scalar>() const {
+    SparseMatrix<Scalar> result(n_, n_);
+    result.setFromTriplets(triplets_.begin(), triplets_.end());
+    return result;
+}
+
+}

simucene/matrix_builder.h

+#ifndef SIMUCENE_MATRIX_BUILDER_H
+#define SIMUCENE_MATRIX_BUILDER_H
+
+#include "scalar.h"
+#include "vector.h"
+#include <vector>
+
+namespace simucene {
+
+//--------------------------------------------------------------------------------------------------
+class MatrixBuilder {
+public:
+    MatrixBuilder(int32_t n): n_(n) {}
+    // k = 0 means on the diagonal, 1 means one above the diagonal, -1 one below, etc.
+    MatrixBuilder& stripe(int32_t k, Scalar value);
+    operator SparseMatrix<Scalar>() const;
+
+private:
+    int32_t n_;
+    std::vector<Eigen::Triplet<Scalar>> triplets_;
+};
+
+}
+
+#endif

simucene/my_class.cpp

-#include "my_class.h"
-
-#include <iostream>
-
-namespace simucene {
-
-//..................................................................................................
-void MyClass::hello() const {
-    std::cout << "Hello, " << name_ << '\n';
-}
-
-}

simucene/my_class.h

-#ifndef SIMUCENE_MY_CLASS_H
-#define SIMUCENE_MY_CLASS_H
-
-#include <string>
-
-namespace simucene {
-
-//--------------------------------------------------------------------------------------------------
-class MyClass {
-public:
-    MyClass(std::string const& name): name_(std::move(name)) {}
-
-    void hello() const;
-
-private:
-    std::string name_;
-};
-
-}
-
-#endif

simucene/scalar.h

+#ifndef SIMUCENE_SCALAR_H
+#define SIMUCENE_SCALAR_H
+
+#include <Eigen/Core>
+
+namespace simucene {
+
+using Scalar = double;
+
+}
+
+#endif

simucene/vector.h

@@ -2,6 +2,7 @@
 #define SIMUCENE_VECTOR_H
 
 #include <Eigen/Core>
+#include <Eigen/Sparse>
 
 namespace simucene {
 
@@ -29,6 +30,14 @@ using Matrix3 = Eigen::Matrix<T, 3, 3>;
 template <typename T>
 using MatrixX = Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>;
 
+//--------------------------------------------------------------------------------------------------
+template <typename T>
+using MatrixX = Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>;
+
+//--------------------------------------------------------------------------------------------------
+template <typename T>
+using SparseMatrix = Eigen::SparseMatrix<T>;
+
 }
 
 #endif