[12jordans]

What type of PR is this? (check all applicable)

• Refactor
• Feature
• Bug Fix
• Optimization
• Example
• Documentation

Description

The current code include modifies 6 files:
comparison_two_way_wave_md_vs_fd_neumann.m:

• Calls mimetic_diff_two_way_wave_eq_neumann.m and finite_diff_two_way_wave_eq_neumann.m to compare walltime, flops, and error. These files model the existing tutorial code, but adds neumann boundary conditions.

mimetic_diff_two_way_wave_eq_neumann.m:

• Uses mimetic operators to solve the equation with BCs.

finite_diff_two_way_wave_eq_neumann.m:

• Uses finite differences to solve the equation with BCs.

comparison_two_way_wave_md_vs_fd_neumann.cpp:

• This file is the C++ equivalent of comparison_two_way_wave_md_vs_fd_neumann.m,

mimetic_diff_two_way_wave_eq_neumann.m, and finite_diff_two_way_wave_eq_neumann.m. Main calls two functions mimetic_diff() and finite_diff() and compares the two.

index.md:

• Added a line for Two_Way_Wave_With_Neumann_BC.md.

Two_Way_Wave_With_Neumann_BC.md:

• Documentation for the new files.

CMakeLists.txt:

Related Issues & Documents

• I have created an Issue that is paired with this PR
• Boundary Condition Example #313

QA Instructions, Screenshots, Recordings

Please replace this line with instructions on how to test your changes, a note
on the devices and browsers this has been tested on, as well as any relevant
images for UI changes.

Added/updated tests?

_We encourage you to test all code included with MOLE, including examples.

• Yes
• No, and this is why: I just looked at the results and made sure it had expected behavior.
• I need help with writing tests

Read Contributing Guide and Code of Conduct

• 📖 I have read the MOLE Contributing Guide: https://github.com/csrc-sdsu/mole/blob/main/CONTRIBUTING.md
• 📖 I have read the MOLE Code of Conduct: https://github.com/csrc-sdsu/mole/blob/main/CODE_OF_CONDUCT.md

[optional] Are there any post deployment tasks we need to perform?

[optional] What gif best describes this PR or how it makes you feel?

[12jordans]

I have finished converting my BC tutorial to C++, but am still working on the documentation. Please let me know if I should make any changes to the Matlab or C++ code.

[valeriabarra]

Hi @12jordans , thank you for your contribution. Please update this branch with the main branch (you can use git merge or git rebase <- recommended, but then it requires a force-push). Please let me know if you need any assistance with that. Thank you

[valeriabarra]

Hey @12jordans, there appears to be a git issue. Your PR should not contain 46 commits, but only the ones you authored for your intended contribution.

Please follow these steps, after verifying that you have set up your origin (your fork's repo) and upstream (the public repo) correctly. You can verify that by running git remote -v

The git pull command above has the --rebase option that might be necessary (without the parentheses)

[valeriabarra]

Hi @12jordans was this closed by accident?

[12jordans]

Yes, I am not sure how I closed it. Is there a way for me to reopen it?

[valeriabarra]

Okay, good. Something seems still strange with the commit history. Now it shows 0 commits

[12jordans]

I tried to revert my branch and must have closed this pr in the process. I'll be adding all of my changes to my branch.

[valeriabarra]

Hi @12jordans thank you. When you think this PR is ready for review, please remove the "draft" state so that it can be reviewed (the review process may take some time and multiple rounds). Thank you

[12jordans]

@jbrzensk I am struggling a bit trying to use addScalerBC1D with an explicit scheme. From the original tutorial code, the scheme is U2_md = (2.0 * U1_md) + (L * U1_md) - U0_md. I can rearrange it to the form U2_md = (2*I+L)U1_md-U0_md, but I am unsure what to do from here. From examples I’ve looked at it seems like you would call addScalerBC1D on L and U1_md then compute U2_md. My confusion lies in what to do about the -U0_md, as I have not found any examples that have an extra U0. I’ve tried applying BCs to L and U1_md then subtracting U0_md, but my code has been unsuccessful. Is this the right way to approach the problem or is there something I am missing?

[jbrzensk]

I think all you want/need is the addScalarBClhs() function. If you look inside addScalarBC1D, the operator gets some lhs and rhs features, and b gets modified. You only need the operator modifications.

Then, you can add the addScalarBClhs right after you make the L operator, then use that throughout.

In fact, I get pretty good results taking my example, and extending the time, and adding these two lines to the operator creation.

Because the wave hits the boundary and reflects back, the method you use for the analytic is no longer valid.

[12jordans]

@jbrzensk I have tried making the changes you suggested, but I still cannot seem to get it to work correctly. I have changed my analytic solution and have added a plot to show the time steps. Once my mimetic version hits the boundary, the wave disappears rather than bounce back. I have pushed the changes I have made to my Matlab code. If you have some time could you review my mimetic_diff_two_way_wave_eq_neumann.m and give me any advice on what I am missing?

[jbrzensk]

@12jordans For sanity, add these lines to the code for mimetics:

U0_md(1)   = U0_md(2);
U0_md(end) = U0_md(end-1);

Ignore the BC parts. The addScalarBC works for a static solve, but the behavior for a wave-type equation is different, probably (@mdumett is the BC expert). You can leave that whole part out if you add the two lines above. (Default BC for a plain Laplacian are Dirichlet).

For the analytic solution, use the starting condition. When the wave reflects, it will end up back where it started. Just calculate the number of time steps needed and use the initial condition as the true solution.

For 2nd order, mimetic should be way better, maybe order 3. The finite difference should be near 2. For fourth order, those lines you add up there are 2nd order, so the scheme will fall back to second order. You can make more complicated fourth order boundary updates, or maybe @mdumett has a way to use the addScalarBC operators.

[github-actions]

This pull request has had no activity for 14 days. It will be marked stale in 14 days if there is no further activity.

[valeriabarra]

It has been a while now that this PR has been open, and with no recent activity, the github-action bot will mark it as stale soon.

@mdumett , are you now able to kindly provide a review for this PR? Thank you!
cc: @jcastillo001 , @Tony-Drummond

[jbrzensk]

Reminder for @cpaolini.

[valeriabarra]

@cpaolini are you still able to review this PR? Just a friendly reminder. Thanks

[valeriabarra]

@cpaolini are you still able to review this PR? Thank you

[github-actions]

This pull request has had no activity for 14 days. It will be marked stale in 14 days if there is no further activity.

[github-actions]

Thanks for the discussion and work here. This pull request has had no activity for 28 days, so it is now marked stale. If there is no further activity, it will be closed in 42 days.

[Copilot]

🟡 Changes recommended

The new tutorial/docs have correctness issues and the C++ example currently doesn’t reliably exercise the Neumann boundary behavior or produce consistent metrics, which undermines the stated purpose.

Once you've addressed the issues Copilot identified, you can request another Copilot review.

Pull request overview

Adds a new tutorial/example variant of the 1D two-way wave equation that introduces Neumann boundary conditions and compares mimetic finite differences vs. standard finite differences across MATLAB/Octave and C++ example code, with accompanying Sphinx documentation.

Changes:

• Introduces MATLAB/Octave solvers and a comparison driver for the two-way wave equation with Neumann BCs.
• Adds a C++ comparison example and wires it into the examples CMake build.
• Adds a new Sphinx tutorial page and links it from the Tutorials index.

File summaries

File	Description
examples/matlab_octave/tutorials/two_way_wave_equation/mimetic_diff_two_way_wave_eq_neumann.m	New mimetic-difference solver with Neumann BC handling.
examples/matlab_octave/tutorials/two_way_wave_equation/finite_diff_two_way_wave_eq_neumann.m	New finite-difference solver with Neumann BC handling.
examples/matlab_octave/tutorials/two_way_wave_equation/comparison_two_way_wave_md_vs_fd_neumann.m	New MATLAB/Octave comparison script to plot error/walltime/FLOPs for MD vs FD.
examples/cpp/comparison_two_way_wave_md_vs_fd_neumann.cpp	New C++ comparison example intended to mirror the MATLAB/Octave Neumann tutorial.
examples/cpp/CMakeLists.txt	Adds the new C++ example to the examples build list.
doc/sphinx/source/examples/Tutorials/Two_Way_Wave_Equation_With_Neumann_BC.md	Adds the Sphinx tutorial describing the Neumann-BC comparison and derivations.
doc/sphinx/source/examples/Tutorials/index.md	Links the new tutorial into the Sphinx toctree.

Review details

Suppressed comments (4)

examples/matlab_octave/tutorials/two_way_wave_equation/finite_diff_two_way_wave_eq_neumann.m:71

• Remove the disp(U0_fd) debug print; it will spam the console for each grid size and makes the example harder to use in batch runs.

    U0_fd = u(grid_fd', 0);
    disp(U0_fd)

doc/sphinx/source/examples/Tutorials/Two_Way_Wave_Equation_With_Neumann_BC.md:85

• The CFL stability condition is dimensionally incorrect in the text (c \Delta t^2 / \Delta x^2). For CTCS on the wave equation it should be expressed as c\Delta t/\Delta x \le 1 (equivalently c^2\Delta t^2/\Delta x^2 \le 1).

### Notes
We are using an explicit two-step scheme (leapfrog). The CFL condition must be satisfied: $ c \Delta t^2 / \Delta x^2 \le 1 $. If you change the stepping, or increase the number of cells, be sure to change this value. 

examples/cpp/comparison_two_way_wave_md_vs_fd_neumann.cpp:216

• Like the mimetic path, the FD path uses t = ceil(1/(c*dt)), which stops before the wave reaches the domain boundaries (time-to-boundary is ~2/c from x=0 to x=±2). This means the Neumann BC logic is unlikely to affect the solution, undermining the stated goal of providing a boundary-condition example and making the MATLAB/C++ versions inconsistent.

        const Real dx = (east - west) / static_cast<Real>(m);       // spatial discretization
        const Real r2_fd = c * c * (dt * dt / (dx * dx));           // Courant-like factor
        const u32 t = static_cast<u32>(std::ceil(1.0 / (c * dt)));  // first step euler, so t is one less

examples/matlab_octave/tutorials/two_way_wave_equation/mimetic_diff_two_way_wave_eq_neumann.m:88

• The comment says the initial displacement is a cosine curve, but the initial condition is currently defined via a Gaussian (f = @(x) exp(-x.^2/0.1)) and the d’Alembert/reflection construction. This is misleading for readers following the tutorial.

    %% Initial Displacement is cos curve
    U0_md = u(grid_md', 0);

• Files reviewed: 7/7 changed files
• Comments generated: 8
• Review effort level: Lite

---

💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.

[jbrzensk]

There are ways to do this using an explicit method, but they all seem way more complicated than this:

    %% Technique number 1, direct solve of boundary
    dc = [0; 0];                % pure Neumann (no Dirichlet component)
    nc = [1; 1];                % Neumann coefficient
    v  = [0; 0];                % zero flux — reflecting walls
    
    % Use addScalarBC1D JUST for the boundaries!!!!
    [L_bc, b_bc] = addScalarBC1D(L, zeros(m+2, 1), k, m, dx, dc, nc, v);

    L = r2_md * L;
    
    U1_md = U0_md + (0.5 * L * U0_md); % First step, Euler

    % Solve boundaries implicitly
    U1_md(1)   = (b_bc(1)   - L_bc(1,   2:end)   * U1_md(2:end))   / L_bc(1,1);
    U1_md(end) = (b_bc(end) - L_bc(end, 1:end-1) * U1_md(1:end-1)) / L_bc(end,end);

For instance, you can pre-solve for those coefficients, and rebuild the Laplacian, and have a "right-hand-side" additive matrix, separate from the Laplacian, which seems way too complicated, and edits pretty some parts of the Laplacian and boundary condition.

   %% Technique #2, over-complicated reformulation
    n = m+2;    % number of centers
    ii = 2:m+1; % non boundary indices
   
    % There is probably a clever way to do this but I do not know other than step by step.
    [Abc, bbc] = addScalarBC1D(sparse(n,n), zeros(n,1), k, m, dx, dc, nc, v);

    % Solve each boundary row for its own face value. Two scalar divisions, not solves
    aL = Abc(1,1);
    wL = -Abc(1,ii)/aL;
    qL = bbc(1)/aL;
    
    aR = Abc(n,n);
    wR = -Abc(n,ii)/aR;
    qR = bbc(n)/aR;

    % Separate into interior and boundary parts
    % u_full = P * u_interior + q
    P = [ wL; speye(m); wR ];
    q = sparse([1 n], [1 1], [qL qR], n, 1);

    %% Mimetic Scheme Matrix, condensed onto the interior
    L    = lap(k, m, dx);           % mimetic Laplacian operator
    S    = L(ii,:) * P;             % m-by-m explicit operator
    cvec = L(ii,:) * q;             % boundary data contribution (0 if v = 0)

    R    = r2_md * S;               % premultiply by dt^2 * c^2
    rhsc = r2_md * cvec;

    %% Initial Displacement
    U0_md   = u(grid_md, 0);
    Ui0_md  = U0_md(ii);                            % interior unknowns only

    % First step, Euler (du/dt = 0 at t = 0)
    Ui1_md  = Ui0_md + 0.5 * (R * Ui0_md + rhsc);
    Ui2_md  = Ui1_md;                               % Increment time step

    %% Loop over time values
    for i = 1 : t
        Ui2_md = (2.0 * Ui1_md) + (R * Ui1_md + rhsc) - Ui0_md;
        % Shift everyone back for leapfrog scheme
        Ui0_md = Ui1_md;
        Ui1_md = Ui2_md;

We are going with this first example (direct solve), and technique, and adding a note to the addScalarBC1d(2d,3d) clearly stating that Neumann boundary conditions are meant for implicit time step methods. Unless anyone wants to say or show some other way of doing it.

[github-actions]

This pull request has had no activity for 14 days. It will be marked stale in 14 days if there is no further activity.

[valeriabarra]

We are going with this first example (direct solve), and technique, and adding a note to the addScalarBC1d(2d,3d) clearly stating that Neumann boundary conditions are meant for implicit time step methods. Unless anyone wants to say or show some other way of doing it.

Hi @jbrzensk , thank you so much for going over this in much details. Who can contribute the idea # 1 so that we can get this PR through the finish line? Thank you

[jbrzensk]

@12jordans I have a fix, would it be alright if I add that to the pull request and merge it?

[mdumett]

The attached code is a 1D Heat equation with Neumann boundary conditions that are enforced at each iteration using the addScalarBC1D feature. This PDE does not have the reflexion on the walls (the issue) but it demonstrates the usage of the functionality in each time iteration for an evolution equation (the first attached file is a plot of the solution; the second is the .m source code (change the extension from txt to m before running it)).

Heat1DNeumann_addScalarBCIter.txt

I updated the branch. The review is the commits.

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.

📢 Thoughts on this report? Let us know!

[jbrzensk]

Miguel, the code you showed still solves using an 'implicit' method.

    % Solve the resulting linear system to advance the field forward in time
    u = A_bc \ b_bc; 

We are concerned here with using an explicit method for the time step with Neumann BCs.

[jbrzensk]

This seems to meet all of the required attributes.

The addScalarBC has been noted to need additional effort when solving using a explicit time step.