Implement the third set of NRBC for 1D Euler solver

For the 1D Euler solver, implement the third set of NRBC documented in [NASA/TM-2003-212495-REV1](http://gltrs.grc.nasa.gov/Citations.aspx?id=498), i.e., NRBC type 3. $\lambda$ is chosen to be 0.
solvcon · Apr 30, 2024 · 92316dd · 92316dd
1 parent bc9e26f
commit 92316dd
Show file tree

Hide file tree

Showing 2 changed files with 21 additions and 17 deletions.
diff --git a/cpp/modmesh/onedim/Euler1DCore.cpp b/cpp/modmesh/onedim/Euler1DCore.cpp
@@ -153,39 +153,43 @@ void Euler1DCore::march_half_so0(bool odd_plane)
 
 void Euler1DCore::treat_boundary_so0()
 {
+    /* Non-reflecting boundary condition (NRBC) type 3 with $\lambda=0$
+     * (the third set in Chang 05) */
     // Set outside value from inside value.
     {
         // Left boundary.
-        size_t const ic = 0;
-        m_so0(ic, 0) = m_so0(ic + 2, 0);
-        m_so0(ic, 1) = m_so0(ic + 2, 1);
-        m_so0(ic, 2) = m_so0(ic + 2, 2);
+        size_t const ic = 1;
+        m_so0(ic, 0) = m_so0(ic + 1, 0);
+        m_so0(ic, 1) = m_so0(ic + 1, 1);
+        m_so0(ic, 2) = m_so0(ic + 1, 2);
     }
     {
         // Right boundary.
-        size_t const ic = ncoord() - 1;
-        m_so0(ic, 0) = m_so0(ic - 2, 0);
-        m_so0(ic, 1) = m_so0(ic - 2, 1);
-        m_so0(ic, 2) = m_so0(ic - 2, 2);
+        size_t const ic = ncoord() - 2;
+        m_so0(ic, 0) = m_so0(ic - 1, 0);
+        m_so0(ic, 1) = m_so0(ic - 1, 1);
+        m_so0(ic, 2) = m_so0(ic - 1, 2);
     }
 }
 
 void Euler1DCore::treat_boundary_so1()
 {
+    /* Non-reflecting boundary condition (NRBC) type 3 with $\lambda=0$
+     * (the third set in Chang 05) */
     // Set outside value from inside value.
     {
         // Left boundary.
-        size_t const ic = 0;
-        m_so1(ic, 0) = m_so1(ic + 2, 0);
-        m_so1(ic, 1) = m_so1(ic + 2, 1);
-        m_so1(ic, 2) = m_so1(ic + 2, 2);
+        size_t const ic = 1;
+        m_so1(ic, 0) = m_so1(ic + 1, 0);
+        m_so1(ic, 1) = m_so1(ic + 1, 1);
+        m_so1(ic, 2) = m_so1(ic + 1, 2);
     }
     {
         // Right boundary.
-        size_t const ic = ncoord() - 1;
-        m_so1(ic, 0) = m_so1(ic - 2, 0);
-        m_so1(ic, 1) = m_so1(ic - 2, 1);
-        m_so1(ic, 2) = m_so1(ic - 2, 2);
+        size_t const ic = ncoord() - 2;
+        m_so1(ic, 0) = m_so1(ic - 1, 0);
+        m_so1(ic, 1) = m_so1(ic - 1, 1);
+        m_so1(ic, 2) = m_so1(ic - 1, 2);
     }
 }
 

diff --git a/cpp/modmesh/onedim/Euler1DCore.hpp b/cpp/modmesh/onedim/Euler1DCore.hpp
@@ -338,8 +338,8 @@ inline void Euler1DCore::march_alpha(size_t steps)
     {
         march_half1_alpha<ALPHA>();
         treat_boundary_so0();
-        march_half2_alpha<ALPHA>();
         treat_boundary_so1();
+        march_half2_alpha<ALPHA>();
     }
 }