Skip to content

Preview: Use cell area and volume for finite volume code #3392

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

Merged
dschwoerer merged 13 commits into from
Jun 26, 2026
Merged

Preview: Use cell area and volume for finite volume code #3392

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
b3ccacb
Merge branch 'cell_area_volume' into move-ops-area
dschwoerer Jun 16, 2026
639808c
Use cell area and volumes for div_par_mod
dschwoerer Jun 16, 2026
5a2a3dc
Use cell area and volumes for Div_par_fvv
dschwoerer Jun 16, 2026
820e221
Use cell area and volumes for Div_par
dschwoerer Jun 16, 2026
3b0117e
Also fill the boundaries to the extend possible
dschwoerer Jun 17, 2026
cee53db
Merge remote-tracking branch 'origin/next' into move-ops-area
dschwoerer Jun 24, 2026
381af95
Merge branch 'move-ops' into move-ops-area
dschwoerer Jun 24, 2026
bc980b1
Ensure calls to _compute_... are guarded
dschwoerer Jun 24, 2026
e68d165
Do not return from critical areas
dschwoerer Jun 24, 2026
1672208
Move critical section to compute functions
dschwoerer Jun 24, 2026
b89867d
Add DDY template for expression
dschwoerer Jun 24, 2026
586c06b
Merge branch 'move-ops' into move-ops-area
dschwoerer Jun 26, 2026
26876b7
Add missing used header
dschwoerer Jun 26, 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
70 changes: 42 additions & 28 deletions include/bout/coordinates.hxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -116,101 +116,115 @@ public:
FieldMetric& g_22_yhigh();
// Cell Areas
const FieldMetric& cell_area_xlow() const {
if (!_cell_area_xlow.has_value()) {
_compute_cell_area_x();
if (_cell_area_xlow.has_value()) {
return *_cell_area_xlow;
}
_compute_cell_area_x();
ASSERT2(_cell_area_xlow.has_value());
return *_cell_area_xlow;
}
const FieldMetric& cell_area_xhigh() const {
if (!_cell_area_xhigh.has_value()) {
_compute_cell_area_x();
if (_cell_area_xhigh.has_value()) {
return *_cell_area_xhigh;
}
_compute_cell_area_x();
ASSERT2(_cell_area_xhigh.has_value());
return *_cell_area_xhigh;
}
const FieldMetric& cell_area_ylow() const {
if (!_cell_area_ylow.has_value()) {
_compute_cell_area_y();
if (_cell_area_ylow.has_value()) {
return *_cell_area_ylow;
}
_compute_cell_area_y();
ASSERT2(_cell_area_ylow.has_value());
return *_cell_area_ylow;
}
const FieldMetric& cell_area_yhigh() const {
if (!_cell_area_yhigh.has_value()) {
_compute_cell_area_y();
if (_cell_area_yhigh.has_value()) {
return *_cell_area_yhigh;
}
_compute_cell_area_y();
ASSERT2(_cell_area_yhigh.has_value());
return *_cell_area_yhigh;
}
const FieldMetric& cell_area_zlow() const {
if (!_cell_area_zlow.has_value()) {
_compute_cell_area_z();
if (_cell_area_zlow.has_value()) {
return *_cell_area_zlow;
}
_compute_cell_area_z();
ASSERT2(_cell_area_zlow.has_value());
return *_cell_area_zlow;
}
const FieldMetric& cell_area_zhigh() const {
if (!_cell_area_zhigh.has_value()) {
_compute_cell_area_z();
if (_cell_area_zhigh.has_value()) {
return *_cell_area_zhigh;
}
_compute_cell_area_z();
ASSERT2(_cell_area_zhigh.has_value());
return *_cell_area_zhigh;
}
FieldMetric& cell_area_xlow() {
if (!_cell_area_xlow.has_value()) {
_compute_cell_area_x();
if (_cell_area_xlow.has_value()) {
return *_cell_area_xlow;
}
_compute_cell_area_x();
ASSERT2(_cell_area_xlow.has_value());
return *_cell_area_xlow;
}
FieldMetric& cell_area_xhigh() {
if (!_cell_area_xhigh.has_value()) {
_compute_cell_area_x();
if (_cell_area_xhigh.has_value()) {
return *_cell_area_xhigh;
}
_compute_cell_area_x();
ASSERT2(_cell_area_xhigh.has_value());
return *_cell_area_xhigh;
}
FieldMetric& cell_area_ylow() {
if (!_cell_area_ylow.has_value()) {
_compute_cell_area_y();
if (_cell_area_ylow.has_value()) {
return *_cell_area_ylow;
}
_compute_cell_area_y();
ASSERT2(_cell_area_ylow.has_value());
return *_cell_area_ylow;
}
FieldMetric& cell_area_yhigh() {
if (!_cell_area_yhigh.has_value()) {
_compute_cell_area_y();
if (_cell_area_yhigh.has_value()) {
return *_cell_area_yhigh;
}
_compute_cell_area_y();
ASSERT2(_cell_area_yhigh.has_value());
return *_cell_area_yhigh;
}
FieldMetric& cell_area_zlow() {
if (!_cell_area_zlow.has_value()) {
_compute_cell_area_z();
if (_cell_area_zlow.has_value()) {
return *_cell_area_zlow;
}
_compute_cell_area_z();
ASSERT2(_cell_area_zlow.has_value());
return *_cell_area_zlow;
}
FieldMetric& cell_area_zhigh() {
if (!_cell_area_zhigh.has_value()) {
_compute_cell_area_z();
if (_cell_area_zhigh.has_value()) {
return *_cell_area_zhigh;
}
_compute_cell_area_z();
ASSERT2(_cell_area_zhigh.has_value());
return *_cell_area_zhigh;
}
// Cell Volume
const FieldMetric& cell_volume() const {
if (!_cell_volume.has_value()) {
_compute_cell_volume();
if (_cell_volume.has_value()) {
return *_cell_volume;
}
_compute_cell_volume();
ASSERT2(_cell_volume.has_value());
return *_cell_volume;
}
FieldMetric& cell_volume() {
if (!_cell_volume.has_value()) {
_compute_cell_volume();
if (_cell_volume.has_value()) {
return *_cell_volume;
}
_compute_cell_volume();
ASSERT2(_cell_volume.has_value());
return *_cell_volume;
}
Expand Down
8 changes: 8 additions & 0 deletions include/bout/derivs.hxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,7 @@
#include "bout/vector3d.hxx"

#include "bout/bout_types.hxx"
#include <type_traits>

////////// FIRST DERIVATIVES //////////

Expand Down Expand Up @@ -86,6 +87,13 @@ Field3D DDY(const Field3DParallel& f, CELL_LOC outloc = CELL_DEFAULT,
const std::string& method = "DEFAULT",
const std::string& region = "RGN_NOBNDRY");

template <typename E>
std::enable_if_t<is_expr_field3d_v<E> && !bout::utils::is_Field3D_v<E>, const Field3D>
Comment thread
dschwoerer marked this conversation as resolved.
DDY(const E& expr, CELL_LOC outloc = CELL_DEFAULT, const std::string& method = "DEFAULT",
const std::string& region = "RGN_NOBNDRY") {
return DDY(Field3D{expr}, outloc, method, region);
}

/// Calculate first partial derivative in Y
///
/// \f$\partial / \partial y\f$
Expand Down
122 changes: 34 additions & 88 deletions include/bout/fv_ops_impl.hxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -334,41 +334,24 @@ Field3D Div_par(const Field3D& f_in, const Field3D& v_in, const Field3D& wave_sp
// Pre-calculate factors which multiply fluxes
#if not(BOUT_USE_METRIC_3D)
// For right cell boundaries
BoutReal common_factor = (coord->J(i, j) + coord->J(i, j + 1))
/ (sqrt(coord->g_22(i, j)) + sqrt(coord->g_22(i, j + 1)));

const BoutReal flux_factor_rc = common_factor / (coord->dy(i, j) * coord->J(i, j));
const BoutReal flux_factor_rp =
common_factor / (coord->dy(i, j + 1) * coord->J(i, j + 1));

const BoutReal area_r = coord->cell_area_yhigh()(i, j);
const BoutReal flux_factor_rc = area_r / coord->cell_volume()(i, j);
const BoutReal flux_factor_rp = area_r / coord->cell_volume()(i, j + 1);
// For left cell boundaries
common_factor = (coord->J(i, j) + coord->J(i, j - 1))
/ (sqrt(coord->g_22(i, j)) + sqrt(coord->g_22(i, j - 1)));

const BoutReal flux_factor_lc = common_factor / (coord->dy(i, j) * coord->J(i, j));
const BoutReal flux_factor_lm =
common_factor / (coord->dy(i, j - 1) * coord->J(i, j - 1));
const BoutReal area_l = coord->cell_area_ylow()(i, j);
const BoutReal flux_factor_lc = area_l / coord->cell_volume()(i, j);
const BoutReal flux_factor_lm = area_l / coord->cell_volume()(i, j - 1);
#endif
for (int k = mesh->zstart; k <= mesh->zend; k++) {
#if BOUT_USE_METRIC_3D
// For right cell boundaries
BoutReal common_factor =
(coord->J(i, j, k) + coord->J(i, j + 1, k))
/ (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j + 1, k)));

BoutReal flux_factor_rc =
common_factor / (coord->dy(i, j, k) * coord->J(i, j, k));
BoutReal flux_factor_rp =
common_factor / (coord->dy(i, j + 1, k) * coord->J(i, j + 1, k));

const BoutReal area_r = coord->cell_area_yhigh()(i, j, k);
const BoutReal flux_factor_rc = area_r / coord->cell_volume()(i, j, k);
const BoutReal flux_factor_rp = area_r / coord->cell_volume()(i, j + 1, k);
// For left cell boundaries
common_factor = (coord->J(i, j, k) + coord->J(i, j - 1, k))
/ (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j - 1, k)));

BoutReal flux_factor_lc =
common_factor / (coord->dy(i, j, k) * coord->J(i, j, k));
BoutReal flux_factor_lm =
common_factor / (coord->dy(i, j - 1, k) * coord->J(i, j - 1, k));
const BoutReal area_l = coord->cell_area_ylow()(i, j, k);
const BoutReal flux_factor_lc = area_l / coord->cell_volume()(i, j, k);
const BoutReal flux_factor_lm = area_l / coord->cell_volume()(i, j - 1, k);
#endif

////////////////////////////////////////////
Expand Down Expand Up @@ -656,12 +639,10 @@ Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in,
const auto iym = i.ym();

result[i] = (0.25 * (f_in[i] + f_up[iyp]) * (v_in[i] + v_up[iyp])
* (coord->J[i] + coord->J.yup()[iyp])
/ (sqrt(coord->g_22[i]) + sqrt(coord->g_22.yup()[iyp]))
* coord->cell_area_yhigh()[i]
- 0.25 * (f_in[i] + f_down[iym]) * (v_in[i] + v_down[iym])
* (coord->J[i] + coord->J.ydown()[iym])
/ (sqrt(coord->g_22[i]) + sqrt(coord->g_22.ydown()[iym])))
/ (coord->dy[i] * coord->J[i]);
* coord->cell_area_ylow()[i])
/ coord->cell_volume()[i];
}
return result;
}
Expand Down Expand Up @@ -702,57 +683,30 @@ Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in,
// Pre-calculate factors which multiply fluxes
#if not(BOUT_USE_METRIC_3D)
// For right cell boundaries
const BoutReal common_factor_r =
(coord->J(i, j) + coord->J(i, j + 1))
/ (sqrt(coord->g_22(i, j)) + sqrt(coord->g_22(i, j + 1)));
const BoutReal area_rp = coord->cell_area_yhigh()(i, j);

const BoutReal flux_factor_rc =
common_factor_r / (coord->dy(i, j) * coord->J(i, j));
const BoutReal flux_factor_rp =
common_factor_r / (coord->dy(i, j + 1) * coord->J(i, j + 1));

const BoutReal area_rp =
common_factor_r * coord->dx(i, j + 1) * coord->dz(i, j + 1);
const BoutReal flux_factor_rc = area_rp / coord->cell_volume()(i, j);
const BoutReal flux_factor_rp = area_rp / coord->cell_volume()(i, j + 1);

// For left cell boundaries
const BoutReal common_factor_l =
(coord->J(i, j) + coord->J(i, j - 1))
/ (sqrt(coord->g_22(i, j)) + sqrt(coord->g_22(i, j - 1)));

const BoutReal flux_factor_lc =
common_factor_l / (coord->dy(i, j) * coord->J(i, j));
const BoutReal flux_factor_lm =
common_factor_l / (coord->dy(i, j - 1) * coord->J(i, j - 1));
const BoutReal area_lc = coord->cell_area_ylow()(i, j);

const BoutReal area_lc = common_factor_l * coord->dx(i, j) * coord->dz(i, j);
const BoutReal flux_factor_lc = area_lc / coord->cell_volume()(i, j);
const BoutReal flux_factor_lm = area_lc / coord->cell_volume()(i, j - 1);
#endif
for (int k = 0; k < mesh->LocalNz; k++) {
#if BOUT_USE_METRIC_3D
// For right cell boundaries
const BoutReal common_factor_r =
(coord->J(i, j, k) + coord->J(i, j + 1, k))
/ (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j + 1, k)));

const BoutReal flux_factor_rc =
common_factor_r / (coord->dy(i, j, k) * coord->J(i, j, k));
const BoutReal flux_factor_rp =
common_factor_r / (coord->dy(i, j + 1, k) * coord->J(i, j + 1, k));
const BoutReal area_rp = coord->cell_area_yhigh()(i, j, k);

const BoutReal area_rp =
common_factor_r * coord->dx(i, j + 1, k) * coord->dz(i, j + 1, k);
const BoutReal flux_factor_rc = area_rp / coord->cell_volume()(i, j, k);
const BoutReal flux_factor_rp = area_rp / coord->cell_volume()(i, j + 1, k);

// For left cell boundaries
const BoutReal common_factor_l =
(coord->J(i, j, k) + coord->J(i, j - 1, k))
/ (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j - 1, k)));
const BoutReal area_lc = coord->cell_area_ylow()(i, j, k);

const BoutReal flux_factor_lc =
common_factor_l / (coord->dy(i, j, k) * coord->J(i, j, k));
const BoutReal flux_factor_lm =
common_factor_l / (coord->dy(i, j - 1, k) * coord->J(i, j - 1, k));

const BoutReal area_lc =
common_factor_l * coord->dx(i, j, k) * coord->dz(i, j, k);
const BoutReal flux_factor_lc = area_lc / coord->cell_volume()(i, j, k);
const BoutReal flux_factor_lm = area_lc / coord->cell_volume()(i, j - 1, k);
#endif

////////////////////////////////////////////
Expand Down Expand Up @@ -937,24 +891,16 @@ Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in,

for (int k = 0; k < mesh->LocalNz; k++) {
// For right cell boundaries
const BoutReal common_factor_r =
(coord->J(i, j, k) + coord->J(i, j + 1, k))
/ (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j + 1, k)));
const BoutReal area_r = coord->cell_area_yhigh()(i, j, k);

const BoutReal flux_factor_rc =
common_factor_r / (coord->dy(i, j, k) * coord->J(i, j, k));
const BoutReal flux_factor_rp =
common_factor_r / (coord->dy(i, j + 1, k) * coord->J(i, j + 1, k));
const BoutReal flux_factor_rc = area_r / coord->cell_volume()(i, j, k);
const BoutReal flux_factor_rp = area_r / coord->cell_volume()(i, j + 1, k);

// For left cell boundaries
const BoutReal common_factor_l =
(coord->J(i, j, k) + coord->J(i, j - 1, k))
/ (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j - 1, k)));

const BoutReal flux_factor_lc =
common_factor_l / (coord->dy(i, j, k) * coord->J(i, j, k));
const BoutReal flux_factor_lm =
common_factor_l / (coord->dy(i, j - 1, k) * coord->J(i, j - 1, k));
const BoutReal area_l = coord->cell_area_ylow()(i, j, k);

const BoutReal flux_factor_lc = area_l / coord->cell_volume()(i, j, k);
const BoutReal flux_factor_lm = area_l / coord->cell_volume()(i, j - 1, k);

////////////////////////////////////////////
// Reconstruct f at the cell faces
Expand Down
Loading
Loading