WIP: Tidying FCI implementation #3092

bendudson · 2025-03-15T22:15:27Z

Fix Div_par(f, v) for FCI

When using FCI yup/down fields, each poloidal plane uses a different coordinate system. Quantities like J therefore can't be averaged between planes. Magnetic field strength is a scalar that can be interpolated, and is used here to calculate the divergence of a parallel flow.

github-actions

clang-tidy made some suggestions

github-actions · 2025-03-15T22:17:30Z

src/mesh/difops.cxx

-            (fluxRight - fluxLeft) / (coord->dy(i, j, k) * coord->J(i, j, k));
-      }
+  Field3D result{emptyFrom(f)};
+  BOUT_FOR(i, f.getRegion("RGN_NOBNDRY")) {


warning: no header providing "BOUT_FOR" is directly included [misc-include-cleaner]

src/mesh/difops.cxx:27:

- #include <bout/assert.hxx> + #include "bout/region.hxx" + #include <bout/assert.hxx>

github-actions · 2025-03-15T22:17:30Z

src/mesh/difops.cxx

+
+#if CHECK > 0
+    if (!std::isfinite(result[i])) {
+      output.write("{} {} {} {}\n", f_up[i], v_up[i], f_down[i], v_down[i]);


warning: no header providing "output" is directly included [misc-include-cleaner]

src/mesh/difops.cxx:27:

- #include <bout/assert.hxx> + #include "bout/output.hxx" + #include <bout/assert.hxx>

github-actions · 2025-03-15T22:17:30Z

src/mesh/difops.cxx

+    if (!std::isfinite(result[i])) {
+      output.write("{} {} {} {}\n", f_up[i], v_up[i], f_down[i], v_down[i]);
+      output.write("{} {} {} {} {}\n", B[i], B_up[i], B_down[i], dy[i], sqrt(g_22[i]));
+      throw BoutException("Non-finite value in Div_");


warning: no header providing "BoutException" is directly included [misc-include-cleaner]

src/mesh/difops.cxx:25:

- #include "bout/build_defines.hxx" + #include "bout/boutexception.hxx" + #include "bout/build_defines.hxx"

Use Bxy rather than J on neighboring slices: B can be compared between slices, but J cannot (different coordinate system).

Each poloidal slice has a different coordinate system, so metrics can't be directly compared or averaged. Only Bxy and the parallel connection length (dy, g_22) make sense to calculate.

Don't calculate parallel slices in a derivative: The result is unlikely to be correct.

If the argument to `sqrt` has parallel slices, then the slices will be operated on and the result will have parallel slices. To avoid unnecessary work, discard slices before calling.

Returns a shallow copy without parallel slices. Enables user to avoid performing calculations on slices if not needed.

Always perform calculations in the yup/down slices if present in BOTH arguments. i.e. automagic always on for arithmetic operators. To avoid unnecessary work, use `withoutParallelSlices` to pass arguments without parallel slices.

dschwoerer

I think we should talk about this, I do not think we want to remove the part about the metric coefficients?

dschwoerer · 2025-03-17T08:27:35Z

include/bout/field.hxx

+      /* Operate on parallel slices */                                                   \
+      result.splitParallelSlicesAndAllocate();                                           \
+      for (size_t i{0}; i != f.numberParallelSlices(); ++i) {                            \
+        BOUT_FOR(d, result.getRegion(rgn)) { result.yup(i)[d] = func(f.yup(i)[d]); }     \


but here we need a shifted region (at least for FCI) ...

My thinking with this was that if we are operating over a given region of a field, then we also want to operate over all the yup/down fields connected to those points. Hence the region indices should be the same for yup/down. I guess you have a different model for how this should work?

dschwoerer · 2025-03-17T08:29:09Z

include/bout/field3d.hxx


+  /// Returns a shallow copy without parallel slices
+  Field3D withoutParallelSlices() const {
+    Field3D result{getMesh(), getLocation(), getDirections(), getRegionID()};


Would that not be more general:

Suggested change

Field3D result{getMesh(), getLocation(), getDirections(), getRegionID()};

auto result{emptyFrom(*this)};

dschwoerer · 2025-03-17T08:31:12Z

src/mesh/coordinates.cxx

-  return setName(coords->J / sqrt(coords->g_22) * Grad_par(f_B, outloc, method),
-                 "Div_par({:s})", f.name);
+  return setName(Bxy * Grad_par(f_B, outloc, method), "Div_par({:s})", f.name);
 }


But this about the metric coefficients, not about the strength about the metric field.
Only for Clebsch coordinates this is the same, but not for FCI.

It works because the area of a flux tube is inversely proportional to B. If we're mapping along flux tubes from one plane to the next, then the ratio of the flux tube areas is (inverse) the ratio of B on the two planes.

dschwoerer · 2025-03-17T08:34:44Z

src/mesh/difops.cxx

+    result[i] = B[i]
+                * ((f_up[i] * v_up[i] / B_up[i]) - (f_down[i] * v_down[i] / B_down[i]))
+                / (dy[i] * sqrt(g_22[i]));


This is only valid in clebsch

Ah no. This is always valid for all coordinate systems, not just Clebsch. It's just a vector identity from Div(B) = 0

Div(b v) = Div(B v / |B|) = Div(B) v/|B| + B dot Grad(v / |B|)

so the divergence of a parallel vector Div_par(v) is always B * Grad_par( v / B )

Because the mesh is locally field aligned - I see 👍

dschwoerer · 2025-03-17T08:36:05Z

src/mesh/difops.cxx

+      output.write("{} {} {} {}\n", f_up[i], v_up[i], f_down[i], v_down[i]);
+      output.write("{} {} {} {} {}\n", B[i], B_up[i], B_down[i], dy[i], sqrt(g_22[i]));


We probably want to move that to the exception message below, and also add i to make it easier to figure out why that is non-finite ...

Added comment to explain why it's ok to use the yup/down slices of the metric components if they are read from the grid file.

Iterated for `i < localmesh->ystart - abs_offset()` so when ystart is 1 and the boundary offset is 1, no boundary is applied.

Reverting some changes back to David's version.

Prefer to have consistent behavior so that code works and can be optimized without breaking. The field.withoutParallelSlices() method can be used to strip slices before calling an operation, if slices are not needed.

If the field has parallel slices, assign the value also to the slices.

applyParallelBoundary can be a no-op for non-fci

bendudson marked this pull request as draft March 15, 2025 22:15

Apply clang-format changes

81ee64e

github-actions bot reviewed Mar 15, 2025

View reviewed changes

bendudson and others added 7 commits March 15, 2025 15:21

Coordinates: Fix Div_par(f) for FCI

3a77437

Use Bxy rather than J on neighboring slices: B can be compared between slices, but J cannot (different coordinate system).

FCItransform: Only load some parallel metrics

1d78559

Each poloidal slice has a different coordinate system, so metrics can't be directly compared or averaged. Only Bxy and the parallel connection length (dy, g_22) make sense to calculate.

DDY: Remove FCI_AUTOMAGIC calcParallelSlices

e5834c1

Don't calculate parallel slices in a derivative: The result is unlikely to be correct.

Field: Extend operations to calculate in parallel slices

e058507

If the argument to `sqrt` has parallel slices, then the slices will be operated on and the result will have parallel slices. To avoid unnecessary work, discard slices before calling.

Field: Add withoutParallelSlices() method

a46a4f7

Returns a shallow copy without parallel slices. Enables user to avoid performing calculations on slices if not needed.

Apply clang-format changes

cedeb2b

field ops: Remove FCI_AUTOMAGIC flag

0b976ec

Always perform calculations in the yup/down slices if present in BOTH arguments. i.e. automagic always on for arithmetic operators. To avoid unnecessary work, use `withoutParallelSlices` to pass arguments without parallel slices.

dschwoerer reviewed Mar 17, 2025

View reviewed changes

bendudson and others added 10 commits March 17, 2025 21:12

FCI parallel: Revert deletion of coordinate slices

fd93d69

Added comment to explain why it's ok to use the yup/down slices of the metric components if they are read from the grid file.

parallel_boundary_region: Fix ITER macro

192d105

Iterated for `i < localmesh->ystart - abs_offset()` so when ystart is 1 and the boundary offset is 1, no boundary is applied.

Apply clang-format changes

8d79e50

Fixes for FCI: Load dy, improve errors, use J in Div_par

e23b104

Reverting some changes back to David's version.

Apply clang-format changes

f169b8c

Remove FCI_AUTOMAGIC compile option

ce70d86

Prefer to have consistent behavior so that code works and can be optimized without breaking. The field.withoutParallelSlices() method can be used to strip slices before calling an operation, if slices are not needed.

Remove stray #endif

3288808

Field3D: Fix operator=(BoutReal)

c840122

If the field has parallel slices, assign the value also to the slices.

applyParallelBoundary can be a no-op for non-fci

568ad42

Merge pull request #3202 from boutproject/fci-higher-order-tidy-noop

90c4e6c

applyParallelBoundary can be a no-op for non-fci

	Field3D result{getMesh(), getLocation(), getDirections(), getRegionID()};
	auto result{emptyFrom(*this)};

		output.write("{} {} {} {}\n", f_up[i], v_up[i], f_down[i], v_down[i]);
		output.write("{} {} {} {} {}\n", B[i], B_up[i], B_down[i], dy[i], sqrt(g_22[i]));

WIP: Tidying FCI implementation #3092

Are you sure you want to change the base?

WIP: Tidying FCI implementation #3092

Uh oh!

Conversation

bendudson commented Mar 15, 2025

Uh oh!

github-actions bot left a comment

Choose a reason for hiding this comment

Uh oh!

github-actions bot Mar 15, 2025

Choose a reason for hiding this comment

Uh oh!

github-actions bot Mar 15, 2025

Choose a reason for hiding this comment

Uh oh!

github-actions bot Mar 15, 2025

Choose a reason for hiding this comment

Uh oh!

dschwoerer left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants