From fe041da7830c322b6b04ba5ad3ca63087a9ee406 Mon Sep 17 00:00:00 2001 From: totork Date: Tue, 27 Oct 2026 14:41:49 +0100 Subject: [PATCH 1/7] Add proper calculation of amax for div_par_mod/div_par_fvv --- include/bout/fv_ops_impl.hxx | 35 ++++++++++++++++++++++++++--------- 1 file changed, 26 insertions(+), 9 deletions(-) diff --git a/include/bout/fv_ops_impl.hxx b/include/bout/fv_ops_impl.hxx index 67b08664c8..68a2918f6d 100644 --- a/include/bout/fv_ops_impl.hxx +++ b/include/bout/fv_ops_impl.hxx @@ -651,6 +651,9 @@ Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in, const auto& v_up = v_in.yup(); const auto& v_down = v_in.ydown(); + const auto& wave_speed_up = wave_speed_in.yup(); + const auto& wave_speed_down = wave_speed_in.ydown(); + Field3D result{emptyFrom(f_in)}; flow_ylow = zeroFrom(f_in); BOUT_FOR(i, f_in.getRegion("RGN_NOBNDRY")) { @@ -664,11 +667,17 @@ Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in, / coord->cell_volume()[i]; if (dissipative) { - const BoutReal amax = BOUTMAX(fabs(wave_speed_in[i]), fabs(v_in[i]), - fabs(v_up[iyp]), fabs(v_down[iym])); - result[i] += (0.5 * amax * (f_in[i] - f_up[iyp]) * coord->cell_area_yhigh()[i] - + 0.5 * amax * (f_in[i] - f_down[iym]) * coord->cell_area_ylow()[i]) - / coord->cell_volume()[i]; + const BoutReal amax_up = BOUTMAX(fabs(wave_speed_in[i]), fabs(wave_speed_up[iyp]), + fabs(v_up[iyp]), fabs(v_in[i])); + + const BoutReal amax_down = + BOUTMAX(fabs(wave_speed_in[i]), fabs(wave_speed_down[iym]), fabs(v_down[iym]), + fabs(v_in[i])); + + result[i] += + (0.5 * amax_up * (f_in[i] - f_up[iyp]) * coord->cell_area_yhigh()[i] + + 0.5 * amax_down * (f_in[i] - f_down[iym]) * coord->cell_area_ylow()[i]) + / coord->cell_volume()[i]; } } return result; @@ -854,6 +863,9 @@ Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, const auto& v_up = v_in.yup(); const auto& v_down = v_in.ydown(); + const auto& wave_speed_up = wave_speed_in.yup(); + const auto& wave_speed_down = wave_speed_in.ydown(); + const auto& g_22 = coord->g_22(); const auto& dy = coord->dy(); @@ -863,8 +875,12 @@ Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, const auto iym = i.ym(); // Maximum local wave speed - const BoutReal amax = BOUTMAX(wave_speed_in[i], std::abs(v_in[i]), - std::abs(v_up[iyp]), std::abs(v_down[iym])); + const BoutReal amax_up = BOUTMAX(fabs(wave_speed_in[i]), fabs(wave_speed_up[iyp]), + fabs(v_up[iyp]), fabs(v_in[i])); + + const BoutReal amax_down = + BOUTMAX(fabs(wave_speed_in[i]), fabs(wave_speed_down[iym]), fabs(v_down[iym]), + fabs(v_in[i])); result[i] = B[i] @@ -872,8 +888,9 @@ Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, - (f_down[iym] * v_down[iym] * v_down[iym] / B_down[iym]) // Penalty terms. This implementation is very dissipative. // Note: This version adds a viscosity that damps gradients of velocity - + amax * (f_in[i] + f_up[iyp]) * (v_in[i] - v_up[iyp]) / (B[i] + B_up[iyp]) - + amax * (f_in[i] + f_down[iym]) * (v_in[i] - v_down[iym]) + + amax_up * (f_in[i] + f_up[iyp]) * (v_in[i] - v_up[iyp]) + / (B[i] + B_up[iyp]) + + amax_down * (f_in[i] + f_down[iym]) * (v_in[i] - v_down[iym]) / (B[i] + B_down[iym])) / (2 * dy[i] * sqrt(g_22[i])); From 8e437ecd2db436791082e81b43118327c743df42 Mon Sep 17 00:00:00 2001 From: totork Date: Tue, 27 Oct 2026 15:13:05 +0100 Subject: [PATCH 2/7] Add asserts for parallel slices of wave_speed --- include/bout/fv_ops_impl.hxx | 2 ++ 1 file changed, 2 insertions(+) diff --git a/include/bout/fv_ops_impl.hxx b/include/bout/fv_ops_impl.hxx index 68a2918f6d..77342b6772 100644 --- a/include/bout/fv_ops_impl.hxx +++ b/include/bout/fv_ops_impl.hxx @@ -644,6 +644,7 @@ Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in, ASSERT1(f_in.hasParallelSlices()); ASSERT1(v_in.hasParallelSlices()); + ASSERT1(wave_speed_in.hasParallelSlices()); const auto& f_up = f_in.yup(); const auto& f_down = f_in.ydown(); @@ -852,6 +853,7 @@ Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, // FCI version, using yup/down fields ASSERT1(f_in.hasParallelSlices()); ASSERT1(v_in.hasParallelSlices()); + ASSERT1(wave_speed_in.hasParallelSlices()); const auto& B = coord->Bxy(); const auto& B_up = coord->Bxy().yup(); From 8f724e01eeed5bbdfc9dca4b310aca4494fb8126 Mon Sep 17 00:00:00 2001 From: totork Date: Tue, 27 Oct 2026 16:00:11 +0100 Subject: [PATCH 3/7] Fix error message --- include/bout/fv_ops_impl.hxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/bout/fv_ops_impl.hxx b/include/bout/fv_ops_impl.hxx index 77342b6772..87069bba42 100644 --- a/include/bout/fv_ops_impl.hxx +++ b/include/bout/fv_ops_impl.hxx @@ -901,8 +901,8 @@ Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, throw BoutException("Non-finite value in Div_par_fvv at {}\n" "fup {} vup {} fdown {} vdown {} amax {}\n", "B {} Bup {} Bdown {} dy {} sqrt(g_22} {}", i, f_up[i], - v_up[i], f_down[i], v_down[i], amax, B[i], B_up[i], B_down[i], - dy[i], sqrt(g_22[i])); + v_up[i], f_down[i], v_down[i], amax_up, amax_down, B[i], + B_up[i], B_down[i], dy[i], sqrt(g_22[i])); } #endif } From f32febf90a4ee37efd0bb68c214238d08142f709 Mon Sep 17 00:00:00 2001 From: totork <116705180+totork@users.noreply.github.com> Date: Wed, 28 Oct 2026 17:09:40 +0100 Subject: [PATCH 4/7] Apply suggestion from @dschwoerer Co-authored-by: David Bold --- include/bout/fv_ops_impl.hxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/bout/fv_ops_impl.hxx b/include/bout/fv_ops_impl.hxx index 87069bba42..67efe6ccd9 100644 --- a/include/bout/fv_ops_impl.hxx +++ b/include/bout/fv_ops_impl.hxx @@ -676,8 +676,8 @@ Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in, fabs(v_in[i])); result[i] += - (0.5 * amax_up * (f_in[i] - f_up[iyp]) * coord->cell_area_yhigh()[i] - + 0.5 * amax_down * (f_in[i] - f_down[iym]) * coord->cell_area_ylow()[i]) + 0.5 * (amax_up * (f_in[i] - f_up[iyp]) * coord->cell_area_yhigh()[i] + + amax_down * (f_in[i] - f_down[iym]) * coord->cell_area_ylow()[i]) / coord->cell_volume()[i]; } } From 0bb1e2e0a7b85abba04e6546cbc60e135f6630f0 Mon Sep 17 00:00:00 2001 From: totork <116705180+totork@users.noreply.github.com> Date: Wed, 28 Oct 2026 17:10:37 +0100 Subject: [PATCH 5/7] Apply suggestion from @dschwoerer Co-authored-by: David Bold --- include/bout/fv_ops_impl.hxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/bout/fv_ops_impl.hxx b/include/bout/fv_ops_impl.hxx index 67efe6ccd9..9f93d01e85 100644 --- a/include/bout/fv_ops_impl.hxx +++ b/include/bout/fv_ops_impl.hxx @@ -899,7 +899,7 @@ Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, #if CHECK > 0 if (!std::isfinite(result[i])) { throw BoutException("Non-finite value in Div_par_fvv at {}\n" - "fup {} vup {} fdown {} vdown {} amax {}\n", + "fup {} vup {} fdown {} vdown {} amax_up {} amax_down {}\n", "B {} Bup {} Bdown {} dy {} sqrt(g_22} {}", i, f_up[i], v_up[i], f_down[i], v_down[i], amax_up, amax_down, B[i], B_up[i], B_down[i], dy[i], sqrt(g_22[i])); From 35c973f17dcf62f8837fca6b82b1acdd9802de42 Mon Sep 17 00:00:00 2001 From: totork Date: Wed, 28 Oct 2026 18:46:38 +0100 Subject: [PATCH 6/7] Formatting --- include/bout/fv_ops_impl.hxx | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/include/bout/fv_ops_impl.hxx b/include/bout/fv_ops_impl.hxx index 9f93d01e85..c8be8a4693 100644 --- a/include/bout/fv_ops_impl.hxx +++ b/include/bout/fv_ops_impl.hxx @@ -676,8 +676,9 @@ Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in, fabs(v_in[i])); result[i] += - 0.5 * (amax_up * (f_in[i] - f_up[iyp]) * coord->cell_area_yhigh()[i] - + amax_down * (f_in[i] - f_down[iym]) * coord->cell_area_ylow()[i]) + 0.5 + * (amax_up * (f_in[i] - f_up[iyp]) * coord->cell_area_yhigh()[i] + + amax_down * (f_in[i] - f_down[iym]) * coord->cell_area_ylow()[i]) / coord->cell_volume()[i]; } } From bb6a1df7bfd5c6fd05686110eabb912fc014460b Mon Sep 17 00:00:00 2001 From: David Bold Date: Thu, 29 Oct 2026 10:08:11 +0100 Subject: [PATCH 7/7] Fix indices in debug message --- include/bout/fv_ops_impl.hxx | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/include/bout/fv_ops_impl.hxx b/include/bout/fv_ops_impl.hxx index c8be8a4693..bdf96d34ec 100644 --- a/include/bout/fv_ops_impl.hxx +++ b/include/bout/fv_ops_impl.hxx @@ -901,9 +901,9 @@ Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, if (!std::isfinite(result[i])) { throw BoutException("Non-finite value in Div_par_fvv at {}\n" "fup {} vup {} fdown {} vdown {} amax_up {} amax_down {}\n", - "B {} Bup {} Bdown {} dy {} sqrt(g_22} {}", i, f_up[i], - v_up[i], f_down[i], v_down[i], amax_up, amax_down, B[i], - B_up[i], B_down[i], dy[i], sqrt(g_22[i])); + "B {} Bup {} Bdown {} dy {} sqrt(g_22} {}", i, f_up[iyp], + v_up[iyp], f_down[iym], v_down[iym], amax_up, amax_down, B[i], + B_up[iyp], B_down[iym], dy[i], sqrt(g_22[i])); } #endif }