Refactor Coordinates with more encapsulation#3425
Conversation
![github-actions[bot]](https://avatars.githubusercontent.com/in/15368?s=60&v=4){kind=small}
| } else { | ||
| // Use full expression with all terms | ||
| result = b0xGrad_dot_Grad(p, f) / coord->Bxy; | ||
| result = b0xGrad_dot_Grad(p, f) / coord->Bxy(); |
There was a problem hiding this comment.
warning: no header providing "b0xGrad_dot_Grad" is directly included [misc-include-cleaner]
result = b0xGrad_dot_Grad(p, f) / coord->Bxy();
^| @@ -331,7 +331,7 @@ class CWM : public PhysicsModel { | |||
| result = VDDZ(-DDX(p), f); | |||
There was a problem hiding this comment.
warning: no matching function for call to 'VDDZ' [clang-diagnostic-error]
result = VDDZ(-DDX(p), f);
^Additional context
include/bout/derivs.hxx:472: candidate function not viable: no known conversion from 'const Field3D' to 'const Field2D' for 2nd argument
^include/bout/derivs.hxx:455: candidate function not viable: no known conversion from 'BinaryExpr<Field2D, Constant, Field2D, bout::op::Mul>' (aka 'BinaryExpr<Field2D, Constant, Field2D, bout::op::Mul>') to 'const Field3D' for 1st argument
^include/bout/derivs.hxx:489: candidate function not viable: no known conversion from 'BinaryExpr<Field2D, Constant, Field2D, bout::op::Mul>' (aka 'BinaryExpr<Field2D, Constant, Field2D, bout::op::Mul>') to 'const Field3D' for 1st argument
^| } else { | ||
| // Use full expression with all terms | ||
| result = b0xGrad_dot_Grad(p, f) / coord->Bxy; | ||
| result = b0xGrad_dot_Grad(p, f) / coord->Bxy(); |
There was a problem hiding this comment.
warning: no header providing "b0xGrad_dot_Grad" is directly included [misc-include-cleaner]
result = b0xGrad_dot_Grad(p, f) / coord->Bxy();
^| } else { | ||
| // Use full expression with all terms | ||
| result = b0xGrad_dot_Grad(p, f) / coord->Bxy; | ||
| result = b0xGrad_dot_Grad(p, f) / coord->Bxy(); |
There was a problem hiding this comment.
warning: no header providing "b0xGrad_dot_Grad" is directly included [misc-include-cleaner]
result = b0xGrad_dot_Grad(p, f) / coord->Bxy();
^|
|
||
| hyper_mu_x = hyperviscos * metric->g_11 * SQ(metric->dx) | ||
| * abs(metric->g11 * D2DX2(U)) / (abs(U) + 1e-3); | ||
| hyper_mu_x = hyperviscos * metric->g_11() * SQ(metric->dx()) |
There was a problem hiding this comment.
warning: no header providing "SQ" is directly included [misc-include-cleaner]
hyper_mu_x = hyperviscos * metric->g_11() * SQ(metric->dx())
^| + metric_x->g_13()[i] * z_at_x[i]; | ||
| y[i] = metric_y->g_22()[i] * y[i] + metric_y->g_12()[i] * x_at_y[i] | ||
| + metric_y->g_23()[i] * z_at_y[i]; | ||
| z[i] = metric_z->g_33()[i] * z[i] + metric_z->g_13()[i] * x_at_z[i] |
There was a problem hiding this comment.
warning: '*' has higher precedence than '+'; add parentheses to explicitly specify the order of operations [readability-math-missing-parentheses]
| z[i] = metric_z->g_33()[i] * z[i] + metric_z->g_13()[i] * x_at_z[i] | |
| z[i] = metric_z->g_33()[i] * z[i] + (metric_z->g_13()[i] * x_at_z[i]) |
| y[i] = metric_y->g_22()[i] * y[i] + metric_y->g_12()[i] * x_at_y[i] | ||
| + metric_y->g_23()[i] * z_at_y[i]; | ||
| z[i] = metric_z->g_33()[i] * z[i] + metric_z->g_13()[i] * x_at_z[i] | ||
| + metric_z->g_23()[i] * y_at_z[i]; |
There was a problem hiding this comment.
warning: '*' has higher precedence than '+'; add parentheses to explicitly specify the order of operations [readability-math-missing-parentheses]
| + metric_z->g_23()[i] * y_at_z[i]; | |
| + (metric_z->g_23()[i] * y_at_z[i]); |
| gx[i] = metric->g_11[i] * x[i] + metric->g_12[i] * y[i] + metric->g_13[i] * z[i]; | ||
| gy[i] = metric->g_22[i] * y[i] + metric->g_12[i] * x[i] + metric->g_23[i] * z[i]; | ||
| gz[i] = metric->g_33[i] * z[i] + metric->g_13[i] * x[i] + metric->g_23[i] * y[i]; | ||
| gx[i] = metric->g_11()[i] * x[i] + metric->g_12()[i] * y[i] |
There was a problem hiding this comment.
warning: '*' has higher precedence than '+'; add parentheses to explicitly specify the order of operations [readability-math-missing-parentheses]
| gx[i] = metric->g_11()[i] * x[i] + metric->g_12()[i] * y[i] | |
| gx[i] = (metric->g_11()[i] * x[i]) + metric->g_12()[i] * y[i] |
| gx[i] = metric->g_11[i] * x[i] + metric->g_12[i] * y[i] + metric->g_13[i] * z[i]; | ||
| gy[i] = metric->g_22[i] * y[i] + metric->g_12[i] * x[i] + metric->g_23[i] * z[i]; | ||
| gz[i] = metric->g_33[i] * z[i] + metric->g_13[i] * x[i] + metric->g_23[i] * y[i]; | ||
| gx[i] = metric->g_11()[i] * x[i] + metric->g_12()[i] * y[i] |
There was a problem hiding this comment.
warning: '*' has higher precedence than '+'; add parentheses to explicitly specify the order of operations [readability-math-missing-parentheses]
| gx[i] = metric->g_11()[i] * x[i] + metric->g_12()[i] * y[i] | |
| gx[i] = metric->g_11()[i] * x[i] + (metric->g_12()[i] * y[i]) |
| gy[i] = metric->g_22[i] * y[i] + metric->g_12[i] * x[i] + metric->g_23[i] * z[i]; | ||
| gz[i] = metric->g_33[i] * z[i] + metric->g_13[i] * x[i] + metric->g_23[i] * y[i]; | ||
| gx[i] = metric->g_11()[i] * x[i] + metric->g_12()[i] * y[i] | ||
| + metric->g_13()[i] * z[i]; |
There was a problem hiding this comment.
warning: '*' has higher precedence than '+'; add parentheses to explicitly specify the order of operations [readability-math-missing-parentheses]
| + metric->g_13()[i] * z[i]; | |
| + (metric->g_13()[i] * z[i]); |
| #define COPY_STRIPE1(symbol) \ | ||
| if (coords->symbol().isAllocated()) { \ | ||
| data[stripe_size * ind.ind + static_cast<int>(Offset::symbol)] = \ | ||
| coords->symbol()[ind]; \ |
There was a problem hiding this comment.
Will the compiler be smart enough to hoist the function call out of the loop? It may not be critical as the coordinates accessor will be cached.
There was a problem hiding this comment.
Yep, with any level of optimisation, all of the new Coordinates getters are completely inlined (godbolt demo, operator[] should also get inlined, I just made it opaque for this demo)
ZedThree
force-pushed
the
refactor-coords-rebase-rebase
branch
from
603ae3a to
b966661
Compare
Refactor `Coordinates` with a more encapsulated interface. Derived quantities are now automatically updated when the metric components are changed. The most impactful change is the conversion of the metric components and grid sizes to `const&` methods: ```diff - coords.dx + coords.dx() ``` Among other things, this change allows us to defer computation of many of the derived quantities (such as Christoffel symbols) and remove the `Coordinates::geometry` method.
Covariant metrics have `_`. The new code is consistent with e.g. metric reads around line 450.
`new_coordinates` should be moved into `coords_map` after it has been recalculated, not before.
`CoordinatesAccessor` and non-uniform geometry quantities depend on cell spacing. Reset cache so they are updated when needed. Bxy is not calculated on demand, so is calculated in `set*MetricTensor` methods.
z index should be in range. Indexed accessors return BoutReal, not reference to field.
ZedThree
force-pushed
the
refactor-coords-rebase-rebase
branch
from
b966661 to
7da06af
Compare
|
Thanks for the fixes @bendudson! |
| #include <bout/fft.hxx> | ||
| #include <bout/field3d.hxx> | ||
| #include <bout/output.hxx> | ||
| #include <bout/sys/timer.hxx> |
There was a problem hiding this comment.
warning: included header output.hxx is not used directly [misc-include-cleaner]
| #include <bout/sys/timer.hxx> | |
| #include <bout/sys/timer.hxx> |
|
|
||
| // NOTE: For now the X-Z terms are omitted, so check that they are small | ||
| ASSERT2(max(abs(coord->g13)) < 1e-5); | ||
| ASSERT2(max(abs(coord->g13())) < 1e-5); |
There was a problem hiding this comment.
warning: no header providing "ASSERT2" is directly included [misc-include-cleaner]
ASSERT2(max(abs(coord->g13())) < 1e-5);
^|
|
||
| // NOTE: For now the X-Z terms are omitted, so check that they are small | ||
| ASSERT2(max(abs(coord->g13)) < 1e-5); | ||
| ASSERT2(max(abs(coord->g13())) < 1e-5); |
There was a problem hiding this comment.
warning: no header providing "abs" is directly included [misc-include-cleaner]
src/invert/laplacexz/impls/cyclic/laplacexz-cyclic.cxx:2:
+ #include <cstdlib>| const Field2D dy = coord->dy; | ||
| const Field2D J = coord->J; | ||
| const Field2D g_22 = coord->g_22; | ||
| const Field2D dy = coord->dy(); |
There was a problem hiding this comment.
warning: no header providing "Field2D" is directly included [misc-include-cleaner]
const Field2D dy = coord->dy();
^| #define COPY_STRIPE1(symbol) \ | ||
| if (coords->symbol.isAllocated()) \ | ||
| data[stripe_size * ind.ind + static_cast<int>(Offset::symbol)] = coords->symbol[ind]; | ||
| #define COPY_STRIPE1(symbol) \ |
There was a problem hiding this comment.
warning: function-like macro 'COPY_STRIPE1' used; consider a 'constexpr' template function [cppcoreguidelines-macro-usage]
#define COPY_STRIPE1(symbol) \
^| dx4 = SQ(SQ(coords->dx)); | ||
| dy4 = SQ(SQ(coords->dy)); | ||
| dz4 = SQ(SQ(coords->dz)); | ||
| dx4 = SQ(SQ(coords->dx())); |
There was a problem hiding this comment.
warning: no header providing "SQ" is directly included [misc-include-cleaner]
dx4 = SQ(SQ(coords->dx()));
^| dz4 = SQ(SQ(coords->dz)); | ||
| dx4 = SQ(SQ(coords->dx())); | ||
| dy4 = SQ(SQ(coords->dy())); | ||
| dz4 = SQ(SQ(coords->dz())); |
There was a problem hiding this comment.
warning: assigning to 'BoutReal' (aka 'double') from incompatible type 'BinaryExpr<Field2D, BinaryExpr<Field2D, Field2D, Field2D, bout::op::Square>, BinaryExpr<Field2D, Field2D, Field2D, bout::op::Square>, bout::op::Square>' [clang-diagnostic-error]
dz4 = SQ(SQ(coords->dz()));
^|
|
||
| SAVE_REPEAT(Ve); | ||
|
|
||
| output.write("dx = {:e}, dy = {:e}, dz = {:e}\n", coords->dx(2, 2), coords->dy(2, 2), |
There was a problem hiding this comment.
warning: call to consteval function 'fmt::fstring<const double &, const double &, const Field2D &>::fstring<33UL>' is not a constant expression [clang-diagnostic-error]
output.write("dx = {:e}, dy = {:e}, dz = {:e}\n", coords->dx(2, 2), coords->dy(2, 2),
^Additional context
externalpackages/fmt/include/fmt/base.h:1766: non-constexpr function 'report_error' cannot be used in a constant expression
report_error(message);
^externalpackages/fmt/include/fmt/base.h:1650: in call to 'handler.on_error(&"unknown format specifier"[0])'
return handler.on_error("unknown format specifier"), end;
^externalpackages/fmt/include/fmt/base.h:1663: in call to 'parse_replacement_field<char, fmt::detail::format_string_checker<char, 3, 0, false> &>(&"dx = {:e}, dy = {:e}, dz = {:e}\n"[29], &"dx = {:e}, dy = {:e}, dz = {:e}\n"[32], checker(s, arg_pack()))'
begin = p = parse_replacement_field(p - 1, end, handler);
^externalpackages/fmt/include/fmt/base.h:2758: in call to 'parse_format_string<char, fmt::detail::format_string_checker<char, 3, 0, false>>({&"dx = {:e}, dy = {:e}, dz = {:e}\n"[0], 32}, checker(s, arg_pack()))'
if (FMT_USE_CONSTEVAL) parse_format_string<char>(s, checker(s, arg_pack()));
^tests/MMS/GBS/gbs.cxx:301: in call to 'fstring<33UL>("dx = {:e}, dy = {:e}, dz = {:e}\n")'
output.write("dx = {:e}, dy = {:e}, dz = {:e}\n", coords->dx(2, 2), coords->dy(2, 2),
^externalpackages/fmt/include/fmt/base.h:662: declared here
FMT_NORETURN FMT_API void report_error(const char* message);
^| if (ionvis) { | ||
| Field3D tau_i = Omega_ci * tau_i0 * pow(Ti, 1.5) / Ne; | ||
| Gi = -(0.96 * Ti * Ne * tau_i) * (2. * Grad_par(Vi) + C(phi) / coords->Bxy); | ||
| Gi = -(0.96 * Ti * Ne * tau_i) * (2. * Grad_par(Vi) + C(phi) / coords->Bxy()); |
There was a problem hiding this comment.
warning: no header providing "Grad_par" is directly included [misc-include-cleaner]
Gi = -(0.96 * Ti * Ne * tau_i) * (2. * Grad_par(Vi) + C(phi) / coords->Bxy());
^| const auto& I = tokamak_coords.I_unnormalised; | ||
|
|
||
| B0 = tokamak_coords.Bxy; | ||
| B0 = tokamak_coords.Bxy(); |
There was a problem hiding this comment.
warning: no matching function for call to object of type 'const FieldMetric' (aka 'const Field2D') [clang-diagnostic-error]
B0 = tokamak_coords.Bxy();
^Additional context
include/bout/field2d.hxx:351: candidate function not viable: requires single argument 'i', but no arguments were provided
BOUT_DEVICE inline BoutReal operator()(int i) { return View()(i); }
^include/bout/field2d.hxx:352: candidate function not viable: requires single argument 'i', but no arguments were provided
BOUT_DEVICE inline BoutReal operator()(int i) const { return View()(i); }
^include/bout/field2d.hxx:235: candidate function not viable: requires 2 arguments, but 0 were provided
inline BoutReal& operator()(int jx, int jy) {
^include/bout/field2d.hxx:249: candidate function not viable: requires 2 arguments, but 0 were provided
inline const BoutReal& operator()(int jx, int jy) const {
^include/bout/field2d.hxx:268: candidate function not viable: requires 3 arguments, but 0 were provided
BoutReal& operator()(int jx, int jy, [[maybe_unused]] int jz) {
^include/bout/field2d.hxx:271: candidate function not viable: requires 3 arguments, but 0 were provided
const BoutReal& operator()(int jx, int jy, [[maybe_unused]] int jz) const {
^Convert input to Field3DParallel before dividing by Bxy, to ensure that yup/ydown fields are retained.
| } | ||
|
|
||
| ddt(U) -= 10 * (SQ(SQ(coords->dx)) * D4DX4(U) + SQ(SQ(coords->dz)) * D4DZ4(U)); | ||
| ddt(U) -= 10 * (SQ(SQ(coords->dx())) * D4DX4(U) + SQ(SQ(coords->dz())) * D4DZ4(U)); |
There was a problem hiding this comment.
warning: no header providing "SQ" is directly included [misc-include-cleaner]
ddt(U) -= 10 * (SQ(SQ(coords->dx())) * D4DX4(U) + SQ(SQ(coords->dz())) * D4DZ4(U));
^| // Central differencing | ||
| ddt(f) = DDX(g, CELL_CENTRE); // + 20*SQ(coord->dx)*D2DX2(f); | ||
| ddt(g) = DDX(f, CELL_XLOW); // + 20*SQ(coord->dx)*D2DX2(g); | ||
| ddt(f) = DDX(g, CELL_CENTRE); // + 20*SQ(coord->dx())*D2DX2(f); |
There was a problem hiding this comment.
warning: no header providing "CELL_CENTRE" is directly included [misc-include-cleaner]
ddt(f) = DDX(g, CELL_CENTRE); // + 20*SQ(coord->dx())*D2DX2(f);
^|
|
||
| // just define a macro for V_E dot Grad | ||
| #define vE_Grad(f, p) (b0xGrad_dot_Grad(p, f) / coord->Bxy) | ||
| #define vE_Grad(f, p) (b0xGrad_dot_Grad(p, f) / coord->Bxy()) |
There was a problem hiding this comment.
warning: function-like macro 'vE_Grad' used; consider a 'constexpr' template function [cppcoreguidelines-macro-usage]
#define vE_Grad(f, p) (b0xGrad_dot_Grad(p, f) / coord->Bxy())
^|
|
||
| // just define a macro for V_E dot Grad | ||
| #define vE_Grad(f, p) (b0xGrad_dot_Grad(p, f) / coord->Bxy) | ||
| #define vE_Grad(f, p) (b0xGrad_dot_Grad(p, f) / coord->Bxy()) |
There was a problem hiding this comment.
warning: no header providing "b0xGrad_dot_Grad" is directly included [misc-include-cleaner]
#define vE_Grad(f, p) (b0xGrad_dot_Grad(p, f) / coord->Bxy())
^|
|
||
| // DENSITY EQUATION | ||
| ddt(Ni) = -b0xGrad_dot_Grad(phi, Ni0) / coord->Bxy; | ||
| ddt(Ni) = -b0xGrad_dot_Grad(phi, Ni0) / coord->Bxy(); |
There was a problem hiding this comment.
warning: no header providing "b0xGrad_dot_Grad" is directly included [misc-include-cleaner]
tests/integrated/test-interchange-instability/2fluid.cxx:5:
- #include <bout/bout.hxx>
+ #include "bout/difops.hxx"
+ #include <bout/bout.hxx>|
|
||
| // DENSITY EQUATION | ||
| ddt(Ni) = -b0xGrad_dot_Grad(phi, Ni0) / coord->Bxy; | ||
| ddt(Ni) = -b0xGrad_dot_Grad(phi, Ni0) / coord->Bxy(); |
There was a problem hiding this comment.
warning: no header providing "ddt" is directly included [misc-include-cleaner]
tests/integrated/test-interchange-instability/2fluid.cxx:5:
- #include <bout/bout.hxx>
+ #include "bout/field3d.hxx"
+ #include <bout/bout.hxx>|
|
||
| // VORTICITY | ||
| ddt(rho) = 2.0 * coord->Bxy * b0xcv * Grad(pei); | ||
| ddt(rho) = 2.0 * coord->Bxy() * b0xcv * Grad(pei); |
There was a problem hiding this comment.
warning: no header providing "Grad" is directly included [misc-include-cleaner]
tests/integrated/test-interchange-instability/2fluid.cxx:5:
- #include <bout/bout.hxx>
+ #include "bout/vecops.hxx"
+ #include <bout/bout.hxx>Should now be consistent with `next` branch. Avoids calculation of parallel slices for derived metric quantities.
Makes lazy getters thread-safe by adding named critical sections. Named so that nested sections (code in critical section calling a function that also contains a critical section) doesn't lead to deadlock.
|
Don't know what's going on here. Tests pass except with OpenMP. |
3 participants
Replaces #2873
Refactor
Coordinateswith a more encapsulated interface. Derived quantitiesare now automatically updated when the metric components are changed.
The most impactful change is the conversion of the metric components and grid
sizes to
const&methods:Among other things, this change allows us to defer computation of many of the
derived quantities (such as Christoffel symbols) and remove the
Coordinates::geometrymethod.We can now also move all of the derivative operators out of
Coordinatesinto free functions (we could've done most of these already, but we neededCoordinates::DD*for the metric derivatives, which we now defer computing)