diff --git a/doc/sphinx/source/examples/Navier-Stokes/Cylinder-Flow-2D.md b/doc/sphinx/source/examples/Navier-Stokes/Cylinder-Flow-2D.md new file mode 100644 index 00000000..ce025ed4 --- /dev/null +++ b/doc/sphinx/source/examples/Navier-Stokes/Cylinder-Flow-2D.md @@ -0,0 +1,132 @@ +# 2D Channel Flow with a Cylinder-Like Obstacle (Masked Region) + +This example simulates a 2D incompressible channel flow past a cylinder-like obstacle (implemented as an axis-aligned masked no-slip region) using MOLE mimetic operators and a fractional-step (projection / pressure-correction) method. At moderate Reynolds number, the solution exhibits a wake and may show vortex shedding depending on grid resolution, time step, and the obstacle mask representation. + +## Governing Equations + +We solve the incompressible Navier–Stokes equations in two dimensions: + +$$ +\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u}\cdot\nabla)\mathbf{u} += -\frac{1}{\rho}\nabla p + \nu\nabla^{2}\mathbf{u} +$$ + +$$ +\nabla\cdot\mathbf{u} = 0 +$$ + +where: +- $$\mathbf{u}=(u,v)$$ is the velocity field +- $$p$$ is the (dynamic) pressure +- $$\rho$$ is the density +- $$\nu$$ is the kinematic viscosity + +> Note: If you prefer using kinematic pressure $$\pi=p/\rho$$, the momentum equation becomes +> $$\partial_t\mathbf{u}+(\mathbf{u}\cdot\nabla)\mathbf{u}=-\nabla \pi + \nu\nabla^2\mathbf{u}.$$ + +## Domain and Initial/Boundary Conditions + +### Spatial/Temporal Domain + +The computational domain is a 2D channel: + +- $$x \in [0,8]$$ +- $$y \in [-1,1]$$ +- $$t \in [0, t_{\text{final}}]$$ with $$t_{\text{final}} = \texttt{tspan}$$ (default in code: `tspan = 32.0`) + +A “square cylinder” is represented by a masked block of cells located near + +$$ +x/L_x = \mathtt{cylin}\_\mathtt{pos} +$$ + +with a size controlled by `cylin_size` (default `1/10`). Inside this mask, velocity is forced to zero (no-slip). The mask is applied as an axis-aligned cell region (not a geometric immersed boundary). + +### Initial Conditions (t = 0) + +At $$t=0$$: +- $$u(x,y,0)=U_0$$ everywhere in the fluid region +- $$v(x,y,0)=0$$ +- the obstacle mask region is set to $$u=v=0$$ + +(Default in code: `U_init = 1.0`.) + +### Boundary Conditions + +Velocity boundary conditions: +- **Inlet (left)**: Dirichlet inflow with $$u = U_0$$ and $$v = 0$$ +- **Outlet (right)**: zero streamwise gradient (Neumann outflow) with $$\partial u/\partial x = 0$$ and $$\partial v/\partial x = 0$$ +- **Top and bottom walls**: no-slip with $$u=0$$ and $$v=0$$ +- **Obstacle mask**: no-slip enforced by setting $$u=v=0$$ inside the mask region after each time step + +Pressure boundary conditions (pressure Poisson step): +- **Outlet (right)**: Dirichlet reference pressure with $$p = 0$$ +- **Other boundaries**: homogeneous Neumann with $$\partial p/\partial n = 0$$ + +## Implementation Details + +### Projection (Fractional-Step) Method + +Each time step advances a predictor–corrector scheme to enforce incompressibility: + +1. **Prediction (intermediate velocity)** + - Compute the intermediate velocity $$\mathbf{u}^*$$ by advancing the momentum equation + - Nonlinear convection is advanced with AB2 (AB1 for the first step) + - Viscous diffusion is treated with Crank–Nicolson, leading to Helmholtz-type solves for $$u^*$$ and $$v^*$$ + +2. **Pressure Poisson solve** + +$$ +\nabla^2 p^{n+1} = \frac{\rho}{\Delta t}\nabla\cdot \mathbf{u}^* +$$ + +3. **Velocity correction** + +$$ +\mathbf{u}^{n+1} = \mathbf{u}^* - \frac{\Delta t}{\rho}\nabla p^{n+1} +$$ + +4. **Re-apply boundary conditions + obstacle mask** + - Re-enforce all velocity boundary conditions + - Set the masked obstacle region to $$u=0$$ and $$v=0$$. + +### Mimetic Spatial Operators (MOLE) + +The spatial discretization uses MOLE operators: +- Divergence operator $$D$$ +- Gradient operator $$G$$ +- Laplacian operator $$L = DG$$ + +along with interpolation operators to map between cell-centered and face-based quantities used in flux evaluations and in the projection step. + +## Output Products + +The solver writes the final cell-centered fields: +- `U_final.csv`, `V_final.csv`, `p_final.csv` + +A plotting script is provided to visualize these CSV outputs using gnuplot: +- `cylinder_flow_2D_plot.gnu` → produces `cylinder_flow_2D_plot_cpp.png` + +## Running the Example + +Assuming you already configured and built MOLE following the tutorial, run the example from the `examples/cpp` directory so that all outputs (CSV + plots) are written next to the source and plotting scripts: + +```bash +cd /examples/cpp +../../build/examples/cpp/cylinder_flow_2D +gnuplot cylinder_flow_2D_plot.gnu +``` + +## MATLAB/Octave Version + +A MATLAB/Octave implementation of the same channel-flow-with-masked-obstacle setup is also provided. It is configured with the same domain, grid resolution, time step, final time, Reynolds number, and boundary conditions as the C++ example, so the final fields are directly comparable. + +The corresponding output figures (generated in this same directory) are: +- `cylinder_flow_2D_plot_cpp.png` +- `cylinder_flow_2D_plot_matlab.png` + +### C++ result +![C++: final U/V/p fields](cylinder_flow_2D_plot_cpp.png) + +### MATLAB result +![MATLAB: final U/V/p fields](cylinder_flow_2D_plot_matlab.png) \ No newline at end of file diff --git a/doc/sphinx/source/examples/Navier-Stokes/cylinder_flow_2D_plot_cpp.png b/doc/sphinx/source/examples/Navier-Stokes/cylinder_flow_2D_plot_cpp.png new file mode 100644 index 00000000..017542f3 Binary files /dev/null and b/doc/sphinx/source/examples/Navier-Stokes/cylinder_flow_2D_plot_cpp.png differ diff --git a/doc/sphinx/source/examples/Navier-Stokes/cylinder_flow_2D_plot_matlab.png b/doc/sphinx/source/examples/Navier-Stokes/cylinder_flow_2D_plot_matlab.png new file mode 100644 index 00000000..d2d0da20 Binary files /dev/null and b/doc/sphinx/source/examples/Navier-Stokes/cylinder_flow_2D_plot_matlab.png differ diff --git a/doc/sphinx/source/examples/Navier-Stokes/index.md b/doc/sphinx/source/examples/Navier-Stokes/index.md index 6e1988f4..0012e3ab 100644 --- a/doc/sphinx/source/examples/Navier-Stokes/index.md +++ b/doc/sphinx/source/examples/Navier-Stokes/index.md @@ -7,4 +7,5 @@ The Navier-Stokes equations describe the motion of viscous fluid substances, for :caption: Contents Lock-Exchange +Cylinder-Flow-2D ``` diff --git a/examples/cpp/cylinder_flow_2D.cpp b/examples/cpp/cylinder_flow_2D.cpp new file mode 100644 index 00000000..14d25401 --- /dev/null +++ b/examples/cpp/cylinder_flow_2D.cpp @@ -0,0 +1,606 @@ +/** + * @file cylinder_flow_2D.cpp + * @brief Solves the 2D incompressible Navier–Stokes equations for channel flow past an obstacle + * using a projection (pressure-correction) method. + * + * This example demonstrates how to assemble and use MOLE mimetic operators (Laplacian, Divergence, + * Gradient, and interpolation operators) to advance the 2D incompressible Navier–Stokes equations + * in time on a structured Cartesian grid. The convective term is advanced with AB2 (AB1 on the first + * step) and viscous diffusion is treated with Crank–Nicolson. Incompressibility is enforced via a + * pressure Poisson solve and velocity correction. + * + * Equation: + * - Momentum: + * \f[ + * \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u}\cdot\nabla)\mathbf{u} + * = -\nabla p + \nu\nabla^{2}\mathbf{u} + * \f] + * - Incompressibility: + * \f[ + * \nabla\cdot\mathbf{u} = 0 + * \f] + * + * Domain: + * - Rectangular channel: \f$ x \in [0,8],\; y \in [-1,1] \f$. + * - A solid obstacle is imposed by masking a block of cells near \f$ x/L_x = 1/8 \f$ with size + * set by `cylin_size` (note: the obstacle is applied as an axis-aligned cell mask). + * + * Boundary Conditions: + * - Velocity \f$\mathbf{u}=(u,v)\f$: + * - Inlet (left): Dirichlet, \f$ u = U_{\mathrm{init}},\; v=0 \f$. + * - Outlet (right): Neumann (zero normal gradient) on velocity. + * - Walls (top/bottom): no-slip, \f$ u=v=0 \f$. + * - Obstacle: no-slip enforced by zeroing masked cells each step. + * - Pressure \f$p\f$: + * - Outlet (right): Dirichlet, \f$ p = 0 \f$ (reference). + * - Elsewhere: Neumann (zero normal gradient). + * + * Run: + * @code + * cd examples/cpp + * ../../build/examples/cpp/cylinder_flow_2D + * @endcode + * + * Outputs: + * - U_final.csv, V_final.csv, p_final.csv + * - cylinder_flow_2D_plot.gnu, cylinder_flow_2D_plot.png (if gnuplot is available) + */ + +#include + +// Optional Eigen (for reusable sparse factorization) +#if defined(__has_include) + #if __has_include() && __has_include() + #define HAS_EIGEN 1 + #include + #include + #else + #define HAS_EIGEN 0 + #endif +#else + #define HAS_EIGEN 0 +#endif + +// Mimetic operator library (in ./cpp) +#include "mole.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using arma::sp_mat; +using arma::vec; +using arma::uvec; +using arma::mat; + +// Reusable sparse linear solver +class SparseSolver { +public: + SparseSolver() = default; + + void factorize(const sp_mat& A) { + n_ = static_cast(A.n_rows); +#if HAS_EIGEN + eigen_A_.resize(static_cast(A.n_rows), static_cast(A.n_cols)); + std::vector> trip; + trip.reserve(A.n_nonzero); + for (auto it = A.begin(); it != A.end(); ++it) { + trip.emplace_back(static_cast(it.row()), static_cast(it.col()), static_cast(*it)); + } + eigen_A_.setFromTriplets(trip.begin(), trip.end()); + solver_.analyzePattern(eigen_A_); + solver_.factorize(eigen_A_); + if (solver_.info() != Eigen::Success) { + throw std::runtime_error("Eigen::SparseLU factorization failed"); + } +#else + A_ = A; // fallback +#endif + } + + vec solve(const vec& b) const { + if (static_cast(b.n_elem) != n_) { + throw std::runtime_error("SparseSolver::solve dimension mismatch"); + } +#if HAS_EIGEN + Eigen::Map rhs(b.memptr(), static_cast(b.n_elem)); + Eigen::VectorXd x = solver_.solve(rhs); + if (solver_.info() != Eigen::Success) { + throw std::runtime_error("Eigen::SparseLU solve failed"); + } + return vec(x.data(), static_cast(x.size())); +#else + vec x; + bool ok = arma::spsolve(x, A_, b); + if (!ok) { + throw std::runtime_error("arma::spsolve failed (no Eigen available)"); + } + return x; +#endif + } + + static bool using_eigen() { +#if HAS_EIGEN + return true; +#else + return false; +#endif + } + +private: + long long n_{0}; +#if HAS_EIGEN + Eigen::SparseMatrix eigen_A_; + Eigen::SparseLU, Eigen::COLAMDOrdering> solver_; +#else + sp_mat A_; +#endif +}; + +// Scalar BC wrapper +struct BCSide { + std::string type; // Dirichlet/Neumann/Robin + std::vector coeffs; // {a} or {b} or {a,b} + bool active{false}; +}; + +static BCSide make_side(double a, double b) { + BCSide s; + s.active = (a != 0.0) || (b != 0.0); + if (!s.active) { + s.type = "Dirichlet"; + s.coeffs = {0.0}; + return s; + } + if (a != 0.0 && b == 0.0) { + s.type = "Dirichlet"; + s.coeffs = {a}; + } else if (a == 0.0 && b != 0.0) { + s.type = "Neumann"; + s.coeffs = {b}; + } else { + s.type = "Robin"; + s.coeffs = {a, b}; + } + return s; +} + +static uvec bc_left_indices(unsigned m, unsigned n) { + uvec idx(n); + const unsigned nx = m + 2; + for (unsigned jj = 0; jj < n; ++jj) { + unsigned j = jj + 1; + idx(jj) = 0 + nx * j; + } + return idx; +} + +static uvec bc_right_indices(unsigned m, unsigned n) { + uvec idx(n); + const unsigned nx = m + 2; + const unsigned i = m + 1; + for (unsigned jj = 0; jj < n; ++jj) { + unsigned j = jj + 1; + idx(jj) = i + nx * j; + } + return idx; +} + +static uvec bc_bottom_indices(unsigned m, unsigned n) { + (void)n; + uvec idx(m + 2); + for (unsigned i = 0; i < m + 2; ++i) idx(i) = i; + return idx; +} + +static uvec bc_top_indices(unsigned m, unsigned n) { + uvec idx(m + 2); + const unsigned nx = m + 2; + const unsigned j = n + 1; + for (unsigned i = 0; i < m + 2; ++i) idx(i) = i + nx * j; + return idx; +} + +struct BCSystem { + sp_mat A_bc; + vec b0; + uvec rowsbc; +}; + +static BCSystem build_scalar_bc_system( + const sp_mat& A, + unsigned k, unsigned m, double dx, unsigned n, double dy, + const std::array& dc, const std::array& nc, + const std::array& v) { + + const arma::uword N = A.n_rows; + if (A.n_rows != A.n_cols) throw std::runtime_error("A must be square"); + + const BCSide left = make_side(dc[0], nc[0]); + const BCSide right = make_side(dc[1], nc[1]); + const BCSide bottom = make_side(dc[2], nc[2]); + const BCSide top = make_side(dc[3], nc[3]); + + MixedBC bc_op(static_cast(k), static_cast(m), dx, + static_cast(n), dy, + left.type, left.coeffs, + right.type, right.coeffs, + bottom.type, bottom.coeffs, + top.type, top.coeffs); + + std::vector pieces; + if (left.active) pieces.push_back(bc_left_indices(m, n)); + if (right.active) pieces.push_back(bc_right_indices(m, n)); + if (bottom.active) pieces.push_back(bc_bottom_indices(m, n)); + if (top.active) pieces.push_back(bc_top_indices(m, n)); + + uvec rowsbc; + if (!pieces.empty()) { + rowsbc = pieces[0]; + for (size_t i = 1; i < pieces.size(); ++i) rowsbc = arma::join_cols(rowsbc, pieces[i]); + rowsbc = arma::unique(rowsbc); + } else { + rowsbc.set_size(0); + } + + sp_mat P = arma::speye(N, N); + for (arma::uword ii = 0; ii < rowsbc.n_elem; ++ii) { + const arma::uword r = rowsbc(ii); + P(r, r) = 0.0; + } + + BCSystem out; + out.A_bc = P * A + static_cast(bc_op); + out.b0 = arma::zeros(N); + + if (left.active) out.b0.elem(bc_left_indices(m, n)) = v[0]; + if (right.active) out.b0.elem(bc_right_indices(m, n)) = v[1]; + if (bottom.active) out.b0.elem(bc_bottom_indices(m, n)) = v[2]; + if (top.active) out.b0.elem(bc_top_indices(m, n)) = v[3]; + + out.rowsbc = rowsbc; + return out; +} + +// enforce inlet/outlet/walls + cylinder +static void applyVelocityBCAndMask( + mat& U, mat& V, + double Uin, + int i1, int i2, int j1, int j2) { + + const arma::uword nx = U.n_rows; + const arma::uword ny = U.n_cols; + + U.row(0).fill(Uin); + V.row(0).zeros(); + + U.row(nx - 1) = U.row(nx - 2); + V.row(nx - 1) = V.row(nx - 2); + + U.submat(1, 0, nx - 1, 0).zeros(); + V.submat(1, 0, nx - 1, 0).zeros(); + U.submat(1, ny - 1, nx - 1, ny - 1).zeros(); + V.submat(1, ny - 1, nx - 1, ny - 1).zeros(); + + U(0, 0) = 0.0; U(0, ny - 1) = 0.0; + V(0, 0) = 0.0; V(0, ny - 1) = 0.0; + + U.submat(i1, j1, i2, j2).zeros(); + V.submat(i1, j1, i2, j2).zeros(); +} + +static void write_gnuplot_script_cylinder_flow(const std::string& filename) { + const char* script = R"GNUPLOT( +set datafile separator "," +set term pngcairo size 1200,1500 +set output "cylinder_flow_2D_plot.png" + +unset key +set view map +set pm3d map +set tics out +set border lw 1 +set size ratio -1 + +set multiplot layout 3,1 rowsfirst + +set lmargin 8 +set rmargin 8 +set tmargin 1 +set bmargin 1 + +tcmd(f) = sprintf("awk -F',' '{line=$0; sub(/,+$/,\"\",line); nf=split(line,a,\",\"); for(i=1;i<=nf;i++) A[NR,i]=a[i]; if(nf>max) max=nf} END{for(i=1;i<=max;i++){for(j=1;j<=NR;j++){printf \"%%s%%s\", A[j,i], (j(m); + const double dy = (y_end - y_start) / static_cast(n); + + // Obstacle (cylinder) geometry parameters + const double cylin_pos = 1.0 / 8.0; + const double cylin_size = 1.0 / 10.0; + + // Physical parameters + const double rho = 1.0; + const double D0 = 2.0 * cylin_size; + const double U_init = 1.0; + const double nu = U_init * D0 / Re; + + const arma::uword nx = m + 2; + const arma::uword ny = n + 2; + const arma::uword Ncell = nx * ny; + const arma::uword Nfaces_x = (m + 1) * n; + const arma::uword Nfaces_y = m * (n + 1); + const arma::uword Nfaces = Nfaces_x + Nfaces_y; + + std::cout << "Ncell=" << Ncell << ", Nfaces=" << Nfaces + << " (" << Nfaces_x << "+" << Nfaces_y << ")\n"; + +#if !HAS_EIGEN + std::cerr << "[warn] Eigen headers not found; using arma::spsolve.\n"; +#endif + + // Construct mimetic operators (MOLE) + Laplacian L(static_cast(k), static_cast(m), static_cast(n), dx, dy); + Divergence D(static_cast(k), static_cast(m), static_cast(n), dx, dy); + Gradient G(static_cast(k), static_cast(m), static_cast(n), dx, dy); + + Interpol I_cf(static_cast(m), static_cast(n), 0.5, 0.5); + Interpol Ix_fc(true, static_cast(m), 0.5); + Interpol Iy_fc(true, static_cast(n), 0.5); + + sp_mat Im(nx, m); + Im.submat(1, 0, m, m - 1) = arma::speye(m, m); + + sp_mat In(ny, n); + In.submat(1, 0, n, n - 1) = arma::speye(n, n); + + sp_mat Sx = Utils::spkron(In, static_cast(Ix_fc)); + sp_mat Sy = Utils::spkron(static_cast(Iy_fc), Im); + + // AB2 + CN diffusion + sp_mat Icell = arma::speye(Ncell, Ncell); + sp_mat M = Icell - (0.5 * dt * nu) * static_cast(L); + sp_mat Mp = Icell + (0.5 * dt * nu) * static_cast(L); + + // Initial conditions / fields + mat U(nx, ny, arma::fill::ones); + U *= U_init; + mat V(nx, ny, arma::fill::zeros); + + const int m_unit = static_cast(std::floor(cylin_pos * static_cast(m))); + const int halfN1 = static_cast(0.5 * static_cast(n + 3)); + const int rad = static_cast(std::floor(cylin_size * static_cast(m_unit))); + + const int i1 = (m_unit - rad) - 1; + const int i2 = (m_unit + rad) - 1; + const int j1 = (halfN1 - rad) - 1; + const int j2 = (halfN1 + rad) - 1; + + U.submat(i1, j1, i2, j2).zeros(); + V.submat(i1, j1, i2, j2).zeros(); + + vec U_flat = arma::vectorise(U); + vec V_flat = arma::vectorise(V); + + vec AdvU_prev(Ncell, arma::fill::zeros); + vec AdvV_prev(Ncell, arma::fill::zeros); + vec p_new_flat(Ncell, arma::fill::zeros); + + // Velocity Helmholtz operator boundary conditions + const std::array dcU = {1, 0, 1, 1}; + const std::array ncU = {0, 1, 0, 0}; + std::array vU = { + vec(n, arma::fill::ones), + vec(n, arma::fill::zeros), + vec(nx, arma::fill::zeros), + vec(nx, arma::fill::zeros) + }; + + const std::array dcV = {1, 0, 1, 1}; + const std::array ncV = {0, 1, 0, 0}; + std::array vV = { + vec(n, arma::fill::zeros), + vec(n, arma::fill::zeros), + vec(nx, arma::fill::zeros), + vec(nx, arma::fill::zeros) + }; + + BCSystem bcU = build_scalar_bc_system(M, k, m, dx, n, dy, dcU, ncU, vU); + BCSystem bcV = build_scalar_bc_system(M, k, m, dx, n, dy, dcV, ncV, vV); + + SparseSolver Au_solver, Av_solver; + Au_solver.factorize(bcU.A_bc); + Av_solver.factorize(bcV.A_bc); + + // Pressure Poisson operator boundary conditions + const std::array dcP = {0, 1, 0, 0}; + const std::array ncP = {1, 0, 1, 1}; + std::array vP = { + vec(n, arma::fill::zeros), + vec(n, arma::fill::zeros), + vec(nx, arma::fill::zeros), + vec(nx, arma::fill::zeros) + }; + + BCSystem bcP = build_scalar_bc_system(static_cast(L), k, m, dx, n, dy, dcP, ncP, vP); + SparseSolver Ap_solver; + Ap_solver.factorize(bcP.A_bc); + + const int nSteps = static_cast(std::llround(tspan / dt)); + const int plotEvery = 100; + + // Time integration loop + for (int step = 1; step <= nSteps; ++step) { + vec U_faces = static_cast(I_cf) * U_flat; + vec V_faces = static_cast(I_cf) * V_flat; + + vec U_on_u = U_faces.head(Nfaces_x); + vec U_on_v = U_faces.tail(Nfaces_y); + + vec V_on_u = V_faces.head(Nfaces_x); + vec V_on_v = V_faces.tail(Nfaces_y); + + vec UU_on_u = U_on_u % U_on_u; + vec UV_on_u = U_on_u % V_on_u; + + vec VV_on_v = V_on_v % V_on_v; + vec UV_on_v = U_on_v % V_on_v; + + vec u_div = arma::join_cols(UU_on_u, UV_on_v); + vec v_div = arma::join_cols(UV_on_u, VV_on_v); + + vec AdvU_n = static_cast(D) * u_div; + vec AdvV_n = static_cast(D) * v_div; + + vec AdvU_ab, AdvV_ab; + if (step == 1) { + AdvU_ab = AdvU_n; + AdvV_ab = AdvV_n; + } else { + AdvU_ab = 1.5 * AdvU_n - 0.5 * AdvU_prev; + AdvV_ab = 1.5 * AdvV_n - 0.5 * AdvV_prev; + } + + vec rhsU = Mp * U_flat - dt * AdvU_ab; + vec rhsV = Mp * V_flat - dt * AdvV_ab; + + rhsU.elem(bcU.rowsbc).zeros(); + rhsV.elem(bcV.rowsbc).zeros(); + rhsU += bcU.b0; + rhsV += bcV.b0; + + vec U_star_flat = Au_solver.solve(rhsU); + vec V_star_flat = Av_solver.solve(rhsV); + + mat U_star = arma::reshape(U_star_flat, nx, ny); + mat V_star = arma::reshape(V_star_flat, nx, ny); + + U_star.submat(i1, j1, i2, j2).zeros(); + V_star.submat(i1, j1, i2, j2).zeros(); + U_star(0, 0) = 0.0; U_star(0, ny - 1) = 0.0; + V_star(0, 0) = 0.0; V_star(0, ny - 1) = 0.0; + + U_star_flat = arma::vectorise(U_star); + V_star_flat = arma::vectorise(V_star); + + vec U_star_faces = static_cast(I_cf) * U_star_flat; + vec V_star_faces = static_cast(I_cf) * V_star_flat; + + vec U_star_on_u = U_star_faces.head(Nfaces_x); + vec V_star_on_v = V_star_faces.tail(Nfaces_y); + + vec UV_star_div = arma::join_cols(U_star_on_u, V_star_on_v); + vec RHS = (rho / dt) * (static_cast(D) * UV_star_div); + + RHS.elem(bcP.rowsbc).zeros(); + RHS += bcP.b0; + + p_new_flat = Ap_solver.solve(RHS); + + vec gradp_faces = static_cast(G) * p_new_flat; + vec gradp_x = gradp_faces.head(Nfaces_x); + vec gradp_y = gradp_faces.tail(Nfaces_y); + + vec u_corr = Sx * gradp_x; + vec v_corr = Sy * gradp_y; + + vec U_new_flat = U_star_flat - (dt / rho) * u_corr; + vec V_new_flat = V_star_flat - (dt / rho) * v_corr; + + mat U_new = arma::reshape(U_new_flat, nx, ny); + mat V_new = arma::reshape(V_new_flat, nx, ny); + + applyVelocityBCAndMask(U_new, V_new, U_init, i1, i2, j1, j2); + + U_flat = arma::vectorise(U_new); + V_flat = arma::vectorise(V_new); + AdvU_prev = AdvU_n; + AdvV_prev = AdvV_n; + + if ((step % plotEvery) == 0 || step == 1 || step == nSteps) { + const double maxU = arma::abs(U_new).max(); + const double maxV = arma::abs(V_new).max(); + const double CFL = dt * (maxU / dx + maxV / dy); + const double inletMean = arma::mean(U_new.row(0)); + + std::cout << "step " << std::setw(6) << step << "/" << std::setw(6) << nSteps + << " | t=" << std::fixed << std::setprecision(6) << (dt * step) + << " | CFL~" << std::setprecision(3) << CFL + << " | max|U|=" << std::scientific << std::setprecision(3) << maxU + << " | max|V|=" << std::scientific << std::setprecision(3) << maxV + << " | mean(U_in)=" << std::fixed << std::setprecision(3) << inletMean + << "\n"; + } + } + + // Output: save fields + plot via gnuplot (if available) + { + mat U_final = arma::reshape(U_flat, nx, ny); + mat V_final = arma::reshape(V_flat, nx, ny); + mat p_final = arma::reshape(p_new_flat, nx, ny); + + U_final.save("U_final.csv", arma::csv_ascii); + V_final.save("V_final.csv", arma::csv_ascii); + p_final.save("p_final.csv", arma::csv_ascii); + + std::cout << "Wrote U_final.csv, V_final.csv, p_final.csv\n"; + + try { + const std::string gnuFile = "cylinder_flow_2D_plot.gnu"; + write_gnuplot_script_cylinder_flow(gnuFile); + run_gnuplot_script(gnuFile); + std::cout << "Wrote cylinder_flow_2D_plot.png\n"; + } catch (const std::exception& e) { + std::cerr << "[warn] Plotting skipped: " << e.what() << "\n"; + } + } + + return 0; +} diff --git a/examples/matlab_octave/cylinder_flow_2D.m b/examples/matlab_octave/cylinder_flow_2D.m new file mode 100644 index 00000000..d0b765dd --- /dev/null +++ b/examples/matlab_octave/cylinder_flow_2D.m @@ -0,0 +1,252 @@ +%% 2D incompressible flow over a square cylinder +% Projection method: AB2 (advection) + CN (diffusion) +% MATLAB version aligned with the C++ example settings/output plotting. + +close all; clear; clc; +addpath('../../src/matlab_octave') + +%% Settings (match C++) +Re = 200; % C++: 200 +k = 2; % C++: 2 +tspan = 32; % C++: 32.0 +dt = 0.005; % C++: 0.005 + +%% Domain and grid (match C++) +x_start = 0; x_end = 8; +y_start = -1; y_end = 1; + +m = 481; +n = 121; + +dx = (x_end - x_start)/m; +dy = (y_end - y_start)/n; + +% Cell-centered "with ghost" coords (length m+2, n+2) +xgrid = [0 dx/2:dx:(x_end-x_start)-dx/2 (x_end-x_start)]; +ygrid = [0 dy/2:dy:(y_end-y_start)-dy/2 (y_end-y_start)]; +[Y, X] = meshgrid(ygrid, xgrid); + +%% Obstacle (match C++) +cylin_pos = 1/8; +cylin_size = 1/10; + +%% Physical parameters (match C++ formula) +rho = 1; +D0 = 2*cylin_size; +U_init = 1; +nu = U_init * D0 / Re; + +%% Mimetic operators (MOLE) +L = lap2D(k, m, dx, n, dy); +D = div2D(k, m, dx, n, dy); +G = grad2D(k, m, dx, n, dy); + +I = interpolCentersToFacesD2D(k, m, n); +II = interpolFacesToCentersG2D(k, m, n); + +Ncell = (m+2)*(n+2); +Icell = speye(Ncell); + +% CN diffusion matrices +M = (Icell - 0.5*dt*nu*L); +Mp = (Icell + 0.5*dt*nu*L); + +%% Initial conditions (match C++) +U = 0.*X + U_init; +V = 0.*X; + +% cylinder mask indices +m_unit = floor(cylin_pos*m); +halfN1 = 0.5*(n+3); % n odd => integer +rad = floor(cylin_size*m_unit); + +i1 = m_unit - rad; +i2 = m_unit + rad; +j1 = halfN1 - rad; +j2 = halfN1 + rad; + +U(i1:i2, j1:j2) = 0; +V(i1:i2, j1:j2) = 0; + +U_flat = U(:); +V_flat = V(:); + +AdvU_prev = zeros(size(U_flat)); +AdvV_prev = zeros(size(V_flat)); + +p_new_flat = zeros(Ncell,1); + +%% Helmholtz BCs for U*, V* +dcU = [1;0;1;1]; ncU = [0;1;0;0]; +vU = {ones(n,1); zeros(n,1); zeros(m+2,1); zeros(m+2,1)}; + +dcV = [1;0;1;1]; ncV = [0;1;0;0]; +vV = {zeros(n,1); zeros(n,1); zeros(m+2,1); zeros(m+2,1)}; + +[Au, bU0] = addScalarBC2D(M, zeros(Ncell,1), k, m, dx, n, dy, dcU, ncU, vU); +[Av, bV0] = addScalarBC2D(M, zeros(Ncell,1), k, m, dx, n, dy, dcV, ncV, vV); + +diffU = Au - M; +diffV = Av - M; +rowsbcU = find(sum(spones(diffU),2) ~= 0); +rowsbcV = find(sum(spones(diffV),2) ~= 0); + +Au_fac = decomposition(Au, 'lu'); +Av_fac = decomposition(Av, 'lu'); + +%% Pressure Poisson BCs +dcP = [0;1;0;0]; ncP = [1;0;1;1]; +vP = {zeros(n,1); zeros(n,1); zeros(m+2,1); zeros(m+2,1)}; + +[Ap, bP0] = addScalarBC2D(L, zeros(Ncell,1), k, m, dx, n, dy, dcP, ncP, vP); +diffP = Ap - L; +rowsbcP = find(sum(spones(diffP),2) ~= 0); + +Ap_fac = decomposition(Ap, 'lu'); + +%% Time integration (match C++ structure) +nSteps = round(tspan/dt); +plotEvery = 100; + +for t_step = 1:nSteps + + % Interpolate cell-centered U/V to faces + U_stag = I * [U_flat; U_flat]; + U_on_u = U_stag(1:(m+1)*n); + U_on_v = U_stag((m+1)*n+1:end); + + V_stag = I * [V_flat; V_flat]; + V_on_u = V_stag(1:(m+1)*n); + V_on_v = V_stag((m+1)*n+1:end); + + % Nonlinear fluxes on faces + UU_on_u = U_on_u .* U_on_u; + UV_on_u = U_on_u .* V_on_u; + + VV_on_v = V_on_v .* V_on_v; + UV_on_v = U_on_v .* V_on_v; + + u_div = [UU_on_u; UV_on_v]; + v_div = [UV_on_u; VV_on_v]; + + AdvU_n = D * u_div; + AdvV_n = D * v_div; + + % AB2 (AB1 on first step) + if t_step == 1 + AdvU_ab = AdvU_n; + AdvV_ab = AdvV_n; + else + AdvU_ab = 1.5*AdvU_n - 0.5*AdvU_prev; + AdvV_ab = 1.5*AdvV_n - 0.5*AdvV_prev; + end + + % Helmholtz solves (CN diffusion) with BCs + rhsU = Mp*U_flat - dt*AdvU_ab; + rhsV = Mp*V_flat - dt*AdvV_ab; + + rhsU_bc = rhsU; rhsU_bc(rowsbcU) = 0; rhsU_bc = rhsU_bc + bU0; + rhsV_bc = rhsV; rhsV_bc(rowsbcV) = 0; rhsV_bc = rhsV_bc + bV0; + + U_star_flat = Au_fac \ rhsU_bc; + V_star_flat = Av_fac \ rhsV_bc; + + U_star = reshape(U_star_flat, m+2, n+2); + V_star = reshape(V_star_flat, m+2, n+2); + + % Mask + corner consistency + U_star(i1:i2, j1:j2) = 0; + V_star(i1:i2, j1:j2) = 0; + + U_star(1,1) = 0; U_star(1,end) = 0; + V_star(1,1) = 0; V_star(1,end) = 0; + + U_star_flat = U_star(:); + V_star_flat = V_star(:); + + % Pressure Poisson RHS from u* + U_star_stag = I * [U_star_flat; U_star_flat]; + U_star_on_u = U_star_stag(1:(m+1)*n); + + V_star_stag = I * [V_star_flat; V_star_flat]; + V_star_on_v = V_star_stag((m+1)*n+1:end); + + UV_star_div = [U_star_on_u; V_star_on_v]; + RHS = (rho/dt) * (D * UV_star_div); + + RHS_bc = RHS; RHS_bc(rowsbcP) = 0; RHS_bc = RHS_bc + bP0; + p_new_flat = Ap_fac \ RHS_bc; + + % Projection correction + U_V_flat = [U_star_flat; V_star_flat] - (dt/rho) * (II * G * p_new_flat); + + U_new = reshape(U_V_flat(1:Ncell), m+2, n+2); + V_new = reshape(U_V_flat(Ncell+1:end), m+2, n+2); + + % CRITICAL: re-apply velocity BCs AFTER projection + [U_new, V_new] = applyVelocityBCAndMask(U_new, V_new, U_init, i1,i2,j1,j2); + + % Update + U_flat = U_new(:); + V_flat = V_new(:); + + AdvU_prev = AdvU_n; + AdvV_prev = AdvV_n; + + if (mod(t_step, plotEvery) == 0) || (t_step == 1) || (t_step == nSteps) + maxU = max(abs(U_new(:))); + maxV = max(abs(V_new(:))); + CFL = dt * (maxU/dx + maxV/dy); + inletMean = mean(U_new(1,:)); + fprintf('step %6d/%6d | t=%.6f | CFL~%.3f | max|U|=%.3e | max|V|=%.3e | mean(U_in)=%.3f\n', ... + t_step, nSteps, dt*t_step, CFL, maxU, maxV, inletMean); + end +end + +%% Final plot: U, V, p in one column (3x1) +U = reshape(U_flat, m+2, n+2); +V = reshape(V_flat, m+2, n+2); +p = reshape(p_new_flat, m+2, n+2); + +figure('Name','Final U/V/p','NumberTitle','off'); +set(gcf, 'WindowState', 'maximized'); + +subplot(3,1,1) +imagesc(U'); axis image; axis tight; colorbar; +xlabel('x (grid index)'); ylabel('y (grid index)'); +title(sprintf('U at t = %.3f', tspan)) + +subplot(3,1,2) +imagesc(V'); axis image; axis tight; colorbar; +xlabel('x (grid index)'); ylabel('y (grid index)'); +title(sprintf('V at t = %.3f', tspan)) + +subplot(3,1,3) +imagesc(p'); axis image; axis tight; colorbar; +xlabel('x (grid index)'); ylabel('y (grid index)'); +title('Pressure p') + +%% -------- local helper: enforce velocity BCs + corner + cylinder mask ----- +function [U, V] = applyVelocityBCAndMask(U, V, Uin, i1,i2,j1,j2) + % inlet Dirichlet + U(1,:) = Uin; + V(1,:) = 0; + + % outlet Neumann (zero-gradient) + U(end,:) = U(end-1,:); + V(end,:) = V(end-1,:); + + % walls no-slip + U(2:end,1) = 0; + V(2:end,1) = 0; + U(2:end,end) = 0; + V(2:end,end) = 0; + + % inlet-wall corners should obey wall + U(1,1) = 0; U(1,end) = 0; + V(1,1) = 0; V(1,end) = 0; + + % cylinder mask + U(i1:i2, j1:j2) = 0; + V(i1:i2, j1:j2) = 0; +end