Skip to content

Adding cpp&MATLAB example - 2D imcompressible square cylinder flow#274

Open
yiyuef wants to merge 19 commits intocsrc-sdsu:mainfrom
yiyuef:cylinder_branch
Open

Adding cpp&MATLAB example - 2D imcompressible square cylinder flow#274
yiyuef wants to merge 19 commits intocsrc-sdsu:mainfrom
yiyuef:cylinder_branch

Conversation

@yiyuef
Copy link

@yiyuef yiyuef commented Feb 7, 2026

Description

This PR adds a new C++ example cylinder_flow_2D that solves 2D incompressible channel flow past a cylinder-like obstacle (implemented as a masked no-slip region) using a projection (pressure-correction) method built on MOLE mimetic operators.

It also adds Sphinx/MyST documentation for the example and includes a representative output figure.

Type of Change

  • New core functionality
  • New example
  • Documentation update
  • Bug fix
  • Performance improvement

Mathematical Details

Not applicable (no new operators added).
The example solves the incompressible Navier–Stokes equations with a fractional-step (projection) method and uses MOLE discrete operators (divergence, gradient, Laplacian) for the Helmholtz/Poisson solves.

Testing

  • Unit tests pass
  • Examples run successfully
  • Convergence studies completed (if applicable)
  • Cross-platform compatibility verified

Documentation

  • Code is well-commented
  • API documentation updated
  • Example documentation added (if applicable)
  • Mathematical background provided

Related Issues

#211

Additional Notes

  • Output image cylinder_flow_2D_output1.png is included alongside the documentation page.

@valeriabarra
Copy link
Collaborator

Thank you, @yiyuef for this terrific contribution!

I ran the example and can see the output figure is quite small. Can you please increase its size? Also, any specific reason why the extension .ppm was chosen for non-MAC-OS users? It is a quite uncommon and unsupported file format.

Finally, I thought we also wanted to include the MATLAB/octave version to the library. I thought you had completed that implementation, first? Thank you!

yiyuef and others added 2 commits February 11, 2026 14:43
Co-authored-by: Valeria Barra <39932030+valeriabarra@users.noreply.github.com>
Signed-off-by: Yiyue Feng <79966731+yiyuef@users.noreply.github.com>
@yiyuef
Copy link
Author

yiyuef commented Feb 12, 2026

@valeriabarra Thanks for the careful review!

  1. You’re right about the inconsistent message and the .ppm choice. So I removed the platform-dependent PPM/PNG output and now the example only writes the three CSV outputs (U_final.csv, V_final.csv, p_final.csv). Visualization is handled via a separate gnuplot script, so the output messaging is now consistent. The reason I chose ppm is that I found cpp can output this format without outside packages.

  2. Regarding the figure size: I updated the gnuplot script to generate a larger PNG.

  3. markdown file also modified, to be consistant with the above change.

  4. For the MATLAB version, I thougth I should open another PR. Because source functions in MATLAB have been changed a little bit so my old code was not good. I have looked into in and now it should work. Matlab exmaple is now also part of this PR.

@yiyuef yiyuef changed the title Adding cpp example - 2D imcompressible cylinder flow Adding cpp&MATLAB example - 2D imcompressible square cylinder flow Feb 12, 2026
@valeriabarra valeriabarra requested a review from mdumett February 12, 2026 01:04
@valeriabarra valeriabarra added the Examples Addition or improvement of examples that showcase the main functionality of the library label Feb 12, 2026
@valeriabarra valeriabarra linked an issue Feb 12, 2026 that may be closed by this pull request
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 umax = maxU;
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

umax is made twice? maxU and umax are the same.

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;
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Split these onto their own separate lines. ( clang-format )

return idx;
}

static uvec bc_bottom_indices(unsigned m, unsigned n) {
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If n is void, then

c_bottom_indices(unsigned m, unsigned /*n*/) 

#include <armadillo>

// Optional Eigen (for reusable sparse factorization)
#if defined(__has_include)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

nit: # directives are not indented ( so says clang-format )

#endif

// Mimetic operator library (in ./cpp)
#include "mole.h"
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I guess this should have armadillo down here, but maybe it doesn't work with EIGEN?

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;
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

needs braces around the idx(i) = i

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;
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

braces around if

uvec rowsbc;
if (!pieces.empty()) {
rowsbc = pieces[0];
for (size_t i = 1; i < pieces.size(); ++i) rowsbc = arma::join_cols(rowsbc, pieces[i]);
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

for loop needs braces {}

@mdumett
Copy link
Collaborator

mdumett commented Mar 2, 2026

This comment refers to the MATLAB implementation.

This is a very interesting and complex example. My main suggestions are mainly with respect to its documentation.

  1. The clear command in line 5 should be removed
  2. This example is complex. To give to the user an easy experience while reading the code, keep in mind that in the first lines of the code, one should describe what PDE is being solved, domain and geometry, initial and boundary conditions, parameter values and the algorithm. In particular de projection method.
  3. Line 24 mentions ghost coordinates but in Mimetic Differences there are no ghost points, like in Finite Differences.
  4. What are CN diffusion matrices?
  5. What are cylinder mask indices?
  6. What are U* and V*?
  7. What are AB2 and AB1?

From the code, it is not clear why the need to apply boundary conditions (applyVelocityBCAndMask) in a finite difference manner. If the code is more documented the user can appreciate the decision of the developer.

This an important example from which many users can learn about how to deal with complicated geometries. Please, provide more documentation.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

Examples Addition or improvement of examples that showcase the main functionality of the library

Projects

None yet

Development

Successfully merging this pull request may close these issues.

adding C++/matlab example, channel flow / cylinder flow

4 participants