From eff6291e0b12b2fe11c68c25cf1c3db03e4452c6 Mon Sep 17 00:00:00 2001 From: Cory Date: Sun, 17 Dec 2023 13:33:06 +0100 Subject: [PATCH] Minimum, messy code --- include/WaveSimulation.hpp | 29 ++++-- include/constants.hpp | 17 ++++ include/literals.hpp | 19 ++++ include/utils.hpp | 10 ++- lib/WaveSimulation.cpp | 180 +++++++++++++++++++++++++++++++++++++ lib/literals.cpp | 17 ++++ lib/utils.cpp | 43 +++++++++ src/main.cpp | 120 +++++++++++++++++++++++++ src/test_suite.cpp | 4 + 9 files changed, 427 insertions(+), 12 deletions(-) create mode 100644 include/constants.hpp create mode 100644 include/literals.hpp create mode 100644 lib/literals.cpp diff --git a/include/WaveSimulation.hpp b/include/WaveSimulation.hpp index 39e51b0..0ef9acf 100644 --- a/include/WaveSimulation.hpp +++ b/include/WaveSimulation.hpp @@ -12,21 +12,34 @@ #ifndef __WAVE_SIMULATION__ #define __WAVE_SIMULATION__ +#include "constants.hpp" +#include "literals.hpp" + #include +#include class WaveSimulation { protected: - int M; - arma::cx_mat U; - arma::cx_mat V; - arma::cx_mat A; - arma::cx_mat B; + uint32_t M; + arma::sp_cx_mat B; double h; double dt; + double T; + public: - virtual void solve() = 0; - void build_A(arma::cx_vec A_vec); - void build_B(arma::cx_vec B_vec); + int32_t N; + arma::cx_mat V; + arma::cx_mat U; + arma::sp_cx_mat A; + WaveSimulation(double h, double dt, double T); + virtual void solve(std::ofstream& ofile); + void build_A(); + void build_B(); + void initialize_U(double x_c, double y_c, double sigma_x, double sigma_y, + double p_x, double p_y); + void write_U(std::ofstream &ofile); + void step(); + void build_V(double thickness, double pos_x, double aperture_sparation, double aperture, uint32_t slits); }; #endif diff --git a/include/constants.hpp b/include/constants.hpp new file mode 100644 index 0000000..4373d64 --- /dev/null +++ b/include/constants.hpp @@ -0,0 +1,17 @@ +/** @file constants.hpp + * + * @author Cory Alexander Balaton (coryab) + * @author Janita Ovidie Sandtrøen Willumsen (janitaws) + * + * @version 1.0 + * + * @brief Library of constants + * + * @bug No known bugs + * */ +#ifndef __CONST__ +#define __CONST__ + +#define I std::complex{0., 1.} + +#endif diff --git a/include/literals.hpp b/include/literals.hpp new file mode 100644 index 0000000..3a29d48 --- /dev/null +++ b/include/literals.hpp @@ -0,0 +1,19 @@ +/** @file literals.hpp + * + * @author Cory Alexander Balaton (coryab) + * @author Janita Ovidie Sandtrøen Willumsen (janitaws) + * + * @version 1.0 + * + * @brief Useful literals + * + * @bug No known bugs + * */ +#ifndef __LITERALS__ +#define __LITERALS__ + +#include + +std::complex operator ""_i(long double magnitude); + +#endif diff --git a/include/utils.hpp b/include/utils.hpp index 0618908..83f98a8 100644 --- a/include/utils.hpp +++ b/include/utils.hpp @@ -117,10 +117,10 @@ std::string dirname(const std::string &path); /** @brief Take 2 strings and concatenate them and make sure there is a * directory separator (/) between them. * - * @details This function doesn't care whether or not the values given as - * parameters are valid path strings. It is the responsibility of the user to make - * sure that the values given are valid path strings. - * The function only guarantees that the output string is a valid path string. + * @details This function doesn't care whether or not the values given as + * parameters are valid path strings. It is the responsibility of the user to + * make sure that the values given are valid path strings. The function only + * guarantees that the output string is a valid path string. * * @param left The left hand side of the result string * @param right The right hand side of the result string @@ -129,6 +129,8 @@ std::string dirname(const std::string &path); * */ std::string concatpath(const std::string &left, const std::string &right); +// A function that prints the structure of a sparse matrix to screen. +void print_sp_matrix_structure(const arma::sp_cx_mat &A); } // namespace utils #endif diff --git a/lib/WaveSimulation.cpp b/lib/WaveSimulation.cpp index e69de29..056cb05 100644 --- a/lib/WaveSimulation.cpp +++ b/lib/WaveSimulation.cpp @@ -0,0 +1,180 @@ +/** @file WaveSimulation.cpp + * + * @author Cory Alexander Balaton (coryab) + * @author Janita Ovidie Sandtrøen Willumsen (janitaws) + * + * @version 0.1 + * + * @brief Implementation of the WaveSimulation class. + * + * @bug No known bugs + * */ +#include "WaveSimulation.hpp" + +#include "utils.hpp" + +#include +#include +#include +#include + +WaveSimulation::WaveSimulation(double h, double dt, double T) +{ + this->dt = dt; + this->h = h; + this->T = T; + this->M = 1. / h; + this->N = M - 2; + this->V.set_size(this->N, this->N); + this->V.fill(0.); + this->U.set_size(this->N, this->N); + this->U.fill(0.); +} + +void WaveSimulation::solve(std::ofstream &ofile) +{ + ofile << this->N << '\n'; + uint32_t iterations = this->T / this->dt; + for (size_t i = 0; i < iterations; i++) { + this->write_U(ofile); + this->step(); + } +} + +void WaveSimulation::step() +{ + DEBUG("Inside step"); + arma::cx_vec tmp = this->B * this->U.as_col(); + arma::spsolve(this->U, this->A, tmp); +} + +void WaveSimulation::build_A() +{ + // Create the diagonal + arma::cx_vec diagonal(this->N * this->N); + + // Set diagonal values + std::complex r = (1._i * this->dt) / (2 * h * h); + for (size_t i = 0; i < diagonal.size(); i++) { + diagonal(i) = 1. + 4. * r + (1._i * this->dt / 2.) * this->V(i); + } + + // Create the submatrix + arma::cx_mat sub_matrix(this->N, this->N, arma::fill::zeros); + sub_matrix.diag(-1).fill(-r); + sub_matrix.diag(1).fill(-r); + + // Set the size of A + this->A.set_size(this->N * this->N, this->N * this->N); + + // Fill in the values in the submatrix diagonal + for (size_t i = 0; i < this->A.n_cols; i += this->N) { + this->A.submat(i, i, i + this->N - 1, i + this->N - 1) = sub_matrix; + } + + // Fill the last sub/sup-diagonals + this->A.diag() = diagonal; + this->A.diag(-this->N).fill(-r); + this->A.diag(this->N).fill(-r); +} + +void WaveSimulation::build_B() +{ + std::complex r = (1._i * this->dt) / (2 * h * h); + + // Create the diagonal + arma::cx_vec diagonal(this->N * this->N); + + for (size_t i = 0; i < diagonal.size(); i++) { + diagonal(i) = 1. - 4. * r - (1._i * this->dt / 2.) * this->V(i); + } + + // Create the submatrix + arma::cx_mat sub_matrix(this->N, this->N, arma::fill::zeros); + sub_matrix.diag(-1).fill(r); + sub_matrix.diag(1).fill(r); + + // Set the size of B + this->B.set_size(this->N * this->N, this->N * this->N); + + // Fill in the values in the submatrix diagonal + for (size_t i = 0; i < this->B.n_cols; i += this->N) { + this->B.submat(i, i, i + this->N - 1, i + this->N - 1) = sub_matrix; + } + + // Fill the last sub/sup-diagonals + this->B.diag() = diagonal; + this->B.diag(-this->N).fill(r); + this->B.diag(this->N).fill(r); +} + +void WaveSimulation::initialize_U(double x_c, double y_c, double sigma_x, + double sigma_y, double p_x, double p_y) +{ + double x, y, diff_x, diff_y; + std::complex sum = 0.; + for (size_t j = 0; j < this->U.n_cols; j++) { + x = j * h; + diff_x = x - x_c; + for (size_t i = 0; i < this->U.n_rows; i++) { + y = i * h; + diff_y = y - y_c; + this->U(i, j) = + std::exp(-(diff_x * diff_x) / (2. * sigma_x * sigma_x) + - (diff_y * diff_y) / (2. * sigma_y * sigma_y) + + p_x * x * 1._i + p_y * y * 1._i); + sum += this->U(i, j) * std::conj(this->U(i, j)); + } + } + if (std::abs(sum.imag()) > 1e-7) { + abort(); + } + double norm = 1. / std::sqrt(sum.real()); + + this->U.for_each([norm](std::complex &el) { el *= norm; }); +} + +void WaveSimulation::write_U(std::ofstream &ofile) +{ + this->U.for_each( + [&ofile](std::complex el) { ofile << el << '\t'; }); + ofile << '\n'; +} + +void WaveSimulation::build_V(double thickness, double pos_x, + double aperture_separation, double aperture, + uint32_t slits) +{ + uint32_t mid_y = this->N / 2 - (this->N % 2 == 0); + arma::cx_vec res; + + if (slits % 2 == 0) { + res = arma::cx_vec(aperture_separation/this->h,arma::fill::value(1e10)); + for (size_t i=0; i < slits; i+=2) { + res = arma::join_cols(res, arma::cx_vec(aperture/this->h,arma::fill::zeros)); + res = arma::join_cols(arma::cx_vec(aperture/this->h,arma::fill::zeros), res); + res = arma::join_cols(res, arma::cx_vec(aperture_separation/this->h,arma::fill::value(1e10))); + res = arma::join_cols(arma::cx_vec(aperture_separation/this->h,arma::fill::value(1e10)), res); + } + } + else { + res = arma::cx_vec(aperture/this->h,arma::fill::value(0)); + for (size_t i=0; i < slits-1; i+=2) { + res = arma::join_cols(res, arma::cx_vec(aperture_separation/this->h,arma::fill::value(1e10))); + res = arma::join_cols(arma::cx_vec(aperture_separation/this->h,arma::fill::value(1e10)), res); + res = arma::join_cols(res, arma::cx_vec(aperture/this->h,arma::fill::zeros)); + res = arma::join_cols(arma::cx_vec(aperture/this->h,arma::fill::zeros), res); + } + } + if (res.size() > this->N) { + abort(); + } + uint32_t fill = (this->N - res.size()) / 2; + res = arma::join_cols(arma::cx_vec(fill, arma::fill::value(1e10)), res); + res = arma::join_cols(res, arma::cx_vec(fill + ((this->N - res.size()) % 2), arma::fill::value(1e10))); + + uint32_t start = pos_x/this->h - thickness/this->h/2; + for (size_t i=0; i < thickness/this->h; i++) { + this->V.col(start+i) = res; + } +} diff --git a/lib/literals.cpp b/lib/literals.cpp new file mode 100644 index 0000000..57daf51 --- /dev/null +++ b/lib/literals.cpp @@ -0,0 +1,17 @@ +/** @file literals.cpp + * + * @author Cory Alexander Balaton (coryab) + * @author Janita Ovidie Sandtrøen Willumsen (janitaws) + * + * @version 1.0 + * + * @brief The implementation of the literals. + * + * @bug No known bugs + * */ +#include "literals.hpp" + +std::complex operator""_i(long double magnitude) +{ + return std::complex(0.,magnitude); +} diff --git a/lib/utils.cpp b/lib/utils.cpp index 554ee31..5dabba2 100644 --- a/lib/utils.cpp +++ b/lib/utils.cpp @@ -70,4 +70,47 @@ std::string concatpath(const std::string &left, const std::string &right) } } +void print_sp_matrix_structure(const arma::sp_cx_mat &A) +{ + using namespace std; + using namespace arma; + + // Declare a C-style 2D array of strings. + string S[A.n_rows][A.n_cols]; + + // Initialise all the strings to " ". + for (int i = 0; i < A.n_rows; i++) { + for (int j = 0; j < A.n_cols; j++) { + S[i][j] = " "; + } + } + + // Next, we want to set the string to a dot at each non-zero element. + // To do this we use the special loop iterator from the sp_cx_mat class + // to help us loop over only the non-zero matrix elements. + sp_cx_mat::const_iterator it = A.begin(); + sp_cx_mat::const_iterator it_end = A.end(); + + int nnz = 0; + for (; it != it_end; ++it) { + S[it.row()][it.col()] = "•"; + nnz++; + } + + // Finally, print the matrix to screen. + cout << endl; + for (int i = 0; i < A.n_rows; i++) { + cout << "| "; + for (int j = 0; j < A.n_cols; j++) { + cout << S[i][j] << " "; + } + cout << "|\n"; + } + + cout << endl; + cout << "matrix size: " << A.n_rows << "x" << A.n_cols << endl; + cout << "non-zero elements: " << nnz << endl; + cout << endl; +} + } // namespace utils diff --git a/src/main.cpp b/src/main.cpp index e69de29..079eae5 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -0,0 +1,120 @@ +/** @file main.cpp + * + * @author Cory Alexander Balaton (coryab) + * @author Janita Ovidie Sandtrøen Willumsen (janitaws) + * + * @version 1.0 + * + * @brief Implementation of the testing library + * + * @bug No known bugs + * */ +#include "WaveSimulation.hpp" +#include "utils.hpp" + +void probability_deviation() +{ + WaveSimulation sim_narrow(.005, 2.5e-5, .008); + sim_narrow.initialize_U(.25, .5, .05, .05, 200., 0.); + sim_narrow.build_A(); + sim_narrow.build_B(); + + WaveSimulation sim_wide(.005, 2.5e-5, .008); + sim_wide.build_V(.02, .5, .05, .05, 2); + sim_wide.initialize_U(.25, .5, .05, .10, 200., 0.); + sim_wide.build_A(); + sim_wide.build_B(); + + std::ofstream ofile; + utils::mkpath("data"); + ofile.open("data/probability_deviation.txt"); + + double sum_narrow, sum_wide; + + for (size_t i = 0; i < 320; i++) { + sum_narrow = 0; + sum_wide = 0; + for (size_t j = 0; j < sim_narrow.U.n_elem; j++) { + sum_narrow += (sim_narrow.U(j) * std::conj(sim_narrow.U(j))).real(); + sum_wide += (sim_wide.U(j) * std::conj(sim_wide.U(j))).real(); + } + + sim_narrow.step(); + sim_wide.step(); + ofile << i << ',' << sum_narrow << ',' << sum_wide << '\n'; + } + + ofile.close(); +} + +void color_map() +{ + WaveSimulation sim(.005, 2.5e-5, .008); + sim.build_V(.02, .5, .05, .05, 2); + sim.initialize_U(.25, .5, .05, .10, 200., 0.); + sim.build_A(); + sim.build_B(); + + std::ofstream ofile; + ofile.open("data/color_map.txt"); + ofile << sim.N << '\n'; + sim.write_U(ofile); + for (size_t i = 0; i < 40; i++) { + sim.step(); + } + sim.write_U(ofile); + for (size_t i = 0; i < 40; i++) { + sim.step(); + } + sim.write_U(ofile); + ofile.close(); +} + +void detector_screen() +{ + WaveSimulation sim(.005, 2.5e-5, .008); + sim.build_V(.02, .5, .05, .05, 1); + sim.initialize_U(.25, .5, .05, .10, 200., 0.); + sim.build_A(); + sim.build_B(); + + std::ofstream ofile; + ofile.open("data/screen/single_slit.txt"); + ofile << sim.N << '\n'; + for (size_t i = 0; i < 80; i++) { + sim.step(); + } + sim.write_U(ofile); + ofile.close(); + + sim.build_V(.02, .5, .05, .05, 2); + sim.initialize_U(.25, .5, .05, .10, 200., 0.); + + ofile.open("data/screen/double_slit.txt"); + ofile << sim.N << '\n'; + for (size_t i = 0; i < 80; i++) { + sim.step(); + } + sim.write_U(ofile); + ofile.close(); + + sim.build_V(.02, .5, .05, .05, 3); + sim.initialize_U(.25, .5, .05, .10, 200., 0.); + + ofile.open("data/screen/triple_slit.txt"); + ofile << sim.N << '\n'; + for (size_t i = 0; i < 80; i++) { + sim.step(); + } + sim.write_U(ofile); + ofile.close(); +} + +int main() +{ + probability_deviation(); + color_map(); + detector_screen(); + + return 0; +} diff --git a/src/test_suite.cpp b/src/test_suite.cpp index e69de29..905869d 100644 --- a/src/test_suite.cpp +++ b/src/test_suite.cpp @@ -0,0 +1,4 @@ +int main() +{ + return 0; +}