Skip to content

Add more accurate calculation of amax for div_par_mod/div_par_fvv #3528

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
totork wants to merge 7 commits into
base: next
Choose a base branch
from
Open

Add more accurate calculation of amax for div_par_mod/div_par_fvv #3528

Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
fe041da
Add proper calculation of amax for div_par_mod/div_par_fvv
totork Oct 27, 2026
8e437ec
Add asserts for parallel slices of wave_speed
totork Oct 27, 2026
8f724e0
Fix error message
totork Oct 27, 2026
f32febf
Apply suggestion from @dschwoerer
totork Oct 28, 2026
0bb1e2e
Apply suggestion from @dschwoerer
totork Oct 28, 2026
35c973f
Formatting
totork Oct 28, 2026
bb6a1df
Fix indices in debug message
dschwoerer Oct 29, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
46 changes: 33 additions & 13 deletions include/bout/fv_ops_impl.hxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -644,13 +644,17 @@ 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();

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")) {
Expand All @@ -664,11 +668,18 @@ 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]
+ amax_down * (f_in[i] - f_down[iym]) * coord->cell_area_ylow()[i])
/ coord->cell_volume()[i];
}
}
return result;
Expand Down Expand Up @@ -843,6 +854,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();
Expand All @@ -854,6 +866,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();

Expand All @@ -863,27 +878,32 @@ 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]
* ((f_up[iyp] * v_up[iyp] * v_up[iyp] / B_up[iyp])
- (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]));

#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",
"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]));
"fup {} vup {} fdown {} vdown {} amax_up {} amax_down {}\n",
"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
}
Expand Down
Loading