Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion doc/sphinx/source/examples/Hyperbolic/1D/Burgers1D.md
Original file line number Diff line number Diff line change
Expand Up @@ -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)
184 changes: 92 additions & 92 deletions examples/cpp/Burgers1D.cpp → examples/cpp/burgers1D.cpp
Original file line number Diff line number Diff line change
@@ -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 <armadillo>
#include <cmath>
#include <cstdlib> // for EXIT_SUCCESS / EXIT_FAILURE
#include <fstream>
#include <iostream>
#include <string>
#include <iomanip>
#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<int>(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 <armadillo>
#include <cmath>
#include <cstdlib> // for EXIT_SUCCESS / EXIT_FAILURE
#include <fstream>
#include <iostream>
#include <string>
#include <iomanip>
#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<int>(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;
}



134 changes: 67 additions & 67 deletions examples/cpp/Poisson2D.cpp → examples/cpp/poisson2D.cpp
Original file line number Diff line number Diff line change
@@ -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 <armadillo>
#include "mole.h"
#include <iomanip>
#include <iostream>
#include <cmath> // 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 <armadillo>
#include "mole.h"
#include <iomanip>
#include <iostream>
#include <cmath> // 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;
}
File renamed without changes.