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

@@ -285,7 +286,7 @@ Field3D Div_par_flux(const Field3D& v, const Field3D& f, CELL_LOC outloc,
  auto Bxy_floc = f.getCoordinates()->Bxy();

  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

@@ -595,7 +597,7 @@ Field3D b0xGrad_dot_Grad(const Field3D& phi, const Field3D& A, CELL_LOC outloc)

  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


        // Calculate at left boundary (y-1/2)
        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.

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

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Coordinates: Replacing sqrt(g_22) with JB #3027

Coordinates: Replacing sqrt(g_22) with JB #3027

bendudson commented Nov 16, 2024

github-actions bot left a comment

ZedThree commented Nov 18, 2024

dschwoerer commented Nov 18, 2024 via email

bendudson commented Nov 18, 2024

ZedThree commented Nov 20, 2024

bendudson commented Nov 20, 2024

dschwoerer left a comment

dschwoerer Nov 21, 2024

dschwoerer Nov 21, 2024

dschwoerer Nov 21, 2024

dschwoerer Nov 21, 2024

dschwoerer left a comment •

edited

Loading

dschwoerer Nov 21, 2024

dschwoerer Nov 21, 2024

dschwoerer Nov 21, 2024

bendudson commented Nov 21, 2024

		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

Conversation

bendudson commented Nov 16, 2024

github-actions bot left a comment

Choose a reason for hiding this comment

ZedThree commented Nov 18, 2024

dschwoerer commented Nov 18, 2024 via email

bendudson commented Nov 18, 2024

ZedThree commented Nov 20, 2024

bendudson commented Nov 20, 2024

dschwoerer left a comment

Choose a reason for hiding this comment

dschwoerer Nov 21, 2024

Choose a reason for hiding this comment

dschwoerer Nov 21, 2024

Choose a reason for hiding this comment

dschwoerer Nov 21, 2024

Choose a reason for hiding this comment

dschwoerer Nov 21, 2024

Choose a reason for hiding this comment

dschwoerer left a comment • edited Loading

Choose a reason for hiding this comment

dschwoerer Nov 21, 2024

Choose a reason for hiding this comment

dschwoerer Nov 21, 2024

Choose a reason for hiding this comment

dschwoerer Nov 21, 2024

Choose a reason for hiding this comment

bendudson commented Nov 21, 2024

dschwoerer left a comment •

edited

Loading