From 4100f8828b990ace616793559d2cc97ff808031b Mon Sep 17 00:00:00 2001 From: totork Date: Mon, 17 Aug 2026 09:31:03 +0200 Subject: [PATCH 1/4] Adjust parallel diffusion operator for fci --- src/mesh/difops.cxx | 21 +++++++-------------- 1 file changed, 7 insertions(+), 14 deletions(-) diff --git a/src/mesh/difops.cxx b/src/mesh/difops.cxx index e6913a67eb..c1c89628b9 100644 --- a/src/mesh/difops.cxx +++ b/src/mesh/difops.cxx @@ -385,32 +385,25 @@ Field3D Div_par_K_Grad_par_mod_impl(const Field3DParallel& Kin, Field3D result{zeroFrom(fin)}; flow_ylow = zeroFrom(fin); - + BOUT_FOR(i, result.getRegion("RGN_NOBNDRY")) { const auto iyp = i.yp(); const auto iym = i.ym(); // Upper cell edge const BoutReal c_up = 0.5 * (Kin[i] + K_up[iyp]); // K at the upper boundary - const BoutReal J_up = - 0.5 * (coord->J()[i] + coord->J().yup()[iyp]); // Jacobian at boundary - const BoutReal g_22_up = 0.5 * (coord->g_22()[i] + coord->g_22().yup()[iyp]); - const BoutReal gradient_up = - 2. * (f_up[iyp] - fin[i]) / (coord->dy()[i] + coord->dy().yup()[iyp]); + const BoutReal gradient_up = (f_up[iyp] - fin[i]) / (coord->dy()[i] * sqrt(g_22_yhigh()[i])); - const BoutReal flux_up = c_up * J_up * gradient_up / g_22_up; + const BoutReal flux_up = c_up * gradient_up * coord->cell_area_yhigh()[i]; // Lower cell edge const BoutReal c_down = 0.5 * (Kin[i] + K_down[iym]); // K at the lower boundary - const BoutReal J_down = - 0.5 * (coord->J()[i] + coord->J().ydown()[iym]); // Jacobian at boundary - const BoutReal g_22_down = 0.5 * (coord->g_22()[i] + coord->g_22().ydown()[iym]); - const BoutReal gradient_down = - 2. * (fin[i] - f_down[iym]) / (coord->dy()[i] + coord->dy().ydown()[iym]); + const BoutReal gradient_down = (fin[i] - f_down[iym]) / (coord->dy()[i] * sqrt(g_22_ylow()[i])); - const BoutReal flux_down = c_down * J_down * gradient_down / g_22_down; + const BoutReal flux_down = c_down * gradient_down * coord->cell_area_ylow()[i]; - result[i] = (flux_up - flux_down) / (coord->dy()[i] * coord->J()[i]); + // Add the fluxes + result[i] = (flux_up - flux_down) / (coord->cell_volume()[i]); } return result; From 6ba4f1d5a6816870a999e4439a9b5f80f056b66a Mon Sep 17 00:00:00 2001 From: totork Date: Mon, 17 Aug 2026 09:38:18 +0200 Subject: [PATCH 2/4] Add missing class --- src/mesh/difops.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/mesh/difops.cxx b/src/mesh/difops.cxx index c1c89628b9..66f388c805 100644 --- a/src/mesh/difops.cxx +++ b/src/mesh/difops.cxx @@ -392,13 +392,13 @@ Field3D Div_par_K_Grad_par_mod_impl(const Field3DParallel& Kin, // Upper cell edge const BoutReal c_up = 0.5 * (Kin[i] + K_up[iyp]); // K at the upper boundary - const BoutReal gradient_up = (f_up[iyp] - fin[i]) / (coord->dy()[i] * sqrt(g_22_yhigh()[i])); + const BoutReal gradient_up = (f_up[iyp] - fin[i]) / (coord->dy()[i] * sqrt(coord->g_22_yhigh()[i])); const BoutReal flux_up = c_up * gradient_up * coord->cell_area_yhigh()[i]; // Lower cell edge const BoutReal c_down = 0.5 * (Kin[i] + K_down[iym]); // K at the lower boundary - const BoutReal gradient_down = (fin[i] - f_down[iym]) / (coord->dy()[i] * sqrt(g_22_ylow()[i])); + const BoutReal gradient_down = (fin[i] - f_down[iym]) / (coord->dy()[i] * sqrt(coord->g_22_ylow()[i])); const BoutReal flux_down = c_down * gradient_down * coord->cell_area_ylow()[i]; From 3f35c51e4ff2ce79dff0280faaca96dc76ef846e Mon Sep 17 00:00:00 2001 From: totork Date: Mon, 17 Aug 2026 09:52:55 +0200 Subject: [PATCH 3/4] Delete brackets --- src/mesh/difops.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/mesh/difops.cxx b/src/mesh/difops.cxx index 66f388c805..4def0328e3 100644 --- a/src/mesh/difops.cxx +++ b/src/mesh/difops.cxx @@ -392,13 +392,13 @@ Field3D Div_par_K_Grad_par_mod_impl(const Field3DParallel& Kin, // Upper cell edge const BoutReal c_up = 0.5 * (Kin[i] + K_up[iyp]); // K at the upper boundary - const BoutReal gradient_up = (f_up[iyp] - fin[i]) / (coord->dy()[i] * sqrt(coord->g_22_yhigh()[i])); + const BoutReal gradient_up = (f_up[iyp] - fin[i]) / (coord->dy[i] * sqrt(coord->g_22_yhigh()[i])); const BoutReal flux_up = c_up * gradient_up * coord->cell_area_yhigh()[i]; // Lower cell edge const BoutReal c_down = 0.5 * (Kin[i] + K_down[iym]); // K at the lower boundary - const BoutReal gradient_down = (fin[i] - f_down[iym]) / (coord->dy()[i] * sqrt(coord->g_22_ylow()[i])); + const BoutReal gradient_down = (fin[i] - f_down[iym]) / (coord->dy[i] * sqrt(coord->g_22_ylow()[i])); const BoutReal flux_down = c_down * gradient_down * coord->cell_area_ylow()[i]; From 2111e68726f6dbe3c937aaf5a71a006d75c0b608 Mon Sep 17 00:00:00 2001 From: totork Date: Mon, 17 Aug 2026 10:02:20 +0200 Subject: [PATCH 4/4] Add missing () for dy --- src/mesh/difops.cxx | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/mesh/difops.cxx b/src/mesh/difops.cxx index 4def0328e3..bcd74a622b 100644 --- a/src/mesh/difops.cxx +++ b/src/mesh/difops.cxx @@ -385,20 +385,22 @@ Field3D Div_par_K_Grad_par_mod_impl(const Field3DParallel& Kin, Field3D result{zeroFrom(fin)}; flow_ylow = zeroFrom(fin); - + BOUT_FOR(i, result.getRegion("RGN_NOBNDRY")) { const auto iyp = i.yp(); const auto iym = i.ym(); // Upper cell edge const BoutReal c_up = 0.5 * (Kin[i] + K_up[iyp]); // K at the upper boundary - const BoutReal gradient_up = (f_up[iyp] - fin[i]) / (coord->dy[i] * sqrt(coord->g_22_yhigh()[i])); + const BoutReal gradient_up = + (f_up[iyp] - fin[i]) / (coord->dy()[i] * sqrt(coord->g_22_yhigh()[i])); const BoutReal flux_up = c_up * gradient_up * coord->cell_area_yhigh()[i]; // Lower cell edge const BoutReal c_down = 0.5 * (Kin[i] + K_down[iym]); // K at the lower boundary - const BoutReal gradient_down = (fin[i] - f_down[iym]) / (coord->dy[i] * sqrt(coord->g_22_ylow()[i])); + const BoutReal gradient_down = + (fin[i] - f_down[iym]) / (coord->dy()[i] * sqrt(coord->g_22_ylow()[i])); const BoutReal flux_down = c_down * gradient_down * coord->cell_area_ylow()[i];