Merge #909

909: Fix differentiability of local element map r=simonbyrne a=simonbyrne

I was experimenting with `IntrinsicMap` and noticed some sign errors in the metric terms due to the non-differentiable extra branches for handling of signed zeros. This is a simpler and better approach.

- [X] Code follows the [style guidelines](https://clima.github.io/ClimateMachine.jl/latest/DevDocs/CodeStyle/) OR N/A.
- [X] Unit tests are included OR N/A.
- [X] Code is exercised in an integration test OR N/A.
- [X] Documentation has been added/updated OR N/A.


Co-authored-by: Simon Byrne <[email protected]>

src/Meshes/cubedsphere.jl

@@ -243,18 +243,12 @@ function coordinates(
     ne = mesh.ne
     x, y, panel = elem.I
     # again, we want to arrange the calculation carefully so that signed zeros are correct
-    # - if isodd(ne) and ϕx == 0, then ξ1 == 0, and should have the same sign
-    ϕx = (ξ1 - FT(ne + 1 - 2 * x)) / ne
-    ϕy = (ξ2 - FT(ne + 1 - 2 * y)) / ne
     # - if iseven(ne) and ϕx == 0, then ξ1 == +/-1, and should have the _opposite_ sign
-    if iseven(ne)
-        if ϕx == 0
-            ϕx = copysign(ϕx, -ξ1)
-        end
-        if ϕy == 0
-            ϕy = copysign(ϕy, -ξ2)
-        end
-    end
+    # - if isodd(ne) and ϕx == 0, then ξ1 == 0, and should have the same sign
+    ox = 2 * x - 1 - ne
+    oy = 2 * y - 1 - ne
+    ϕx = (ox <= 0 ? -(-FT(ox) - ξ1) : ox + ξ1) / ne
+    ϕy = (oy <= 0 ? -(-FT(oy) - ξ2) : oy + ξ2) / ne
     return _coordinates(mesh, ϕx, ϕy, panel)
 end
 function coordinates(

test/Meshes/cubedsphere.jl

@@ -1,6 +1,6 @@
 using ClimaCore: Geometry, Domains, Meshes
 using Test
-using StaticArrays, SparseArrays, LinearAlgebra
+using StaticArrays, SparseArrays, LinearAlgebra, ForwardDiff
 
 
 @testset "opposing face" begin
@@ -182,6 +182,30 @@ end
             SVector(0.5, -1.0),
         ).x3 === +0.0
 
+        # derivative handling
+        M = ForwardDiff.jacobian(SVector(-1.0, 1.0)) do ξ
+            Geometry.components(
+                Meshes.coordinates(mesh, CartesianIndex(1, 2, 3), ξ),
+            )
+        end
+        @test M[1, 1] < 0
+        @test M[2, 1] == 0
+        @test M[3, 1] > 0
+        @test M[1, 2] == 0
+        @test M[2, 2] < 0
+        @test M[3, 2] == 0
+
+        M = ForwardDiff.jacobian(SVector(-1.0, -1.0)) do ξ
+            Geometry.components(
+                Meshes.coordinates(mesh, CartesianIndex(1, 3, 3), ξ),
+            )
+        end
+        @test M[1, 1] < 0
+        @test M[2, 1] == 0
+        @test M[3, 1] > 0
+        @test M[1, 2] == 0
+        @test M[2, 2] < 0
+        @test M[3, 2] == 0
     end
 end
 
@@ -214,6 +238,30 @@ end
             CartesianIndex(1, 2, 1),
             SVector(0.5, +0.0),
         ).x3 === +0.0
+
+        # derivative handling
+        M = ForwardDiff.jacobian(SVector(-1.0, -0.0)) do ξ
+            Geometry.components(
+                Meshes.coordinates(mesh, CartesianIndex(1, 2, 3), ξ),
+            )
+        end
+        @test M[1, 1] < 0
+        @test M[2, 1] == 0
+        @test M[3, 1] > 0
+        @test M[1, 2] == 0
+        @test M[2, 2] < 0
+        @test M[3, 2] == 0
+        M = ForwardDiff.jacobian(SVector(-1.0, +0.0)) do ξ
+            Geometry.components(
+                Meshes.coordinates(mesh, CartesianIndex(1, 2, 3), ξ),
+            )
+        end
+        @test M[1, 1] < 0
+        @test M[2, 1] == 0
+        @test M[3, 1] > 0
+        @test M[1, 2] == 0
+        @test M[2, 2] < 0
+        @test M[3, 2] == 0
     end
 end