Coordinates: Replacing sqrt(g_22) with JB #3027

bendudson · 2024-11-16T00:52:54Z

Aiming to consistently handle left-handed coordinate systems, where J is negative. These happen in the standard BOUT++ field-aligned coordinates when the poloidal field is negative.

Implemented on top of refactor-coordinates #2873

Aiming to consistently handle left-handed coordinate systems, where J is negative. These happen in the standard BOUT++ field-aligned coordinates when the poloidal field is negative.

github-actions

clang-tidy made some suggestions

src/mesh/difops.cxx

ZedThree · 2024-11-18T15:23:06Z

Ah, this currently doesn't work because we've not added the overload for Bxy(x, y, z)

dschwoerer · 2024-11-18T15:50:04Z

Does this fix things for FCI where JB is not equal to sqrt(g22)? I wanted to discuss this with Ben on Wednesday.

bendudson · 2024-11-18T21:38:46Z

Does this fix things for FCI where JB is not equal to sqrt(g22)? I wanted to discuss this with Ben on Wednesday.

Hey @dschwoerer I hope this will help with FCI too: I'm trying to remove all use of sqrt(g_22) from the code, and assumptions that J is positive. B is always positive (as it's the magnitude of B), but J can be positive or negative. In the tokamak case J < 0 when Bp < 0.

I haven't yet looked at the boundary conditions, or how to handle parallel components of vector quantities: When J < 0 it means that the y coordinate is in the opposite direction to B. Hence the 'y' component of a vector v is in the opposite direction to b dot v. It'll take some time to work out the implications for the operators e.g. Grad_par() takes a scalar argument, whereas Div_par() takes a vector component (b dot v).

ZedThree · 2024-11-20T14:04:52Z

We also have some checks that J is positive that I guess need to be loosened?

bendudson · 2024-11-20T16:30:57Z

We also have some checks that J is positive that I guess need to be loosened?

Ah yes. I think we should check that all elements of J have the same sign: All positive or all negative.

dschwoerer

It seems some of the cases are (still) only valid for Clebsch.
I have not looked in too much detail, so if you disagree I can look in more detail ...

dschwoerer · 2024-11-21T07:03:28Z

src/mesh/difops.cxx


  if (!f.hasParallelSlices()) {
-    return metric->Bxy() * FDDY(v, f / Bxy_floc, outloc, method) / sqrt(metric->g_22());
+    return FDDY(v, f / Bxy_floc, outloc, method) / metric->J();


Should it really be Bxy_floc? I thought this term comes from the curvi-linear geometry, and thus should be independent of B (for non clebsch)

dschwoerer · 2024-11-21T07:03:49Z

src/mesh/difops.cxx

  f_B.yup() = f.yup() / Bxy_floc;
  f_B.ydown() = f.ydown() / Bxy_floc;
-  return metric->Bxy() * FDDY(v, f_B, outloc, method) / sqrt(metric->g_22());
+  return FDDY(v, f_B, outloc, method) / metric->J();


dschwoerer · 2024-11-21T07:04:32Z

src/mesh/difops.cxx

+  result /= SQ(metric->J())
+            * metric->Bxy(); // J^2 B same sign for left & right handed coordinates


move comment above to avoid strange formatting?

dschwoerer · 2024-11-21T07:14:39Z

src/mesh/difops.cxx

  Field3D result = VDDX(vx, A, outloc) + VDDY(vy, A, outloc) + VDDZ(vz, A, outloc);

-  result /= (metric->J() * sqrt(metric->g_22()));
+  result /= SQ(metric->J()) * metric->Bxy(); // J^2 B


Is that not only true for Clebsch?

dschwoerer

Here are some fixes for the compilation errors.

Alternative is to add the function to metric, like J and g_22 ...

dschwoerer · 2024-11-21T08:45:31Z

src/mesh/difops.cxx

-        BoutReal by = 1. / sqrt(metric->g_22(x, y, z));
+
+        // Note: by is negative in a left-handed coordinate system
+        BoutReal by = 1. / (metric->J(x, y, z) * metric->Bxy(x, y, z));


Suggested change

BoutReal by = 1. / (metric->J(x, y, z) * metric->Bxy(x, y, z));

BoutReal by = 1. / (metric->J(x, y, z) * metric->Bxy()(x, y, z));

dschwoerer · 2024-11-21T08:45:55Z

src/mesh/difops.cxx

        BoutReal const fluxRight =
-            fR * vR * (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)));
+            fR * vR * 2 / (coord->Bxy(i, j, k) + coord->Bxy(i, j + 1, k));


Suggested change

fR * vR * 2 / (coord->Bxy(i, j, k) + coord->Bxy(i, j + 1, k));

fR * vR * 2 / (coord->Bxy()(i, j, k) + coord->Bxy()(i, j + 1, k));

dschwoerer · 2024-11-21T08:46:10Z

src/mesh/difops.cxx

        BoutReal const fluxLeft =
-            fL * vL * (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)));
+            fL * vL * 2 / (coord->Bxy(i, j, k) + coord->Bxy(i, j - 1, k));


Suggested change

fL * vL * 2 / (coord->Bxy(i, j, k) + coord->Bxy(i, j - 1, k));

fL * vL * 2 / (coord->Bxy()(i, j, k) + coord->Bxy()(i, j - 1, k));

bendudson · 2024-11-21T16:44:39Z

It seems some of the cases are (still) only valid for Clebsch. I have not looked in too much detail, so if you disagree I can look in more detail ...

@dschwoerer I agree: The identity sqrt(g_22) = JB is only for Clebsch coordinates. In general how should the 'parallel' direction be defined? Perhaps we use sign(J) * sqrt(g_22) as the parallel metric, which would be consistent with the Clebsch case but I think more general. That would mean that if J > 0 then b is in the same direction as e_y, but if J < 0 then they are in opposite directions.

bendudson · 2025-01-29T20:25:50Z

Note: Mesh checks that J > 0. Change to check that J has the same sign (+ or -) in all cells.

Coordinates: Replacing sqrt(g_22) with JB

14e7eda

Aiming to consistently handle left-handed coordinate systems, where J is negative. These happen in the standard BOUT++ field-aligned coordinates when the poloidal field is negative.

github-actions bot reviewed Nov 16, 2024

View reviewed changes

Apply clang-format changes

19ee7bf

ZedThree approved these changes Nov 18, 2024

View reviewed changes

bendudson marked this pull request as draft November 18, 2024 21:38

dschwoerer reviewed Nov 21, 2024

View reviewed changes

		result /= SQ(metric->J())
		* metric->Bxy(); // J^2 B same sign for left & right handed coordinates

	BoutReal by = 1. / (metric->J(x, y, z) * metric->Bxy(x, y, z));
	BoutReal by = 1. / (metric->J(x, y, z) * metric->Bxy()(x, y, z));

	fR * vR * 2 / (coord->Bxy(i, j, k) + coord->Bxy(i, j + 1, k));
	fR * vR * 2 / (coord->Bxy()(i, j, k) + coord->Bxy()(i, j + 1, k));

	fL * vL * 2 / (coord->Bxy(i, j, k) + coord->Bxy(i, j - 1, k));
	fL * vL * 2 / (coord->Bxy()(i, j, k) + coord->Bxy()(i, j - 1, k));

Coordinates: Replacing sqrt(g_22) with JB #3027

Are you sure you want to change the base?

Coordinates: Replacing sqrt(g_22) with JB #3027

Uh oh!

Conversation

bendudson commented Nov 16, 2024

Uh oh!

github-actions bot left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

ZedThree commented Nov 18, 2024

Uh oh!

dschwoerer commented Nov 18, 2024 via email

Uh oh!

bendudson commented Nov 18, 2024

Uh oh!

ZedThree commented Nov 20, 2024

Uh oh!

bendudson commented Nov 20, 2024

Uh oh!

dschwoerer left a comment

Choose a reason for hiding this comment

Uh oh!

dschwoerer Nov 21, 2024

Choose a reason for hiding this comment

Uh oh!

dschwoerer Nov 21, 2024

Choose a reason for hiding this comment

Uh oh!

dschwoerer Nov 21, 2024

Choose a reason for hiding this comment

Uh oh!

dschwoerer Nov 21, 2024

Choose a reason for hiding this comment

Uh oh!

dschwoerer left a comment • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

dschwoerer Nov 21, 2024

Choose a reason for hiding this comment

Uh oh!

dschwoerer Nov 21, 2024

Choose a reason for hiding this comment

Uh oh!

dschwoerer Nov 21, 2024

Choose a reason for hiding this comment

Uh oh!

bendudson commented Nov 21, 2024

Uh oh!

bendudson commented Jan 29, 2025

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

dschwoerer left a comment •

edited

Loading