diff --git a/doc/sphinx/source/examples/Hyperbolic/1D/Burgers1D.md b/doc/sphinx/source/examples/Hyperbolic/1D/Burgers1D.md index 181a656f..574f826a 100644 --- a/doc/sphinx/source/examples/Hyperbolic/1D/Burgers1D.md +++ b/doc/sphinx/source/examples/Hyperbolic/1D/Burgers1D.md @@ -18,4 +18,4 @@ The wave is allowed to propagate across the domain while the area under the curv This example is implemented in: - [MATLAB/ OCTAVE](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/burgers1D.m) -- [C++](https://github.com/csrc-sdsu/mole/blob/main/examples/cpp/Burgers1D.cpp) +- [C++](https://github.com/csrc-sdsu/mole/blob/main/examples/cpp/burgers1D.cpp) diff --git a/examples/cpp/Burgers1D.cpp b/examples/cpp/burgers1D.cpp similarity index 96% rename from examples/cpp/Burgers1D.cpp rename to examples/cpp/burgers1D.cpp index cb8b77d9..f0007793 100644 --- a/examples/cpp/Burgers1D.cpp +++ b/examples/cpp/burgers1D.cpp @@ -1,92 +1,92 @@ - -/** - * Solving the 1D Advection Equation using a Mimetic Finite Difference Scheme - * - * Equation: ∂U/∂t + ∂(U²)/∂x = 0 (Nonlinear Burgers' Equation in conservative form) - * Domain: x ∈ [-15, 15] with m = 300 grid cells - * Time: Simulated until t = 10.0 with time step dt = dx (CFL condition) - * Initial Condition: U(x,0) = exp(-x² / 50) - * Boundary Conditions: Mimetic divergence and interpolation operators applied (implicit treatment) - * - * Solution is computed using a staggered grid approach, explicit time-stepping, - * and mimetic finite difference operators for divergence and interpolation. - */ -#include -#include -#include // for EXIT_SUCCESS / EXIT_FAILURE -#include -#include -#include -#include -#include "mole.h" -#include "utils.h" - -int main() { - constexpr double west = -15.0; - constexpr double east = 15.0; - constexpr int k = 2; - constexpr int m = 300; - constexpr double t = 10.0; - - const double dx = (east - west) / m; - const double dt = dx; - - Divergence D(k, m, dx); - Interpol I(m, 1.0); - - // Spatial grid (including ghost cells) - arma::vec xgrid(m + 2); - xgrid(0) = west; - xgrid(m + 1) = east; - for (int i = 1; i <= m; ++i) { - xgrid(i) = west + (i - 0.5) * dx; - } - - // Initial condition - arma::vec U = arma::exp(-arma::square(xgrid) / 50.0); - - // Sanity check: matrix dimensions - if (D.n_cols != I.n_rows || I.n_cols != U.n_rows) { - std::cerr << "Error: Incompatible matrix dimensions!" << std::endl; - return EXIT_FAILURE; - } - - int total_steps = static_cast(t / dt); - int plot_interval = total_steps / 5; - - for (int step = 0; step <= total_steps; ++step) { - double time = step * dt; - - // Explicit update - U += (-dt / 2.0) * (D * (I * arma::square(U))); - - if (step % plot_interval == 0) { - double area = Utils::trapz(xgrid, U); - std::cout << "Time step: " << step - << ", Time: " << time - << ", Trapz Area: " << area - << ", U_min: " << U.min() - << ", U_max: " << U.max() - << ", U_center: " << U(U.n_elem / 2) - << std::endl; - - std::string filename = "output_step_" + std::to_string(step) + ".dat"; - std::ofstream outfile(filename); - if (!outfile) { - std::cerr << "Error: Could not open file for writing: " << filename << std::endl; - return EXIT_FAILURE; - } - - outfile << "# x U(x)\n"; - for (arma::uword i = 0; i < xgrid.n_elem; ++i) { - outfile << std::setw(12) << xgrid(i) << " " - << std::setw(12) << U(i) << "\n"; - } - } - } - - return EXIT_SUCCESS; -} - - - + +/** + * Solving the 1D Advection Equation using a Mimetic Finite Difference Scheme + * + * Equation: ∂U/∂t + ∂(U²)/∂x = 0 (Nonlinear Burgers' Equation in conservative form) + * Domain: x ∈ [-15, 15] with m = 300 grid cells + * Time: Simulated until t = 10.0 with time step dt = dx (CFL condition) + * Initial Condition: U(x,0) = exp(-x² / 50) + * Boundary Conditions: Mimetic divergence and interpolation operators applied (implicit treatment) + * + * Solution is computed using a staggered grid approach, explicit time-stepping, + * and mimetic finite difference operators for divergence and interpolation. + */ +#include +#include +#include // for EXIT_SUCCESS / EXIT_FAILURE +#include +#include +#include +#include +#include "mole.h" +#include "utils.h" + +int main() { + constexpr double west = -15.0; + constexpr double east = 15.0; + constexpr int k = 2; + constexpr int m = 300; + constexpr double t = 10.0; + + const double dx = (east - west) / m; + const double dt = dx; + + Divergence D(k, m, dx); + Interpol I(m, 1.0); + + // Spatial grid (including ghost cells) + arma::vec xgrid(m + 2); + xgrid(0) = west; + xgrid(m + 1) = east; + for (int i = 1; i <= m; ++i) { + xgrid(i) = west + (i - 0.5) * dx; + } + + // Initial condition + arma::vec U = arma::exp(-arma::square(xgrid) / 50.0); + + // Sanity check: matrix dimensions + if (D.n_cols != I.n_rows || I.n_cols != U.n_rows) { + std::cerr << "Error: Incompatible matrix dimensions!" << std::endl; + return EXIT_FAILURE; + } + + int total_steps = static_cast(t / dt); + int plot_interval = total_steps / 5; + + for (int step = 0; step <= total_steps; ++step) { + double time = step * dt; + + // Explicit update + U += (-dt / 2.0) * (D * (I * arma::square(U))); + + if (step % plot_interval == 0) { + double area = Utils::trapz(xgrid, U); + std::cout << "Time step: " << step + << ", Time: " << time + << ", Trapz Area: " << area + << ", U_min: " << U.min() + << ", U_max: " << U.max() + << ", U_center: " << U(U.n_elem / 2) + << std::endl; + + std::string filename = "output_step_" + std::to_string(step) + ".dat"; + std::ofstream outfile(filename); + if (!outfile) { + std::cerr << "Error: Could not open file for writing: " << filename << std::endl; + return EXIT_FAILURE; + } + + outfile << "# x U(x)\n"; + for (arma::uword i = 0; i < xgrid.n_elem; ++i) { + outfile << std::setw(12) << xgrid(i) << " " + << std::setw(12) << U(i) << "\n"; + } + } + } + + return EXIT_SUCCESS; +} + + + diff --git a/examples/cpp/Poisson2D.cpp b/examples/cpp/poisson2D.cpp similarity index 97% rename from examples/cpp/Poisson2D.cpp rename to examples/cpp/poisson2D.cpp index 6945b0cd..60d669c9 100644 --- a/examples/cpp/Poisson2D.cpp +++ b/examples/cpp/poisson2D.cpp @@ -1,67 +1,67 @@ -/** - * Solving the 2D Poisson Equation with Robin Boundary Conditions - * - * Equation: ∇²u = f(x, y) (Poisson Equation) - * Domain: Defined on a (m+2) x (n+2) grid with spacing dx, dy - * Boundary Conditions: - * - Bottom boundary (y = 0) has a Dirichlet condition: u = 100 - * - Other boundaries are subject to Robin conditions as defined in RobinBC - * - * Solution is computed using a Mimetic Finite Difference Laplacian and solved via Armadillo's sparse solver. - */ - -#include -#include "mole.h" -#include -#include -#include // For std::abs - -using namespace arma; - -int main() { - constexpr uint16_t k = 2; // Order of accuracy - constexpr uint32_t m = 5; // Vertical resolution - constexpr uint32_t n = 6; // Horizontal resolution - constexpr double dx = 1.0, dy = 1.0; // Grid spacing - - // Construct the 2D Mimetic Laplacian - Laplacian L(k, m, n, dx, dy); - RobinBC BC(k, m, dx, n, dy, 1.0, 0.0); - L = L + BC; - - // Define RHS matrix and apply boundary conditions - mat RHS = zeros(m + 2, n + 2); - RHS.row(0).fill(100.0); // Known value at the bottom boundary - - // Convert RHS to a column vector - vec rhs = vectorise(RHS); - - // Solve the system - vec SOL = spsolve(L, rhs); - - // Reshape solution back to 2D form - mat SOL2D = reshape(SOL, m + 2, n + 2); - - // Display solution without negative zeros or excessive decimal places - std::cout << "2D Poisson Solution:\n"; - for (uint32_t i = 0; i < SOL2D.n_rows; ++i) { - for (uint32_t j = 0; j < SOL2D.n_cols; ++j) { - double value = SOL2D(i, j); - if (std::abs(value) < 1e-10) { // If value is very close to zero, set it to exactly 0.0 - value = 0.0; - } - - // Adjust precision and remove unnecessary decimal places - if (std::abs(value - std::round(value)) < 1e-4) { // If value is close to an integer - std::cout << std::fixed << std::setprecision(0) << value; // No decimals - } else { - std::cout << std::fixed << std::setprecision(4) << value; // Four decimal places - } - - std::cout << "\t"; // Tab separation - } - std::cout << "\n"; - } - - return 0; -} +/** + * Solving the 2D Poisson Equation with Robin Boundary Conditions + * + * Equation: ∇²u = f(x, y) (Poisson Equation) + * Domain: Defined on a (m+2) x (n+2) grid with spacing dx, dy + * Boundary Conditions: + * - Bottom boundary (y = 0) has a Dirichlet condition: u = 100 + * - Other boundaries are subject to Robin conditions as defined in RobinBC + * + * Solution is computed using a Mimetic Finite Difference Laplacian and solved via Armadillo's sparse solver. + */ + +#include +#include "mole.h" +#include +#include +#include // For std::abs + +using namespace arma; + +int main() { + constexpr uint16_t k = 2; // Order of accuracy + constexpr uint32_t m = 5; // Vertical resolution + constexpr uint32_t n = 6; // Horizontal resolution + constexpr double dx = 1.0, dy = 1.0; // Grid spacing + + // Construct the 2D Mimetic Laplacian + Laplacian L(k, m, n, dx, dy); + RobinBC BC(k, m, dx, n, dy, 1.0, 0.0); + L = L + BC; + + // Define RHS matrix and apply boundary conditions + mat RHS = zeros(m + 2, n + 2); + RHS.row(0).fill(100.0); // Known value at the bottom boundary + + // Convert RHS to a column vector + vec rhs = vectorise(RHS); + + // Solve the system + vec SOL = spsolve(L, rhs); + + // Reshape solution back to 2D form + mat SOL2D = reshape(SOL, m + 2, n + 2); + + // Display solution without negative zeros or excessive decimal places + std::cout << "2D Poisson Solution:\n"; + for (uint32_t i = 0; i < SOL2D.n_rows; ++i) { + for (uint32_t j = 0; j < SOL2D.n_cols; ++j) { + double value = SOL2D(i, j); + if (std::abs(value) < 1e-10) { // If value is very close to zero, set it to exactly 0.0 + value = 0.0; + } + + // Adjust precision and remove unnecessary decimal places + if (std::abs(value - std::round(value)) < 1e-4) { // If value is close to an integer + std::cout << std::fixed << std::setprecision(0) << value; // No decimals + } else { + std::cout << std::fixed << std::setprecision(4) << value; // Four decimal places + } + + std::cout << "\t"; // Tab separation + } + std::cout << "\n"; + } + + return 0; +} diff --git a/examples/cpp/Schrodinger2D.cpp b/examples/cpp/schrodinger2D.cpp similarity index 100% rename from examples/cpp/Schrodinger2D.cpp rename to examples/cpp/schrodinger2D.cpp