Skip to content

Commit abc23be

Browse files
dhyan1272dhyan1272
committed
New regularization and matrix inversion
1 parent c0fbe6a commit abc23be

File tree

4 files changed

+276
-31
lines changed

4 files changed

+276
-31
lines changed

‎src/pmpo_MPMesh.cpp‎

Lines changed: 21 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -7,14 +7,18 @@ namespace polyMPO{
77

88
void printVTP_mesh(MPMesh& mpMesh, int printVTPIndex=-1);
99

10-
void MPMesh::calcBasis() {
10+
void MPMesh::calcBasis(bool use3DArea) {
1111
assert(p_mesh->getGeomType() == geom_spherical_surf);
1212
auto MPsPosition = p_MPs->getPositions();
1313
auto mp_basis_field = p_MPs->getData<MPF_Basis_Vals>(); // we can implement MPs->getBasisVals() like MPs->getPositions()
14+
1415
auto elm2VtxConn = p_mesh->getElm2VtxConn();
1516
auto vtxCoords = p_mesh->getMeshField<MeshF_VtxCoords>();
1617
double radius = p_mesh->getSphereRadius();
17-
18+
//For Gnomonic Projection
19+
auto gnomProjVtx = p_mesh->getMeshField<polyMPO::MeshF_VtxGnomProj>();
20+
auto gnomProjElmCenter = p_mesh->getMeshField<polyMPO::MeshF_ElmCenterGnomProj>();
21+
1822
auto calcbasis = PS_LAMBDA(const int& elm, const int& mp, const int& mask) {
1923
if(mask) { //if material point is 'active'/'enabled'
2024
Vec3d position3d(MPsPosition(mp,0),MPsPosition(mp,1),MPsPosition(mp,2));
@@ -30,15 +34,22 @@ void MPMesh::calcBasis() {
3034
v3d[numVtx][1] = vtxCoords(elm2VtxConn(elm,1)-1,1);
3135
v3d[numVtx][2] = vtxCoords(elm2VtxConn(elm,1)-1,2);
3236

33-
double basisByArea3d[maxVtxsPerElm] = {0.0};
34-
initArray(basisByArea3d,maxVtxsPerElm,0.0);
35-
36-
// calc basis
37-
getBasisByAreaGblFormSpherical(position3d, numVtx, v3d, radius, basisByArea3d);
38-
37+
double basisByArea[maxVtxsPerElm] = {0.0};
38+
initArray(basisByArea,maxVtxsPerElm,0.0);
39+
40+
if(!use3DArea){
41+
double mpProjX, mpProjY;
42+
auto gnomProjElmCenter_sub = Kokkos::subview(gnomProjElmCenter, elm, Kokkos::ALL);
43+
computeGnomonicProjectionAtPoint(position3d, gnomProjElmCenter_sub, mpProjX, mpProjY);
44+
auto gnom_vtx_subview = Kokkos::subview(gnomProjVtx, elm, Kokkos::ALL, Kokkos::ALL);
45+
compute2DplanarTriangleArea(numVtx, gnom_vtx_subview, mpProjX, mpProjY, basisByArea);
46+
}
47+
else{
48+
getBasisByAreaGblFormSpherical(position3d, numVtx, v3d, radius, basisByArea);
49+
}
3950
// fill step
4051
for(int i=0; i<= numVtx; i++){
41-
mp_basis_field(mp,i) = basisByArea3d[i];
52+
mp_basis_field(mp,i) = basisByArea[i];
4253
}
4354
}
4455
};
@@ -302,7 +313,7 @@ void MPMesh::T2LTracking(Vec2dView dx){
302313
void MPMesh::reconstructSlices() {
303314
if (reconstructSlice.size() == 0) return;
304315
Kokkos::Timer timer;
305-
calcBasis();
316+
calcBasis(true);
306317
resetPreComputeFlag();
307318
for (auto const& [index, reconstruct] : reconstructSlice) {
308319
if (reconstruct) reconstruct();

‎src/pmpo_MPMesh.hpp‎

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -28,6 +28,7 @@ class MPMesh{
2828
void computeMatricesAndSolve();
2929
void resetPreComputeFlag();
3030
Kokkos::View<double*[vec4d_nEntries]> precomputedVtxCoeffs;
31+
Kokkos::View<double*[vec3d_nEntries][vec4d_nEntries]> precomputedVtxCoeffs_new;
3132

3233
Mesh* p_mesh;
3334
MaterialPoints* p_MPs;
@@ -70,7 +71,7 @@ class MPMesh{
7071
void push_swap();
7172
void push_swap_pos();
7273
void push();
73-
void calcBasis();
74+
void calcBasis(bool use3DArea);
7475

7576
DoubleView assemblyV0();
7677
template <MaterialPointSlice index>
@@ -100,7 +101,8 @@ class MPMesh{
100101
template <MeshFieldIndex meshFieldIndex>
101102
void reconstruct_full();
102103
void reconstruct_coeff_full();
103-
void solveMatrix(const Kokkos::View<double**>& vtxMatrices, double& radius, bool scaling);
104+
void solveMatrix(const Kokkos::View<double**>& vtxMatrices, double& radius, bool scaling);
105+
void invertMatrix(const Kokkos::View<double**>& vtxMatrices, const double& radius);
104106

105107
template<MeshFieldIndex meshFieldIndex>
106108
void setReconstructSlice(int order, MeshFieldType type);

‎src/pmpo_MPMesh_assembly.hpp‎

Lines changed: 165 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -297,7 +297,7 @@ void MPMesh::subAssemblyCoeffs(int vtxPerElm, int nCells, double* m11, double* m
297297
auto dual_triangle_area=p_mesh->getMeshField<MeshF_DualTriangleArea>();
298298

299299
//Material Points
300-
calcBasis();
300+
calcBasis(true);
301301
auto weight = p_MPs->getData<MPF_Basis_Vals>();
302302
auto mpPos = p_MPs->getData<MPF_Cur_Pos_XYZ>();
303303
auto mpAppID = p_MPs->getData<polyMPO::MPF_MP_APP_ID>();
@@ -746,7 +746,9 @@ void MPMesh::reconstruct_coeff_full(){
746746
auto dual_triangle_area=p_mesh->getMeshField<MeshF_DualTriangleArea>();
747747

748748
//Material Points
749-
calcBasis();
749+
bool use3DArea=false;
750+
calcBasis(use3DArea);
751+
750752
auto weight = p_MPs->getData<MPF_Basis_Vals>();
751753
auto mpPos = p_MPs->getData<MPF_Cur_Pos_XYZ>();
752754

@@ -760,7 +762,7 @@ void MPMesh::reconstruct_coeff_full(){
760762
if(p_mesh->getGeomType() == geom_spherical_surf)
761763
radius=p_mesh->getSphereRadius();
762764

763-
bool scaling=true;
765+
bool scaling=false;
764766

765767
//Assemble matrix for each vertex
766768
auto assemble = PS_LAMBDA(const int& elm, const int& mp, const int& mask) {
@@ -774,15 +776,15 @@ void MPMesh::reconstruct_coeff_full(){
774776
mScale=sqrt(dual_triangle_area(vID,0))/radius;
775777

776778
Kokkos::atomic_add(&vtxMatrices(vID,0), w_vtx*mScale*mScale);
777-
Kokkos::atomic_add(&vtxMatrices(vID,1), w_vtx*mScale*(vtxCoords(vID,0)-mpPos(mp,0))/radius);
778-
Kokkos::atomic_add(&vtxMatrices(vID,2), w_vtx*mScale*(vtxCoords(vID,1)-mpPos(mp,1))/radius);
779-
Kokkos::atomic_add(&vtxMatrices(vID,3), w_vtx*mScale*(vtxCoords(vID,2)-mpPos(mp,2))/radius);
780-
Kokkos::atomic_add(&vtxMatrices(vID,4), w_vtx*mScale*(vtxCoords(vID,0)-mpPos(mp,0))*(vtxCoords(vID,0)-mpPos(mp,0))/(radius*radius));
781-
Kokkos::atomic_add(&vtxMatrices(vID,5), w_vtx*mScale*(vtxCoords(vID,0)-mpPos(mp,0))*(vtxCoords(vID,1)-mpPos(mp,1))/(radius*radius));
782-
Kokkos::atomic_add(&vtxMatrices(vID,6), w_vtx*mScale*(vtxCoords(vID,0)-mpPos(mp,0))*(vtxCoords(vID,2)-mpPos(mp,2))/(radius*radius));
783-
Kokkos::atomic_add(&vtxMatrices(vID,7), w_vtx*mScale*(vtxCoords(vID,1)-mpPos(mp,1))*(vtxCoords(vID,1)-mpPos(mp,1))/(radius*radius));
784-
Kokkos::atomic_add(&vtxMatrices(vID,8), w_vtx*mScale*(vtxCoords(vID,1)-mpPos(mp,1))*(vtxCoords(vID,2)-mpPos(mp,2))/(radius*radius));
785-
Kokkos::atomic_add(&vtxMatrices(vID,9), w_vtx*mScale*(vtxCoords(vID,2)-mpPos(mp,2))*(vtxCoords(vID,2)-mpPos(mp,2))/(radius*radius));
779+
Kokkos::atomic_add(&vtxMatrices(vID,1), w_vtx*mScale*(-vtxCoords(vID,0)+mpPos(mp,0))/radius);
780+
Kokkos::atomic_add(&vtxMatrices(vID,2), w_vtx*mScale*(-vtxCoords(vID,1)+mpPos(mp,1))/radius);
781+
Kokkos::atomic_add(&vtxMatrices(vID,3), w_vtx*mScale*(-vtxCoords(vID,2)+mpPos(mp,2))/radius);
782+
Kokkos::atomic_add(&vtxMatrices(vID,4), w_vtx*mScale*(-vtxCoords(vID,0)+mpPos(mp,0))*(-vtxCoords(vID,0)+mpPos(mp,0))/(radius*radius));
783+
Kokkos::atomic_add(&vtxMatrices(vID,5), w_vtx*mScale*(-vtxCoords(vID,0)+mpPos(mp,0))*(-vtxCoords(vID,1)+mpPos(mp,1))/(radius*radius));
784+
Kokkos::atomic_add(&vtxMatrices(vID,6), w_vtx*mScale*(-vtxCoords(vID,0)+mpPos(mp,0))*(-vtxCoords(vID,2)+mpPos(mp,2))/(radius*radius));
785+
Kokkos::atomic_add(&vtxMatrices(vID,7), w_vtx*mScale*(-vtxCoords(vID,1)+mpPos(mp,1))*(-vtxCoords(vID,1)+mpPos(mp,1))/(radius*radius));
786+
Kokkos::atomic_add(&vtxMatrices(vID,8), w_vtx*mScale*(-vtxCoords(vID,1)+mpPos(mp,1))*(-vtxCoords(vID,2)+mpPos(mp,2))/(radius*radius));
787+
Kokkos::atomic_add(&vtxMatrices(vID,9), w_vtx*mScale*(-vtxCoords(vID,2)+mpPos(mp,2))*(-vtxCoords(vID,2)+mpPos(mp,2))/(radius*radius));
786788
}
787789
}
788790
};
@@ -805,6 +807,32 @@ void MPMesh::reconstruct_coeff_full(){
805807
pumipic::RecordTime("Communicate Matrix Values" + std::to_string(self), timer.seconds());
806808

807809
solveMatrix(vtxMatrices, radius, scaling);
810+
811+
invertMatrix(vtxMatrices, radius);
812+
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+
808836
}
809837

810838
void MPMesh::solveMatrix(const Kokkos::View<double**>& vtxMatrices, double& radius, bool scaling){
@@ -849,6 +877,111 @@ void MPMesh::solveMatrix(const Kokkos::View<double**>& vtxMatrices, double& radi
849877
pumipic::RecordTime("SolveMatrix" + std::to_string(self), timer.seconds());
850878
}
851879

880+
void MPMesh::invertMatrix(const Kokkos::View<double**>& vtxMatrices, const double& radius){
881+
882+
int nVertices = p_mesh->getNumVertices();
883+
auto vtxCoords = p_mesh->getMeshField<polyMPO::MeshF_VtxCoords>();
884+
auto dual_triangle_area = p_mesh->getMeshField<MeshF_DualTriangleArea>();
885+
bool isRotated=false;
886+
887+
Kokkos::View<double*[3][vec4d_nEntries]> VtxCoeffs("VtxCoeffs", nVertices);
888+
Kokkos::parallel_for("invertMatrix", nVertices, KOKKOS_LAMBDA(const int vtx){
889+
auto relativeScale = vtxMatrices(vtx, 0)*dual_triangle_area(vtx, 0)/(radius*radius);
890+
if(isRotated){
891+
892+
}
893+
auto cosLat = sqrt(pow(vtxCoords(vtx,0)/radius, 2) + pow(vtxCoords(vtx,1)/radius, 2));
894+
auto vtx_area_sqrt = sqrt(dual_triangle_area(vtx,0)/(radius*radius));
895+
Vec3d v0 = { -vtxCoords(vtx, 1)/(radius*cosLat), -vtxCoords(vtx, 2)*vtxCoords(vtx, 0)/(radius*radius*cosLat), vtxCoords(vtx, 0)/(radius*vtx_area_sqrt) };
896+
Vec3d v1 = { vtxCoords(vtx, 0)/(radius*cosLat), -vtxCoords(vtx, 1)*vtxCoords(vtx, 2)/(radius*radius*cosLat), vtxCoords(vtx, 1)/(radius*vtx_area_sqrt) };
897+
Vec3d v2 = { 0, cosLat, vtxCoords(vtx, 2)/(radius*vtx_area_sqrt) };
898+
899+
Matrix3d rotateScaleM = {v0, v1, v2};
900+
if(isRotated){
901+
902+
}
903+
904+
double invM11 = 1.0 / vtxMatrices(vtx, 0);
905+
Matrix3d subM;
906+
subM(0, 0) = vtxMatrices(vtx, 4) - invM11 * vtxMatrices(vtx, 1) * vtxMatrices(vtx, 1);
907+
subM(0, 1) = vtxMatrices(vtx, 5) - invM11 * vtxMatrices(vtx, 1) * vtxMatrices(vtx, 2);
908+
subM(0, 2) = vtxMatrices(vtx, 6) - invM11 * vtxMatrices(vtx, 1) * vtxMatrices(vtx, 3);
909+
subM(1, 1) = vtxMatrices(vtx, 7) - invM11 * vtxMatrices(vtx, 2) * vtxMatrices(vtx, 2);
910+
subM(1, 2) = vtxMatrices(vtx, 8) - invM11 * vtxMatrices(vtx, 2) * vtxMatrices(vtx, 3);
911+
subM(2, 2) = vtxMatrices(vtx, 9) - invM11 * vtxMatrices(vtx, 3) * vtxMatrices(vtx, 3);
912+
subM(1, 0) = subM(0, 1);
913+
subM(2, 0) = subM(0, 2);
914+
subM(2, 1) = subM(1, 2);
915+
916+
auto subM1 = (rotateScaleM.transpose())*(subM*rotateScaleM);
917+
918+
auto trG = subM1(0, 0) + subM1(1, 1);
919+
if((trG < 1e-4 * relativeScale) || (vtxMatrices(vtx, 0)<0.01)){
920+
VtxCoeffs(vtx, 0, 0) = invM11;
921+
return;
922+
}
923+
924+
Vec3d mVec = {vtxMatrices(vtx, 1), vtxMatrices(vtx, 2), vtxMatrices(vtx, 3)};
925+
auto blockC = rotateScaleM * mVec;
926+
927+
auto delG = sqrt(4.0 * pow(subM1(0, 1), 2) + pow( subM1(0, 0)-subM1(1, 1), 2));
928+
auto minEig = 0.5 * (trG-delG);
929+
double regularization = max(0.1 * relativeScale - minEig, 0.0);
930+
subM1(0, 0) = subM1(0, 0) + regularization;
931+
subM1(1, 1) = subM1(1, 1) + regularization;
932+
933+
double minZ2 = max(pow(subM1(0, 2), 2)/subM1(0, 0), pow(subM1(1, 2), 2)/subM1(1, 1)) + 0.1 * minEig;
934+
subM1(2, 2) = max(subM1(2, 2), minZ2);
935+
936+
double invM2D[6]={0.0};
937+
invM2D[0] = subM1(2, 2) * subM1(1, 1) - subM1(1, 2) * subM1(1, 2);
938+
invM2D[1] = subM1(0, 2) * subM1(1, 2) - subM1(2, 2) * subM1(0, 1);
939+
invM2D[2] = subM1(0, 1) * subM1(1, 2) - subM1(0, 2) * subM1(1, 1);
940+
invM2D[3] = subM1(2, 2) * subM1(0, 0) - subM1(0, 2) * subM1(0, 2);
941+
invM2D[4] = subM1(0, 1) * subM1(0, 2) - subM1(0, 0) * subM1(1, 2);
942+
invM2D[5] = subM1(0, 0) * subM1(1, 1) - subM1(0, 1) * subM1(0, 1);
943+
double det = subM1(0, 0) * invM2D[0] + subM1(0, 1) * invM2D[1] + subM1(0, 2) * invM2D[2];
944+
for (int i=0; i<6; i++) invM2D[i] = invM2D[i]/det;
945+
946+
Vec3d iBlockC(0, 0, 0);
947+
iBlockC[0] = blockC[0] * invM2D[0] + blockC[1] * invM2D[1] + blockC[2] * invM2D[2];
948+
iBlockC[1] = blockC[0] * invM2D[1] + blockC[1] * invM2D[3] + blockC[2] * invM2D[4];
949+
iBlockC[2] = blockC[0] * invM2D[2] + blockC[1] * invM2D[4] + blockC[2] * invM2D[5];
950+
iBlockC = iBlockC*invM11;
951+
952+
VtxCoeffs(vtx, 0, 0) = invM11 + invM11*iBlockC.dot(blockC);
953+
auto temp = -rotateScaleM.rightMultiply(iBlockC);
954+
VtxCoeffs(vtx, 0, 1) = temp[0];
955+
VtxCoeffs(vtx, 0, 2) = temp[1];
956+
VtxCoeffs(vtx, 0, 3) = temp[2];
957+
958+
if(vtx == 182){
959+
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+
979+
}
980+
});
981+
this->precomputedVtxCoeffs_new = VtxCoeffs;
982+
983+
}
984+
852985
template <MeshFieldIndex meshFieldIndex>
853986
void MPMesh::reconstruct_full() {
854987
Kokkos::Timer timer;
@@ -859,6 +992,7 @@ void MPMesh::reconstruct_full() {
859992
MPI_Comm_size(comm, &numProcsTot);
860993

861994
auto VtxCoeffs=this->precomputedVtxCoeffs;
995+
auto VtxCoeffs_new=this->precomputedVtxCoeffs_new;
862996

863997
//Mesh Information
864998
auto elm2VtxConn = p_mesh->getElm2VtxConn();
@@ -889,14 +1023,19 @@ void MPMesh::reconstruct_full() {
8891023
for(int i=0; i<nVtxE; i++){
8901024
int vID = elm2VtxConn(elm,i+1)-1;
8911025
double w_vtx=weight(mp,i);
892-
double CoordDiffs[vec4d_nEntries] = {1, (vtxCoords(vID,0) - mpPositions(mp,0))/radius,
893-
(vtxCoords(vID,1) - mpPositions(mp,1))/radius,
894-
(vtxCoords(vID,2) - mpPositions(mp,2))/radius};
1026+
double CoordDiffs[vec4d_nEntries] = {1, (-vtxCoords(vID,0) + mpPositions(mp,0))/radius,
1027+
(-vtxCoords(vID,1) + mpPositions(mp,1))/radius,
1028+
(-vtxCoords(vID,2) + mpPositions(mp,2))/radius};
8951029

896-
auto factor = w_vtx*(VtxCoeffs(vID,0) + VtxCoeffs(vID,1)*CoordDiffs[1] +
1030+
/*auto factor = w_vtx*(VtxCoeffs(vID,0) + VtxCoeffs(vID,1)*CoordDiffs[1] +
8971031
VtxCoeffs(vID,2)*CoordDiffs[2] +
8981032
VtxCoeffs(vID,3)*CoordDiffs[3]);
1033+
*/
1034+
auto factor = w_vtx*(VtxCoeffs_new(vID,0, 0) + VtxCoeffs_new(vID,0, 1)*CoordDiffs[1] +
1035+
VtxCoeffs_new(vID,0, 2)*CoordDiffs[2] +
1036+
VtxCoeffs_new(vID,0, 3)*CoordDiffs[3]);
8991037

1038+
if(vID == 182) printf("Weight W %.15e\n", w_vtx);
9001039
for (int k=0; k<numEntries; k++){
9011040
auto val = factor*mpData(mp,k);
9021041
Kokkos::atomic_add(&meshField(vID,k), val);
@@ -912,6 +1051,16 @@ void MPMesh::reconstruct_full() {
9121051
if(numProcsTot>1)
9131052
communicate_and_take_halo_contributions(meshField, numVertices, numEntries, 0, 0);
9141053
pumipic::RecordTime("Communicate Field Values" + std::to_string(self), timer.seconds());
1054+
1055+
//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));
1061+
}
1062+
});
1063+
9151064
}
9161065

9171066
void MPMesh::communicate_and_take_halo_contributions(const Kokkos::View<double**>& meshField, int nEntities, int numEntries, int mode, int op){

0 commit comments

Comments
 (0)