PR Change 3

Author: dhyan1272

src/pmpo_MPMesh_assembly.hpp

@@ -157,11 +157,11 @@ void MPMesh::assemblyVtx1() {
 
     //double f_norm=A.frobeniusNorm();
     double A_trace = A.trace();
-    Matrix4d A_new = {v0, v1, v2, v3};
-    A_new.regularize(A_trace*1e-8);
+    Matrix4d A_regularized = {v0, v1, v2, v3};
+    A_regularized.addToDiag(A_trace*1e-8);
 
     double coeff[vec4d_nEntries]={0.0, 0.0, 0.0, 0.0};
-    CholeskySolve4d(A_new, coeff);
+    CholeskySolve4d_UnitRHS(A_regularized, coeff);
     for (int i=0; i<vec4d_nEntries; i++) 
       VtxCoeffs(vtx,i)=coeff[i];
   });
@@ -174,12 +174,15 @@ void MPMesh::assemblyVtx1() {
         int vID = elm2VtxConn(elm,i+1)-1;
         double w_vtx=weight(mp,i); 
         double CoordDiffs[vec4d_nEntries] = {1, (vtxCoords(vID,0) - mpPositions(mp,0))/radius,
-                                   (vtxCoords(vID,1) - mpPositions(mp,1))/radius, 
-				   (vtxCoords(vID,2) - mpPositions(mp,2))/radius};
+                                                (vtxCoords(vID,1) - mpPositions(mp,1))/radius, 
+                                                (vtxCoords(vID,2) - mpPositions(mp,2))/radius};
+        
+        auto factor = w_vtx*(VtxCoeffs(vID,0) + VtxCoeffs(vID,1)*CoordDiffs[1] + 
+                                                VtxCoeffs(vID,2)*CoordDiffs[2] + 
+                                                VtxCoeffs(vID,3)*CoordDiffs[3]);
+  
         for (int k=0; k<numEntries; k++){
-          auto val = w_vtx*(VtxCoeffs(vID,0) + VtxCoeffs(vID,1)*CoordDiffs[1] + 
-                                               VtxCoeffs(vID,2)*CoordDiffs[2] + 
-                                               VtxCoeffs(vID,3)*CoordDiffs[3])*mpData(mp,k);
+          auto val = factor*mpData(mp,k);
           Kokkos::atomic_add(&reconVals(vID,k), val);
         }
       }

src/pmpo_utils.hpp

@@ -285,7 +285,7 @@ class Matrix4d {
 
     //Function to regularize the matrix
     KOKKOS_INLINE_FUNCTION
-    void regularize(double eps) {
+    void addToDiag(double eps) {
       for (int i = 1; i < rows_; i++) {
         data_[i][i]=data_[i][i]+eps;
       }
@@ -303,7 +303,7 @@ void initArray(Vec2d* arr, int n, Vec2d fill){
 
 
 KOKKOS_INLINE_FUNCTION
-void CholeskySolve4d(Matrix4d& A, double* x){
+void CholeskySolve4d_UnitRHS(Matrix4d& A, double* x){
 
     double a_00=A(0,0);
     if (A(0,0)==0){
@@ -355,20 +355,19 @@ void CholeskySolve4d(Matrix4d& A, double* x){
     }
     else{
        x[0]= 1.0/A(0,0);
-       x[1]= -A(0,1)*x[0]/A(1,1);                //- m12 % array(iVertex)*a0(iVertex)/m22 % array(iVertex)
-       x[2]= -(A(0,2)*x[0]+A(1,2)*x[1])/A(2,2);  //-(m13 % array(iVertex)*a0(iVertex) + m23 % array(iVertex)*a1(iVertex))/m33 % array(iVertex)
-       x[3]= -(A(0,3)*x[0]+A(1,3)*x[1]+A(2,3)*x[2])/A(3,3); //-(m14 % array(iVertex)*a0(iVertex) + m24 % array(iVertex)*a1(iVertex) + m34 % array(iVertex)*a2(iVertex))/m44 % array(iVertex)
+       x[1]= -A(0,1)*x[0]/A(1,1);                
+       x[2]= -(A(0,2)*x[0]+A(1,2)*x[1])/A(2,2);  
+       x[3]= -(A(0,3)*x[0]+A(1,3)*x[1]+A(2,3)*x[2])/A(3,3); 
 
        x[3] = x[3]/A(3,3);
        x[2] = ( x[2] - A(2,3)*x[3] )/A(2,2);
        x[1] = ( x[1] - A(1,2)*x[2] - A(1,3)*x[3])/A(1,1);
        x[0] = ( x[0] - A(0,1)*x[1] - A(0,2)*x[2] - A(0,3)*x[3])/A(0,0);
-
     } 
 }
 
 KOKKOS_INLINE_FUNCTION
-void QR_decomp4d(Matrix4d& A, Matrix4d& Q, Matrix4d& R){
+void QRDecomp4d(Matrix4d& A, Matrix4d& Q, Matrix4d& R){
   int m_size=4;
   for (int i=0; i<m_size; i++){
       Vec4d A_column;
@@ -462,6 +461,7 @@ Vec3d xyz_from_lat_lon(double lat, double lon, double r){
   return xyz;
 }
 
+//Rotate 90 degrees about Y axis
 KOKKOS_INLINE_FUNCTION
 Vec3d grid_rotation_backward(Vec3d& xyz_input){
   Vec3d xyz_output;

test/testFortranMPAdvection.f90

@@ -67,13 +67,11 @@ subroutine runAdvectionTest2(mpMesh, numPush, latVertex, lonVertex, nEdgesOnCell
     integer, dimension(:,:), pointer :: verticesOnCell
     real(kind=MPAS_RKIND) :: sphereRadius
 
-    PRINT *, "Foward: "
     do i = 1, numPush
       call calcSurfDispIncr(mpMesh, latVertex, lonVertex, nEdgesOnCell, verticesOnCell, nVertices, sphereRadius)
       call polympo_push(mpMesh)
     end do
 
-    PRINT *, "Backward: "
     call calcSurfDispIncr(mpMesh, latVertex, lonVertex, nEdgesOnCell, verticesOnCell, nVertices, sphereRadius, -numPush)
     call polympo_push(mpMesh)
 

test/testFortranMPReconstruction.f90

@@ -101,9 +101,9 @@ program main
     j = verticesOnCell(1,i)
     mpLatLon(1,i) = latVertex(j)
     mpLatLon(2,i) = lonVertex(j) 
-    mpPosition(1,i) = xCell(i) !xVertex(j)
-    mpPosition(2,i) = yCell(i) !yVertex(j)
-    mpPosition(3,i) = zCell(i) !zVertex(j)
+    mpPosition(1,i) = xCell(i)
+    mpPosition(2,i) = yCell(i)
+    mpPosition(3,i) = zCell(i)
   end do
 
   call polympo_createMPs(mpMesh,nCells,numMPs,c_loc(mpsPerElm),c_loc(mp2Elm),c_loc(isMPActive))