@@ -785,6 +785,9 @@ void MPMesh::reconstruct_coeff_full(){
785785 Kokkos::atomic_add (&vtxMatrices (vID,7 ), w_vtx*mScale *(-vtxCoords (vID,1 )+mpPos (mp,1 ))*(-vtxCoords (vID,1 )+mpPos (mp,1 ))/(radius*radius));
786786 Kokkos::atomic_add (&vtxMatrices (vID,8 ), w_vtx*mScale *(-vtxCoords (vID,1 )+mpPos (mp,1 ))*(-vtxCoords (vID,2 )+mpPos (mp,2 ))/(radius*radius));
787787 Kokkos::atomic_add (&vtxMatrices (vID,9 ), w_vtx*mScale *(-vtxCoords (vID,2 )+mpPos (mp,2 ))*(-vtxCoords (vID,2 )+mpPos (mp,2 ))/(radius*radius));
788+ if (vID==315459 ){
789+ printf (" MP %d in Cell %d positions %.15e %.15e %.15e \n " , mp, elm, mpPos (mp,0 ), mpPos (mp,1 ), mpPos (mp,2 ));
790+ }
788791 }
789792 }
790793 };
@@ -810,29 +813,6 @@ void MPMesh::reconstruct_coeff_full(){
810813
811814 invertMatrix (vtxMatrices, radius);
812815
813- Kokkos::parallel_for (" printSymmetricBlock" , 185 +1 , KOKKOS_LAMBDA (const int vtx){
814- if (vtx ==182 ) {
815- printf (" Vertex %d Matrix\n " , vtx);
816-
817- // Map indices to symmetric 4x4 entries
818- double m11 = vtxMatrices (vtx,0 );
819- double m12 = vtxMatrices (vtx,1 );
820- double m13 = vtxMatrices (vtx,2 );
821- double m14 = vtxMatrices (vtx,3 );
822- double m22 = vtxMatrices (vtx,4 );
823- double m23 = vtxMatrices (vtx,5 );
824- double m24 = vtxMatrices (vtx,6 );
825- double m33 = vtxMatrices (vtx,7 );
826- double m34 = vtxMatrices (vtx,8 );
827- double m44 = vtxMatrices (vtx,9 );
828- // Print symmetric 4x4 block
829- printf (" %25.15e %25.15e %25.15e %25.15e\n " , m11, m12, m13, m14);
830- printf (" %25.15e %25.15e %25.15e %25.15e\n " , m12, m22, m23, m24);
831- printf (" %25.15e %25.15e %25.15e %25.15e\n " , m13, m23, m33, m34);
832- printf (" %25.15e %25.15e %25.15e %25.15e\n\n " , m14, m24, m34, m44);
833- }
834- });
835-
836816}
837817
838818void MPMesh::solveMatrix (const Kokkos::View<double **>& vtxMatrices, double & radius, bool scaling){
@@ -879,6 +859,9 @@ void MPMesh::solveMatrix(const Kokkos::View<double**>& vtxMatrices, double& radi
879859
880860void MPMesh::invertMatrix (const Kokkos::View<double **>& vtxMatrices, const double & radius){
881861
862+ static int count_deb = 1 ;
863+ std::cout<<__FUNCTION__<<count_deb<<std::endl;
864+
882865 int nVertices = p_mesh->getNumVertices ();
883866 auto vtxCoords = p_mesh->getMeshField <polyMPO::MeshF_VtxCoords>();
884867 auto dual_triangle_area = p_mesh->getMeshField <MeshF_DualTriangleArea>();
@@ -954,32 +937,21 @@ void MPMesh::invertMatrix(const Kokkos::View<double**>& vtxMatrices, const doubl
954937 VtxCoeffs (vtx, 0 , 1 ) = temp[0 ];
955938 VtxCoeffs (vtx, 0 , 2 ) = temp[1 ];
956939 VtxCoeffs (vtx, 0 , 3 ) = temp[2 ];
957-
958- if (vtx == 182 ){
940+
941+ // Debugging
942+ if (vtx == 315459 ){
943+ printf (" Matrices %d vtx \n " , vtx);
944+ printf (" [ %.15e %.15e %.15e %.15e ]\n " , vtxMatrices (vtx,0 ), vtxMatrices (vtx,1 ), vtxMatrices (vtx,2 ), vtxMatrices (vtx,3 ));
945+ printf (" [ %.15e %.15e %.15e %.15e ]\n " , vtxMatrices (vtx,1 ), vtxMatrices (vtx,4 ), vtxMatrices (vtx,5 ), vtxMatrices (vtx,6 ));
946+ printf (" [ %.15e %.15e %.15e %.15e ]\n " , vtxMatrices (vtx,2 ), vtxMatrices (vtx,5 ), vtxMatrices (vtx,7 ), vtxMatrices (vtx,8 ));
947+ printf (" [ %.15e %.15e %.15e %.15e ]\n " , vtxMatrices (vtx,3 ), vtxMatrices (vtx,6 ), vtxMatrices (vtx,8 ), vtxMatrices (vtx,9 ));
959948 printf (" Coefficients %d vtx \n " , vtx);
960- // printf("XYZ %.15e %.15e %.15e \n", vtxCoords(vtx, 0), vtxCoords(vtx, 1), vtxCoords(vtx, 2) );
961- // printf("XYZ/R %.15e %.15e %.15e %.15e %.15e \n", vtxCoords(vtx, 0)/radius, vtxCoords(vtx, 1)/radius, vtxCoords(vtx, 2)/radius, vtx_area_sqrt, cosLat);
962- // printf("%.15e %.15e %.15e \n", rotateScaleM(0, 0), rotateScaleM(0, 1), rotateScaleM(0, 2));
963- // printf("%.15e %.15e %.15e \n", rotateScaleM(1, 0), rotateScaleM(1, 1), rotateScaleM(1, 2));
964- // printf("%.15e %.15e %.15e \n", rotateScaleM(2, 0), rotateScaleM(2, 1), rotateScaleM(2, 2));
965- // printf("%.15e %.15e %.15e \n", subM(0, 0), subM(0, 1), subM(0, 2));
966- // printf("%.15e %.15e %.15e \n", subM(1, 0), subM(1, 1), subM(1, 2));
967- // printf("%.15e %.15e %.15e \n", subM(2, 0), subM(2, 1), subM(2, 2));
968- // printf("%.15e %.15e %.15e \n", subM1(0, 0), subM1(0, 1), subM1(0, 2));
969- // printf("%.15e %.15e %.15e \n", subM1(1, 0), subM1(1, 1), subM1(1, 2));
970- // printf("%.15e %.15e %.15e \n", subM1(2, 0), subM1(2, 1), subM1(2, 2));
971- // printf("%.15e %.15e %.15e \n", mVec[0], mVec[1], mVec[2]);
972- // printf("%.15e %.15e %.15e \n", blockC[0], blockC[1], blockC[2]);
973- // printf("%.15e %.15e %.15e %.15e \n", relativeScale, trG, delG, minEig);
974- // printf("%.15e %.15e %.15e \n", regularization, minZ2, det);
975- // for (int i=0; i<6; i++) printf("%.15e \n", invM2D[i]);
976- // for (int i=0; i<3; i++) printf("%.15e \n", iBlockC[i]);
977- for (int i=0 ; i<4 ; i++) printf (" %.15e \n " , VtxCoeffs (vtx, 0 , i));
978-
949+ for (int i=0 ; i<4 ; i++) printf (" %.15e " , VtxCoeffs (vtx, 0 , i));
950+ printf (" \n " );
979951 }
980952 });
981953 this ->precomputedVtxCoeffs_new = VtxCoeffs;
982-
954+ count_deb++;
983955}
984956
985957template <MeshFieldIndex meshFieldIndex>
@@ -1035,7 +1007,6 @@ void MPMesh::reconstruct_full() {
10351007 VtxCoeffs_new (vID,0 , 2 )*CoordDiffs[2 ] +
10361008 VtxCoeffs_new (vID,0 , 3 )*CoordDiffs[3 ]);
10371009
1038- if (vID == 182 ) printf (" Weight W %.15e\n " , w_vtx);
10391010 for (int k=0 ; k<numEntries; k++){
10401011 auto val = factor*mpData (mp,k);
10411012 Kokkos::atomic_add (&meshField (vID,k), val);
@@ -1053,11 +1024,10 @@ void MPMesh::reconstruct_full() {
10531024 pumipic::RecordTime (" Communicate Field Values" + std::to_string (self), timer.seconds ());
10541025
10551026 // Debug Delete
1056- Kokkos::parallel_for (" printSymmetricBlock" , 185 +1 , KOKKOS_LAMBDA (const int vtx){
1057- if (vtx ==182 ) {
1058- printf (" IceArea in %d \n " , vtx);
1059- printf (" %.15e %.15e %.15e %.15e \n " , VtxCoeffs_new (vtx, 0 , 0 ), VtxCoeffs_new (vtx,0 , 1 ), VtxCoeffs_new (vtx,0 , 2 ), VtxCoeffs_new (vtx, 0 , 3 ));
1060- printf (" %25.15e \n " , meshField (vtx, 0 ));
1027+ Kokkos::parallel_for (" printSymmetricBlock" , numVertices, KOKKOS_LAMBDA (const int vtx){
1028+ if (vtx ==315459 ) {
1029+ printf (" IceArea in %d: " , vtx);
1030+ printf (" %25.15e \n " , meshField (vtx, 0 ));
10611031 }
10621032 });
10631033
0 commit comments