Enabling rotated co-ordinates

SCOREC · Sep 17, 2024 · 8fbb43b · 8fbb43b
1 parent 58d0bbe
commit 8fbb43b
Show file tree

Hide file tree

Showing 7 changed files with 89 additions and 108 deletions.
diff --git a/src/pmpo_materialPoints.hpp b/src/pmpo_materialPoints.hpp
@@ -206,26 +206,31 @@ class MaterialPoints {
         auto tgtPosRotLatLon = MPs->get<MPF_Tgt_Pos_Rot_Lat_Lon>();
         auto tgtPosXYZ = MPs->get<MPF_Tgt_Pos_XYZ>();
         auto rotLatLonIncr = MPs->get<MPF_Rot_Lat_Lon_Incr>();
-
+        auto is_rotated=getRotatedFlag();
+
         auto updateRotLatLon = PS_LAMBDA(const int& elm, const int& mp, const int& mask){
             if(mask){
                 auto rotLat = curPosRotLatLon(mp,0) + rotLatLonIncr(mp,0); // phi
                 auto rotLon = curPosRotLatLon(mp,1) + rotLatLonIncr(mp,1); // lambda
-                auto geoLat = rotLat;
+
+		tgtPosRotLatLon(mp,0) = rotLat;
+		tgtPosRotLatLon(mp,1) = rotLon;
+
+		//These are updated if rotation enabled
+		auto geoLat = rotLat;
                 auto geoLon = rotLon;
-                tgtPosRotLatLon(mp,0) = geoLat;
-                tgtPosRotLatLon(mp,1) = geoLon;
-                // x = cosLon cosLat, y = sinLon cosLat, z = sinLat
+                if(is_rotated){
+		    auto xyz_rot = xyz_from_lat_lon(rotLat, rotLon, radius);
+		    auto xyz_geo = grid_rotation_backward(xyz_rot);
+		    lat_lon_from_xyz(geoLat, geoLon, xyz_geo, radius);
+		}
+
+		// x = cosLon cosLat, y = sinLon cosLat, z = sinLat
                 tgtPosXYZ(mp,0) = radius * std::cos(geoLon) * std::cos(geoLat);
                 tgtPosXYZ(mp,1) = radius * std::sin(geoLon) * std::cos(geoLat);
                 tgtPosXYZ(mp,2) = radius * std::sin(geoLat); 
             } 
         };
-        if(isRotatedFlag){
-            //TODO rotation lat lon calc
-            fprintf(stderr, "rotational lat lon in MP is not support yet!");
-            PMT_ALWAYS_ASSERT(false);
-        } 
         ps::parallel_for(MPs, updateRotLatLon,"updateRotationalLatitudeLongitude"); 
         pumipic::RecordTime("PolyMPO_updateRotLatLonAndXYZ2Tgt", timer.seconds());
     } 

diff --git a/src/pmpo_utils.hpp b/src/pmpo_utils.hpp
@@ -210,6 +210,34 @@ double sphericalTriangleArea2(Vec3d &a, Vec3d &b, Vec3d &c, double radius){
 //this is a lazy comparison and shouldn't be relied on beyond simple testing
 bool isEqual(double a, double b, double tol=1e-9);
 
+
+KOKKOS_INLINE_FUNCTION
+Vec3d xyz_from_lat_lon(double lat, double lon, double r){
+    Vec3d xyz;
+    xyz[0] = r*cos(lat)*cos(lon);
+    xyz[1] = r*cos(lat)*sin(lon);
+    xyz[2] = r*sin(lat);
+
+    return xyz;
+}
+
+KOKKOS_INLINE_FUNCTION
+Vec3d grid_rotation_backward(Vec3d& xyz_input){
+    Vec3d xyz_output;
+    xyz_output[0] =  xyz_input[2];
+    xyz_output[1] =  xyz_input[1];
+    xyz_output[2] = -xyz_input[0];
+
+    return xyz_output;
+}
+
+KOKKOS_INLINE_FUNCTION
+void lat_lon_from_xyz(double& lat, double& lon, Vec3d& xyz, double r){
+    lon = Kokkos::atan2(xyz[1], xyz[0]);
+    lat = Kokkos::asin(xyz[2]/r);
+} 
+
+
 }//namespace polyMPO end
 
 #endif

diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
@@ -69,7 +69,6 @@ pmpo_add_fortran_exe(testFortranMPAppIDs testFortranMPAppIDs.f90)
 pmpo_add_fortran_exe(testFortranMPAdvection testFortranMPAdvection.f90)
 pmpo_add_fortran_exe(testFortranCreateRebuildMPs testFortranCreateRebuildMPs.f90)
 pmpo_add_fortran_exe(testFortranMPReconstruction testFortranMPReconstruction.f90)
-pmpo_add_fortran_exe(testFortraAdjacency testFortraAdjacency.f90)
 
 #add test
 find_program(MPIRUN_EXECUTABLE NAMES mpirun) 

diff --git a/test/calculateDisplacement.f90 b/test/calculateDisplacement.f90
@@ -1,4 +1,5 @@
-subroutine calcSurfDispIncr(mpMesh, latVertex, lonVertex, nEdgesOnCell, verticesOnCell, nVertices, sphereRadius)
+subroutine calcSurfDispIncr(mpMesh, latVertex, lonVertex, nEdgesOnCell, verticesOnCell, nVertices, sphereRadius, &
+                rotated_coords)
   use :: polympo
   use :: readMPAS
   use :: iso_c_binding
@@ -11,7 +12,8 @@ subroutine calcSurfDispIncr(mpMesh, latVertex, lonVertex, nEdgesOnCell, vertices
   integer :: i, j, nVertices, nCompsDisp
   real(kind=MPAS_RKIND), dimension(:,:), pointer :: dispIncr
   type(c_ptr) :: mpMesh
-
+  LOGICAL, INTENT(IN) :: rotated_coords 
+
   nCompsDisp = 2
   allocate(dispIncr(nCompsDisp,nVertices))
 
@@ -29,12 +31,18 @@ subroutine calcSurfDispIncr(mpMesh, latVertex, lonVertex, nEdgesOnCell, vertices
   end do
 
   deltaLon = maxlon - minlon
-
-  do i = 1,nVertices
-    dispIncr(1,i) = sphereRadius*cos(latVertex(i))*deltaLon
-    dispIncr(2,i) = 0.0_MPAS_RKIND
-  end do
+  IF (rotated_coords) THEN
+      do i = 1,nVertices
+          dispIncr(2,i) = sphereRadius*cos(latVertex(i))*deltaLon
+          dispIncr(1,i) = 0.0_MPAS_RKIND
+      end do
+  ELSE 
+      do i = 1,nVertices
+          dispIncr(1,i) = sphereRadius*cos(latVertex(i))*deltaLon
+          dispIncr(2,i) = 0.0_MPAS_RKIND
+      end do
+  END IF
   call polympo_setMeshVtxOnSurfDispIncr(mpMesh,nCompsDisp,nVertices,c_loc(dispIncr))
 
   deallocate(dispIncr)
-end subroutine
+end subroutine
diff --git a/test/testFortraAdjacency.f90 b/test/testFortraAdjacency.f90
diff --git a/test/testFortranMPAdvection.f90 b/test/testFortranMPAdvection.f90
@@ -2,7 +2,8 @@ module advectionTests
   contains
   include "calculateDisplacement.f90"
 
-  subroutine runAdvectionTest(mpMesh, numPush, latVertex, lonVertex, nEdgesOnCell, verticesOnCell, nVertices, sphereRadius)
+  subroutine runAdvectionTest(mpMesh, numPush, latVertex, lonVertex, nEdgesOnCell, verticesOnCell, nVertices, sphereRadius, &
+                  rotated_coords)
     use :: polympo
     use :: readMPAS
     use :: iso_c_binding
@@ -14,17 +15,18 @@ subroutine runAdvectionTest(mpMesh, numPush, latVertex, lonVertex, nEdgesOnCell,
     integer, dimension(:), pointer :: nEdgesOnCell
     integer, dimension(:,:), pointer :: verticesOnCell
     real(kind=MPAS_RKIND) :: sphereRadius
-
+    LOGICAL, INTENT(IN) :: rotated_coords
 
     do i = 1, numPush
-      call calcSurfDispIncr(mpMesh, latVertex, lonVertex, nEdgesOnCell, verticesOnCell, nVertices, sphereRadius)
+      call calcSurfDispIncr(mpMesh, latVertex, lonVertex, nEdgesOnCell, verticesOnCell, nVertices, sphereRadius, &
+              rotated_coords)
       call polympo_push(mpMesh)
     end do
 
   end subroutine
 
   subroutine runReconstructionTest(mpMesh, numMPs, numPush, nCells, nVertices, mp2Elm, &
-                                  latVertex, lonVertex, nEdgesOnCell, verticesOnCell, sphereRadius)
+                                  latVertex, lonVertex, nEdgesOnCell, verticesOnCell, sphereRadius, rotated_coords)
     use :: polympo
     use :: readMPAS
     use :: iso_c_binding
@@ -41,6 +43,7 @@ subroutine runReconstructionTest(mpMesh, numMPs, numPush, nCells, nVertices, mp2
     real(kind=MPAS_RKIND) :: sphereRadius
     real(kind=MPAS_RKIND) :: TEST_VAL = 1.1_MPAS_RKIND
     real(kind=MPAS_RKIND) :: TOLERANCE = 0.0001_MPAS_RKIND
+    LOGICAL, INTENT(IN) :: rotated_coords
 
     allocate(mpMass(1,numMPs))
     allocate(mpVel(2,numMPs))
@@ -66,7 +69,8 @@ subroutine runReconstructionTest(mpMesh, numMPs, numPush, nCells, nVertices, mp2
     ! Test push reconstruction
 
     do j = 1, numPush
-      call calcSurfDispIncr(mpMesh, latVertex, lonVertex, nEdgesOnCell, verticesOnCell, nVertices, sphereRadius)
+      call calcSurfDispIncr(mpMesh, latVertex, lonVertex, nEdgesOnCell, verticesOnCell, nVertices, sphereRadius, &
+              rotated_coords)
       call polympo_setReconstructionOfMass(mpMesh,0,polympo_getMeshFElmType())
       call polympo_setReconstructionOfMass(mpMesh,0,polympo_getMeshFVtxType())
       call polympo_setReconstructionOfVel(mpMesh,0,polympo_getMeshFVtxType())
@@ -90,7 +94,8 @@ subroutine runReconstructionTest(mpMesh, numMPs, numPush, nCells, nVertices, mp2
     deallocate(meshElmMass)
   end subroutine
 
-  subroutine runApiTest(mpMesh, numMPs, nVertices, nCells, numPush, mpLatLon, mpPosition, xVertex, yVertex, zVertex, latVertex)
+  subroutine runApiTest(mpMesh, numMPs, nVertices, nCells, numPush, mpLatLon, mpPosition, xVertex, yVertex, zVertex, latVertex, & 
+                  rotated_coords)
     use :: polympo
     use :: readMPAS
     use :: iso_c_binding
@@ -102,6 +107,7 @@ subroutine runApiTest(mpMesh, numMPs, nVertices, nCells, numPush, mpLatLon, mpPo
     real(kind=MPAS_RKIND), dimension(:), pointer :: meshVtxMass, meshElmMass, meshVtxVel
     real(kind=MPAS_RKIND), dimension(:), pointer :: latVertex
     real(kind=MPAS_RKIND) :: TEST_VAL = 1.1_MPAS_RKIND
+    LOGICAL, INTENT(IN) :: rotated_coords
 
     nCompsDisp = 2
     allocate(dispIncr(nCompsDisp,nVertices))
@@ -176,6 +182,7 @@ program main
   integer, parameter :: MP_ACTIVE = 1
   integer, parameter :: MP_INACTIVE = 0
   integer, parameter :: INVALID_ELM_ID = -1
+  LOGICAL :: rotated_coords = .TRUE.
 
   call mpi_init(ierr)
   call mpi_comm_rank(mpi_comm_handle, self, ierr)
@@ -292,15 +299,25 @@ program main
   call assert(numMPsCount == numMPs, "num mps miscounted")
 
   call polympo_createMPs(mpMesh,nCells,numMPs,c_loc(mpsPerElm),c_loc(mp2Elm),c_loc(isMPActive))
+
+  if (rotated_coords) then
+      call polympo_setMPLatLonRotatedFlag(mpMesh, 1)
+  else
+      call polympo_setMPLatLonRotatedFlag(mpMesh, 0)
+  endif
+
+
   call polympo_setMPRotLatLon(mpMesh,2,numMPs,c_loc(mpLatLon))
   call polympo_setMPPositions(mpMesh,3,numMPs,c_loc(mpPosition))
 
-  ! call runAdvectionTest(mpMesh, numPush, latVertex, lonVertex, nEdgesOnCell, verticesOnCell, nVertices, sphereRadius)
+  call runAdvectionTest(mpMesh, numPush, latVertex, lonVertex, nEdgesOnCell, verticesOnCell, nVertices, sphereRadius, &
+          rotated_coords)
 
-  ! call runReconstructionTest(mpMesh, numMPs, numPush, nCells, nVertices, mp2Elm, &
-  !                                 latVertex, lonVertex, nEdgesOnCell, verticesOnCell, sphereRadius)
+  call runReconstructionTest(mpMesh, numMPs, numPush, nCells, nVertices, mp2Elm, &
+                                   latVertex, lonVertex, nEdgesOnCell, verticesOnCell, sphereRadius, rotated_coords)
 
-  call runApiTest(mpMesh, numMPs, nVertices, nCells, numPush, mpLatLon, mpPosition, xVertex, yVertex, zVertex, latVertex)
+  call runApiTest(mpMesh, numMPs, nVertices, nCells, numPush, mpLatLon, mpPosition, xVertex, yVertex, zVertex, latVertex, & 
+          rotated_coords)
 
   call polympo_summarizeTime();
 

diff --git a/test/testFortranMPReconstruction.f90 b/test/testFortranMPReconstruction.f90
@@ -36,6 +36,7 @@ program main
   integer, parameter :: MP_ACTIVE = 1
   integer, parameter :: MP_INACTIVE = 0
   integer, parameter :: INVALID_ELM_ID = -1
+  LOGICAL :: rotated_coords = .FALSE.
 
   call mpi_init(ierr)
   call mpi_comm_rank(mpi_comm_handle, self, ierr)
@@ -121,7 +122,8 @@ program main
   ! Test push reconstruction
 
   do j = 1, 5
-    call calcSurfDispIncr(mpMesh, latVertex, lonVertex, nEdgesOnCell, verticesOnCell, nVertices, sphereRadius)
+    call calcSurfDispIncr(mpMesh, latVertex, lonVertex, nEdgesOnCell, verticesOnCell, nVertices, sphereRadius, &
+            rotated_coords)
     call polympo_setReconstructionOfMass(mpMesh,0,polympo_getMeshFElmType())
     call polympo_setReconstructionOfMass(mpMesh,0,polympo_getMeshFVtxType())
     call polympo_push(mpMesh)